Epidemiological studies show that exposure to the organochlorine pesticide dieldrin is associated with increased risk of Parkinson’s disease (PD). Animal studies support a link between developmental dieldrin exposure and increased neuronal susceptibility in the α-synuclein preformed fibril (α-syn PFF) and MPTP models in adult male C57BL/6 mice. In a previous study, we showed that developmental dieldrin exposure was associated with sex-specific changes in DNA modifications within genes related to dopaminergic neuron development and maintenance at 12 weeks of age. Here, we used capture hybridization-sequencing with custom baits to interrogate DNA modifications across the entire genetic loci of the previously identified genes at multiple time points – birth, 6 weeks, 12 weeks, and 36 weeks old. We identified largely sex-specific dieldrin-induced changes in DNA modifications at each time point that annotated to pathways important for neurodevelopment, potentially related to critical steps in early neurodevelopment, dopaminergic neuron differentiation, synaptogenesis, synaptic plasticity, and glial-neuron interactions. Despite large numbers of age-specific DNA modifications, longitudinal analysis identified a small number of DMCs with dieldrin-induced deflection of epigenetic aging. The sex-specificity of these results adds to evidence that sex-specific responses to PD-related exposures may underly sex-specific differences in disease. Overall, these data support the idea that developmental dieldrin exposure leads to changes in epigenetic patterns that persist after the exposure period and disrupt critical neurodevelopmental pathways, thereby impacting risk of late life diseases, including PD.  

Animals: Adult female and male C57BL/6 mice (RRID:MGI:2159769) were purchased from Jackson Laboratory (Bar Harbor, Maine). Seven-week-old female mice were allowed to habituate after arrival for 1 week prior to beginning the developmental dieldrin exposure. Male mice were 11 weeks old upon arrival and were also allowed 1 week to habituate prior to mating. Mice were maintained on a 12-hour:12-hour reverse light/dark cycle. Mice were housed in Thoren ventilated caging systems with automatic water and 1/8-inch Bed-O-Cobs bedding with Enviro-Dri for enrichment. Food and water were available ad libitum. Mice were maintained on Teklad 8940 rodent diet (Envigo). After arrival, females were separated and individually housed during dieldrin dosing, except during the mating phase. After birth, F1 pups were group-housed by sex, with no more than 5 animals housed in each cage. No singly housed F1 animals were used in this study; as animals aged, all cages included a “buddy” littermate to ensure that even animals used for the last time point were never individually housed. All procedures were conducted in accordance with the National Institutes of Health Guide for Care and Use of Laboratory Animals and approved by the Institutional Animal Care and Use Committee at Michigan State University.

Developmental dieldrin exposure: Dieldrin exposure was carried out as previously described and is summarized in Figure 1 (Richardson et al., 2006,Gezer et al., 2020,Kochmanski et al., 2019,Boyd et al., 2023). Adult (8-week-old) C57Bl/6 female mice (n=20 per group) were treated throughout breeding, gestation, and lactation. Following 3 days of habituation to peanut butter feeding, mice were administered 0.3 mg/kg dieldrin (ChemService) dissolved in corn oil vehicle and mixed with peanut butter pellets every 3 days (Richardson et al., 2006,Gezer et al., 2020,Kochmanski et al., 2019,Boyd et al., 2023). Control mice received an equivalent amount of corn oil vehicle mixed in peanut butter. Mice were exposed via oral ingestion by the dam because the most likely route of exposure to dieldrin in humans is through ingestion of contaminated foods and ingestion of the resulting contaminated breast milk (ATSDR, 2022).

The dieldrin dose was based on previous results showing low toxicity, but clear effects on the epigenome and neuronal susceptibility to neurotoxic insults (Richardson et al., 2006,Gezer et al., 2020,Kochmanski et al., 2019,Boyd et al., 2023). In addition, because dieldrin is a persistent compound, even though this daily dose is likely higher than human doses, this dose was chosen to produce a similar body burden in animals compared to known levels in humans. There is a very wide range of values reported for both brain (~4 ppb – 1 ppm) and adipose (~50 ppb – 50 ppm) tissue levels (U.S. Environmental Protection Agency, 1994,Corrigan et al., 2000,Agency for Toxic Substances and Disease Registry, 2022,Fleming, L. et al., 1994) Our measurements in this cohort in dam brain and adipose tissue are within the reported ranges of multiple studies (Table 1). In this study, we measured these as a quality control step on a subset of animals, but the range of levels between animals suggests that there are variable exposure levels in the pups and that measuring dieldrin in dam and pup tissue in future studies is warranted.

After four weeks of exposure, unexposed C57BL/6 males (12 weeks old) were introduced for breeding for 48 hours such that male mice were not present for any peanut butter pellet feedings. To ensure adequate animal numbers for each sex at all four time points, this was carried out in two cohorts: one to generate animals for the birth time point, and one to generate animals for the 6-week, 12-week, and 36-week time points. For the birth cohort, pups were sacrificed at PND0, tissue was collected for sex determination by genotyping, and whole brains were dissected and frozen. For the remaining cohort, F1 pups were weaned and separated by litter and by sex at 3 weeks of age, with 2-5 animals per cage. At each time point, male and female littermates from independent litters were selected. Animals were assigned to cages such that for each experimental group and sex, all animals were from independent litters and no animals were singly housed.

Group sizes were determined based on the previously published RRBS data that this analysis was based  (Kochmanski et al., 2019). Based on that data, we determined prior to starting the exposure cohorts that a sample size of at least 7 per sex per treatment group per time point was sufficient to detect >80% of all true differences with an effect size >1.4 (among the smallest effect-sizes from our previous data). Samples sizes of 8-10 are powered to detect >95% of effects >1.8 (the median effect-size from the previous data). At the PND0 timepoint, we only generated 5 female and 4 male control pups and 7 dieldrin exposed pups for each sex from independent litters, limiting our statistical power at birth to 80% power to detect >80% of effects >1.8. We have reported this birth analysis for completeness, but it is important to acknowledge that the statistical power at birth was lower than for other time points for smaller effect sizes. Sample sizes are shown in Table 2.

The 12-week time point was selected based on previous results demonstrating increased neuronal susceptibility to MPTP and α-synuclein pre-formed fibrils (PFFs) administered at 12 weeks of age (Richardson et al., 2006,Gezer et al., 2020). The 9-month time point is equivalent to the 6-month post-PFF injection time point where we observed dieldrin-induced exacerbation of PFF-induced deficits in motor behavior and DA handling (Gezer et al., 2020). Birth and six weeks were selected to address the question of whether dieldrin-associated changes reflect a change in establishment or maintenance of epigenetic marks across time.

Sex determination genotyping at birth: At the birth time point, sex determination was not possible by visual inspection and was determined using PCR for Rmb31x/y gametologs using an established primer pair (Tunster, 2017). DNA was isolated from tail clips taken from neonatal mice during sacrifice using the Kapa Express Extract Kit (Kapa Biosystems, Cat. #KK7100) with one minor modification – sample lysis was assisted via physical homogenization using a 1.5 mL tube plastic pestle. After lysis, samples were briefly centrifuged to pellet cellular debris, and the DNA extract was diluted 10-fold with 10 mM Tris-HCl (pH 8.0–8.5) to dilute cellular debris and digested proteins to prevent inhibition of downstream PCR. PCR was performed using 1 µL of input DNA as described (Richardson et al., 2006,Gezer et al., 2020,Kochmanski et al., 2019,Boyd et al., 2023). PCR products were visualized on a 1.2% Agarose TBE gel with ethidium bromide. During visualization, two bands of DNA were produced for male samples, and a single band of DNA was produced for female samples, allowing for accurate determination of sex at birth.

Mass spectrometry: Dieldrin levels were measured from frozen neonatal brain, dam brain and dam adipose tissue by the RTSF Mass Spectrometry and Metabolomics core at MSU using established methods (Hong et al., 2004,US EPA, 2014).

SN microdissections: Frozen brains were mounted on a freezing cryostat (Leica, Model CM3050S) and sliced to the midbrain. Unilateral substantia nigra (SN) punches were collected using a chilled 1.0 mm micropunch and immediately placed in a frozen 1.5 mL tube on dry ice.

DNA isolation: DNA was isolated from unilateral SN punches using Qiagen QIAamp DNA Micro Kits (Qiagen, Cat. #56304) according to the included protocols with the following minor modifications. First, prior to adding ATL buffer to SN punches, 80 µL PBS was added to each sample and a 1.5 mL tube pestle was used to break up the sample. After that, 100 µL ATL buffer was added to each homogenized sample, which were allowed to equilibrate to room temperature. Second, for the proteinase K digestion step, the 56°C incubation time was extended to overnight. Third, carrier RNA was added to Buffer AL (this is an optional step). Fourth, the incubation time after addition of 100% EtOH was increased to 10 minutes, and the incubation time for the elution step was increased to 5 minutes. Lastly, to increase yield, the elution step was repeated by reapplying elution buffer to QIAamp MinElute column. DNA was eluted in 54 µL of 10 mM Tris-HCl, pH 8.0. DNA yield and purity were determined using a Qubit 3 fluorometer (ThermoFisher) and a NanoDrop spectrophotometer (ThermoFisher). Isolated DNA was stored at -80°C prior to sequencing.

Selection of candidate regions: To select candidate regions, we targeted all intragenic genomic features and known enhancers at genes annotated to DMCs or DMRs in our previous study (Kochmanski et al., 2019). Targeted regions include 255 male-specific annotations and 1043 female-specific annotations and represent a total of ~11.4 Mb of genomic space (Supplementary File 1). We included regions identified for both male- and female-specific changes to determine if these are truly sex-specific. Intergenic regions annotated to our candidate genes account for 151 Mb of genomic space, making inclusion cost-prohibitive with the SeqCapEpi platform. In addition, these intergenic regions remain largely unexplored, making data generated from these regions of limited interpretability. SeqCapEpi custom bait probes for regions of interest were designed using the Roche NimbleDesign software (Supplementary File 2). Baits covered 91.6% of target bases for a total capture space of ~10.4 Mb.

Capture hybridization sequencing: Capture hybridization-sequencing libraries were prepared by the Van Andel Genomics Core from 100 ng of high molecular weight DNA using the Accel-NGS Methyl-Seq DNA Library kit (v3.0) (Swift Biosciences, Cat. #30024). DNA was sheared following manufacturer’s protocol to an average size of 250bp, and sheared DNA was bisulfite converted using the EZ DNA Methylation-Gold kit (Zymo Research, Cat. #D5005) with an elution volume of 15 µL. Following adapter ligation, 8 cycles of library amplification were performed. Amplified libraries were pooled in batches of 8 and targeted enrichment of a custom 11Mb region was performed using Roche SeqCapEpi developer probes and SeqCapHyperEpi workflow starting at step 4.0 with the following modification: the capture was performed using IDT xGen Universal blockers to replace the SeqCap HE Universal Oligo and SeqCap HE Index Oligo. The post-capture amplification was also adjusted to 11 cycles of amplification and the final extension changed from 30 seconds to 1 minute. Quality and quantity of the finished library pools were assessed using a combination of Agilent DNA High Sensitivity chip (Agilent Technologies, Inc.) and QuantiFluor® dsDNA System (Promega Corp.). Sequencing (100 bp, paired end) was performed on an Illumina NovaSeq6000 sequencer using an S4, 200 bp sequencing kit (Illumina Inc.) with 10% PhiX included to improve base diversity.  Base calling was done by Illumina Real Time Analysis 3 and output of NextSeq Control Software was demultiplexed and converted to FastQ format with the Illumina Bcl2fastq software (version 1.9.0).

Data processing: FastQ files were processed using a slightly modified form of the bioinformatics pipeline previously established by our group to analyze reduced representation bisulfite sequencing (RRBS) data (Richardson et al., 2006,Gezer et al., 2020,Kochmanski et al., 2019,Boyd et al., 2023). Command line tools and the open-source statistical software R (version 4.1.2) were used for all analyses. For all sequencing data, the FastQC tool (version 0.11.7) was used for data quality control, and the trim_galore tool (version 0.4.5) was used for adapter trimming (Andrews, 2016,Krueger, 2017). During adapter trimming, we used the default minimum quality score and added a stringency value of 6, thereby requiring a minimum overlap of 6 bp. Trimmed CapHyb-seq reads were aligned to the mm10 reference genome using bismark (version 0.19.1) (Krueger and Andrews, 2011). Methylation data were extracted from the aligned reads in bismark using a minimum threshold of 5 reads to include a CpG site in analysis.

Verification of sequencing coverage: To determine the sequencing coverage across the targeted regions, we used BedTools Basic Protocol 3 for measuring coverage in targeted DNA sequencing experiments as described (Quinlan, 2014). BAM files for each sample obtained from bismark were compared with targeted bases (Supplementary File 2) using the coverage and multicov functions from bedtools (version 2.31.0) to determine the fraction of target bases that were captured and their sequencing depths (Quinlan and Hall, 2010). In differential modification analysis, only sites with coverage ≥ 5 in all samples were included in analysis. Of 6604 target regions, 1359 remained (21%), for ~2.8 Mb of captured sequence.

Differential modification analysis: The DSS (version 2.48.0) and DMRcate (version 2.14.1) R packages were used to test CapHyb-seq data for differential methylation (Feng et al., 2014,Peters et al., 2015). Given that dieldrin exposure has shown sex-specific effects on the dopaminergic system, all differential modification models were stratified by sex (Kochmanski et al., 2019,Boyd et al., 2023,Richardson et al., 2006). All pups included in modeling were from independent litters. Of note, we refer to DNA modifications, rather than DNA methylation, since our BS-based method does not differentiate between DNA methylation and hydroxymethylation. While DNA hydroxymethylation plays a critical role in gene expression in the brain and is particularly sensitive to environmental factors, methods for differentiating these marks remain limited and cost-prohibitive, especially across multiple time points (Kochmanski and Bernstein, 2020).  

Cross-sectional differential DNA modification analysis: To test for differentially modified cytosines (DMCs) by dieldrin exposure at each cross-sectional timepoint, we used the DMLtest function in DSS to perform two-group Wald tests. For DMLtest modeling, the equal dispersion parameter was set to FALSE and smoothing was set to TRUE. DMCs were considered significant at false discovery rate (FDR) <0.05. To test for differentially modified regions (DMRs) at each timepoint, we combined outputs from the callDMR function in DSS and the dmrcate function in DMRcate. For callDMR modeling, the p-value threshold was set to 0.05, minimum length was set to 50 base pairs, and minimum CpGs was set to 3. Meanwhile, for dmrcate modeling, the lambda value was set to 500, the C value was set to 4, and minimum CpGs was set to 3. DMR significance for the dmrcate output was set to a minimum smoothed FDR < 0.05.

Longitudinal differential DNA modification analysis: To test for DMCs by age, as well as simultaneous age and exposure in a multivariate model, we used the DMLfit.multiFactor function in DSS to perform linear models using a general experimental design. Given that age was an ordered variable, we coded each age group as a number for modeling – 6 weeks old = 1, 12 weeks old = 2, and 36 weeks old = 3. The birth timepoint was not included in longitudinal models because the F1 offspring came from a separate cohort of exposed animals, meaning they were not matched littermates like the later three time points. In the model for age only, coded age was included as the only independent variable. Meanwhile, in the multivariate model, exposure was included as a two-group categorical variable (“dieldrin”, “control”), and age:exposure was included as an interaction term. For the age alone and age:exposure models, DMCs were considered significant at false discovery rate (FDR) <0.05.

Genomic annotation: After differential methylation testing, the annotatr R package (version 1.26.0) was used to annotate identified DMCs and DMRs to the reference mm10 genome (Cavalcante and Sartor, 2017). Within annotatr, the annotate_regions function was used to generate CpG context, gene body, and regulatory feature annotations. Annotations for miRNA (miRbase), ENCODE predicted mouse midbrain enhancers (Accession: ENCSR114ZIJ), and custom full stack ChromHMM chromatin states for mm10 databases were added to the annotation cache in annotatr (Vu and Ernst, 2023,Kozomara et al., 2019,Luo et al., 2020,Kundaje et al., 2012).

Data visualization: Raw CapHyb-seq beta values were extracted using the bsseq R package (1.36.0). Volcano plots were generated using the ggplot2 R package (version 3.4.4). Euler diagrams were generated with the eulerr R package (version 7.0.1) (Larsson, 2024). UpSet plots were generated using the UpSetR package (version 1.4.0) (Conway et al., 2017). The ComplexHeatmap R package (version 2.18.0) was used to visualize beta value differences over time for all DMCs significant at one or more of the four time points. Specific genes of interest were visualized using the WASHU Epigenome Browser to determine genomic location and context of differentially modified (DM) regions; chromHMM annotations were loaded to compare DM loci with chromatin state annotations (Conway et al., 2017,Li et al., 2022). The R packages ggplot2, ggpubr (version 0.6.0), and ggeasy (version 0.1.4) were used to generate violin plots to display raw beta values for premature aging DMCs (Carroll et al., 2023,Kassambara, 2023).  

Pathway and Network Analysis: Gene ontology (GO) term enrichment testing and pathway analysis was performed on genes annotated to male and female DMCs and DMRs using the ClueGO application in Cytoscape (version 3.10.1) (Bindea et al., 2009,Smoot et al., 2011).

For ClueGO testing, genes annotated to DMCs and DMRs stratified by sex and timepoint were input as separate gene lists, “groups” was selected as the visual style, and the GO- biological process (GOBP) term was included for enrichment testing. Network specificity was set to “Medium”, such that the GO Tree Interval minimum was equal to 3 and the maximum was equal to 8. Only terms with at least 3 genes and a Bonferroni-corrected p-value < 0.05 were included in pathway visualizations. The connectivity score (Kappa) was set at 0.4, and default GO Term Grouping settings were used in all analyses. The genes found in enriched GOBP terms by ClueGO were used for STRING network analysis and to create UpSet plots for each sex displaying the overlap at different time points.

Protein-protein interaction network analysis was performed using STRING (version 12.0) with the genes found in enriched GOBP terms by ClueGO as input (Szklarczyk et al., 2017,Szklarczyk et al., 2015). STRING network analysis was performed using default parameters, including a minimum required interaction score = 0.4 and all interaction sources activated.