Pharmacological approach: effects of BPA, BPS and the GPER-specific agonist (G-1) and antagonist (G-15) on ovine GC steroidogenesis

The treatments (48 h) had no effect on cell viability regardless of the method used (Supplementary Figure S1B): Live/Dead staining (Supplementary Figure S1B.a), dehydrogenase activity assay (Supplementary Figure S1B.b) and lactate dehydrogenase activity assay in the medium (Supplementary Figure S1B.c).

Progesterone secretion (Fig. 1A) was increased (p < 0.001) by 53% (58.22 ± 15.37 ng progesterone/mg protein) after G-1 treatment at 1 µM compared with the control (37.07 ± 4.68 ng progesterone/mg protein) (Supplementary Table 1). Conversely, G-15 at 10 μM decreased (p < 0.001) progesterone secretion (26.39 ± 2.92 ng progesterone/mg protein) by 26% compared with the control. Similarly, BPA and BPS also decreased (p < 0.001) progesterone secretion by 16% and 14%, respectively, at 10 µM (30.31 ± 2.95 and 30.33 ± 2.60 ng progesterone/mg protein, respectively), and by 38% and 20%, respectively, at 50 µM (22.59 ± 2.52 and 27.79 ± 2.52 ng progesterone/mg protein, respectively) compared with the control. G-15 at 10 µM impaired progesterone secretion of ovine GC similarly to BPA or BPS. When G-15 was combined with BPA or BPS at 10 µM, progesterone secretion was decreased (p < 0.001) by 29% and 30%, respectively (25.45 ± 2.84 and 24.34 ± 2.23 ng progesterone/mg protein, respectively) compared with the control, and similarly to G-15 alone. When G-15 was combined with BPA or BPS at 50 µM, progesterone secretion was further decreased (p < 0.001) by 52% and 47%, respectively (17.08 ± 1.70 and 18.71 ± 1.83 ng progesterone/mg protein, respectively) compared with the control. This was also the case when compared with G-15 alone: progesterone was reduced by 35% (p < 0.001) after BPA 50 µM treatment, while it was reduced by 29% (p < 0.001) after BPS 50 µM treatment. When G-15 was combined with BPA or BPS at 50 µM, progesterone secretion was reduced by 28% (p = 0.003) or 34% (p < 0.001) compared with BPA or BPS at 50 µM.

Fig. 1

The effects of the GPER-specific agonist (G-1) or antagonist (G-15) and/or bisphenol A (BPA) and bisphenol (BPS) on ovine granulosa cells (GC) steroidogenesis. The progesterone (A) and oestradiol (B) concentrations were determined in culture medium after 48 h of culture in complemented serum-free McCoy’s 5A media in the presence or absence (control) of BPA or BPS at 10 or 50 μM, and/or G-15 at 10 µM or G-1 at 1 μM. The results are expressed as the mean ± standard error of the mean of 10 independent cultures. Each condition was performed in duplicate and normalised to the control condition of each culture experiment. Bars with different superscripts indicate a significant difference (p ≤ 0.05). The actual control values were 37.07 ± 4.68 ng progesterone/mg protein (A) and 52.04 ± 2.21 pg oestradiol/mg protein B

Oestradiol secretion (Fig. 1B) was not impacted significantly by G-1 at 1 µM or BPA at 10 or 50 µM (52.66 ± 7.78, 59.26 ± 10.36 and 54.67 ± 6.12 pg oestradiol/mg protein, respectively) compared with control (52.04 ± 5.21 pg oestradiol/mg protein) (Supplementary Table 1). G-15 at 10 μM decreased (p < 0.001) oestradiol secretion (23.02 ± 3.55 pg oestradiol/mg protein) by 56% compared with the control. BPS at 10 or 50 µM also decreased (p < 0.001) oestradiol secretion by 35% and 47%, respectively (32.90 ± 6.44 and 23.18 ± 2.87 pg oestradiol/mg protein, respectively) compared with the control. BPS at 10 or 50 µM alone impaired oestradiol secretion of ovine GC, similarly to G-15 at 10 µM. When G-15 was combined with BPS at 10 or 50 µM, oestradiol secretion was decreased (p = 0.006 and p = 0.019, respectively) by 73% and 71%, respectively (17.26 ± 3.09 and 16.68 ± 2.59 pg oestradiol/mg protein, respectively) compared with the control. This was also the case when compared with G-15: oestradiol was reduced by 38% (p = 0.006) after BPS 10 µM treatment, while it was reduced by 35% (p = 0.019) after BPS 50 µM treatment. When G-15 was combined with BPS at 10 or 50 µM, oestradiol secretion was reduced by 58% (p < 0.001) or 46% (p = 0.012) compared with BPS at 10 or 50 µM, respectively. On the other hand, the association of BPA at 10 or 50 µM with G-15 at 10 µM did not affect oestradiol secretion (32.25 ± 2.92 and 43.50 ± 4.91 pg oestradiol/mg protein, respectively) compared with the control.

Transcriptomic approach: effects of BPA, BPS and the GPER-specific agonist (G-1) on the ovine GC transcriptome

No effect on cell viability was observed after 1 h of treatment with G1 10 μM, BPA 50 µM or BPS 50 μM, using two complementary methods (Supplementary Figure S1A): Live/Dead staining (Supplementary Figure S1A.a) and the dehydrogenase activity assay (Supplementary Figure S1A.b).

To investigate the early mechanisms of action of BPA and BPS on the steroidogenesis of ovine GC, 24 ovine GC samples corresponding to four experimental conditions (control, 50 µM BPA, 50 µM BPS and 10 µM G-1) with six biological replicates for each condition were used for RNA-seq. The analysis generated an average of 23.7 ± 0.5 million reads per sample (Supplementary Table 2). Around 73.7% ± 0.2% of the reads were aligned uniquely with the sheep reference genome. This analysis was performed by using the genetic profiles of 14,883 genes generated by the RNA-seq experiment. Among them, 346 DEG were detected in the four comparisons: control versus BPA, control versus BPS, BPA versus BPS and control versus G1, with 259, 12, 69 and 75 DEG, respectively (Fig. 2A and Supplementary Table 3). Only three DEG, SREBF1, SPRY2 and TXNIP, were common between the control versus BPA, control versus BPS and control versus G-1 comparisons. The control versus BPA and control versus G-1 comparisons revealed 11 common DEG: TXNIP, JUN, NR4A1, SPRY2, C15orf39, ENSOARG00000018728 novel gene, TUBB, PPP1R3C, GADD45B, SREBF1 and STAR. However, the majority of DEG induced by BPA (200/259, 77%) and G-1 (69/75, 92%) were specific modifications, not found in other comparisons. Moreover, nine DEG (TXNIP, ENSOARG00000022521 novel gene, SPRY2, ENSOARG00000024133 novel gene, FAM217A, ANKRD1, DDIT4, MAPS2, SREBF1) were common between the control versus BPA and control versus BPS comparisons. Thus, the majority of DEG found in the control versus BPS comparison (9/12, 75%) were with the same DEG found in the control versus BPA comparison (Fig. 2A). A hierarchical clustering showed that BPA modified the expression of the genes in an opposite direction compared with G1. Indeed, cluster 1 regrouped the genes, which were upregulated by BPA, whereas cluster 2 showed the genes that were upregulated by G-1 (Supplementary Table 3 and Supplementary Table 4). In addition, after 1 h of treatment, BPS 50 μM only mildly modified the transcriptomic profile of the ovine GC compared with the control (Fig. 2B and Supplementary Table 4), as confirmed by a Principal Component Analysis (Fig. 2C).

Fig. 2

Differentially expressed genes (DEG) in ovine granulosa cells (GC) treated with bisphenol A (BPA), bisphenol S (BPS) or the GPER-specific agonist (G-1). After 1 h culture in complemented serum-free McCoy’s 5A media in the presence or absence (control) of BPA or BPS at 50 μM, or G-1 at 10 μM, six biological replicates were analysed with RNA-sequencing, and DEG were identified (padj ≤ 0.05). The Venn diagram (A) shows common and specific DEG from four comparisons: control versus BPA, control versus BPS, BPA versus BPS, and Control versus G-1. The heatmap (B) shows hierarchical clustering of DEG. The ordered list of heatmap genes have been annotated in Supplementary Table S3. Principle component analysis of gene expression in GC treated or not with BPA, BPS or G1 performed on expression values of 24 differential genes C

These results suggest that different genes were affected by BPA or BPS treatment; moreover, treatment with G-1, the specific agonist of GPER, affected a different set of genes. Therefore, the following analyses focused on the 283 DEG identified from three differential analyses: control versus BPA, control versus BPS and BPA versus BPS (Supplementary Table 3). From these DEG, 150 genes were downregulated by BPA (ARRDC4, KLF10, LIF, KLF11, EGR3, KLF9, ENC1, NOTCH1, and others) whereas 109 genes were upregulated by this molecule (CHAC1, JUN, LSMEM1, NR4A1, and others). Regarding BPS, only four genes were downregulated (RHOB, ANKRD1, SPRY2 and TXNIP) and eight were upregulated (RASL11B, PDK4, SREBF1, MASP2, DDIT4, FAM217A, ENSOARG00000024133 novel gene and ENSOARG00000022521 novel gene) compared with the control. In addition, TXNIP and SPRY2 were downregulated and six genes were upregulated (FAM217A, SREBF1, DDIT4, MAPS2, ENSOARG00000024133 novel gene, ENSOARG00000022521 novel gene) by both BPA and BPS compared with the control.

Effects of BPA and BPS on biological processes

A functional analysis based on GO analysis was performed with the ViSEAGO R package [48] by hierarchical clustering of enriched GO terms. There were 348 enriched GO terms among the three lists of DEG (Supplementary Fig. 2 and Supplementary Table 5). These terms were grouped into 18 clusters that could be further categorised into three major processes: metabolic processes (8 clusters and 137 GO terms), cellular communication (4 clusters and 80 GO terms) and structure and regulation of development (6 clusters and 131 GO terms) (Fig. 3 and Supplementary Table 5).

Fig. 3

Functional analysis of differentially expressed genes (DEG) from RNA-sequencing of ovine granulosa cells (GC). After 1 h of treatment in the presence or absence (control) of bisphenol A (BPA) or bisphenol S (BPS) at 50 μM, GC (six replicates per condition) were analysed with RNA-sequencing to obtain a list of DEG (p ≤ 0.05). The global clustering heatmap plot of functional sets of gene ontology (GO; p ≤ 0.01) terms was obtained by using ViSEAGO. From left to right are: the major processes, the cluster name, a heatmap of GO term counts from functional enrichment tests and a dendrogram based on Wang’s semantic similarity distance and Ward’s clustering criterion

The top 20 GO terms specifically altered by BPA correspond to three biological processes: metabolic processes (10 GO terms and 76/259 [29%] BPA DEG), transcription regulation (9 GO terms and 53/259 [20%] BPA DEG) and structure development (1 GO term and 25/259 [10%] BPA DEG). Regarding BPS, two others biological processes were impaired: cellular communication (18 GO terms) and cellular process (2 GO terms), representing 75% of BPS DEG (Fig. 2B). These five biological processes were also highlighted in the BPA versus BPS comparison. However, 18/20 GO (90%) corresponded to the effect of BPA while only 2/20 GO (10%) corresponded to the effect of BPS (Supplementary Table 6).

Interactions of DEG regulated by BPA

STRING analysis [49] was performed on 214 out of 259 (83%) BPA-regulated DEG to investigate the potential interaction network of the corresponding proteins (Supplementary Figure S3). Of note, the remaining DEG have not yet been annotated for the O. aries genome. This analysis highlighted the presence of several proteins involved in metabolic pathways (NDL4, MVD, SQLE, MAT2A, CYTB, ODC1, SETMAR, SGMS2, CYP2J, PIGH, DHCR24 and CHAC1) or in mitogen-activated protein kinase (MAPK) signalling pathways (SRF, ATF4, IL1A, MAP3K14, GADD45A, NR4A1 and GADD45B). These results confirm that metabolism is the primary biological process altered by BPA, therefore corroborating the GO analysis of DEG at the transcript level (Fig. 3). Using 12 DEG regulated by BPS, no common signalling pathways were identified.

Effects of BPA and BPS on the expression of most DEG after 12, 24 and 48 h

Among the list of 283 DEG, obtained from the analysis of RNA-seq, the four most upregulated genes (CHAC1, JUN, LSMEM1 and NR4A1) and the four most downregulated genes (TXNIP, ARRDC4, KLF10 and LIF) in BPA compared to control were selected (Supplementary Table 3), with a fold-change of 1.5 (|FC|> 1.5), an adjusted p-value ≤ 0.05 and a sufficient number of reads (> 118) to perform qPCR. qPCR analyses of expression of these eight genes were performed in GC 12, 24 or 48 h after treatment with BPA or BPS 50 µM (Fig. 4 and Supplementary Table 7).

Fig. 4

The effects of the GPER-specific antagonist (G-15) and/or bisphenol A (BPA) and bisphenol S (BPS) on gene expression of ovine granulosa cells (GC), according to changes over time. The expression of eight genes (CHAC1, JUN, LSMEM1, NR4A1, TXNIP, ARRDC4, KLF9 and LIF) were assessed in ovine GC. Gene expression was determined after 12, 24 and 48 h of culture in complemented serum-free McCoy’s 5A media in the presence or absence (control) of BPA or BPS at 50 μM. Total messenger RNA (mRNA) was extracted and reverse transcribed from ovine GC, then real-time polymerase chain reaction (qPCR) was performed. To normalise gene expression, the geometric mean of two housekeeping genes (β-actin [ACTB] and ribosomal protein L19 [RPL19]) was used. The results are expressed as the mean ± standard error of the mean of eight independent cultures and normalised to the mean of the control condition at 12 h. Statistical analysis was performed between the conditions of a same time (12, 24 or 48 h). Bars with different superscripts indicate a significant difference (p ≤ 0.05)

Similarly to 1 h of treatment, after 12 h of treatment, CHAC1 expression was significantly increased after both BPA (+ 93%, p < 0.001) and BPS (+ 66%, p < 0.001) treatment compared with the control. After 24 h and 48 h, no difference in CHAC1 expression was observed among the treatments (Fig. 4). Nevertheless, while in control condition CHAC1 expression increased as the culture time increased, the same pattern was not observed after BPA or BPS treatment (Supplementary Table 7).

After 12, 24 or 48 h of treatment, no significant effect on JUN expression was observed in any condition, except BPS that increased JUN expression after 48 h (+ 43%, p = 0.002) compared with the control.

LSMEM1 expression was still increased after 12 h and 24 h treatment with BPA 50 μM (+ 66%, p < 0.001 and + 40%, p = 0.005, respectively) or BPS 50 μM (+ 50%, p = 0.007 and + 38%, p < 0.001, respectively) compared with the control. There was no difference in LSMEM1 gene expression after 48 h of treatment (Fig. 4).

While NR4A1 gene expression was upregulated after 1 h of BPA 50 µM treatment compared with control (Supplementary Table 3), there was no difference after 12 h of treatment in any condition. NR4A1 gene expression was reduced with BPA (-33%, p = 0.046) and BPS 50 µM (-33%, p = 0.030),after 24 h treatment. Similarly, NR4A1 gene expression was reduced with BPA 50 µM (-43%, p = 0.013) and BPS (-42%, p = 0.013).

Thus, among CHAC1, JUN, LSMEM1 and NR4A1 – all upregulated after 1 h of BPA 50 µM treatment according to RNseq analysis (Supplementary Table 3) – only LSMEM1 expression was still increased after 12 h and 24 h of treatment according to qPCR analysis (Fig. 4).

TXNIP was downregulated after 1 h of BPA 50 µM (Supplementary Table 3). This downregulation was also observed after 12 h of BPS 50 µM treatment (-22%, p = 0.014) (Fig. 4). No effect was observed after 24 h culture. After 48 h culture, both BPA 50 µM (+ 26%, p = 0.018) and BPS 50 µM (+ 24%, p = 0.002) increased TXNIP expression compared with the control. An increase of TXNIP expression is also observed over time for both BPA and BPS 50 µM, but not for the control (Supplementary Table 7).

ARRDC4 gene expression was not reduced anymore after 12 h and 24 h of BPA or BPS 50 µM treatments compared with control. No difference was observed in ARRDC4 gene expression after 48 h, regardless of the treatment (Fig. 4). ARRDC4 gene expression expression is also increased over time in all conditions (Supplementary Table 7).

KLF10 expression did not change after 12 h, 24 h or 48 h of treatment compared to the control. Two other members of the KLF family, KLF9 and KLF11, were also BPA-induced DEG. After 12 and 24 h of treatment, KLF9 expression was reduced after treatment with BPA 50 µM (-46%, p < 0.001, and -44%, p < 0.001, respectively) and BPS 50 µM (-43%, p < 0.001, and -33%, p = 0.070, respectively) compared with the control. After 48 h, KLF9 expression was still reduced after BPS 50 µM treatment (-38%, p = 0.008). KLF9 expression increased over time (Supplementary Table 7). KLF11 expression did not differ between conditions at any time (Supplementary Table 7).

Finally, LIF expression was still reduced after 24 h (-54%, p < 0.001) and 48 h (-54%, p < 0.001) BPA 50 µM treatment compared with the control. The same effect was observed after 24 h (-57%, p < 0.001) and 48 h (-53%, p < 0.001) BPS 50 µM treatment (Fig. 4). LIF expression increased over time (Supplementary Table 7).

Among the four genes (TXNIP, ARRDC4, KLF10 and LIF) that were downregulated after 1 h of BPA 50 µM treatment according to RNAseq analysis (Supplementary Table 3), only TXNIP and LIF were also downregulated after either 12 h or 24 h or 48 h according to qPCR analysis (Fig. 4).

Two other downregulated genes after BPA treatment compared to control were also analysed by qPCR (EGR3 and NOTCH1). No change in EGR3 expression was observed between conditions. NOTCH1 expression was decreased after 1 h of 50 µM BPA treatment (FC = 0.729, p = 0.036) compared with the control (Supplementary Table 3). On the contrary, after 12 h of treatment with 50 µM BPA or BPS, NOTCH1 expression was increased (+ 52%, p < 0.001, and + 64%, p < 0.001, respectively).