Repeat number variation in adolescent males

The experimental design is provided in the Methods section as well as in Additional file 1: Fig. S1. To assess the effects of BPA on genetic stability, we first analysed repeat numbers at variable regions of the human genome: specifically, we looked at the relative repeat numbers of three satellite DNA regions (pericentromeric SAT2, centromeric and pericentromeric SATA), at LINE1 retroelements, which encode for the viral protein ORF2, and at Alu elements, as Alu is the most abundant transposable element. Two groups of LINEs were analysed. The first group (ORF2-1 region) matches 4560 site in the human genome, while the second group (ORF2-2) matches 2918 regions in the genome [27]. We also assessed the telomere repeats. Genomic DNA was extracted from human blood samples, and each DNA sample was adjusted to 0.2 ng/µL. Given the worldwide distribution of BPA, virtually every human has been exposed to some extent. As a consequence, it is not possible to measure the genetic and epigenetic “baseline” in an unexposed control group. Therefore, our analysis was performed by comparing groups with low and high BPA levels. The samples from each group are presented in Additional file 1: Table S1. Repeat numbers were normalized to the RPLP0 gene. Normalized DNA copy number variations were plotted and compared (Fig. 1). We observed a significant increase in the repeat number of centromeric SATA (fold change (FC) = 1.6, p = 0.05) and a tendency to increase for Alu repeat number (FC = 1.2, p = 0.06) in the high BPA group (Fig. 1A) (n = 15, BPA range 4.85–42.64 µg/g) compared to the low BPA group (n = 16, BPA range 0.66–4.67 µg/g, (Additional file 1: Table S1)). We also observed a tendency to an increased telomere length (TL) in the high BPA group (Fig. 1A).

Fig. 1

DNA analysis in the blood of adolescent males: A Repeat number variations identified in the low BPA and high BPA groups. B DNA methylation analysis of imprinted genes and repeat regions. C DNA methylation analysis of genes regulated by ESR1. *p < 0.05, Mann‒Whitney test

Analysis of DNA methylation in adolescent males

Next, differences in DNA methylation in the CpG regions of various promoters were analysed in 12 samples from the low BPA group and 9 from the high BPA group. We chose the available samples with the most extreme BPA values for the analysis to reveal the most profound effects of BPA, in low group, n = 12, BPA range: 0.66–3.94 µg/g, on high group, n = 9, BPA range: 9.44–42.64 µg/g, (Additional file 1: Table S2).

The MEDIP technique was used, and several promoter gene targets were selected for the analysis. To design the primers, the CpG island map of the reference hg38 genome was downloaded from UCSC using the Table Browser tool, and the Primer-blast tool from PubMed was used. Several control regions were chosen for the MEDIP assay, and a set of imprinted genes (H19, KCNQ1OT1, IGF2) was also analysed, as the imprinting of these genes requires DNA methylation. Satellite DNAs (SAT2) and retroelements such as SINEB1 and the 5′UTR of LINE1 were also chosen due to the high levels of DNA methylation in these regions. The ratio of precipitated methylated DNA to input DNA was calculated. Most targets showed strong MEDIP to input ratios, indicating that these targets are highly methylated. The data are shown in Fig. 1B, C. Adolescents in the high BPA group showed lower DNA methylation than those in the low BPA group for all sequences tested, and these differences were statistically significant—approximately 0.5-fold—for the H19, KCNQ10T1, IGF2 and LINE1 promoters.

On the other hand, since BPA exhibits oestrogenic activity, we wondered whether DNA methylation would be altered at the promoters of genes that have indications of regulation by ESR1 and have ESR1 binding sites in the vicinity. For this purpose, genes regulated by ESR1 were chosen, and we used the Chea 2022 database of ESR1 binding ChIP-seq data from Enrichr [28]. Thus, we chose to analyse DNA methylation at the promoters of the ARID2, ESRRA, EGFR, PGR, ESR2, SUV39H1, VEGFA and TERT genes. The results showed a significant decrease in the DNA methylation of ESRRA (FC = 0.6, p = 0.0096) and EGFR (FC = 0.5, p = 0.038) and a decreasing trend in ARID2 (FC = 0.7, p = 0.15) and SUV39H1 (FC = 0.63, p = 0.13) in the high BPA group compared to the low BPA group (Fig. 1C). We also analysed DNA methylation at the telomerase (TERT), since it is known that oestrogenic compounds affect telomere length (TL); which was corroborated by our results showing a significant decrease of TERT DNA methylation in the high BPA group (FC = 0.16, p = 0.003), Fig. 1C.

Overall, the above results indicate that high BPA exposure generally correlates with lower DNA methylation. This decrease in DNA methylation was observed at the promoters of genes that regulated by ESR1/ERa, imprinted genes and LINE1 sequences.

Genome-wide analysis of H3K4me3

To investigate the potential new role of BPA in the regulation of the epigenetic landscape at the genome level, genome-wide analysis of H3K4me3 occupancy was performed using 100 µL of whole blood samples from 10 adolescents in the low BPA group, BPA range: 0.66–2.55 µg/g, and 7 adolescents in the high BPA group, BPA range: 10.07–42.64 µg/g (Additional file 1: Table S3). Chromatin immunoprecipitation was performed using anti-H3K4me3 antibody followed by high-throughput sequencing (all details are provided in the Methods section). This histone mark was chosen because it is associated with transcription start sites and its presence correlates well with gene expression. Sequencing reads were analysed, and the differential peaks were detected. Using cut-off values of FC >  = 1.5 and FDR < 0.05, 5707 peaks were found to be differential between the two exposure groups (Additional file 1). The vast majority of differential peaks that we found had increased H3K4me3 occupancy in the high BPA group compared to the low BPA group.

Differential H3K4me3 occupancy regions were then screened for genes with shared biological functions. To this end, differential peaks were assigned to genes using GREAT (details and parameters of these analyses are given in the Methods section). A total of 666 genes related to different “biological process” (Additional file 1: Fig. S2A), “molecular function” (Fig. 2A), or “cellular component” (Additional file 1: Fig. S2B) terms were assigned. Interestingly, the “molecular function” analysis revealed H3K4me3 differential occupancy in gene clusters related to histone protein transporter activity, histone H4 acetyltransferase activity, and telomeric DNA binding (Fig. 2A). For example, the enriched cluster ‘H4 histone acetyltransferase activity’ includes genes such as KAT2A, KAT8, and OG, and the enriched cluster ‘telomeric DNA binding’ contains 22 genes (Additional file 1: Table S4) important for telomere function, such as PURA, RPA1, TEN1, TERT, TERF2IP and ZBTB48. An example of differential peaks for TERF2IP and ZBTB48 is illustrated in Fig. 2B.

Fig. 2

Differential H3K4me3 occupancy in the low and high BPA groups. Functional annotation analysis of differential regions. A Gene ontology annotation with “molecular function” terms was performed by GREAT. Terms were sorted by p value. B H3K4me3 regions in TERF2IP and ZBTB48. Reads were mapped and plotted on a graph by IGV; blue represents the low BPA group, and red represents the high BPA group. C Comparison of genes located in differential peaks and protein‒protein interaction (PPI) analysis; the second column shows the number of genes identified for each transcription factor compared to the total number of known interacting proteins. D Motif found by MEME. Parts of this motif are similar to the binding sites of E the methyl binding protein MBD2 and F the methyl binding protein MECP2. The E-value shows the statistical significance of the motif. G ESR1 binding sites. The E-value shows the statistical significance of the motif presence in differential peaks

The list of genes identified in differential regions assigned by GREAT was then compared with the Enrichr protein‒protein interaction (PPI) database [28]. This database combines the data of protein‒protein interactions previously described in experimental studies.

This analysis revealed that targets of 9 proteins, including ESR1/ERa, DNA repair factors (BRCA1, TP53), transcriptional factors (MYC, FOXP3, POU5F1, NOTCH1) and chromatin regulators (EP300, TRIM28), were all enriched in our dataset.

Our analysis revealed that genes encoding ESR1-interacting proteins were strongly enriched for H3K4me3 in the high BPA group (Fig. 2C). According to the Enrichr PPI database, ESR1 could interact with 871 proteins, of which 510 were encoded by genes located within differential H3K4me3 peaks (adj. p value = 6.95E−53). Among these genes were the AHR, ESRRA, NCOR1, NCOR2, RARA and RXRA receptors (Additional file 1: Table S5). Furthermore, a PPI database search showed that 236 genes were located in differential peaks corresponding to genes encoding proteins that interact with BRCA1 (BRCA1 could interact with up to 418 proteins). Many BRCA1-interacting proteins encoded by genes located in differential peaks were DNA repair proteins. The AMIGO database includes 519 “DNA repair” genes. Of them, 39 genes were detected in our differential peaks, including ATM, ATR, BLM, FANCA, MDC1, MLH1, MSH2, MSH3 and RAD51 (Additional file 1: Table S6). These genes encode proteins critical for various DNA repair processes, suggesting that exposure to BPA could activate DNA repair pathways.

Overall, whole-genome H3K4me3 occupancy analysis revealed that the methylation at genes involved in protecting genetic stability, such as DNA repair genes and telomeric DNA binding protein-encoding genes, was altered by BPA exposure, suggesting the potential toxicity of BPA in these crucial processes.

Motif analysis in differential H3K4me3 peaks

To address whether regulatory DNA motifs were present within the regions containing the altered peaks, MEME-CHIP was used [29]. We investigated whether there was enrichment for transcription factor-binding sites in the differential peaks. For this, sequences located in differential peaks were extracted in fasta format, and repeated masking was performed to remove repetitive sequences from the analysis. The resulting data were processed in MEME-CHIP to reveal the significantly enriched motifs in the differential peaks. A recurrent 15-mer motif was detected by MEME-ChIP (Fig. 2D). Parts of this motif were similar to the MBD2 binding site (E-value = 1.3−e15) (Fig. 2E) and the MECP2 binding sites (Fig. 2F).

We also looked for the presence of ESR1 binding sites in the differential H3K4me3 peaks. We used the FIMO tool for ESR1 motif scanning. The ESR1 binding motif was downloaded from HOCOMO (https://hocomoco11.autosome.org). FIMO detected 12,839 sites in the human reference genome hg38 (p value less than 1e−05). Next, we scanned the sequences from differential H3K4me3 peaks for the presence of potential ESR1 binding sites. There were 2279 occurrences of ESR1 motifs with a p value less than 1e−05. A comparison of the ESR1 sites found in differential peaks to a total ESR1 sites in the hg38 reference genome showed a significant increase in ESR1 binding sites in the high BPA group (p value = 1e−16, Chi-square test, Fig. 2G) compared to the low BPA group, suggesting a significant impact of BPA exposure on the epigenetic regulation of ESR1 targets.

Comparison of alterations in histone marks in human samples with gene expression changes in the HeLa cell line

To confirm the observed changes and further mechanistically assess the impacts of BPA on human cells, several in vitro assays were performed in HeLa cells exposed to a low dose of BPA (10 nM), a dose that has not revealed general toxicity [30, 31] and coincides with urinary BPA concentrations showed by the general population (around 2 ng/mL [32]).

The important observed effects of BPA exposure were on histone H4 acetylation enzymes, so we chose to analyse the levels of histone H4 acetylation which is abundant at open chromatin. We also looked at H3K9 trimethylation, which is enriched at heterochromatin regions. For this analysis, cells were fixed, permeabilized and immunostained against the selected histone marks. Signals were analysed by fluorescence microscopy. To quantitatively analyse the histone mark distribution, z-stack analysis was performed. The z-stack signals were merged and normalized to the merged DAPI signals. We observed that H4 acetylation appeared everywhere throughout the nucleus (Fig. 3A). Quantitative analysis showed a slight but significant 1.2-fold increase in H4 intensity in BPA-exposed cells (Fig. 3B). H3K9me3 was preferentially localized in the cell periphery (Fig. 3C), and in contrast to H4 acetylation, H3K9me3 global amount appeared to be decreased in BPA-exposed cell samples (Fig. 3D).

Fig. 3

Global histone modification analysis in HeLa cells. A H4 acetylation in control cells (top) and in BPA-exposed cells (bottom); B Quantitative analysis of H4 acetylation; C H3K9me3 in control cells (top) and in BPA-exposed cells (bottom); D Quantitative analysis of H3K9me3. The sum of fluorescence for all slices of z-stacks was calculated for each channel. The histone modification signal was normalized to DAPI and presented as combined corrected fluorescence compared to the control, *p < 0.05, Mann‒Whitney test, n = 4 experiments for each group

Next, we evaluated whether the changes observed in histone H3K4me3 occupancy in human blood samples could lead to gene expression changes in cultured cells. To this end, the expression of genes encoding DNA repair, cell cycle, chromatin remodelling and transcription factors was measured in HeLa cells exposed or not exposed to BPA. The analysis showed 1.7- 2.7–1.5-fold increases in the expression of the ATM, ARID2, TP53BP1 DNA repair genes, respectively, in BPA-exposed compared to nonexposed cells (Fig. 4A). Strikingly, gene expression was found to be increased in almost all of the telomeric-DNA binding factors tested, including 1.5-, 2.1-, 1.7-, 1.3-, 1.7- and twofold increases in the SMG7, TERT, UPF1, UPF2, TEN1 and ZBTB48 genes, respectively (Fig. 4B). In addition, we observed a 1.7-fold increase in STAT2 expression. STAT2 acts as a transcriptional activator in inflammatory processes [33]. We did not observe changes in genes encoding oestrogen receptors, except for ESRRA, which increased in 1.5 times.

Fig. 4

Gene expression analysis by RT‒qPCR in BPA-exposed HeLa cells. A Gene expression analysis of DNA repair and chromatin factors. B Genes encoding transcription, signalling and telomeric DNA binding proteins, *p < 0.05, **p < 0.01, Mann‒Whitney test, n = 6 experiments for each group

Therefore, these data confirmed that BPA exposure of human cell lines leads to alteration of histone marks and changes in gene expression programs. Modifications of the expression of DNA repair, inflammation and telomeric-DNA binding proteins encoding genes correlate with the differences in DNA methylation and H3K4 occupancy previously observed in human blood samples.

H3K4me3 and H3K9me3 occupancies and DNA methylation state at selected promoters in HeLa cells

Next, we asked whether the changes in gene expression were a consequence of changes in epigenetic regulation. We chose to analyse the H3K4me3 marker, as we had analysed it in human samples. We performed chromatin immunoprecipitation against H3K4me3 in exposed and unexposed cells. qPCR analysis was then performed. The qPCR analysis of individual genes showed a global tendency toward increased H3K4me3 occupancy after BPA treatment. The observed increase in H3K4me3 levels was especially significant for the DNA repair genes, including 2.4-, 1.4-, 1.3- and 1.3-fold increases in the promoters of ARID2, ATM, OGG1 and H2AX, respectively (Fig. 5A); and also, for the telomeric-DNA binding factors, including 1.3- and 1.3-fold increases in the promoters of UPF1 and UPF2 genes, respectively. Finally, we determined 1.2-fold increase in H3K4me3 at the promoter of KAT5 (Fig. 5B) and 0.5-fold decrease in H3K4me3 at the promoter of ESR1 (Fig. 5C).

Fig. 5

H3K4me3 and H3K9me3 promoter occupancies analysis in HeLa cells. A H3K4me3 occupancy at the promoters of DNA repair genes. B H3K4me3 occupancy at the promoters of oestrogen signalling genes. C H3K4me3 occupancy at the promoters of the telomere maintenance and chromatin factors encoding genes. D H3K9me3 occupancy at the promoters of DNA repair and transcription factors genes and at repeated sequences. E H3K9me3 occupancy at the promoters of telomere maintenance genes

We also analysed the H3K9me3 mark as a proxy for silencing, as a decrease in the occupancy of this mark would lead to enhanced gene expression [34]. In addition, this mark is particularly vulnerable to BPA exposure, and alterations in H3K9me3 could be inherited; e.g., BPA exposure causes a reduction in the repressive marks H3K9me3 and H3K27me3 in whole worms and germline nuclei in third-generation offspring [35]. To this end, we performed chromatin immunoprecipitation against H3K9me3 in cells exposed or not exposed to BPA. qPCR analysis was then performed for selected targets. In agreement with our immunofluorescence data (Fig. 3C, D), the analysis of individual genes showed a global tendency towards a decrease in H3K9me3 occupancy following BPA treatment. This H3K9me3 decrease was especially significant for the STAT2 gene and the 5’UTR of LINE1 and slightly less significant for the SAT2 and SATA loci and PGR loci (Fig. 5D). We also detected decreases in H3K9me3 at the promoters of genes encoding factors involved in telomere function or maintenance, with SMG7 and XRCC5 being the most significantly affected (Fig. 5E). Thus, our data show that BPA treatment leads to decreased H3K9me3 occupancy at a large number of loci, suggesting that modifications of H3K9me3 occupancy could play a pivotal role in gene expression dysregulation following BPA exposure.

To reveal whether there was a change in DNA methylation, we extracted DNA from exposed and unexposed HeLa cells to BPA and performed MEDIP analysis at targets analysed in human blood. We observed an overall decrease in DNA methylation, although only a few genes showed a statically significant alteration, such as ARID2, UPF2, and a trend toward a decrease in PGR, BRCA1, ESR2, VEGF, LINE1 and IGF2 (Additional file 1: Fig. S3).

Overall, our analysis of H3K4me3 and DNA methylation showed similar changes in BPA-exposed human blood and HeLa cells.

ESR1/ERα binding to its targets

To test whether ESR1/ERα could be involved in the BPA-induced dysregulation of specific genes by directly binding to the promoters of KAT5, RARA and some telomere-associated genes, chromatin immunoprecipitation using anti-ESR1 antibody was performed, along with qPCR analysis of the targets. The design of the ChIP‒qPCR primers is detailed in the Methods section. The analysis showed a 1.7-, 1.8- and 1.9-fold increase in ESR1 binding at the KAT5, KMT2E and STAT2 loci, respectively (Fig. 6A). For the five telomeric DNA-binding genes, an increase in ESR1 binding was identified only for TERF2IP (1.9-fold), suggesting that the alterations in telomeric DNA binding factor expression could be mediated by other factors.

Fig. 6

ESR1 binding analysis. A ESR1 ChIP‒qPCR analysis of ESR1 target genes; differential peak regions and those with ESR1 binding motifs were chosen for the analysis. *p < 0.05, **p < 0.01, Mann‒Whitney test, n = 6 experiments for each group. B ESR1 motif conservation analysis in the promoters of UPF2 and TERF2IP

To further reveal whether at least some gene promoters located in differential peaks have evolutionarily conserved ESR1 binding motifs, we used the ECR browser tool [36]. This tool is designed to highlight candidate functional elements in the genome by comparing sequences from several species and finding conserved blocks. Notably, both UPF2 and TERF2IP gene regulatory regions were found to contain an ESR1 binding motif that was conserved between the mouse and human reference genome sequences (Fig. 6B), suggesting that these could be functional elements in both species.

Analysis of ESR1 and p53/TP53 proteins in cells

Our analysis in both human blood and HeLa cells revealed that some of the ESR1 pathway targets were altered at epigenetic and transcriptional levels, we did not observe any changes in the ESR1 gene (which encodes ESR1/ERα) expression. We hypothesized that the observed changes seen in ESR1-target genes expression could be mediated by the changes in ESR1/ERα protein levels. We therefore performed analysis in HeLa cells with preserved nuclei and z-stack analysis of cells immunolabeled with an anti-ESR1 antibody. ESR1 staining appeared as a dot in both the nucleus and the cytoplasm (Fig. 7A); in exposed cells, the staining was mostly nuclear (Fig. 7A). Quantitative analysis revealed that there was a 3.4-fold increase in nuclear ESR1 staining in exposed cells compared to unexposed cells (Fig. 7B), suggesting that there is activation of ESR1 at the protein level.

Fig. 7

ESR1 and p53 immunofluorescence analysis in HeLa cells. A ESR1 staining in control cells (top) and in BPA-exposed cells (bottom); B Quantitative analysis of ESR1 fluorescence intensity; C p53 localisation in control cells (top) and in BPA-exposed cells (bottom); D Quantitative analysis of p53 foci per nucleus

Thus, our data suggest that BPA exposure probably stabilizes the ERa protein.

Since BPA exposure could lead to genotoxic stress, we decided to check whether cell exposure to 10 nM BPA could induce p53, a key player of the DNA damage response pathway. Using an immunoassay against p53, we observed that p53 also appeared as a focus, localised in both the nucleus and the cytoplasm (Fig. 7C). Quantitative analysis of foci showed an increased amount of p53 in exposed cells (Fig. 7D), mainly in the nucleus, 1.5-fold higher than that observed in unexposed cells (Fig. 7D).

Our immunofluorescence data therefore show that ESR1 and p53 localise more in the nucleus in treated cells than in unexposed cells.