Behavioral assessment of adult offspring after prenatal exposure to BPA

Timed-pregnant females were randomly assigned to the control (Ctrl) and BPA groups. Corn oil or BPA with a dosage of 40 µg/kg bw/day (body weight per day) was administered via oral gavage from gestational day (GD) 0.5 to GD 13.5, respectively (Fig. S1A). We collected 12–14 litters from each group and found no significant difference in sex ratio between the Ctrl and BPA groups (Fig. S1B). We then evaluated behavioral performance of adult offspring in multiple paradigms, including open field, elevated plus maze, tail suspension, forced swim. Only mild abnormality, if any, was observed in the BPA compared to the Ctrl group at this low dosage of exposure (Fig. S1C-D). Given the widely reported sex-effects in BPA studies (Rebuli and Patisaul 2016), we analyzed data from both sexes separately. In the male offspring, we observed a significant increase in total distance in open field test (n = 26–29/group), but no significant change was detected in other paradigms (Fig. S1C). No significant behavioral changes were observed in females except increased time in the center from open filed test (n = 15–30/group) (Fig. S1D). We then collected brains at the age of 4–5-month-old and performed transcriptomic (RNA-seq) and epigenomic profiling for DNA methylation (Reduced Representation Bisulfite Sequencing, RRBS), H3K27me3 and H3K9me3 (ChIP-seq) (Fig. S1A). For the molecular study, we selected animals from different litters to minimize any potential litter effect.

Downregulation of genes in energy metabolic pathway in adult cortex after prenatal BPA exposure

RNA-seq was performed on the prefrontal cortex of adult offspring from both sexes (Fig. 1). Despite the mild behavioral changes observed in male offspring (Fig. S1C), significant transcriptional changes were observed in response to prenatal BPA exposure, with 1,182 significant differentially expressed genes (DEGs) identified between the BPA and Ctrl groups (P < 0.05, n = 3/group) (Fig. 1A, left, Table S1). However, the impact of prenatal BPA exposure on females was considerably less pronounced, with only 145 genes showing significant changes (P < 0.05, n = 3/group) (Fig. 1A, right, Table S2). Furthermore, there was minimal overlap in DEGs between males and females, with only 9 upregulated and 17 downregulated genes common to both (Fig. 1B). These results, together with findings from the literatures (Farabollini et al. 2002; Harley et al. 2013), suggest that males are more vulnerable to prenatal BPA exposure. This also corroborates reported sex differences in response to adverse environmental events during brain development (Pérez-Cerezales et al. 2018).

Fig. 1

Impaired mitochondrial function in adult cortex after prenatal low-dose BPA exposure. (A) Volcano plot of DEGs in adult prefrontal cortex of male offspring and female offspring. Red, up-regulated genes; blue, down-regulated genes. BPA vs. Ctrl, P < 0.05, n = 3/group. (B) Overlapped down-regulated genes (Left) and up-regulated genes (Right) between male and female. (C) ShinyGO analysis (“Biological Process”) of up- (Top) and down-regulated (Bottom) genes (P < 0.05) of male offspring. (D) qRT-PCR validation in male cortex, n = 6/group. Mean ± SEM. Two-way ANOVA. BPA effect, *P < 0.05 (E) ATP levels in the cortex of adult male offspring measured by ATP Detection Assay Kit. n = 6/group. Mean ± SEM. *P < 0.05, Unpaired t test, two-tailed. (F) Bar plots show the gene set enrichment score (-log10Pvalue) of “Down_genes” and “Up_genes” in each cell types. The eight cell-types specific marker genes were identified using the published scRNA-seq (GSE124952). Notably, only “Down_genes” was significantly enriched in neurons. The dashed line represents the significance cutoff threshold at P = 0.05. (G) Transmission electron microscopy (TEM) image of prefrontal cortex in adult male offspring. (Left) Neurons are depicted by red dashed line. Red arrows indicate mitochondria. Scale bar, 2 μm. (Right) Statistical analysis of the number of mitochondria in individual neuron. N = 9/group. Mean ± SEM. P = 0.054, unpaired t test, one-tailed

We focused on males for subsequent analysis. Gene Ontology (GO) analysis of the male DEGs revealed that upregulated genes were enriched in the pathways associated with neuronal functions, including neurogenesis, neuron differentiation and chemical synaptic transmission (Fig. 1C, top). Meanwhile, downregulated genes were primarily enriched in the oxidative phosphorylation (OXPHOS) pathway and other energy metabolic pathways (Fig. 1C, bottom). This finding was also confirmed by Gene Set Enrichment Analysis (GSEA), which assesses the general trend of all detected genes without applying a cutoff threshold for DEGs. The OXPHOS pathway was at the top of the list of significantly downregulated pathways associated with the BPA group (FDR < 0.05) (Fig. S2). Furthermore, we performed qRT-PCR and validated the decreased expression of genes in mitochondrial Complex I in the BPA group, compared to the Ctrl group, using different batch of samples (P < 0.05, n = 6/group) (Fig. 1D).

In addition, we performed an integrated analysis of our bulk RNA-seq, referencing the published scRNA-seq data (Bhattacherjee et al. 2019). When visualizing the average expression of DEGs on the t-SNE plots of the published scRNA-seq, both upregulated (Fig. S3B, D) and downregulated (Fig. S3C, E) genes exhibited relatively high expression in excitatory neurons. However, no significant enrichment was calculated for the upregulated DEGs, while the downregulated DEGs were significantly enriched in both excitatory and inhibitory neurons (Fig. 1F). We performed additional analysis using a different set of published scRNA-seq data (GSE211099). Despite some variations, we consistently observed enrichment of downregulated genes in excitatory neurons (Fig. S4). Given the robust nature of energy metabolism in neurons, and the observed downregulation of genes associated to energy metabolic pathway (Fig. 1C-D, Fig. S2), we speculated that prenatal BPA exposure may potentially influence the mitochondrial function of neurons in adult offspring. To explore this further, we measured ATP levels and observed a significant decrease in the cortex of the BPA compared to the Ctrl group (P < 0.05, n = 6/group) (Fig. 1E). Additionally, we examined mitochondrial morphology under electron microscopy. We observed a trend (P = 0.054, n = 9/group) towards fewer mitochondria in neurons in the BPA group (Fig. 1G). Collectively, these findings suggested potentially compromised energy metabolism in the adult brain following prenatal exposure to BPA.

Differential regulation patterns of DEGs via chromatin interactions in adult cortex after prenatal BPA exposure

Transcriptomic analysis reveals that prenatal BPA exposure leads to alterations in numerous functional pathways in the adult cortex (Fig. 1). These changes encompass the collective expression of a vast number of genes, suggesting a generalized regulatory mechanism as opposed to one restricted to individual genes. Previous research has shown that chromatin interactions, within the scope of spatial chromatin organization, are essential drivers of the coordinated regulation of multiple genes (Rajarajan et al. 2018). Therefore, we utilized published Hi-C data to gain deep insight into the transcriptional regulation patterns of DEGs in our current study (Jiang et al. 2017). We first reconstruct the chromatin contact map, retrieving significant contacts at a resolution of 20 Kb. We then selected the significant contacts associated with all genes detected in our bulk RNA-seq and computed the cumulative interaction score for each gene. Intriguingly, when we arranged the genes according to their interaction scores, the upregulated genes were significantly overrepresented among the top 5% (Fig. 2A, B left). In contrast, downregulated genes were underrepresented in this population (Fig. 2A, B right). Violin plots depict the distribution of interaction scores for upregulated (Up_genes), downregulated (Down_genes), and unchanged (Non_genes) genes. Consistently, the interaction scores of Up_genes were significantly higher compared to those of Non_genes, while Down-genes exhibited the opposite pattern (Fig. 2C). The number of significant contacts associated with Up_genes was larger compared to Down_genes (Fig. 2D), indicating more interactions with Up_genes. O/E value reflets the intensity for each contact, and the contact intensity for Up_genes was higher compared to Down_genes (Fig. 2E), indicating stronger interactions for Up_genes. Moreover, empirical cumulative distribution function (ECDF) plot demonstrated the distribution of interaction distance, and revealed higher proportion of Up_genes contacts with long distances as compared to Down_genes (Fig. 2F). Collectively, these data indicated that there were more and stronger chromatin interactions with longer distances associated with Up_genes compared to Down_genes. Notably, there was a much higher proportion of interactions associated with promoters for Down_genes (80%), compared with 20% for Up_genes (Fig. 2G). We also performed additional analysis with a different set of published Hi-C data (Chandrasekaran et al. 2021) and obtained consistent results (Fig. S5A-G).

Fig. 2

Distinct 3D epigenomic signature on the upregulated and downregulated genes. Published Hi-C data (GSE99363) were used to obtain gene-associated contacts. (A) Gene-associated contacts were extracted out and the cumulative interaction scores for each gene were calculated. Genes were ranked according to the cumulative interaction scores. The dashed line separates the 5% top genes with highest cumulative interaction score with other genes. (B) Permutation plot of number of overlapped genes between the randomly sampled background genes and the 5% top ranked genes. Yellow and blue line indicates the overlap observed between the up- (Left) or down-regulated (Right) genes and the 5% top ranked genes. (C) Violin plot shows distribution of cumulative interaction score of each gene. *P < 0.05, Dunn's test. (D) Bar plot shows the number of significant contacts associated with up- and down-regulated genes. Gene-associated significant contacts were defined as either anchor of the contacts overlapping with genes. (E) Violin plot shows the distribution of O/E value of each contact. (F) ECDF plot of interaction distance shows more proportion of long-range interactions for Up_genes. (G) Bar plot shows the proportion of promoter associated with significant contacts. (H) Chromatin states were categorized by ChromHMM using our ATAC-seq, H3K27me3, H3K9me3 ChIP-seq data and published H3K27ac (GSE99363), CTCF (GSE99363), H3K4me1 (GSE90020), H3K4me3 (GSE90020) ChIP-seq datasets. The enrichment score for the categories is shown by color. The red box highlights the most significant enrichment of categories. (I) Representative map tracks show chromatin contacts and enhancer marks (H3K27ac, H3K4me1) associated with Up_genes (Left) and Down_genes (Right)

We further applied ChromHMM analysis to assess the chromatin states for both genes and their corresponding anchors of chromatin interactions associated with upregulated and downregulated genes. Using our ChIP-seq (H3K27me3, H3K9me3) and ATAC-seq data, together with published ChIP-seq (H3K27ac, H3K4me1, H3K4me3 and CTCF) (Jiang et al. 2017; Walsh et al. 2017) datasets, we categorized chromatin into six states for the adult mouse cortex (Fig. 2H left). We then calculated the enrichment of upregulated and downregulated genes, along with their corresponding anchors, in each of the chromatin state respectively. Our results revealed that upregulated genes and their corresponding anchors were primarily enriched in the enhancer regions (2, Enh; 3, EnhLo), characterized by H3K27ac, H3K4me1, and ATAC-seq signals (Fig. 2H right). Meanwhile, downregulated genes and their corresponding anchors were predominantly enriched in active gene promoters (1, Tss), featured with H3K4me3 and other open chromatin signals (Fig. 2H right). Consistent results were revealed (Fig. S5H) when we performed additional ChromHMM analysis using different sets of published ChIP-seq (ENCODE Project Consortium 2012; Hagelkruys et al. 2022; Mo et al. 2015) and ATAC-seq data (Stroud et al. 2020).

In summary, we discovered that the patterns of chromatin spatial interaction associated with upregulated genes markedly differed from those associated with downregulated genes. The upregulated genes, which exhibit more and stronger chromatin interactions, are primarily under the remote control of enhancers (Fig. 2I left, Fig. S5I left). Conversely, downregulated genes appear to be regulated by chromatin interactions among relative adjacent genes (Fig. 2I right, Fig. S5I right). These distinct regulatory patterns provide insight into the unique gene regulation signatures in adult cortex following prenatal exposure to BPA.

Alterations of DNA methylation in adult cortex after prenatal BPA exposure

DNA methylation is one of the most extensively studied repressive epigenetic marks. It has been reported that BPA can influence DNA methylation thus affect the expression of neuronal genes (Alavian-Ghavanini et al. 2018; Kundakovic et al. 2015). To study whether DNA methylation contributed to gene dysregulation, we performed RRBS to profile genome-wide alterations of DNA methylation in adult cortex from both the BPA and Ctrl groups. We identified a total of 10,136 differentially methylated sites (DMSs), comprising 5,329 hypermethylated and 4,807 hypomethylated sites (q < 0.05, difference > 10%, n = 3/group) (Fig. 3A, Table S3). The alterations on the hypomethylation sites were more pronounced (Fig. 3A). Notably, genome annotation of DMSs revealed that 61% hypomethylated and 64% hypermethylated sites were gene-associated, located in the promoter, intron, and exon regions (Fig. 3B). GO analysis showed that hypomethylated genes were enriched in pathways of glutamatergic synapse and cell signaling (Fig. 3C), while hypermethylated genes were enriched in the pathways of neurogenesis, neuron differentiation and cell projection (Fig. 3D).

Fig. 3

Alterations of DNA methylation in adult cortex after prenatal BPA exposure. (A) Volcano plot of DMSs in cortex of adult male offspring. Red, hypermethylated sites; blue, hypomethylated sites. BPA vs. Ctrl, q-value < 0.05 and methylation difference > 10%, n = 3/group. (B) Genome annotation pie-plot of the hypo- (Top) and hyper- (Bottom) methylated sites. ShinyGO analysis (“Biological Process”) of hypo- (C) and hyper- (D) methylated genes. (E) Overlap between down-regulated genes and hypermethylated genes (Left). Overlap between up-regulated genes and hypomethylated genes (Right). (F) ShinyGO analysis (“Biological Process”) of the overlapping hypermethylated and downregulated genes. (G) Permutation test shows mean distance (Left) and number of overlaps (Right) between DNA hypermethylated sites and anchors interacting with the promoters of downregulated genes

Despite the modest overlap between genes associated with DMSs and DEGs from RNA-seq (Fig. 3E), GO analysis of the 61 overlapping hypermethylated and downregulated genes showed significantly enrichment in pathways associated with mitochondrial functions, notably the energy respiratory chain (Fig. 3F). Meanwhile, 73 up-regulated genes were overlapped with hypomethylated genes, but no enriched pathway were identified in this group. Further analysis revealed that, compared to chromatin interaction anchors for all gene promoters, a significantly higher number of overlaps and a shorter mean distance were observed between DNA hypermethylation sites and regions (anchors) interacted with the promoters of downregulated genes (Fig. 3G). These findings suggest that prenatal BPA exposure leads to alterations in DNA methylation within the cerebral cortex of adult offspring, and DNA hypermethylation might partially account for the observed gene downregulation.

Alterations of repressive histone modifications in adult cortex after prenatal BPA exposure

In addition to DNA methylation, repressive histone modifications such as H3K27me3 and H3K9me3 also mediate long-term effects of BPA (Fatma Karaman et al. 2019; Senyildiz et al. 2017). We, therefore, performed H3K27me3 ChIP-seq on adult cortex from both BPA and Ctrl groups. On average, we detected 10,581 H3K27me3 peaks per sample. however, only a small percentage (1.6%) of these peaks were altered in the BPA compared to the Ctrl group. Among them, 125 peaks showed increased H3K27me3 occupancy in the BPA group (P < 0.05, n = 4/group) (Fig. 4A, Table S4), with 86% of these peaks located on gene promoters (Fig. 4B, left). Only 41 peaks were downregulated in the BPA group (P < 0.05, n = 4/group), and 32% on the gene promoters (Fig. 4B, right). Subsequently, we investigated the association between differential H3K27me3 peaks and DEGs. Permutation tests revealed that the location of upregulated H3K27me3 peaks were significantly closer to the TSS of downregulated genes (Fig. 4C). Conversely, downregulated H3K27me3 peaks were significantly further away from the TSS of upregulated genes (Fig. 4D). Moreover, H3K27me3 signals peaked at the TSS of downregulated genes, and these signals were more pronounced in the BPA compared to the Ctrl group (Fig. 4E, left), especially for the OXPHOS genes (Fig. 4E, right). These findings suggest that H3K27me3 might regulate gene expression by directly repressing the promoter activity of downregulated genes in adult cortex after BPA exposure.

Fig. 4

Alterations of repressive histone modifications in adult cortex after prenatal BPA exposure. (A) Volcano plot of differential H3K27me3 modification sites. Red, enhanced modification; blue, decreased modification. BPA vs. Ctrl, P < 0.05, n = 4/group. (B) Genome annotation pie-plot of the increased (Left) and decreased (Right) H3K27me3 modification sites. (C) Permutation test shows mean distance between increased H3K27me3 peaks and the TSS of downregulated genes. (D) Permutation test shows mean distance between decreased H3K27me3 peaks and the TSS of upregulated genes. (E) Heatmap of H3K27me3 signals on ± 3 Kb of Down_genes (Left) and OXPHOS gene promoters (Right). In each group, the signal of three biological replicates were merged. The gradient color bar indicates the signal strength. (F) Neurons of adult male cortex were enriched by FANS sorting for H3K9me3 ChIP-seq. Volcano plot of differential H3K9me3 modification sites were shown. Red, enhanced modification; blue, decreased modification. BPA vs. Ctrl, for differential H3K9me3 sites, cut-off was set as: Padj < 0.05, basemean > 50, n = 4 Ctrl/5 BPA. (G) Genome annotation pie-plot of the increased (Left) and decreased (Right) H3K9me3 sites. (H) Permutation test shows mean distance (Left) and number of overlaps (Right) between decreased H3K9me3 loci and TSS of upregulated genes. (I) Permutation test shows mean distance (Left) and number of overlaps (Right) between decreased H3K9me3 loci and anchors interacting with the promoters of upregulated genes. (J) H3K9me3 signals at the overlapping regions between H3K9me3 decreased loci and anchors interacting with the promoters of upregulated genes

To improve detection sensitivity and minimize interference from various cell types in the brain, we applied fluorescence activated cell nuclei sorting (FANS) to isolate NeuN + nuclei from the adult cortex for H3K9me3 ChIP-seq. Differential analysis, performed using diffReps with a window size of 1 Kb, identified significant alterations in H3K9me3 occupancy in response to prenatal BPA exposure (Padj < 0.05, n = 4 Ctrl/5 BPA) (Fig. 4F, Table S5). Unlike H3K27me3, genome annotation showed that the majority of differential H3K9me3 sites were located in gene deserts and other intergenic regions. Only about 25% of these sites were gene-related for both upregulated and downregulated sites (Fig. 4G), implying H3K9me3 might primarily function on distal regulatory elements. Indeed, when we examined the association between downregulated H3K9me3 loci and upregulated genes, we only detected mild enrichment in terms of distance or overlapping number compared to the background (Fig. 4H). However, permutation tests revealed a significantly higher number and a shorter mean distance between downregulated H3K9me3 loci and regions (anchors) interacting with upregulated genes (Fig. 4I). Moreover, H3K9me3 signals were notably lower in the BPA group compared to the Ctrl group at these overlapping regions (Fig. 4J). These findings suggest that the loss of H3K9me3 occupancy on distal regulatory regions may contribute to the upregulation of gene expression in adult cortex after prenatal BPA exposure.