Longitudinal Effects of Developmental Bisphenol A Exposure on Epigenome-Wide DNA Hydroxymethylation at Imprinted Loci in Mouse Blood.

BACKGROUND: Epigenetic machinery plays an important role in genomic imprinting, a developmental process that establishes parent-of-origin-specific monoallelic gene expression. Although a number of studies have investigated the role of 5-methylcytosine in imprinting control, the contribution of 5-hydroxymethylcytosine (5-hmC) to this epigenetic phenomenon remains unclear.
OBJECTIVES: Using matched mouse blood samples (from mice at 2, 4, and 10 months of age), our objective was to examine the effects of perinatal bisphenol A (BPA) exposure (50 μg/kg diet) on longitudinal 5-hmC patterns at imprinted regions. We also aimed to test the hypothesis that 5-hmC would show defined patterns at imprinted genes that persist across the life course.
METHODS: Genome-wide 5-hmC levels were measured using hydroxymethylated DNA immunoprecipitation sequencing (HMeDIP-seq). Modeling of differential hydroxymethylation by BPA exposure was performed using a pipeline of bioinformatics tools, including the csaw R package.
RESULTS: Based on BPA exposure, we identified 5,950 differentially hydroxymethylated regions (DHMRs), including 12 DHMRs that were annotated to murine imprinted genes—Gnas, Grb10, Plagl1, Klf14, Pde10a, Snrpn, Airn, Cmah, Ppp1r9a, Kcnq1, Phactr2, and Pde4d. When visualized, these imprinted gene DHMRs showed clear, consistent patterns of differential 5-hmC by developmental BPA exposure that persisted throughout adulthood.
CONCLUSIONS: These data show long-term establishment of 5-hmC marks at imprinted loci during development. Further, the effect of perinatal BPA exposure on 5-hmC at specific imprinted loci indicates that developmental exposure to environmental toxicants may alter long-term imprinted gene regulation via an epigenetic mechanism. https://doi.org/10.1289/EHP3441.

Introduction

DNA methylation is an epigenetic mark that typically occurs at cytosines (5-methylcytosine; 5-mC) in cytosine-phosphate-guanine (CpG) dinucleotides. Research has shown that 5-mC plays a role in X chromosome inactivation (Cotton et al. 2015), regulation of gene transcription (Medvedeva et al. 2014), and genomic imprinting (Bartolomei and Ferguson-Smith 2011). In addition to 5-mC, recent studies have shown that the oxidized form of 5-mC—5-hydroxymethylcytosine (5-hmC)—is a stable epigenetic mark present in a variety of mammalian tissues (Globisch et al. 2010; Hahn et al. 2014; Wu et al. 2011). Active processing of 5-mC to 5-hmC occurs via a Ten-Eleven Translocation (TET) methylcytosine dioxygenase–mediated oxidative pathway (Shen et al. 2014), and previous studies have shown that exposure-induced oxidative stress can alter both TET enzyme expression (Coulter et al. 2013) and global 5-hmC levels (Delatte et al. 2015). In addition to these characteristics, 5-hmC has a complex role as both a positive and negative regulator of transcription (Hahn et al. 2014; Wu et al. 2011), suggesting that it may represent an important secondary genomic regulator in the methylome.

During early embryonic development, DNA methylation undergoes a distinct wave of demethylation and de novo methylation (Reik et al. 2001; Smallwood and Kelsey 2012), processes that assist in the regulation of stem cell proliferation and differentiation (Messerschmidt et al. 2014). After differentiation, there is a second phase of DNA methylation reprograming that occurs to establish a baseline methylome in primordial germ cells prior to birth (Smallwood and Kelsey 2012). Given these multiple waves of reprograming during development, DNA methylation has been proposed as a mechanism driving the developmental origins of health and disease (DOHaD) hypothesis. The DOHaD posits that exposure to environmental factors (e.g., diet, chemicals) during critical periods of development influences developmental plasticity, thereby altering disease susceptibility later in life (Bateson et al. 2004; Heindel et al. 2015; Jirtle and Skinner 2007). Of particular interest in the DOHaD field is genomic imprinting, an epigenetic phenomenon by which parent-of-origin–specific monoallelic gene expression is established during development (Bartolomei and Ferguson-Smith 2011; Das et al. 2009). Genomic imprinting is critical for proper placental maturation, embryonic growth, and development, and it also plays important roles in postnatal development—especially in parent–offspring interactions (e.g., milk suckling, maternal care) (Plasschaert and Bartolomei 2014). A vast amount of research on genomic imprinting has shown that differential DNA methylation at imprinting control regions (ICRs) plays an important role in establishing genomic imprinting during development, and that these regions are not demethylated during postfertilization reprograming (Arnaud 2010; Kelsey and Feil 2013; Kim et al. 2017; Pidsley et al. 2012; Smallwood and Kelsey 2012). Given that imprinted gene expression is controlled via developmentally programed epigenetic marks, imprinted genes have been investigated as potential targets of developmental environmental exposures in both human and rodent studies (Haycock and Ramsay 2009; Heijmans et al. 2008; Nye et al. 2015; Susiarjo et al. 2013).

However, virtually all of the existing studies on genomic imprinting have relied upon sodium bisulfite treatment for the generation of DNA methylation levels from biological samples. Although bisulfite treatment is a useful quantitative tool for measuring the methylome, it also has an inherent weakness: It does not distinguish between 5-mC and 5-hmC. This is because both 5-mC and 5-hmC are resistant to sodium bisulfite-induced cytosine deamination (Huang et al. 2010). This means that bisulfite treatment–based measurements of the “methylome” actually represent combined levels. In recent years, methods to specifically measure 5-hmC, including hydroxymethylated DNA immunoprecipitation sequencing (HMeDIP-seq) (Tan et al. 2013), have been developed. Despite these recent advances, however, little work has been done to identify the role of 5-hmC in genomic imprinting. One recent study showed enrichment of 5-hmC in imprinted regions in the human brain and placenta (Hernandez Mora et al. 2018), but no study has yet investigated the contribution of 5-hmC to imprinting control in animal models. Furthermore, existing animal studies have not investigated effects of the environment on 5-hmC levels at imprinted genes.

Developmental exposure to a number of environmental factors, including endocrine disrupting chemicals (EDCs), has been linked to changes in DNA methylation. Bisphenol A (BPA) is a commercial EDC that is found in a variety of consumer products (e.g., receipt paper, metal can linings) and has near ubiquitous human exposure across the world (Calafat et al. 2008). BPA can activate growth-related transcription factors, bind nuclear receptors involved in cell growth and maturation, and also alter DNA methylation levels across the epigenome (Krüger et al. 2008; Manikkam et al. 2013; Singh and Li 2012; Sui et al. 2012; Watson et al. 2007; Wolstenholme et al. 2011; Zhang et al. 2012). Developmental BPA exposure in mouse models has been shown to alter the methylome (Anderson et al. 2012; Kim et al. 2014; Singh and Li 2012), including specific effects on DNA methylation at imprinted loci (Susiarjo et al. 2013). Despite these results, however, it remains unclear whether developmental BPA exposure can specifically affect DNA hydroxymethylation.

Here, we used the HMeDIP-seq method to measure epigenome-wide 5-hmC from longitudinal mouse blood samples in an effort to examine DNA hydroxymethylation at imprinted loci throughout murine adulthood. HMeDIP-seq is an antibody-based high-throughput sequencing method that measures the genome-wide distribution of 5-hmC. Although this method has some inherent antibody inefficiency and fails to provide base-pair resolution sequencing data, it does provide cost-effective amplification of high-frequency, low-signal 5-hmC marks within gene bodies (Skvortsova et al. 2017; Tan et al. 2013). Given the discovery nature of this novel study and the lack of a priori knowledge regarding 5-hmC levels in mouse blood or in response to toxicant exposures, HMeDIP-seq was chosen to measure genome-wide 5-hmC from longitudinal blood samples. From this data, we tested the hypothesis that perinatal BPA exposure alters longitudinal 5-hmC levels at imprinted genes. We also tested the hypothesis that imprinted genes would show defined 5-hmC levels that persist across the life course.

Materials and Methods

Study Animals and Blood Collection

Mice included in the longitudinal analysis were offspring sourced from a genetically invariant viable yellow agouti mouse colony maintained by sibling mating and forced heterozygosity for more than 220 generations (Waterland and Jirtle 2003). Within this colony, the allele is passed through the heterozygous male line, which has a genetically constant background 93% identical to the C57BL/6J strain (Waterland and Jirtle 2003; Weinhouse et al. 2014). Two weeks prior to mate-pairing with males, 8- to 10-wk-old wild-type dams (Control: ; BPA: ) were placed on one of two experimental diet groups: a) Control (modified, phytoestrogen-free 7% corn oil AIN-93G), or b) Control diet (see Figure S1) (Anderson et al. 2012; Weinhouse et al. 2014). For each mouse, the experimental diet was assigned using an online random number generator. After randomization, any dam litter mates assigned to the same exposure group were switched to different diets. All dams in a diet group were co-housed for the first 2 wk of exposure, prior to mate-pairing. During mate-pairing, dams were individually housed and a sire was immediately added to those cages; both the dam and sire had access to the experimental diet until pregnancy was confirmed, and then the sire was removed. In an effort to limit cage bias, all exposure groups were in the same room with the same environmental conditions, and cages were rotated on the shelves so that the same mice were not always closest to the lights or the floor. However, during preconception, gestation, and weaning, different colored chow pellets were used for the Control and BPA chow to ensure that the correct diet was placed in the correct cages. BPA dose was chosen based on a previous dosing study from our lab that showed intake of diet produced, on average, liver (Anderson et al. 2012). This exposure level is within the range of human exposure levels measured in human fetal liver samples (range: nondetect to liver) (Anderson et al. 2012). Dietary exposure was continued through pregnancy and lactation. At postnatal day (PND) 21, BPA-exposed offspring (total ; used for these analyses) were weaned to the modified AIN-93G Control diet and followed along with Control offspring (total ; used for these analyses) until 10 months of age. All diets were provided by Envigo. At 2 and 4 months of age, tail-vein blood samples were collected from all mice (see Figure S1). The mice were sacrificed at 10 months of age, and blood samples were again collected, this time using cardiac puncture. Six mouse blood samples from each exposure group were selected for inclusion in next-generation sequencing analyses. In accordance with methods established by the National Institute of Environmental Health Sciences (NIEHS) Toxicant Exposures and Responses by Genomic and Epigenomic Regulators of Transcription (TaRGET) II consortium (Wang et al. 2018), we selected male () and female () blood samples collected from six separate litters of Control and BPA-exposed mice. All animals were housed in polycarbonate-free cages with ad libitum access to food and drinking water, and were maintained in accordance with Institute for Laboratory Animal Research (ILAR) guidelines (National Research Council (US) Committee for the Update of the Guide for the Care and Use of Laboratory Animals 2011). The study protocol was approved by the University of Michigan Committee on Use and Care of Animals (UCUCA PRO00004797).

DNA and RNA Isolation

To allow for matched analyses of DNA hydroxymethylation and gene expression from the same set of mice ( per group), genomic DNA and RNA were co-isolated from frozen mouse blood that was collected at 2, 4, and 10 months of age using the Qiagen Allprep DNA/RNA/miRNA Universal Kit (Catalog #80,224; Qiagen). Yield and purity of all DNA and RNA was measured using a NanoDrop spectrophotometer. All samples were stored at prior to DNA and RNA isolations.

Real-Time Quantitative PCR

The Bio-Rad iScript cDNA Synthesis Kit (Catalog #1,708,890) was used to reverse transcribe complementary DNA (cDNA) from of RNA for each sample. In preparation for real-time quantitative polymerase chain reaction (RT-qPCR), cDNA samples were diluted 1:2.5 in RNase-free water, then mixed with forward/reverse primers, nuclease-free water, and Bio-Rad iQ SYBR Green Supermix (Catalog #1,708,880). RT-qPCR was then performed using a Bio-Rad CFX96 Real-Time System C1000 Thermal Cycler (Bio-Rad). The preprogramed 2-step curve protocol was used for all qPCR reactions: 95°C for 3 min, [95°C for 10 s, 55°C for 30 s, plate read] X 40, 95°C for 10 s. The melt curve for each plate was 65–95°C with a 0.5°C increment for 5 s and a plate read at each temperature. RT-qPCR analyses were performed for the Gnas gene in triplicate for each cDNA sample. Three housekeeping genes—Actb, 18S, and Gapdh—were included as internal controls in all RT-qPCR runs. In addition to the housekeeping genes, an inter-plate calibrator control of brain cDNA was included for calculation of relative gene expression across multiple plates. Expression levels were calculated following the method (Livak and Schmittgen 2001). RT-qPCR primers for the Gnas and Gapdh genes (see Table S1) were designed using the online Genscript Real-time PCR Primer Design software (https://www.genscript.com/tools/real-time-pcr-tagman-primer-design-tool). The and 18S gene primer pairs were sourced from the literature (Dolinoy et al. 2010; La Salle et al. 2004). Primer pair specificity for all designed primers was checked using the National Center for Biotechnology Information (NCBI) Primer-BLAST online tool (https://www.ncbi.nlm.nih.gov/tools/primer-blast/).

Next-Generation Sequencing of 5-hmC

Epigenome-wide hydroxymethylation levels were quantified in blood from mice using HMeDIP-seq (Tan et al. 2013). Unlike some other methods, HMeDIP-seq is an antibody-based approach that does not provide base-pair resolution data. This method is also subject to the inherent bias of a multistep, immunoprecipitation-based sequencing method, including potential for antibody binding inefficiency, as well as library preparation bias, read mapping bias, or PCR amplification bias (Meyer and Liu 2014). Despite these weaknesses, HMeDIP-seq is a cost-effective method to measure genome-wide 5-hmC, including at regions of high-frequency weak 5-hmC signal (Tan et al. 2013). Sequencing data were generated by C. Lalancette at the University of Michigan Epigenomics Core Facility for a subset of Control () and BPA-exposed () mice (Figure 1). Blood samples from the six mice (3 male, 3 female) in each exposure group were sequenced at three time points across the life course (2, 4, and 10 months of age), for a total sequencing data sample size of 36 samples. The six mice for each exposure group were selected from different litters to minimize litter effects. Sample quality and quantitation were assessed using the Agilent TapeStation genomic DNA kit (Catalog #G2991AA; Agilent) and Qubit broad range dsDNA (Catalog #Q32850; Invitrogen), respectively. Indexed adapters and PCR primers were synthesized by Integrated DNA Technologies (IDT). Enzymes used for library preparation were sourced from New England BioLabs (NEB).

Figure 1.

Sequencing data collection and analysis workflow. Two weeks prior to mate-pairing with males, 8- to 10-wk-old wild-type dams were placed on one of two experimental diet groups: (1) Control (modified, phytoestrogen-free 7% corn oil AIN-93G), or (2) Control diet. Dietary exposure was continued through gestation and lactation until weaning at postnatal day 21. Genomic DNA was isolated from matched wild-type offspring blood samples at 2, 4, and 10 months of age. DNA was isolated from a subset of Control ( per age group) and BPA-exposed ( per age group) mice, then processed in preparation for HsMeDIP-seq. All processed samples were amplified and sequenced on an Illumina HiSeq 4000 sequencer using single-end, reads. BPA-related differentially hydroxymethylated regions (DHMRs) were identified and annotated using a bioinformatics pipeline. Annotated DHMRs were then visualized in the genome browser. The target gene region—Gnas—was then validated using RT-qPCR on available RNA from the blood samples.

A total of of genomic DNA (gDNA) was sheared by adaptive focused acoustics, using the Covaris S220 (Catalog #4,465,653; Covaris). This sheared DNA was next blunt-ended and phosphorylated. A single adenine nucleotide was then added to the 3′ end of the fragments in preparation for the ligation of the adapter duplex with a thymine overhang. The ligated fragments were cleaned using Qiagen’s MinElute PCR purification columns (Catalog #28,004; Qiagen). DNA standards for HMeDIP (Diagenode 5-hmC, 5-mC, and cytosine DNA standard pack for HMeDIP; Catalog #AF-107-0040) were added to each sample before denaturation. Resuspension was then performed in ice-cold immunoprecipitation buffer ( sodium phosphate pH 7.0, sodium chloride (NaCl), 0.05% Triton X-100). A 10% volume (input) was retrieved before of a 5-hmC–specific antibody (Active Motif, Catalog # 39,791) was added for immunoprecipitation overnight at with rotation. Dynabeads Protein-G (Catalog #10003D; Invitrogen) was used to pull down 5-hmC–enriched fragments. The 5-hmC-enriched DNA fragments (from immunoprecipitation) were then released from the antibody by digestion with Proteinase K (Catalog #AM2548; Ambion). After cleanup with AMPure XP beads (Product #A63880; Beckman Coulter), the percentage input enrichment (%input) in the IP was evaluated by qPCR, using hydroxymethylated, methylated, and unmethylated primers for spike-ins. Samples with high %input for the 5-hmC spike-in—typical inclusion threshold was —were then PCR amplified for the final library production, cleaned using AMPure XP beads, and quantified using the Qubit assay (Catalog #Q32850; Invitrogen) and TapeStation High Sensitivity D1000 kit (Catalog #G2991AA; Agilent). Single-end, reads were obtained for each library by sequencing on the HiSeq 4000 system (Illumina). Each HMeDIP-seq sample was run on a single sequencing lane.

Next-Generation Sequencing of 5-mC

Enhanced reduced representation bisulfite (ERRBS) was performed at the University of Michigan Epigenomics Core as described previously (Akalin et al. 2012; Garrett-Bakelman et al. 2015). Briefly, of genomic DNA was digested using MspI, a restriction enzyme that preferentially cuts CG-rich sites. The digested DNA was then purified using phenol:chloroform extraction and ethanol precipitation in the presence of glycogen, before blunt-ending and phosphorylation. A single adenine nucleotide was next added to the 3′ end of the fragments in preparation for the ligation of the adapter duplex with a thymine overhang. The ligated fragments were cleaned, then processed for size selection on agarose gel. Selected fragments were treated with sodium bisulfite to convert unmethylated cytosines to uracils, which are then replaced with thymines during PCR amplification. After cleanup with AMPure XP beads (Product #A63880; Beckman Coulter), libraries were quantified using the Agilent TapeStation genomic DNA kit (Catalog #G2991AA; Agilent) and Qubit broad range dsDNA (Catalog #Q32850; Invitrogen). Single-end, reads were obtained for each library by sequencing on the HiSeq 4000 system (Illumina). ERRBS samples were multiplexed, with three samples per sequencing lane.

Sequencing Data Quality Control, Trimming, and Alignment

HMeDIP-seq data for all Control and BPA-exposed mice were compared using a suite of bioinformatics tools (Figure 1). First, the Sartor lab mint pipeline was used for data quality control (FastQC and MultiQC), adapter trimming (trim_galore), and alignment (bowtie2) (Cavalcante et al. 2017). For quality and adapter trimming, we required a minimum overlap of and minimum quality score of 20, along with special ERRBS trimming of from the 3′ end. For bismark methylation extractor (ERRBS data), the minimum threshold to consider a CpG site in analysis was five reads. Although the default minimum overlap in trim_galore is , we selected a less stringent minimum overlap of in an effort to include all legitimate genomic sequences and improve read depth. Despite these benefits, the less stringent minimum overlap length increases the possibility for adapter contamination in our data (Krueger 2017). The quality score cutoff was selected based on previous research showing an optimal tradeoff between correct read mapping and read survival at quality scores between 20 and 30 (Del Fabbro et al. 2013). The lower end of this range was selected to maximize read survival after trimming. Default parameters were used for bowtie2.

Bioinformatics Pipeline for Differential 5-hmC

After data quality check, trimming, and alignment, the csaw package was used to test for differential 5-hmC by BPA exposure (Lun and Smyth 2016). Aligned HMeDIP-seq data were read into R (version 3.4.0; R Development Core Team) using the window Counts function in csaw. When reading in the data, extension was set to 52, window width was set to 100, and sex chromosomes were removed. The data were then filtered twice—by count (average count ), and by local enrichment using the filter Windows function ()—to remove regions of negligible binding. After filtering, normalization factors were calculated for each sample using the normOffsets function on binned data (). Normalization factors were linked to filtered data using the asDGEList function, and then the estimateDisp function was used to generate dispersion factors based on a multifactorial design matrix. The glmQLFit function was used to fit a model [] for differential 5-hmC binding; an empirical Bayesian method used to stabilize the QL dispersion estimates. Contrast statements were used to extract modeling results for the variable of interest: BPA exposure. To correct for multiple testing, filtered data was clustered using a window length of ; the combineTests function was then used to compute a combined p-value for each cluster. Multiple testing correction was performed using the Benjamini-Hochberg method (Benjamini and Hochberg 1995, 1997). Differential hydroxymethylated region (DHMR) length cutoff was set at , and significance was set at a false discovery rate (FDR) of .

Bioinformatics Pipeline for Differential 5-mC

The DSS R package was used to test for differential 5-mC in ERRBS data by BPA exposure (Feng et al. 2014; Park and Wu 2016; Wu et al. 2013, 2015). Within the DSS package, the DMLfit.multiFactor function was used to test for differential methylation by BPA exposure using a multifactorial modeling approach according to the following formula: . The DMLtest.multiFactor function was used to test the null hypothesis that the BPA exposure coefficient was equal to 0. Differentially methylated CpG sites (DMCs) by BPA exposure were then sorted and filtered according to a FDR cutoff of .

After testing for differential hydroxymethylation using csaw, the annotatr R package was used to annotate all DHMRs to the mm10 genome (Cavalcante and Sartor 2017). The annotate_regions function was used to generate genomic annotations, which include classifications by region class (e.g., intron, exon, promoter) and annotated gene identifiers (IDs). A list of mouse imprinted loci was then sourced from the MouseBook online database (Williamson et al. 2013). All available imprinted genes were manually cross-checked with the generated list of BPA-related DHMR-annotated gene IDs.

Sequencing Data Visualization

5-hmC levels were visualized using the csaw and GViz R packages (Hahne and Ivanek 2016; Lun and Smyth 2016). Using csaw, data was first read into R as a GRanges object using the extractReads function. Next, the GeneRegionTrack in the Gviz R package was used to define regions of interest for visualization. Separate blue and red genome tracks were used to represent the forward and reverse strand reads, respectively. The plotTracks function was then used to plot 5-hmC levels for these regions. Reads-per-million was used for the y-axis in all graphs, and scale was adjusted to individual region coverage. 5-mC levels were visualized using the RnBeads R package (Assenov et al. 2014). Within this package, the rnb.execute.import function was used to import the aligned ERRBS data. Next, the rnb.sample.groups function was used to group samples by age, and the rnb.plot.locus.profile function was used to plot relative 5-mC level in a heat map output for defined regions of the genome. A custom color palette was chosen from RColorBrewer for the heat map gradient.

BPA-related DHMRs at imprinted loci were visualized using the GViz R package to determine directionality and magnitude of differential 5-hmC. Plots generated in R (version 3.4.0; R Development Core Team) were formatted for publication in Adobe Illustrator CS6 (version 16.0.5).

Pathway Analysis

Pathway analysis for DHMRs was performed using ChIP-enrich (Welch et al. 2014). Within the ChIP-enrich online interface (http://chip-enrich.med.umich.edu/), the gene set filter was set to 2,000, peak threshold was set to 1, and adjustment for mappability of gene locus regions was set to false. The genome used for pathway analyses was mm10, and the ChIP-enrich method was used for enrichment testing. DHMRs were split into hypo- and hyper-hydroxymethylated regions prior to separate analyses. Only regions from TSS were included in the ChIP-enrich analysis. All pathway analyses included Gene Ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways. After correcting for multiple testing, only pathways with an were considered significant.

Results

Epigenome-Wide Differential 5-hmC and 5-mC

Using the csaw R package, multivariate models were constructed to test for differential 5-hmC by perinatal BPA exposure ( diet). Based on the csaw models, we identified 5,950 DHMRs by BPA exposure (Table 1). Comparing the directionality of DHMRs by BPA exposure, we found more hypo-hydroxymethylated regions (; 71.4%) than hyper-hydroxymethylated regions (; 26.2%) (Table 1). We also identified a small fraction of BPA-related DHMRs (; 2.4%) that showed both hypo- and hyper-hydroxymethylation. Taken together, these epigenome-wide results show that the directionality of differential hydroxymethylation varies by region, with a skew toward decreased 5-hmC by BPA exposure. Using the annotatr R package, BPA-related DHMRs were annotated to the mm10 genome. Compared with a random distribution, BPA-related DHMRs were depleted at CpG shores, CpG shelves, promoters, 3′-UTRs, 5′-UTRs, exons, and introns, but slightly enriched at regions greater than from a CpG island (interCGI) (see Figure S2).

Table 1

Differential 5-hmC in mouse blood by BPA exposure.

Δ5-hmC	BPA
	DHMRsa
Hyper-hydroxymethylated	1,559
Hypo-hydroxymethylated	4,247
Both	144
Total	5,950

Note: The csaw R package was used to examine the effect of BPA exposure ( diet) on 5-hmC levels across the epigenome. Models included a paired mouse ID variable to account for within-individual effects, as well as a sex variable. Directionality of DHMRs (relative to Control) is indicated in separate rows. BPA, bisphenol A; DHMR, differentially hydroxymethylated region; ID, identification; , direction of DNA hydroxymethylation change (BPA exposed relative to Control).

DHMR-calling significance threshold was set at a false discovery rate (FDR) of .

To investigate differential DNA methylation by BPA exposure, ERRBS sequencing data were generated from the same blood DNA samples used for HMeDIP-seq data. Unfortunately, due to biased and limited genomic coverage inherent in the ERRBS method, overlapping DNA methylation data were available at only one of the identified imprinted gene DHMRs: Kruppel-like factor 14 (Klf14). At the Klf14 DHMR, there were no significant changes in ERRBS data by BPA exposure (see Figure S3). However, in our differential methylation analyses of the ERRBS data by BPA exposure, we identified a number of DMCs annotated to each of the top six imprinted loci that had DHMRs (see Table S2). Of note, we identified three CpG sites within the guanine nucleotide binding protein, alpha stimulating, complex locus and the neuroendocrine secretory protein antisense (Gnasxl/Nespas) ICR that demonstrated significant hypomethylation with BPA exposure (see Table S2). None of the identified DMCs at the other top imprinted genes—growth factor receptor bound protein 10 (Grb10), pleiomorphic adenoma gene-like 1(Plagl1), phosphodiesterase 10A (Pde10a), Klf14, and small nuclear ribonucleoprotein polypeptide N (Snrpn)—fell in known murine ICRs.

BPA-Related DHMRs at Imprinted Genes

Previous research has shown effects of developmental exposures on DNA methylation at imprinted loci (Gallou-Kabani et al. 2010; Susiarjo et al. 2013), but it remains unclear whether exposures can also alter DNA hydroxymethylation at imprinted genes. Further complicating this question, despite recent research in humans showing enrichment of 5-hmC at imprinted loci (Hernandez Mora et al. 2018), the distribution and potential function of 5-hmC at murine imprinted genes is unknown. In an effort to broadly identify BPA-sensitive regions of 5-hmC at imprinted loci, we cross-checked the BPA-related DHMR annotations with a database of known murine imprinted genes (Williamson et al. 2013). In total, 12 of the 151 known imprinted genes had annotated BPA-related DHMRs (Table 2). Five of these 12 genes had increased 5-hmC peaks by BPA exposure—Grb10, Pde10a, phosphodiesterase 4D (Pde4d), Plagl1, and Gnas. The remaining 7 imprinted loci—protein phosphatase 1, regulatory subunit 9A (Ppp1r9a), phosphatase and actin regulator 2 (Phactr2), Klf14, potassium voltage-gated channel, subfamily Q, member 1 (Kcnq1), cytidine monophospho-N-acetylneuraminic acid hydroxylase (Cmah), antisense Igf2r RNA (Airn), and Snrpn—had decreased 5-hmC peaks by BPA exposure. The number of CpG sites within these DHMRs varied by region, with 5 of the DHMRs showing zero CpG sites (see Table S3). In an effort to confirm the DHMRs, we visualized the 5 most significant imprinted gene DHMRs and 1 DHMR at the well-characterized Snrpn gene, showing marked changes in 5-hmC peaks by BPA exposure that persisted across all of the three matched sample ages (Figure 2). We also found striking changes in 5-hmC peaks at the remaining 6 imprinted gene DHMRs, despite their decreased significance (see Figure S4). In an effort to determine whether this was a shared trait at imprinted genes, we visualized 5-hmC peaks at insulin like growth factor 2 (Igf2) and H19, imprinted maternally expressed transcript (H19), 2 well-characterized imprinted genes that did not have any annotated DHMRs with . Despite their lack of significant DHMRs, Igf2 and H19 still had specific 5-hmC peaks that showed nonsignificant changes by BPA exposure (see Figure S5).

Table 2

Imprinted genetic loci with BPA-related DHMR.

Chr	Start	End	Gene ID	Gene Name	Location	Direction with BPA exposure	csaw score	FDR
chr2	174315451	174315750	14683	Gnas	Intronic	+	113.4	<0.001
chr11	12004801	12005100	14783	Grb10	Intronic	+	70.5	<0.001
chr10	13130651	13131050	22634	Plagl1	Exonic	+	60.7	<0.001
chr6	30957751	30958050	619665	Klf14	Exonic	−	55.9	<0.001
chr17	8746901	8747050	23984	Pde10a	intronic	+	50.1	<0.001
chr17	12828001	12828550	104103	Airn	Intronic	−	16.2	0.02
chr13	24335001	24335150	12763	Cmah	Intronic	−	13.1	0.05
chr7	60207401	60207600	20646	Snrpn	Intronic	−	13.0	0.05
chr6	5079501	5079850	243725	Ppp1r9a	Intronic	−	12.8	0.05
chr7	143348051	143348350	16535	Kcnq1	Intronic	−	11.5	0.07
chr10	13402051	13402250	215789	Phactr2	Intronic	−	11.3	0.07
chr13	108963951	108964200	238871	Pde4d	Intronic	+	10.1	0.10

Note: List of imprinted loci with annotated, significant () DHMRs generated from csaw modeling results. For csaw results, a transformed FDR is used for the score, such that ; higher score indicates a larger degree of significance (lower FDR value). Additionally, direction of effect by BPA exposure is shown, with a () indicating increased 5-hmC with BPA exposure and a () indicating decreased 5-hmC with BPA exposure. BPA, bisphenol-A; Chr, chromosome; DHMR, differentially hydroxymethylated region; FDR, false discovery rate; ID, identification.

Figure 2.

Differential imprinted gene 5-hmC peaks by BPA exposure. 5-hmC coverage was visualized at six imprinted loci with significant BPA-related DHMRs: (A) Gnas; (B) Grb10; (C) Plagl1; (D) Klf14; (E) Pde10a; and (F) Snrpn. 5-hmC levels are shown for matched 2-, 4-, and 10-month blood samples, as indicated by y-axis labels. Blue and red peaks represent forward and reverse strand 5-hmC enrichment, respectively. The Gnas, Plagl1, and Pde10a DHMRs occurred on the forward strand, whereas the Grb10, Klf14, and Snrpn DHMRs occurred on the reverse strand.

In examining the imprinted loci visually, two of the DHMRs—Plagl1 and Gnas—had 5-hmC peaks that were cut off by the DHMR boundaries, suggesting larger scale 5-hmC patterns at these genes. To provide a complete picture for these two loci, 5-hmC levels across the entire gene were visualized. Both Plagl1 and Gnas had widespread 5-hmC peaks along the entire gene-coding region, with only the annotated DHMRs showing apparent modifications by BPA exposure (Figure 3). Of note, 5-hmC peaks across these two genes were consistent between mice and across time, a pattern that was also observed at the Igf2 and H19 imprinted genes (see Figure S5). Together, these data show that imprinted genes have predictable 5-hmC peaks along the length of their gene bodies.

Figure 3.

Genomic context of Gnas and Plagl1 DHMRs. The DHMRs annotated to the Gnas and Plagl1 genes were visualized in the context of their respective imprinted genes using the csaw and Gviz R packages. 5-hmC levels across the complete Gnas and Plagl1 imprinted loci are presented, showing distinct 5-hmC peaks along the length of both genes. 5-hmC levels are shown for matched 2-, 4-, and 10-month blood samples, as indicated by y-axis labels. Boxes indicate significant DHMRs that were identified using csaw differential hydroxymethylation models. Blue and red peaks represent forward and reverse strand 5-hmC enrichment, respectively. Bars below University of California, Santa Cruz genome browser gene tracks represent CpG islands.

RT-qPCR Gene Expression Data

To follow up on the DHMR annotated to the Gnas gene, we performed RT-qPCR on RNA from matched mouse blood samples to examine longitudinal mRNA expression. At the Gnas locus, there was a significant increase in mean mRNA expression from 2 to 10 months of age in mice exposed to BPA () (Figure 4). This increase was not significant in Control samples. Additionally, mean Gnas expression was significantly lower in Control mice than BPA-exposed blood at 10 months of age ().

Figure 4.

Gnas expression by BPA exposure and age. Based on BPA-related DHMR annotated to the Gnas gene, real-time quantitative polymerase chain reaction (RT-qPCR) was used to investigate longitudinal blood Gnas mRNA expression levels. RNA was isolated from the same longitudinal Control ( per age group) and BPA-exposed ( per age group) mouse blood samples used for DNA hydroxymethylation analyses. RT-qPCR was performed on the Gnas gene in triplicate. Three housekeeping genes—, 18S, and Gapdh—were included as internal controls in all RT-qPCR runs. In addition to housekeeping genes, an inter-plate calibrator control of brain cDNA was included for calculation of relative gene expression across multiple plates; all expression values are shown relative to this inter-plate calibrator. Expression levels were calculated following the method. aMean Gnas expression in BPA-exposed mouse blood showed a significant increase between 2 and 10 months of age (); this pattern was not found in Control samples. bMean Gnas expression was significantly lower in Control blood than BPA-exposed blood at 10 months of age ().

Pathway Analysis for DHMRs

The ChIP-enrich pathway analysis was performed on separate lists of hyper- and hypo-methylated DHMRs. In an effort to maximize biological relevance and limit number of comparisons, DHMR datasets were restricted to regions from the transcription start site. After correction for multiple testing, a small number of gene ontology (GO) terms were enriched in hypo- and hyper-hydroxymethylated DHMRs (see Table S4). Of note, enriched pathways showed no overlap between hyper- and hypo-hydroxymethylated DHMRs. In the BPA-related hypo-hydroxymethylated DHMRs, the only two significantly enriched pathways were hippocampus development and spinal cord association neuron differentiation. In the BPA-related hyper-hydroxymethylated DHMRs, the two most significant enriched pathways were l-ascorbic acid binding and oxidoreductase activity. These pathway analysis results indicate that only a few GO terms were enriched in the significant DHMRs, and that enrichment was specific to directionality of BPA-related DHMRs.

Discussion

Epigenome-Wide Differential 5-hmC

We found statistically significant effects of perinatal BPA exposure on DNA hydroxymethylation across the epigenome. Although the exact mechanism driving this relationship is unclear, it is possible that BPA-related DHMRs are a result of BPA-induced oxidative stress (OS) (Gassman 2017). Based on the available literature, free BPA is thought to induce OS via enzymatic formation of BPA phenoxyl radicals (Sakuma et al. 2010), which can then be further processed to other reactive oxygen species (ROS), including superoxides and peroxides (Babu et al. 2013). In addition to potential ROS-induced cytotoxicity (Gassman 2017), the pro-oxidant activity of BPA also has the potential to modify 5-hmC levels across the epigenome. Research has shown that Tet enzyme activity is activated under oxidative conditions (Chia et al. 2011; Coulter et al. 2013; Zhao et al. 2014), indicating that active processing of 5-mC to 5-hmC may be affected by ROS production. Further supporting this idea, recent research has shown that OS-inducing exogenous factors—buthionine sulfoximine and fine particulate matter with an aerodynamic diameter of —can alter DNA hydroxymethylation levels (Delatte et al. 2015; Wei et al. 2017). These prior results support the hypothesis that BPA exposure could be triggering changes in 5-hmC via the induction of OS. This idea should be further explored in future studies examining the effects of endocrine disrupting chemicals on 5-hmC levels across the epigenome.

In addition to examining the effects of BPA on broad-scale DNA hydroxymethylation, we also specifically investigated whether BPA exposure altered 5-hmC levels at murine imprinted genes. Of the 151 interrogated imprinted loci, 12 had annotated DHMRs (7.95%) (Table 2). As a broad-scale comparison, of the estimated 24,360 known protein-coding genes in mice (Blake et al. 2017), 1,616 unique gene IDs had annotated DHMRs (6.63%). These results indicate that imprinted genes show slight enrichment for DHMRs compared with all known genes combined. Although none of the 12 imprinted gene BPA-related DHMRs overlapped with known ICRs (Figure 5), expanded visualization of Gnas and Plagl1 showed widespread 5-hmC peaks across imprinted gene regions. The directionality and magnitude of the imprinted gene DHMRs was also gene specific. Of note, 5 of the imprinted gene DHMRs had zero CpG sites (see Table S3), a result that could be driven by several possible scenarios: a) DNA hydroxymethylation at noncanonical methylation contexts (i.e., CHG sites, where , C, or T); b) inexact DHMR definitions due to inherent lack of specificity in pulldown method (HMeDIP-seq); or c) false positives in our data. Supporting the first option, recent work in mice has shown intragenic enrichment of non-CpG methylation at particular domains in clusters of genes related to embryonic development (He et al. 2017). Further fitting with this idea, 8 of the 12 imprinted gene DHMRs had at least 10 CHG sites within their chromosomal ranges (see Table S3). Nevertheless, it remains difficult to distinguish between the three outlined scenarios using the data available in this project. The most significant imprinted gene DHMR was annotated to the Gnas locus, which encodes the G-protein alpha subunit protein, a key component of G-protein coupled signal transduction (Plagge and Kelsey 2006). Gnas is an imprinted gene that has a complex expression pattern, four alternative promoters, a number of isoforms, and may be involved in energy homeostasis (Peters and Williamson 2007; Plagge and Kelsey 2006). Given the complexity of transcriptional control at this locus, the identified intronic Gnas DHMR may represent a long-range regulatory region. This hypothesis is supported by research demonstrating that 5-hmC is enriched at enhancers (Ehrlich and Ehrlich 2014; Stroud et al. 2011; Sun et al. 2013; Wen et al. 2014), regulatory regions that can be quite distant from the gene promoters they activate (Pennacchio et al. 2013; Shlyueva et al. 2014). Based on this idea of distant regulatory regions, we hypothesized that the documented increase in intronic Gnas 5-hmC would show a corresponding BPA-related change in Gnas expression. Using RT-qPCR, we found increased Gnas expression by BPA exposure across all three time points, with the magnitude of this increase reaching significance at 10 months of age (Figure 4). Given the complex regulation of this locus, the age-specific effect of BPA exposure on Gnas expression may be a result of changes in alternative splicing, which is dynamic during the aging process (Li et al. 2017). During aging, Gnas splicing could shift toward specific isoforms that are controlled by 5-hmC at the BPA-related DHMR, leading to a late adulthood effect of developmental BPA exposure. Future functional work should explore this idea, as well as test for potential effects of BPA exposure at the human GNAS locus.

Figure 5.

Organization of the Gnas, Grb10, Plagl1, Klf14, Pde10a, and Snrpn imprinted loci. The thick black midline represents the gene sequence. Maternally expressed exons are indicated by black boxes, and paternally expressed exons are indicated by gray boxes. Exon locations for each gene are relative, but not to scale. Exon 1 of the Gnas gene is both maternally and paternally expressed and is therefore indicated by a white box. Arrows show directionality and start sites for transcription. Parental origin-specific expression patterns are represented by the location of the arrow above or below the midline, and context-specific expression is indicated by dotted arrows. Filled and empty circles represent methylated or unmethylated alleles, respectively. Imprinting control regions (ICRs) are indicated by dotted line boxes, and BPA-related DHMRs are indicated by boxes bordered by dotted vertical lines. None of the identified BPA-related DHMRs are located in ICRs. Figure for Gnas adapted from Tibbit et al. (2015); Grb10 from Hikichi et al. (2003); Plagl1 from Iglesias-Platas et al. (2012) and Smith et al. (2002); Klf14 from Parker-Katiraee et al. (2007); Pde10a from Wang et al. (2011); and Snrpn from Sanli and Feil (2015) and Shemer et al. (1997).

In addition to Gnas, three other imprinted genes—Grb10, Plagl1, and Pde10a—showed significant increases in 5-hmC by BPA exposure. The first of these additional hyper-hydroxymethylated genes, Grb10, encodes the growth factor receptor-bound protein 10 (Grb10). Grb10 is involved in a number of biological processes, including cellular proliferation, apoptosis, and metabolism (Holt and Siddle 2005; Kabir and Kazi 2014; Plasschaert and Bartolomei 2015; Riedel 2004). Grb10 is maternally expressed in most tissues, but it is paternally expressed in the brain; this complex imprinting pattern has been established through tissue-specific alternate promoters (Sanz et al. 2008). The tissue-specific maternal and paternal expression patterns of Grb10 have been linked to fetal growth and adult social behavior in mice, respectively (Garfield et al. 2011). This complex expression patterning is thought to be controlled via epigenetic modifications (Sanz et al. 2008). Like Gnas, the Grb10 DHMR is intronic, meaning any functional effects of differential 5-hmC at this region would have to be through long-distance contacts. The second additional gene, Plagl1, encodes the pleomorphic adenoma of the salivary gland gene like 1 (Plagl1) protein. Plagl1 is a zinc finger transcription factor that regulates other imprinted loci involved in embryonic growth, including H19 and Igf2 (Varrault et al. 2006, 2017). As such, BPA-related alterations in 5-hmC at Plagl1 have the potential to affect an entire network of imprinted genes. Like Gnas and Grb10, Plagl1 has several alternative transcripts in the mouse that are controlled by alternate promoters (Iglesias-Platas et al. 2013). The Plagl1 DHMR overlaps Exon 8 in some Plagl1 isoforms, but not others, suggesting that this region may play a role in alternative splicing for this gene. The third additional hyper-hydroxymethylated imprinted gene, Pde10a, encodes phosphodiesterase 10A (Pde10a). Pde10a is a member of the cyclic nucleotide phosphodiesterases (PDEs), a family of enzymes that regulate intracellular levels of cyclic AMP (cAMP) and cyclic GMP (cGMP), endogenous molecules involved in signal transduction (Bender and Beavo 2006; Soderling et al. 1999). Pde10a shows its highest expression levels in the brain, with minimal expression in peripheral tissues (Bender and Beavo 2006; Soderling et al. 1999). Recent research shows that pharmaceutical inhibition of Pde10a leads to increased energy expenditure, decreased food intake, reduced adiposity, and improved insulin sensitivity in mice with high-fat diet–induced obesity (Nawrocki et al. 2014). As such, altered regulation of the Pde10a gene has the potential to modify murine energy homeostasis. Although we found a Pde10a DHMR in blood tissue, where expression is minimal, our developmental BPA exposure occurred throughout tissue differentiation, so it is possible that the BPA-related DHMR annotated to Pde10a is present in multiple tissues, including the brain.

We also visualized two hypo-hydroxymethylated DHMRs annotated to imprinted genes—Klf14 and Snrpn. The first of these genes, Klf14, encodes the maternally expressed Krüppel-like factor 14 (Klf14), a member of the Cys2/His2 zinc finger transcription factors (Small et al. 2011; Wang et al. 2017). Research shows that KLF14 acts as master regulator of adipose gene expression in humans (Small et al. 2011), and may be involved in metabolic disease risk. In addition, recent results in mice suggest that Klf14 may interact with peroxisome proliferator-activated coactivator (), a transcription coactivator that regulates a number of metabolic pathways, including hepatic gluconeogenesis (Wang et al. 2017). Based on these prior studies, changes in 5-hmC at Klf14 have the potential to not only alter murine energy homeostasis, but also to modify the risk of metabolic disorders. The second imprinted gene to show BPA-related decrease in 5-hmC was Snrpn, a paternally expressed gene that encodes the SmN protein, a key component of the spliceosome in the brain (Shemer et al. 1997). Snrpn has multiple alternate promoters and splice variants, is directly related to neurological function, and is located in the human Prader-Willi and Angelman syndrome (PWS/AS) imprinted domain (Wu et al. 2012). The PWS/AS imprinted domain is highly complex, containing a large number of C/D box small nucleolar RNAs (snoRNAs) and two ICRs: Prader-Willi Syndrome Imprinting Control (PWS-IC) and Angelman Syndrome Imprinting Control (AS-IC) (Wu et al. 2012). Research suggests the snoRNAs play a role in the etiology of the murine PWS phenotype (Relkovic and Isles 2013; Skryabin et al. 2007), which is often characterized by weight gain, decreased activity, and impaired attention (Relkovic and Isles 2013). Here, we identified a BPA-related hypo-hydroxymethylated region in an intron of Snrpn. Although the functional relevance of this region remains unknown, previous work has shown that developmental BPA exposure alters DNA methylation at SNRPN/Snrpn (Faulk et al. 2015; Susiarjo et al. 2013), providing support for the idea that BPA can affect the epigenome at this locus. Additionally, the human SNRPN domain is enriched for 5-hmC across the expressed allele in brain tissue, indicating that 5-hmC may have a functional role in control of SNRPN gene expression (Hernandez Mora et al. 2018). Therefore, it is possible that BPA-related changes in Snrpn 5-hmC could alter gene transcription, thereby modifying risk for PWS, a neurobehavioral disorder.

Remarkably, BPA-related DHMRs at the described imprinted loci persisted throughout adulthood despite BPA exposure ending at PND21. These data indicate that perinatal exposure to BPA can have long-lasting effects on 5-hmC. Given the complex role of 5-hmC in regulating transcription activation (Hahn et al. 2014; Wu et al. 2011), the identified BPA-related DHMRs may reflect programed changes in gene regulation. Supporting this idea, a number of recent studies have shown that 5-hmC is a stable epigenetic mark that has an important role in gene regulation (reviewed by López et al. (2017). Additionally, previous work has shown that there is a unique binding protein (Mbd3) that recognizes 5-hmC (Yildirim et al. 2011), suggesting that 5-hmC has its own epigenetic function. Expanding on this idea, a review of the 5-hmC literature suggested that 5-hmC may be maintained across DNA replication via a complex of the DNMT1, Tet, and UHRF1 proteins (Shen and Zhang 2013), although the validity of this theory has not yet been determined. Building on these ideas, our results show differential 5-hmC at imprinted genes with complex regulatory patterns, indicating that 5-hmC marks at these genes may play a role in the establishment of alternative splicing in response to BPA exposure. Though limited, our RNA expression data support this idea, showing that BPA exposure has long-term functional consequences at the imprinted Gnas locus. Based on these results, a new hypothesis has emerged—that differential 5-hmC at imprinted genes could be an important mechanism driving the developmental origins of adult disease.

In addition to the BPA-specific DHMRs, our data separately show that widespread 5-hmC peaks are established across entire imprinted genes during development, and that these patterns persist throughout life. This was apparent at genes with annotated DHMRs—Gnas and Plagl1—and well-characterized genes without annotated DHMRs—Igf2 and H19. These results match recent research in humans, which showed 5-hmC enrichment at the H19-IGF2 locus in brain and at GNAS A/B in placenta (Hernandez Mora et al. 2018). Combined, the available data indicate that 5-hmC may be involved in imprinting control. As technologies for measuring genome-wide 5-hmC continue to advance, efforts should be made to further examine the regulatory role of 5-hmC in genomic imprinting.

Due to the biased and limited genomic coverage of the ERRBS method, we were not able to compare 5-hmC and 5-mC data at most of our identified DHMRs. At the Klf14 imprinted locus, we showed that 5-mC levels at the identified DHMR did not significantly differ by BPA exposure (see Figure S3). This lack of significance may be a result of inconsistent coverage of ERRBS data across samples at this region; alternatively, 5-hmC at this region may be uniquely sensitive to BPA exposure. Future studies should utilize alternative methods—such as oxidative bisulfite treatment—to coinvestigate 5-hmC and 5-mC at identified imprinted gene DHMRs.

Given that the identified DHMRs had almost no coverage in our ERRBS data, we also looked for differential methylation at known imprinted gene ICRs. Of particular note, we identified three CpG sites within the Gnasxl/Nespas ICR that showed significant BPA-related hypomethylation (see Table S2). In previous studies, hypomethylation at this ICR has been linked to increased expression of the Gnas locus (Coombes et al. 2003; Williamson et al. 2006), a pattern that matches our expression data. These results, combined with the previously identified DHMR, suggest that developmental BPA exposure has a dynamic effect on the epigenetic landscape at the Gnas imprinted gene. Future studies should utilize chromosomal conformation capture techniques to examine whether our identified DHMR interfaces with the ICR to assist in gene regulation.

GO terms showed specific enrichment based on directionality of differential hydroxymethylation. This suggests that BPA exposure is associated with both up- and down-regulation of several biological processes. In the BPA-related hypo-hydroxymethylated DHMRs, the only two significantly enriched pathways were hippocampus development and spinal cord association neuron differentiation. Although it is difficult to interpret these pathways from blood samples, 5-hmC has its highest levels in brain tissue, where it is suspected to play a role in neuron development (Kinde et al. 2015). As such, it is possible that these enriched pathways reflect a predifferentiation effect of developmental BPA exposure on 5-hmC levels in stem cells. Previous work has shown that epigenetic marks at imprinted genes are maintained during postfertilization reprograming (Plasschaert and Bartolomei 2014), suggesting that the epigenetic effects of developmental exposure could be maintained at imprinted loci during cellular differentiation. Building on this idea, our results indicate that 5-hmC establishment at imprinted genes could be involved in regulating neural plasticity. Additional work in predifferentiated cells and brain tissue could help determine whether BPA actually alters neuronal differentiation via 5-hmC.

Limitations

Little is known about the stability of 5-hmC at imprinted loci across blood cell type proportions, but previous studies have shown that epigenetic marks can vary across blood cell types (Houseman et al. 2012; Skinner 2016). As such, it is possible that the documented effects of BPA on genome-wide 5-hmC in blood are simply a reflection of shifts in blood cell type proportions. Counter to this idea, however, previous studies have shown that DNA methylation is often conserved across different cell types at imprinted loci (Skinner 2016; Talens et al. 2010), suggesting that investigating differential methylation from whole blood should still be valid at imprinted genes. Although this evidence supports the validity our BPA-related DHMRs at imprinted loci, the dynamics of 5-hmC at these imprinted genes across blood cell types remains to be elucidated. To further explore the idea that altered 5-hmC in blood impacts imprinted gene regulation, future studies should also investigate allele-specific expression in individual cell types using new single-cell RNA-seq methods (Santoni et al. 2017). Despite lingering blood cell type questions, the use of matched blood samples allowed for direct measurement of 5-hmC from the same mice over time, reducing the potential confounding of interindividual variability. Additionally, matched blood samples provide greater translatability to human epigenetics studies, which typically rely on peripheral blood samples.

Due to low RNA yields from the longitudinal blood samples, it was not possible to examine gene expression at all identified imprinted loci with annotated DHMRs. RNA yields were diminished due to the small amount of blood collected during in vivo tail vein collection at 2 and 4 months of age. Despite the low amounts of blood RNA available, the use of longitudinal samples allowed for direct measurement of expression at the Gnas imprinted locus from the same longitudinal samples used for HMeDIP-seq data generation. As such, BPA-related changes in Gnas expression directly coincide with BPA-related alterations in DNA hydroxymethylation at this locus.

DNA hydroxymethylation has only recently become recognized as an important consideration in the field of genomic imprinting, meaning its role in imprinting control is poorly defined. As a result, the exact functional effects of the BPA-related changes shown in this paper remain undetermined. So far, the available research suggests that 5-hmC can have very different regulatory effects depending upon its genomic context, so the relationship between BPA-related DHMRs and gene expression could vary in a gene-specific manner. Similarly, the phenotypic effects of BPA-related DHMRs remain unclear without additional measures of murine biology throughout the life course. Despite these limitations in interpretation, we identified a large number of BPA-related DHMRs, indicating that DNA hydroxymethylation is sensitive to environmental factors during development.

Conclusions

We measured 5-hmC in matched blood samples collected from isogenic mice at 2, 4, and 10 months of age, and then examined the effects of BPA on the longitudinal DNA hydroxymethylation. Across the epigenome, we identified a number of exposure-related DHMRs, suggesting that perinatal BPA exposure can alter this DNA modification throughout life. At 12 imprinted loci, including the Gnas locus, developmental BPA exposure significantly altered 5-hmC levels in blood across all three measured ages. Echoing this result, we showed that BPA exposure modified Gnas expression throughout murine adulthood. These results suggest that BPA-related increases in Gnas hydroxymethylation may have long-lasting effects on gene expression at this complex imprinted locus, possibly through shifts in alternative splicing. In addition to BPA-related DHMRs, we also found stable patterns of 5-hmC along a number of imprinted loci, including Igf2 and H19. Combined, our data indicate that 5-hmC may play an important regulatory role in imprinting control, and that establishment of DNA hydroxymethylation at imprinted loci may be sensitive to developmental BPA exposure. Future studies should examine the contribution of DNA hydroxymethylation to the methylome at imprinted loci, as well as the impact of additional environmental exposures on this recently rediscovered epigenetic mark.

Click here for additional data file.