Low-dose exposure to bisphenols A, F and S of human primary adipocyte impacts coding and non-coding RNA profiles

Bisphenol A (BPA) exposure has been suspected to be associated with deleterious effects on health including obesity and metabolically-linked diseases. Although bisphenols F (BPF) and S (BPS) are BPA structural analogs commonly used in many marketed products as a replacement for BPA, only sparse toxicological data are available yet. Our objective was to comprehensively characterize bisphenols gene targets in a human primary adipocyte model, in order to determine whether they may induce cellular dysfunction, using chronic exposure at two concentrations: a “low-dose” similar to the dose usually encountered in human biological fluids and a higher dose. Therefore, BPA, BPF and BPS have been added at 10 nM or 10 μM during the differentiation of human primary adipocytes from subcutaneous fat of three non-diabetic Caucasian female patients. Gene expression (mRNA/lncRNA) arrays and microRNA arrays, have been used to assess coding and non-coding RNA changes. We detected significantly deregulated mRNA/lncRNA and miRNA at low and high doses. Enrichment in “cancer” and “organismal injury and abnormalities” related pathways was found in response to the three products. Some long intergenic non-coding RNAs and small nucleolar RNAs were differentially expressed suggesting that bisphenols may also activate multiple cellular processes and epigenetic modifications. The analysis of upstream regulators of deregulated genes highlighted hormones or hormone-like chemicals suggesting that BPS and BPF can be suspected to interfere, just like BPA, with hormonal regulation and have to be considered as endocrine disruptors. All these results suggest that as BPA, its substitutes BPS and BPF should be used with the same restrictions.

Introduction usually encountered for BPA in biological fluids in a general population [18]. A higher dose (10 μM), frequently used in literature despite poorly relevant to human exposure, was also assayed to detect potential dose-specific effects. In addition, we noted that non-monotonic dose-response relationships were demonstrated for bisphenol A [19], which dissuaded us from trying to get dose-response information and using a wider range of concentrations.
We therefore compared the effects of BPA, BPS and BPF added at 10 nM and at 10 μM in the differentiation medium of primary human pre-adipocytes, on different RNA levels. We used cell cultures derived from three different patients in order to avoid overestimating effects resulting from individual susceptibility and to focus on general mechanisms. In addition, our strategy allowed the study of the broad family of RNAs: mRNA, microRNAs (miRNAs) and other non-coding RNAs in order to determine a potential role of epigenetic mechanisms in bisphenols effects. A microRNA (miRNA) is a small non-coding RNA molecule (constituted of about 22 nucleotides) that plays a role in RNA silencing and post-transcriptional regulation of gene expression [20]. Small nucleolar RNAs (snoRNAs) are a class of small RNA molecules that primarily guide chemical modifications of other RNAs, mainly ribosomal RNAs, transfer RNAs and small nuclear RNAs or may have a role in the regulation of alternative splicing [21]. Long non-coding RNAs (lncRNA) are non-protein coding transcripts longer than 200 nucleotides. They play a role in the regulation of gene transcription and also control various aspects of post-transcriptional mRNA processing [22]. These different non-coding RNA classes are detected with our set of microarrays.
Using modern hypothesis-free genome-wide technologies may contribute to make significant progress in endocrine disruptor mechanistic knowledge and open the way to unsuspected responsive pathways of interest, either for basic science or for public health.

Materials and methods
Chemicals BPA, BPS, and BPF from Sigma-Aldrich (Darmstadt, Germany; S1 Fig) were diluted in dimethyl sulfoxide (DMSO) and added to the adipocyte differentiation medium. They were used at two concentrations (10 nM and 10 μM) for the microarray experiment.

Adipocyte cell culture
Primary pre-adipocytes coming from subcutaneous fat of non-diabetic Caucasian female patients were used (Lonza, Basel, Switzerland). Pre-adipocytes are precursor cells that develop into adipocytes when fully differentiated. They were supplied with the differentiation kit and the culture medium. Primary human white subcutaneous pre-adipocytes were cultured in PBM-2 medium. Confluent pre-adipocytes were induced to differentiate with PBM-2 medium supplemented with insulin, dexamethasone, isobutylmethylxanthine and indomethacin (all supplied by Lonza) for ten days according to the instructions of the manufacturer. They were induced to differentiate into adipocytes with either the selected chemical or vehicle (DMSO; Fig 1). After ten days, cells were harvested for total RNA extraction.

Nucleic acid extraction
Specific kits (Kit Small & Large RNA; Macherey Nagel) were used to extract all species of RNAs. The RNA concentration was determined by absorbance at 260 nm (A260), and the purity was estimated by determining the A260/A230 ratio with a Nanodrop spectrophotometer (Nanodrop Technologies, Wilmington, DE). RNA Integrity Number (RIN) was assessed with a Bioanalyzer (Agilent). This kit allowed the extraction of small and large RNAs with no degradation (RIN>9).

Microarrays
To measure different species of RNAs, the mRNA and lncRNA microarray SurePrint G3 Human V2 (Agilent Technologies) and miRNA microarray SurePrint G3 miRNA release 19.0 (Agilent Technologies) were used. The microarray experiments from Agilent Technologies were performed and analyzed by the Transcriptomics and Applied Genomics (TAG) team in the Institut Pasteur de Lille (Dr David Hot). This combination of microarrays allowed the detection of mRNAs, but also of long intergenic non-protein coding RNAs, small nucleolar RNAs and miRNAs.

Data normalization
Henceforth the term probe will be used to refer to a mRNA, lncRNA, snoRNA or miRNA probe. For mRNAs and other long non-coding RNAs, raw data were composed of 62,976 RNA expressions (probes). A first filter consisted in excluding non-expressed probes. To do so, a detection p-value, provided by Agilent which indicated how well the expressions of the gene were detected in a given sample, was used to filter out undetected probes. Probes, which were expressed in a certain proportion of the samples of interest, were kept for further analysis, leading to 28,244 probes (for the study of differentiation) and 22,935 probes (for the study of BPs on differentiation) in the final datasets. In addition, several normalization steps were applied, such as taking into account housekeeping genes or taking the log 2 of the measures. The normalization was performed with the R package limma [23].
For miRNAs, raw data were composed of a total of 62,344 miRNA probes which represented 2,027 miRNAs (30 probes per miRNA + normalization probes). After a preliminary between-array normalization, the first step of the normalization consisted in estimating a summary measure of the miRNA expression using a linear model to take into account the probe affinity over the 30 probes characterizing each one of the miRNA. Secondly, as for the mRNA data, a detection p-value was used to filter out non-expressed miRNA leading to consider 325 (differentiation) or 483 probes (BPs) as expressed. The normalization was performed with the R package AgiMicroRna [24].

Statistical analyses
For each one of the product (BPA, BPF, BPS), 12 assays of cell cultures differentiated without any product (negative control, 4 per patient), 12 assays of cell cultures differentiated with 10 nM (2 per patient) or 10 μM of the given product (2 per patient) were compared. For each product, differentially expressed genes were identified with a mixed model, where the variations of a given probe were explained by a fixed effect of the 10 nM concentration and/or the 10 μM concentration in comparison with the control condition, plus the patients which were not of direct interest, were considered as a random effect. P-values were corrected for multiple testing using Benjamini and Hochberg's method to control the false discovery rate with a 5% threshold. All the analyses were performed with the R package limma.

Biological analyses
Lists of genes significantly induced or repressed after exposure to BPA were uploaded into Ingenuity Pathway Analysis Software (IPA, Ingenuity Systems, www.ingenuity.com) for biological analysis by comparison with the Ingenuity Knowledge Database. These lists of altered genes were then processed to investigate the functional distribution of the genes, as defined by Gene Ontology. Datasets and known canonical pathway associations were measured by IPA by using a ratio of the number of genes from a dataset that map a specific pathway divided by the total number of genes which map this canonical pathway. Fisher's exact test was used to determine a p-value representing the significance of these associations.

Taqman RT-PCR
The relative quantification of mRNA was performed with a quantitative RT-PCR assay. 400 ng of total RNA was transcribed into cDNA using the cDNA Archive Kit (Life Technologies, Foster City, CA, USA). Each cDNA sample was analyzed by quantitative real-time PCR (qPCR) using the fluorescent TaqMan 5'-nuclease assays (Life Technologies). The analysis was performed with ViiA 7 detection system and software (Life Technologies). Gene expression was standardized to POLR2A expression (2 -ΔCt ). miRNAs were reverse transcribed using TaqMan Micro RNA Reverse Transcription Kit (Life Technologies). Taqman microRNA assays (Life Technologies) were used for miRNA quantification.

Results
All microarray data results were deposited on Gene Expression Omnibus website, (http:// www.ncbi.nlm.nih.gov/geo/) under the accession number GSE98682 which includes both series of mRNA/lncRNA data (GPL16699) and miRNA data (GPL19730). The samples series providing both the raw and normalized data is accessible through the accession numbers from GSM2609729 to GSM2609851.
Adipocyte differentiation: mRNA/lncRNA and miRNA expression Pre-adipocytes coming from subcutaneous fat of three non-diabetic obese females of European origin were used (age: 37, 52 and 38 years old and BMI: 30.3, 39.5 and 32.9 respectively). After ten days of differentiation, as expected, an accumulation of refringent lipid droplets was observed in the cytoplasm of the adipocytes of all three patients, evidencing their ability to fully differentiate.
The heatmap representation of the gene expression in the pre-adipocyte and adipocyte demonstrates a clear-cut set of upregulated and downregulated probes (S2 Fig). Among up-regulated genes classical adipocyte differentiation markers are encountered (LPL, ADIPOQ, FABP4, PLIN4,. . .) which confirms the successful and complete differentiation of the pre-adipocytes.
For mRNA/lncRNA, 28,244 probes were detected in undifferentiated or differentiated adipocytes. Among these, 122 were significantly differentially up-regulated with a log 2 FC>3 and 275 with a log 2 FC>2; 23 were down-regulated with a log 2 FC<-3 and 114 with a log 2 FC<-2 in differentiated adipocytes compared to undifferentiated adipocytes (S1A Table).
For miRNA, 325 probes were detected in undifferentiated or differentiated adipocytes. Among these, 4 were significantly differentially up-regulated with a log 2 FC>2 and 36 with a log 2 FC>1 and 2 were down-regulated with a log 2 FC<-2 and 23 with a log 2 FC<-1 in differentiated adipocytes compared to undifferentiated adipocytes (S1B Table).
The analysis of probes differentially expressed with a log 2 FC>2 for each subject demonstrated that 56% of differentially expressed probes were shared by at least two patients and 36% by the three patients. Nevertheless, this disparity can be turned to our advantage to help us to decipher the shared differentially expressed probes relevant to a general mechanism and not to an individual susceptibility.
For miRNAs, the up or down-regulation during adipocyte differentiation observed was less obvious with lower fold changes and a very high variability between patients. With a log 2 FC>2, no miRNA probe was shared between the three patients (S4A Fig). Shared miRNA probes were detected only with a log 2 FC>1 (S4B Fig). Only 7% of probes unregulated with a log 2 FC>1 in the first patient were also up-regulated with a log 2 FC>1 in the other two patients and 42% of probes up-regulated with a log 2 FC>1 were down-regulated in at least one of the other two patients. In opposition to the PCA of the mRNA/lncRNA probes, the three patients and the differentiation stage could not be clearly distinguished using the first two principal components for the miRNA probes (S4C and S4D Fig).
On the mRNA/lncRNA microarray, the expression of a specific gene can be measured with several probes. The probes refereeing to the same gene were not aggregated, but analyzed separately since they might capture alternative splicing. Among the 20 most differentially expressed mRNA probes in differentiated adipocytes according to p-values, only up-regulated genes were found (S1A Table). Classical genes of adipocyte differentiation, such as LPL, PLIN1/4, ADIPOQ or FABP4 were found within this top list. For miRNA probes, the observed fold changes were lower along with less significant p-values. In contrast to mRNAs, both up and down-regulated miRNAs were found (S1B Table). Two long intergenic non-coding RNAs (LINC01140, LINC01088) were up-regulated in differentiated adipocytes and one (LINC01048) was down-regulated compared to pre-adipocytes (S1A Table).
Alteration of the global profiles by exposure to bisphenols mRNA and lncRNA. After an exposition to 10 nM of BPA (or BPS or BPF) during adipocyte differentiation, using a mixed model, 846 probes coding for mRNA and lncRNA (respectively 1,173 and 216) were found to be significantly down-regulated and 417 (respectively 367 and 621) up-regulated when compared to the control cells without bisphenol (S2 Table). After the exposition to the higher dose of 10 μM of BPA (or BPS or BPF) during adipocyte differentiation, using a mixed model, 774 probes (respectively 748 and 1,228) were found to be significantly down-regulated and 1106 (respectively 323 and 989) up-regulated when compared to the control cells without bisphenol (S2 Table). Among the probes detected using mRNA/ lncRNA arrays, some long intergenic non-coding RNAs and some small nucleolar RNAs were differentially expressed after exposition to the bisphenols during differentiation (S2 Table).
Considering all the deregulated probes, several of these deregulated probes were shared between concentrations (Table 1) and between bisphenols with the majority in the same direction (Table 2). For example, A_33_P3321657 tagging HSPG2 is in the top 20 of most differentially expressed probes and down-regulated for BPS and BPF at both concentrations.
To complement the previous results, we represented the global expression profiles of the significantly differentially expressed genes according to the three bisphenols in a heatmap ( Fig  2). What is sticking is the great homogeneity within the profiles. Indeed when a gene is over or underexpressed in one concentration of a bisphenol, it tends to have a similar response to the others, both in terms of direction and magnitude of effect. It is worth noting that BPF at 10 nM seems to behave slightly differently and shows an overall reduced magnitude of effect, which might suggest that a higher dose of BPF is required to have a similar effect as BPA and S in terms of expression profile.
miRNA. After an exposition to 10 nM of BPA (or BPS or BPF) during adipocyte differentiation, 39 probes coding for miRNA (respectively 61 and 0) were found to be significantly down-regulated and 33 (respectively 40 and 12) up-regulated when compared to the control cells without bisphenol (S2 Table). After the exposition to the higher dose of 10 μM of BPA (or BPS or BPF) during adipocyte differentiation, 18 probes (respectively 11 and 50) were found to be significantly down-regulated and 13 (respectively 8 and 48) up-regulated when compared to the control cells without bisphenol (S2 Table). miRNA probes shared between concentrations and between bisphenols are displayed in Table 3 and Table 4. As for mRNA, shared probes are dysregulated in the same direction. The most significant fold-changes values for miRNA are higher than fold changes observed for mRNA/lncRNA.
In addition, we represented the global expression profiles of the significantly differentially expressed miRNAs according to the three bisphenols in a heatmap (Fig 3). We observed, just like for mRNAs/lncRNAs, an overall homogeneity within the expression profiles, although we clearly see a separation in the clustering. As previously stated for the mRNAs/lncRNAs, BPF at 10 nM shows an overall reduced expression and the response to BPF at 10 μM is extremely similar to the response to BPA and S at 10nM.
Analysis of biological functions. Similarities between the alterations of the profile were assessed using Ingenuity Pathway Analysis software. The mRNA differentially expressed in response to 10 nM and 10 μM exposure of all three bisphenols (BPs) were analyzed: top canonical associated pathways, upstream regulators as well as related diseases and disorders are displayed in Table 5. We found shared enriched pathways for all three BPs. For instance, eIF2 signaling was among the top canonical pathways both for BPA and BPF; ESR1, MYCN or MYC were found as upstream regulators; and "cancer" as well as "organismal injury and abnormalities" were shared between the three products. The same analysis performed on the lists of differentially expressed mRNA probes of all three bisphenols (BPs) for the two concentrations gave similar results, suggesting general dose-independent mechanisms.
Probes shared by the three bisphenols. The selection of the probes differentially expressed at a concentration of 10 nM for all three BPs led to highlight 40 (37) mRNA probes (genes), but no miRNA. Among these 40 probes, 16 were up-regulated and 24 were down-regulated (Table 6). Cytoskeleton and ribosomal proteins were found within this set of genes suggesting that crucial mechanisms of cell regulation and of protein traduction could be disrupted in presence of low dose bisphenols. The selection of the probes differentially expressed at a concentration of 10 μM for all three BPs led to highlight 161 (134) mRNA probes (genes), and only one miRNA (hsa-miR-4655-3p). Among these 161 probes, 33 were up-regulated and 128  were down-regulated (S3 Table). If only the probes differentially expressed in both concentrations of all three BPs were selected, we highlighted 12 (10) mRNA probes (genes). Among these ten genes, two were up-regulated and eight were down-regulated; they are highlighted with a grey background in Table 6 (10 nM) and with yellow filling in S3 Table (10 μM). From our analyses, only one miRNA (hsa mir-4655-3p) was differentially expressed in all three BPs, but only at a high concentration of 10 μM for BPA and at both concentrations for BPS and BPF.
In addition, we searched with IPA for potential upstream regulators of the latter set of ten genes differentiated for the three bisphenols at both concentrations (Table 7). Among the upstream regulators, beta-estradiol and ESR1 already known to be targeted by BPA were found as potential targets for all three BPs.
In order to validate the microarray data by quantitative PCR, we first screened classical housekeeping gene expressions. POLR2A (Polymerase (RNA) II Subunit A) demonstrated no significant differential expression for the three bisphenols at both doses and was chosen as housekeeping gene (S4 Table). Taqman qPCR analysis was performed using "best coverage" probes recommended by the manufacturer (S5 Table).
The direction and values of the expression fold-changes obtained by qPCR were consistent with the data from the microarray analysis for eight of the ten genes studied (S6 Table). SET and SORD showed an up-regulation when analyzed by Taqman qPCR which is opposite to the microarray data. In the microarray analysis, only one probe among three annotated for SET and SORD was differentially expressed therefore different transcripts were evaluated by microarray probes or by qPCR probes: this explains the observed discrepancies. Without SET and SORD probes, the correlation coefficient obtained by Spearman analysis was r = 0.708 (p = 2.56 x 10 −8 ) in expression between the data from the microarray analysis and from Taqman qPCR (S6 Table). mir-4655-3p expression analyzed by Taqman was very low and was only detectable above the threshold of 35 cycles.
We integrated data from the Comparative ToxicoGenomics Database (CTD, http://ctdbase. org/) which tells us if BPA has already been shown to modulate the expression of a given gene. CTD provides manually curated information about chemical-gene/protein interactions, coming from published microarray data. BPA is highly represented but no association has been curated for BPS and BPF yet. Among the 40 probes deregulated for the three bisphenols at 10 nM, all the probes, with the exception of the probes SNORD100 and ANKRD11, were found to interact with BPA in this database.

Discussion
We have intended here to make progress in toxicology analyses of common chemicals present in the human environment, by the use of human primary cells chronically exposed to these  Bisphenols A, S and F impacts on human primary adipocyte RNA profiles chemicals, at a concentration similar to the range found in human fluids in many studies [18], with a hypothesis-free genomic approach analyzing both coding and non-coding RNA levels. We considered this unbiased strategy important to analyze the differential effects of BPA and of its substitutes, BPS and BPF, on human cellular genomic profiles. We first validated the use of primary cell cultures derived from three different patients. Indeed, mRNA/lncRNA expression induced by pre-adipocyte differentiation was very consistent between subjects. Despite some expected variability, we were indeed able to highlight classical markers of differentiation (LPL, FABP4, ADIPOQ, PLIN1/4 . . .) demonstrating both the achievement of differentiation of the cells and the validity of our design. Concerning the miRNA study, we showed that the miRNA profile for adipocyte differentiation is very heterogeneous between individuals. Nevertheless, among the top probes (according to p-values) differentially expressed during pre-adipocyte differentiation, some are consistent with existing literature. miR-26a, which was down-regulated in our study, induces characteristics of brown adipocytes during human adipocyte differentiation [25]. miR-30c up-regulated in our study, was shown to promote human adipocyte differentiation [26]. miR-136, down-regulated in our study, was up-regulated during chondrogenesis of human adipose-derived stem cells suggesting that it can be considered as a player in the switch from white adipocyte to chondrocyte [27].
More importantly, these three cell cultures were sufficient to detect significantly differentially expressed RNA species in response to bisphenols, demonstrating an effect of bisphenols even at "low doses" on cellular metabolism. We showed here that both BPA and its substitutes BPS and BPF, induced a dysregulation of major regulators in cell homoeostasis, which is illustrated by highlighting cancer and translation pathways, whatever the concentration used. Among the mRNA/lncRNA differentially expressed probes for BPA, 46% were shared with BPS and only 10% with BPF. For miRNA probes, 80% were shared between BPA and BPS while none were shared between BPA and BPF. Interestingly, among the mRNA/lncRNA differentially expressed probes for BPA at 10 nM, 46% were shared with probes differentially expressed for BPF, when BPF was used at 10 μM. These results suggest that BPS effect is closer than BPF to BPA effect and that a higher concentration of BPF could be required for an effect on the same probes. Nevertheless, when focusing on the probes deregulated for the three bisphenols, we highlighted extracellular matrix and cytoskeleton genes, transcription regulators, cyclins: all genes implied in cell metabolism regulation. In addition, the analysis of upstream regulators of the set of ten genes differentiated for the three bisphenols at both concentrations highlighted hormones or hormone-like chemicals (beta-estradiol, dihydrosterone) as well as ESR1 (estrogen receptor 1). This suggests that BPS and BPF, can be suspected to interfere, just like BPA, with hormonal regulation, and therefore deserve to undergo the same restrictions from a regulation standpoint.
In this regard, we highlighted an impairment of basic pathways of cell metabolism by bisphenols, such as eIF2 (Eukaryotic Initiation Factor 2) signaling; eIF2 is a eukaryotic initiation factor required in the initiation of translation. Moreover, among the more represented disease categories, corresponding to differentially expressed species, we found "cancer" and "organismal injuries and abnormalities". It is worth noting that these oncogenic features cannot be attributed to the cellular model since we have used primary cells. Surprisingly, in our conditions, neither adipogenesis nor inflammation probes were overrepresented as it was observed in other cellular adipocyte models [28] [29]. BPA was demonstrated to significantly enhanced adipogenesis via an ER-mediated pathway [29]. BPA has been first considered a weak estrogen, based on its low binding affinity to the nuclear estrogen receptors, however there is also evidence of a potent BPA-estrogenic action through non-classical estrogen pathways e.g. via ERalpha located outside the nucleus [30]. In beta cells 1 nM BPA regulates insulin gene transcription via ERalpha in an estrogen response element independent manner [31]. Extranuclearly located ERbeta mediates rapid actions of BPA (1 nM) as well [32]. The analysis of upstream regulators of the set of ten genes differentiated for the three bisphenols at both concentrations highlighted ESR1 (estrogen receptor 1 or ERalpha). Therefore, our results suggest that BPA, and its substitutes BPS and BPF, are potent estrogens acting at nanomolar concentrations via ERalpha. Further studies are required to determine whether this could occur through the activation of ER-extranuclear signaling pathways, as demonstrated in other models [33].
The originality of our study resides in the origin of our cells: widely used 3T3-L1 [28] [34] are murine pre-adipocytes, in addition many studies have demonstrated an effect only in a micromolar range [35] [36]. More importantly, the differentiation stage is very crucial: in a recent study from Boucher et al. [37], it was nicely demonstrated that 10 nM BPS induces a significant increase of 1.6-fold in PPAR-gamma expression at day 4 in primary human pre-adipocytes, but this over-expression was no longer detectable at day 10. In addition, in this study from Boucher et al., other adipogenic markers such as FABP4, LPL, PLIN were significantly increased only by the use of a 25 μM concentration of BPS.
It is worth mentioning that some long intergenic non-coding RNAs and some small nucleolar RNAs were differentially expressed after the exposition to the bisphenols during differentiation. To our knowledge, there is few data reporting the induction of long intergenic noncoding RNAs or small nucleolar RNAs by any bisphenol. The long non-coding RNA HOTAIR was induced in breast cancer cells by BPA [38]. Fifteen small nucleolar RNAs with C/D motif (SNORD) were shown to be reduced in cultures of primary prostate epithelial cells treated with BPA [39]. SNORD are noncoding, small nucleolar RNAs known to regulate ribosomal RNA assembly and function, suggesting that bisphenol A and its substitutes may also act through epigenetic modifications, as already reported [40].
Interestingly, two references were frequently found along comparison of our data with the Comparative ToxicoGenomics Database. The first one is from Ali et al. [41] and describes a microarray experiment on a rat seminiferous tubule culture model exposed for several days with low dose bisphenol A (1 nM and 10 nM) and intersects 31 of our 40 probes. The culture conditions are very similar to ours and it is striking that in a different model undergoing DNA reprogramming such as meiosis, the same genes can be affected by BPA. The second reference was found for 8 of our 40 probes, and describes the ovarian transcriptome of two fish species exposed to BPA [42].
Among the genes identified by this study as bisphenol targets, inhibin beta B subunit (INHBB) has drawn our attention because it has already been described as a target of estrogens [43]. The beta B subunit forms a homodimer activin B, and in combination with the beta A subunit forms a heterodimer activin AB, both of which stimulate FSH secretion [44]. INHBB produced in adipocytes may play a role in the metabolic syndrome, since it was demonstrated that its expression is high in human adipocytes, reduced by weight loss, and correlates with factors implicated in metabolic disease [45]. In addition, INHBB was shown to inhibit lipolysis in adipocytes [46] and decrease thermogenesis in primary cultures of mouse white adipocytes [47]. Considering that INHBB may be used as a marker of spermatogenesis function and male infertility [48], our results strengthen the need of regulations that guarantee a high level of protection of human health.
In conclusion, our results suggest that human primary adipocytes undergoing chronic exposure to BPA or its substitutes BPS and BPF, during differentiation, present deleterious effects on their transcriptome even at a "low" dose (e.g. frequently encountered in a general population). The upstream regulators of the deregulated genes belong to hormonal chemicals and the most highlighted pathways are related to oncogenesis. We noticed an interaction with different families of non-coding RNAs (miRNA, long non-coding and small nucleolar RNAs) suggesting an effect of the three bisphenols on the regulation of gene transcription and also with various aspects of post-transcriptional mRNA processing. Altogether these data suggest that more caution should be taken for the use of BPA and its so-called "substitutes" as previously unsuspected impaired cellular event could be initiated even at a low dose.  Variability between patients (mRNA/lncRNA). A) Venn diagram of differentially expressed probes for each patient before and after ten days of differentiation, a probe is considered as differentially expressed if log 2 FC>2. B) Individual representation of the first two axes of the principal component analysis on all probes for each patient before and after ten days of differentiation (colored according to patient). C) Individual representation of the first two axes of the principal component analysis on all probes for each patient before and after ten days of differentiation (colored according to differentiation stage). (TIF) S4 Fig. Variability between patients (miRNA). A) Venn diagram of differentially expressed probes for each patient before and after ten days of differentiation, a probe is considered as differentially expressed if log 2 FC>2. B) Venn diagram of differentially expressed probes for each patient before and after ten days of differentiation, a probe is considered as differentially expressed if log 2 FC>1. C) Individual representation of the first two axes of the principal component analysis on all probes for each patient before and after ten days of differentiation (colored according to patient). D) Individual representation of the first two axes of the principal component analysis on all probes for each patient before and after ten days of differentiation (colored according to differentiation stage). (TIF) S1 Table. Probes with up-regulated or down-regulated expression for mRNA/lncRNA (Spreadsheet A) or miRNA (Spreadsheet B) in differentiated adipocytes. Fold changes are log 2 transformed fold changes (log 2 FC) of the expression in differentiated adipocytes compared to undifferentiated pre-adipocytes (positive when the probe is over-expressed; negative when the probe is under-expressed). P-values are corrected for multiple testing using Benjamini and Hochberg's method (false discovery rate<5%). Only probes with log 2 FC>2 or log 2 FC<-2 for mRNA (A) and log 2 FC>1 or log 2 FC<-1 for miRNA (B) are listed. Gene names (A) or miRNA symbols (B) have been determined with Ingenuity Pathway. D indicates the presence of a duplication: two different probes are tagging the same gene. Long non-coding RNA data are displayed with a blue filing. ID corresponds to the probe id. (XLSX) S2 Table. List of RNA (mRNA and lncRNA) and miRNA probes differentially expressed for the three bisphenols (A, F and S). Gene names (A) or miRNA symbols (B) have been determined with IPA. D indicates the presence of a duplication: two different probes are tagging the same gene. Fold changes (FC) are actually log 2 transformed fold changes (log 2 FC), therefore if a FC is positive the probe is over-expressed in the tested condition compared to the control, whereas it is under-expressed when a FC is negative. In addition, when a FC is equal to +1 or -1, the probe is twice as over or under-expressed. P-values are corrected for multiple testing using Benjamini and Hochberg's method to control the false discovery rate with a 5% threshold. Long intergenic non-coding RNAs data are displayed with a blue filing and small nucleolar RNAs data are displayed with a pink filing. ID corresponds to the probe id. (XLSX) S3 Table. List of differentially expressed probes deregulated for the three bisphenols A, F and S at 10 μM. Fold changes (FC) are actually log 2 transformed fold changes, therefore if a FC is the probe is over-expressed in the tested condition compared to the control, whereas it is underexpressed when a FC is negative. Du column: D indicates the presence of a duplication: two different probes are tagging the same gene. Long intergenic non-coding RNAs data are displayed with a blue filling and small nucleolar RNAs data are displayed with a pink filling. Probes differentially expressed for both concentrations (10 nM and 10 μM) are displayed with a yellow filling. (XLSX) S4 Table. Differential analysis of housekeeping genes for the three bisphenols (A, F and S) at both concentrations (10 nM and 10 μM). The significant log 2 fold changes along with the p-values are displayed. NS stands for non-significant. The highlighted housekeeping genes (pink background) show no differential expression in any condition compared to the control. (XLSX) S5 Table. Taqman probes (mRNA and mirRNA) used in the qPCR. (XLSX) S6 Table. Comparison of the microarray analysis with the qPCR analysis. The ten genes that were differentially expressed in the microarray analysis for the three bisphenols (A, F and S) at both concentrations (10 nM and 10 μM) were analyzed. The log 2 fold changes (positive when the probe is over-expressed; negative when the probe is under-expressed), along with the p-values are displayed for each technology. (XLSX)