Introduction

Exercise has powerful effects on neurobiological processes such as neurotransmission, cognition, and reward-dependent behaviors (Rhodes et al. 2005; Caetano-Anolles et al. 2016; Saul et al. 2017). While motivation and physiological ability are critical determinants of an individual’s physical activity levels (Garland et al. 2017), long-term genetic selection experiments in laboratory mice have also illuminated genes associated with elevated levels of voluntary exercise (Kelly et al. 2012; Kostrzewa and Kas 2014; Vellers et al. 2018; Aasdahl et al. 2021). RNA (Zhang et al. 2018) and whole genome (Xu and Garland 2017; Hillis et al. 2020) sequencing from mice selected for elevated physical activity levels have revealed differential expression of genes associated with behavior, motivation, and athletic ability, supporting a genetic basis for elevated physical activity. In addition to genetic underpinnings, non-genomic, developmental programming effects during early life can also substantially impact exercise behavior and physical activity (Li et al. 2013; Baker et al. 2015; Eclarinal et al. 2016; Garland et al. 2017).

One such developmental programming mechanism involves genomic imprinting. Unlike standard Mendelian inheritance patterns where genes are biallelically expressed, genomically imprinted genes display parent-of-origin and monoallelic expression (Bartolomei and Ferguson-Smith 2011; Barlow and Bartolomei 2014). Approximately 1% of all mammalian genes are thought to be imprinted, with about 150 imprinted genes identified to date (Barlow and Bartolomei 2014; Kalish et al. 2014). Imprinted genes are critical for growth and development (Bouschet et al. 2017; Monk et al. 2019), cognition (Isles and Wilkinson 2000; Franklin and Mansuy 2010; Zamarbide et al. 2018), and exercise (Kelly et al. 2010). Parent-of-origin effects occur through several molecular mechanisms. The most well studied mechanism is differential DNA methylation (Bartolomei and Ferguson-Smith 2011; Barlow and Bartolomei 2014). The parental origin frequently determines the methylation status, with the paternal or maternal genomes typically exerting counteracting influences on gene expression (Azzi et al. 2014; Barlow and Bartolomei 2014; Kalish et al. 2014). For example, paternally expressed genes (PEG) such as Igf2 and Mest generally enhance growth of the fetus, while maternally expressed genes (MEG) such as H19 and Igf2r inhibit fetal growth (DeChiara et al. 1990; Leighton et al. 1995; Thorvaldsen et al. 1998). Moreover, differentially methylated regions exhibit variable methylation patterns following exposure to various early-life conditions such as diet (Waterland et al. 2006; Kovacheva et al. 2007; Heijmans et al. 2008; Tobi et al. 2009; Hoyo et al. 2011) and chemical exposure (Robles-Matos et al. 2021). Therefore, the methylation profile of imprinted genes is an appealing target to investigate in the context of exercise behavior because the expression or silencing of these monoallelic genes is critically involved in growth regulation and energy balance during development (Charalambous et al. 2007; Tunster et al. 2013; Waterland 2014; Millership et al. 2019) which can influence brain function (Perez et al. 2016; Kravitz and Gregg 2019).

The developing brain is a known hotspot for allele-specific gene expression (Perez et al. 2016; Huang et al. 2017, 2018; Kravitz and Gregg 2019). Imprinted genes in the brain encode various proteins and non-coding RNAs, from ubiquitination-related proteins to neurotransmitter receptor subunits. Similar to their counteracting effects on body size, imprinted genes have reciprocal effects on brain size, with MEGs enhancing and PEGs reducing brain size, supporting a role for imprinted genes in neurodevelopment (Keverne et al. 1996, Bouschet et al. 2017). The most substantial evidence for a neurodevelopmental role in humans emerged when Angelman syndrome (AS) and Prader-Willi syndrome (PWS) were linked to maternally and paternally transmitted mutations, respectively, on human 15q11-13. In the case of AS, UBE3A is a MEG located among a cluster of imprinted genes on chromosome 15 (Albrecht et al. 1997; LaSalle et al. 2015), and maternally inherited mutations or deletions in UBE3A lead to AS (Cassidy et al. 2000). UBE3A imprinting is regulated by the expression of an antisense transcript, UBE3A-ATS, which arises from a nearby imprinting control region that is methylated only on the maternal allele. The UBE3A-ATS transcript partially overlaps the UBE3A gene, blocking UBE3A expression from the paternal allele (Albrecht et al. 1997; Yamasaki et al. 2003; Chamberlain and Lalande 2010). Reciprocal mutations similar to those of AS cause PWS due to the lack of PEG products in the 15q11-13 chromosomal region (Cassidy et al. 2000). A neurodevelopmental role for genomic imprinting is supported further by the widespread expression of imprinted genes in the mouse brain and by the extensive neural and behavioral phenotypes of mutants of these imprinted genes (Babak et al. 2008; Wang et al. 2008; Gregg et al. 2010; DeVeale et al. 2012; Bouschet et al. 2017) including Cdkn1c (Imaizumi et al. 2020; Laukoter et al. 2020), Gnas (Mouallem et al. 2008; Turan and Bastepe 2013), Sgce (Zimprich et al. 2001; Asmus et al. 2002; Peall et al. 2013), and Trappc9 (Mochida et al. 2009; Babak et al. 2015).

Given the prominent roles of monoallelically expressed genes in growth, development, and energy metabolism in the brain, this study aimed to distinguish DNA methylation from genetic contributions for elevated physical activity in brains of mice that have been genetically selected for increased wheel-running behavior. In a long-term genetic selection experiment using outbred Hsd:ICR mice, four replicate high runner (HR) lines were bred for high voluntary wheel-running behavior and four control (C) lines were bred without regard to their wheel-running (C mice; Swallow et al. 1998, 1999). Compared with C lines, HR lines run ~ threefold more daily revolutions on running wheels. In addition, HR mice also demonstrate physiological and biochemical adaptations that are consistent with increased physical activity, including smaller body size, decreased fat composition, altered lipid profiles, and increased aerobic capacity (Swallow et al. 1999; Malisch et al. 2007; Acosta et al. 2015). More recently, early-life studies on HR mice demonstrated differential responses to early-life factors of altered juvenile diet and exercise opportunity (Acosta et al. 2015; Hiramatsu et al. 2017; Cadney et al. 2021a; McNamara et al. 2021).

In the current experiment, four groups were created wherein the offspring of non-selected C mice were fostered at birth and were reared by either a C or HR dam. Similarly, the offspring of selected HR mice were fostered at birth such that they were raised either by a C or HR dam (Cadney et al. 2021b). Using this cross-fostering approach, the current study sought to address two questions. First, do mice from a line genetically selected for elevated voluntary physical activity have altered DNA methylation profiles of imprinted genes compared to non-selected C mice? Second, does maternal upbringing further modify the DNA methylation status of these imprinted genes? Cross-fostering experiments are a powerful method to determine maternal effects and possible gene × environment interactions (Fish et al. 2004; Weaver et al. 2004; Kessler et al. 2011; Cohen et al. 2015; McCarty 2017). Therefore, to address these questions, we used cross-fostered HR and C male and female offspring from generation 90 of the HR selection experiment (see Cadney et al. 2021b for further methodological details concerning focal mice and their birth and foster parents). Using bisulfite sequencing of 16 known imprinted genes, we demonstrate altered DNA methylation patterns in the cortex and hippocampus for imprinted genes with known roles in fetal growth and energy metabolism. Our results show for the first time that, in addition to genetic underpinnings, differential methylation patterns of imprinted genes may also contribute to increased wheel-running behavior.

Methods

Experimental mice

The artificial selection experiment began in 1993 with a population of 224 mice from the outbred Hsd:ICR strain, which was randomly mated for two generations before being randomly partitioned into eight lines. Four replicate high runner (HR) lines of house mice were bred in the ongoing selection experiment for high voluntary wheel-running, based on average wheel revolutions per day on days five and six of a six-day running period as young adults. Four additional lines were bred randomly as Control (C) lines to the four HR lines (Swallow et al. 1998; Cadney et al. 2021b). The previous experiment used a subset of virgin male and female mice from generation 89 to produce experimental mice of generation 90 in the prior study. All mice were fed standard mouse chow (Teklad Rodent Diet W-8604) and regular drinking water. Pregnant dams were given a breeder diet (Teklad S-2235 Mouse Breeder Sterilizable Diet 7004) through weaning. All experiments involving animal handling and care were approved by the University of California, Riverside IACUC.

For the previous study, one C line (Line 4) and one HR line (Line 7) were used because they represented extremes in body mass among their respective linetypes (Table 2 in Cadney et al. 2021b). As described, it was hypothesized that differences in dam body size would lead to differential cross-fostering effects, including wheel-running behavior of the pups at weaning. Thus, using lines with different body masses could serve as a positive control for offspring body masses at weaning. The wheel-running behavior of these lines was representative of their respective linetypes (Cadney et al. 2021b). The dams used in the previous study were typical of other line 4 and 7 breeders of the same generation in terms of both wheel-running and body mass (Cadney et al. 2021b).

Experimental design

The experimental design has been previously published (Cadney et al. 2021b). Briefly, mice from generation 89 were sampled randomly to create a total of 60 C line 4 and 60 HR line 7 mating pairs, with the constraint that the sample size per foster group was equal (Table 1; Cadney et al. 2021b). A large number of pairings were made because only pups born on the same day could be used for fostering (Cadney et al. 2021b).

Table 1 Litter parameters for bisulfite sequencing analyses

At birth, litters were standardized to eight pups to avoid litter effects. As sex could not be determined at birth, the litter sex ratio could not be controlled. Cross-fostering only occurred between litters born within 24 h of one another. During the 48 h after cross-fostering, fostered pups were checked regularly, and none were rejected by their foster mother.

As births occurred, entire litters were fostered to another dam (no pup was returned to its biological mother). Thus, the previous study did not include a “control” group for the effects of fostering per se. This design was chosen to maximize the sample size in experimental groups sufficient to address their specific hypotheses (i.e., the effects of cross-fostering HR and C mice), given logistical constraints on the total sample size. In the current study, we wanted to determine whether rearing by an HR dam might be necessary for some proportion of high-running variance in DNA methylation profiles. This factorial experimental design produced four groups (Fig. 1 and Cadney et al. 2021b): C offspring raised by cross-fostered C dams (CC), C offspring raised by cross-fostered HR dams (CHR), HR offspring raised by cross-fostered C dams (HRC), and HR offspring raised by cross-fostered HR dams (HRHR). The number of male and female mice from each group used in this study is listed in Table 1.

Fig. 1
figure 1

Experimental design. Four experimental groups were created by cross-fostering between families of C (line 4) and HR (line 7) mice

Mice were removed from wheel access for one day, then sacrificed at seven weeks of age via decapitation without anesthesia. After whole brains were extracted, the hippocampus and cortex were dissected and snap-frozen on dry ice. Samples were stored at −80 °C until further sequencing steps were performed. Because of the reduced number of brain samples (N = 46), separate sex models could not be analyzed.

In-silico assay design

Gene sequences were acquired from the Ensembl genome browser and annotated. The assay target sequences were then re-evaluated against the UCSC mouse GRCm38 genome browser for repeat sequences, including LINE, SINE, LTR elements, and other DNA repeats. Sequences containing repetitive elements, low sequence complexity, high thymidine content, and overall CpG density were excluded from the in-silico design process. Forty-two assays were designed to cover 293 CpG sites across 16 genes, and the percentage methylation of each CpG site was determined in each sample. A list of the genes, their respective coordinates, and the number of CpG sites analyzed in this study are provided in Table 2.

Table 2 Coordinates, genomic context, and number of CpG sites analyzed for 16 imprinted genes analyzed by bisulfite sequencing

Bisulfite sequencing and data analysis

Bisulfite sequencing was performed on 92 samples (46 hippocampus and 46 cortex) by EpigenDx, Inc. (Hopkinton, MA). Tissue samples were digested using 500 μL of ZymoResearch M-digestion Buffer (Zymo, Irvine, CA) and 5–10 μL of protease K (20 mg/mL) with a final M-digestion concentration of 2X. The samples were incubated at 65 °C for a minimum of 2 h. 20 µL of the supernatant from the sample extracts were bisulfite modified using the ZymoResearch EZ-96 DNA Methylation-Direct Kit™ (cat# D5023) kit per the manufacturer’s protocol with minor modification. The bisulfite-modified DNA samples were eluted using M-elution buffer in 46 µL.

All bisulfite-modified DNA samples were amplified using separate multiplex or simplex PCRs with Qiagen (Gaithersburg, MD) HotStar Taq. All PCR products were verified and quantified using the QIAxcel Advanced System. Prior to library preparation, PCR products from the same sample were pooled and purified using QIAquick PCR Purification Kit columns (Qiagen).

Libraries were prepared using a custom Library Preparation method created by EpigenDx. Library molecules were then purified using Agencourt AMPure XP beads (Beckman Coulter) and quantified using the Qiagen QIAxcel Advanced System. Barcoded samples were then pooled in an equimolar fashion before template preparation and enrichment were performed on the Ion Chef™ system (Thermo Fisher) using Ion 520™ & Ion 530™ ExT Chef reagents. Following this, enriched, template-positive library molecules were sequenced on the Ion S5™ sequencer using an Ion 530™ sequencing chip (Thermo Fisher).

FASTQ files from the Ion Torrent S5 server were aligned to the local reference database using open-source Bismark Bisulfite Read Mapper with the Bowtie2 alignment algorithm (https://www.bioinformatics.babraham.ac.uk/projects/bismark/; Krueger and Andrews 2011). Methylation levels were calculated in Bismark by dividing the number of methylated reads by the total number of reads, considering all CpG sites covered by a minimum of 30 total reads. CpG sites with fewer than 30 reads were excluded from analyses.

Statistical analysis

Data were analyzed as mixed models in SAS 9.1.3 (SAS Institute, Cary, NC) Procedure Mixed, with REML estimation and Type III Tests of Fixed Effects. Line (selected line 7 vs. non-selected line 4), foster-line, and sex were fixed effects, while dam ID (n = 28) was a random effect nested within linetype × foster-line. Separate mixed models were created to analyze the percent methylation levels for CpG sites across the entire genomic region for each gene (Supplementary Table 1) as well as the percent methylation levels for CpG sites within distinct genomic regions (i.e., introns, exons, 5’ upstream regions, and 5’ untranslated regions; Supplementary Table 2).

In all tables, we present least squares means (L.S. mean) and standard errors (S.E.) for each imprinted gene for both brain regions (cortex and hippocampus) and the results of each F-statistic. Supplementary Tables 1 and 2 also present the p-values for the differences in the L.S. means between the in-fostered and cross-fostered groups (i.e., the effect of cross-fostering between the HR and C lines by sex). Supplemental material can also be referenced for main effects, as well as interactions (line × foster-line, line × sex, foster-line × sex, line × foster-line × sex).

In all analyses, outliers were iteratively removed when the standardized residuals exceeded ~ 3. Statistical significance was set at p < 0.05. Effect sizes, presented as Hedges’ g values, were also calculated for group comparisons (CC vs. HRHR, CC vs. CHR, and HRHR vs. HRC; Tables 3 and 4 and Supplemental Tables 3 and 4). Hedges’ g values greater than  + 0.8 and less than − 0.8 were viewed as a large effect size, while values between + 0.8 and + 0.5 and between − 0.8 and − 0.5 were considered to be a medium effect size.

Table 3 Type 3 tests of fixed effects for genes in the cortex with at least one significant main effect and/or interaction (Color table online)
Table 4 Type 3 tests of fixed effects for genes in the hippocampus with at least one significant main effect and/or interaction (Color table online)

Results

To address whether mice from a line genetically selected for elevated voluntary physical activity have altered DNA methylation profiles of imprinted genes compared with non-selected C mice, we investigated the methylation profile of C offspring that were cross-fostered at birth and reared by C dams and HR offspring that were cross-fostered at birth and raised by HR dams. We also investigated the methylation profiles of C offspring cross-fostered at birth and reared by HR dams and HR offspring cross-fostered at birth and raised by C dams to investigate additional early-life programming effects (i.e., maternal effects) on genomic imprinting (Fig. 1). This factorial experimental design allowed us to examine the genomic and non-genomic contributions to elevated physical activity. Because the expression of imprinted genes varies with tissue type and male and female mice differ in wheel-running behaviors, we examined the DNA methylation profiles in the hippocampus and cortex of both sexes. The percent methylation of 16 imprinted genes totaling to 293 CpG sites were analyzed by bisulfite sequencing and are reported in Supplemental Tables 1 and 2. The F-statistic and accompanying p-values for those that had p < 0.05 for a main effect or interaction are further reported in Tables 3 and 4 for the cortex and hippocampus, respectively, as well as Supplemental Tables 3 and 4. Hedges’ g value of effect sizes between key groups (CC vs. HRHR, CC vs. CHR, and HRHR vs. HRC) for these genes are also reported.

Individual main effects of line and maternal upbringing on DNA methylation

When analyzed across the entire genomic region, our statistical model revealed several significant individual main effects of line (physical activity) and foster-line (maternal upbringing) on the DNA methylation profiles of the 16 imprinted genes analyzed in this study (Supplemental Table 1). A significant main effect of line for Rasgrf1 in the cortex (p = 0.0120; Table 3) and Zdbf2 in the hippocampus (p = 0.0306; Table 4) was observed, as well as additional trends towards statistical significance (p < 0.10) for foster-line for Zdbf2, (p = 0.0727; Table 4) in the hippocampus.

When the CpG sites were analyzed within distinct genomic regions such as introns, exons, and promoter regions, no significant individual main effects of line and foster-line on the DNA methylation of the 16 imprinted genes were revealed (Supplemental Tables 3 and 4). These results suggest that differences in genomic regions do not significantly influence the overall methylation profile of any of the imprinted genes analyzed.

Interaction of line and maternal upbringing on DNA methylation

When assessed for line and foster-line interaction across the entire gene, there were significant main effects for Ig-DMR (p = 0.0045) and Mest (p = 0.0216) in the cortex (Table 3). The hippocampus had a significant main effect for line and foster-line for Peg3 (p = 0.0126; Table 4). When the CpG sites were analyzed by genomic context, there were additional line × foster-line interactions (Supplemental Tables 3 and 4). In the cortex, there was a significant line × foster-line interaction for the exon (p = 0.0063) and promoter (p = 0.0483) regions of Mest (Supplemental Tables 3). In the hippocampus, there was a significant line × foster-line interaction for the intron region of Peg3 (p = 0.0062; Supplemental Tables 4).

Influence of sex and the interaction of sex, line, and maternal upbringing on DNA methylation

When analyzed across the entire genomic region, sex produced an individual main effect on Snrpn in the hippocampus (p = 0.0461; Table 4). A trend for a significant main effect of sex for Rasgrf1 (p = 0.0857) and Ig-DMR (p = 0.0746) in the cortex (Table 3) and Igf2 (p = 0.0671) in the hippocampus (Table 4) were also observed. When assessed for a three-way line, foster-line, and sex interaction, there were significant main effects for Rasgrf1 (p = 0.0397) in the cortex (Table 3) and Igf2 (p = 0.0212) and Impact (p = 0.0175) in the hippocampus (Table 4). The Sgce intron region within the hippocampus also revealed an additional three-way interaction (p = 0.0168; Supplemental Tables 4).

Analyses of effect size between CC vs. CHR, CC vs. HRHR, and HRHR vs. HRC groups revealed large effect sizes with Hedges’ g values of greater than + 0.8 and/or less than − 0.8 for Rasgrf1, Ig-DMR, Mest, Gnas, Grb10, and Trappc9 in the cortex (Table 3 and Supplemental Tables 3) and Zdbf2, Peg3, Igf2, and Impact in the hippocampus (Table 4 and Supplemental Tables 4), indicating notable effects of line and foster-line, despite non-significant line × foster-line interactions for these genes.

Discussion

The DNA methylation status of genes is an attractive conduit for early-life effects that can persist into adulthood. Moreover, alterations in DNA methylation can interact with other early-life influences and amplify across successive generations (Yin et al. 2013; Short et al. 2017; Yeshurun et al. 2017; McGreevy et al. 2019). Here we analyzed DNA methylation patterns for 16 imprinted genes in a genetically selected line of mice (HR) and a non-selected line (C) whose offspring were cross-fostered in a full factorial experimental design. Our results shed light on non-genomic contributions to elevated physical activity. The results from our statistical model, which accounted for the influence of sex, line, and foster-line, showed that lines differed in DNA methylation patterns of Rasgrf1 in the cortex (Table 3) and Zdbf2 in the hippocampus (Table 4). These results support the hypothesis that the genetic inheritance of elevated physical activity can potentially lead to altered expression of imprinted genes in the brain. Furthermore, early-life effects through cross-fostering modified the methylation status of other imprinted genes involved in fetal growth and energy metabolism, including Ig-DMR, Mest, and Peg3. Finally, the sexes differed in DNA methylation for the paternally expressed genes Snrpn, Rasgrf1, Igf2, and Impact.

Below, we first discuss the significance of the imprinted genes, Rasgrf1 and Zdbf2, given that the methylation patterns of these genes were significantly affected by wheel-running. We then discuss the methylation status of additional imprinted genes modified with cross-fostering and by sex and the implications of these genes on neurodevelopment.

The HR line and genomic imprinting of Rasgrf1 and Zdbf2 in the brain

Our most notable finding was alterations in DNA methylation in the paternally imprinted genes, Rasgrf1 and Zdbf2, in the cortex and hippocampus, respectively, in the HR line (Tables 3 and 4). In mice, Rasgrf1 (Ras protein-specific guanine nucleotide releasing factor 1) is a paternally methylated and paternally expressed imprinted gene located on chromosome 9 (Plass et al. 1996; Yoon et al. 2002; Dockery et al. 2009). All 16 CpG sites assayed for this gene were located approximately 30 kb to 40 kb upstream of the transcription start site, where the Rasgrf1 differentially methylated region (DMR) is located (Shibata et al. 1998; Dockery et al. 2009). The methylation status of the Rasgrf1 DMR directly influences Rasgrf1 expression in the brain. The Rasgrf1 DMR contains CTCF binding sites, which have been shown to function as enhancer blockers at maternally imprinted loci, including H19, Igf2, and KvDMR1 (Bell and Felsenfeld 2000; Hark et al. 2000; Yoon et al. 2005; Holmes et al. 2006; Fitzpatrick et al. 2007). CTCF binds to the unmethylated Rasgrf1 DMR, which functions in cis to silence the maternal allele (Yoon et al. 2005; Fitzpatrick et al. 2007). Experimental conditions that prevent methylation of the paternal DMR result in silencing the normally expressed paternal allele in the brain, consistent with the model that the unmethylated DMR binds CTCF, blocking expression of the Rasgrf1 allele in cis (Yoon et al. 2002, 2005). Adding an enhancer between the DMR and Rasgrf1 overrides the enhancer-blocking activity of the CTCF-bound DMR, implicating enhancers located upstream of the DMR in the regulation of Rasgrf1 in the brain (Yoon et al. 2005; Holmes et al. 2006).

Rasgrf1 is detected mainly in the brain (Plass et al. 1996; Yoon et al. 2002; Holmes et al. 2006; Fernandez-Medarde and Santos 2011), and its expression changes from exclusively paternal before weaning to preferentially paternal after weaning (Drake et al. 2009). Mice with Rasgrf1 paternal deletions are normal at birth but display retarded growth by weaning (Itier et al. 1998), and loss of Rasgrf1 imprinting increases growth and weight after birth (Drake et al. 2009). Synchronization of growth in the developing brain relies on the secretion of growth hormone (GH)-releasing hormone and somatostatin, which stimulate and inhibit GH and IGF1 secretion, respectively. Levels of both GH and IGF1 decrease or increase when Rasgrf1 is deleted or overexpressed, respectively.

Rasgrf1 also mediates oxidative stress and inflammatory responses and can be upregulated with chronic exercise (Table 3; Moreland et al. 2020) and stress (Cheng et al. 2020). Rasgrf1-deficient mice display lower levels of reactive oxygen species, protection against oxidative stress, reduced body size, and improved motor coordination with age (Borras et al. 2011). It is possible that the alteration of Rasgrf1 methylation in HR mice results in a shift in oxidant and anti-oxidant factors and could contribute to the neurobiological alterations observed in these mice. Supporting this, Rasgrf1 is expressed in the postsynaptic densities of neurons and is an intrinsic mediator of brain derived neurotrophic factor-induced small GTPase R-Ras activation, R-Ras-mediated axonal morphological regulation (Umeda et al. 2019), and is involved in memory formation (Brambilla et al. 1997).

We also observed a significant main effect of selective breeding for elevated physical activity for 12 CpG sites located 10 kb upstream of the Zdbf2 transcription start site (Table 4). The Zdbf2 (DBF-type zinc finger-containing protein 2) gene encodes for a protein containing DBF4-type zinc finger domains but has no known function. Like Rasgrf1, Zdbf2 is also paternally methylated and paternally expressed (Kobayashi et al. 2009) and controls the paternal-specific expression of Zdbf2 isoforms during early development (Duffie et al. 2014). The Zdbf2 locus is paternally imprinted, with three intergenic paternally methylated DMRs located between 8.5 kb and 16 kb upstream of the Zdbf2 transcription start site. However, additional work suggests that the Zdbf2 DMR locus may also be under transient maternal imprinting control during specific developmental periods (Kobayashi et al. 2012; Proudhon et al. 2012; Duffie et al. 2014), likely due to multiple promoters (Duffie et al. 2014). Because the function of Zdbf2 is unknown, it is difficult to estimate the functional consequences of altered Zdbf2 methylation in the hippocampus of our HR mouse model. However, studies have shown that aberrant hypermethylation of the maternal Zdbf2-DMR has been observed in some patients with Beckwith-Wiedemann Syndrome (Maeda et al. 2014), an imprinting disorder characterized by body overgrowth and enlargement of internal organs. Investigating the functions of the Zdbf2 gene in exercise behavior may provide further information about how Zdbf2 expression may control fetal growth and metabolism.

Because the two genes that were significantly altered by selective breeding are paternally methylated and paternally expressed genes, our results support a hypothesis that paternal imprinting of genes in the brain may be a conduit through which exercise behavior influences nervous system development and function (Radford et al. 2011; Li et al. 2013; Baker et al. 2015; Eclarinal et al. 2016; Garland et al. 2017). In addition, given the well-characterized changes in allele frequencies for genes associated with elevated levels of exercise (Xu and Garland 2017; Hillis et al. 2020), it is also plausible that genetic and epigenetic factors interact with each other to modify DMR methylation profiles, as has been noted for human IGF2R (Xu et al. 1997; Sandovici et al. 2003) and IGF2 (Sandovici et al. 2003) genes. This possibility remains to be tested in future studies.

Influence of maternal upbringing on genomic imprinting in the brain

The early family life environment shapes human neurodevelopment and behavioral outcomes well into adulthood (Heim et al. 2002; Apter-Levy et al. 2013). In rodents, naturally occurring variations in maternal care shape the offspring’s neuroendocrine stress reactivity, anxiety, and depression-like behavior (Meaney 2001; Schroeder and Weller 2010; Kessler et al. 2011; Cohen et al. 2015). Dams exhibiting different social behaviors also demonstrate variations in maternal care (Clinton et al. 2007, 2010), and cross-fostering their offspring shifts their adult phenotype, leading to greater social interaction and reduced anxiety-like anxiety behavior. Anxiety-like behavior effects are accompanied by gene expression changes in the amygdala that emerge in early postnatal life and persist through adulthood, revealing how an early-life manipulation such as cross-fostering can change brain development and ultimately impact adult behavior.

Variations in maternal care can also alter the expression of imprinted genes. It has been suggested that paternally expressed imprinted genes increase offspring size and survival. Mutations affecting the paternally expressed Mest (Lefebvre et al. 1998) and Peg3 (Li et al. 1999; Curley et al. 2004) genes result in a striking lack of maternal care, supporting the notion that genomic imprinting is associated with the evolution of social behaviors. Targeted mutation of paternally expressed imprinted genes influences several aspects of fetal and postnatal development. Mutation of Peg3 in offspring reared with wild-type mothers, and vice versa, significantly increases offspring mortality, with the combined mutation in mother and offspring exhibiting near 100% lethality (Curley et al. 2004). Moreover, Peg3 mutant mothers exhibit decreased nurturing behaviors, including nest-building, crouching, and pup retrieval, as well as reduced numbers of oxytocin-producing neurons in the hypothalamus (Li et al. 1999). Inactivation of Mest in embryonic cells also resulted in abnormal maternal behavior directly impacting the feeding behavior of newborns (Lefebvre et al. 1998). These data suggest that paternally imprinted genes are critical for maternal care. Consistent with these data, we observed alterations in the methylation status of four paternally expressed genes, Mest, Peg3, Ig-DMR, and Impact, in cross-fostered HR mice (Tables 3 and 4). These genes have been shown to promote the growth of body size but inhibit the growth of the brain (Keverne et al. 1996, Bouschet et al. 2017). Given that these genes play critical roles in fetal growth and metabolism, these data support the hypothesis that DNA methylation may be a mechanism to explain the persistence of developmental programming influences on exercise behavior.

Although the observed changes in DNA methylation levels in the current study may directly influence wheel-running behavior, this hypothesis is challenging to assess. In Cadney et al. (2021b), HR pups raised by HR dams weighed less than C pups raised by C dams and continued to have reduced body masses as adults. As expected, adult HR mice also ran approximately threefold more than their C counterparts, and females ran more than males. However, cross-fostering had no statistical consequence on any aspect of wheel-running behavior, including running distance, duration, or maximum speed. In addition, with body mass as a covariate, HR mice had higher VO2max than C mice, and males had higher VO2max than females, but cross-fostering had no effect. Related to the brain, female HR pups raised by C dams had significantly larger brains than female HR pups raised by HR dams (Cadney et al. 2021b). Because cross-fostering had no effect on adult exercise levels, future studies focusing on the methylation profiles of other tissues relevant to exercise physiology (e.g., skeletal muscle, heart, liver) are not likely to reveal important effects of cross-fostering on voluntary exercise behavior in these lines of mice. However, potential epigenetic effects induced by other early-life environmental factors, such as altered diet and/or exercise (Li et al. 2013; Desai et al. 2014a, 2014b; Meek et al. 2014; Acosta et al. 2015; Hiramatsu et al. 2017; Cadney et al. 2021a, 2022), remain an important area of focus for additional study.

The methylation profiles at most of the imprinted gene DMRs analyzed in the current study had a baseline level of methylation close to the theoretical 50% expected for loci whose methylation status is differentially established on parental alleles. However, two genes, Igf2 and Ig-DMR, had almost 100% methylation in the cortex and hippocampus (Supplemental Table 1 and Tables 3 and 4). These unexpected findings prompted us to consider the identity and nature of such early-life effects that may be capable of shifting gene methylation profiles. Many developmental exposures have been demonstrated to shift the DMR methylation profile for Igf2, including maternal cigarette smoking (Murphy et al. 2012a), maternal stress (Vangeel et al. 2015), anti-depressant medications taken during pregnancy (Soubry et al. 2011), prenatal nutrition (Kovacheva et al. 2007; Heijmans et al. 2008; Tobi et al. 2009; Hoyo et al. 2011), postnatal nutrition (Waterland et al. 2006; Hoyo et al. 2011), and chemical exposure (Robles-Matos et al. 2021). Such deviations in DMR methylation are likely functionally significant, given that imprinted genes are critical for cellular differentiation (Zhang et al. 1997), prenatal and postnatal growth (Bouschet et al. 2017), neurobiological function (Isles and Wilkinson 2000; Perez et al. 2016; Huang et al. 2017; Kravitz and Gregg 2019), and maternal behaviors (Lefebvre et al. 1998; Li et al. 1999).

The imprinting of Igf2 (Gregg et al. 2010; Harper et al. 2014; Ye et al. 2015), Dlk1 (Croteau et al. 2005), Grb10 (Monk et al. 2009; Garfield et al. 2011), and Peg1/Mest (Shi et al. 2004) DMRs in the brain may also not be as strictly maintained in other tissues. Such epigenetic heterogeneity may be partially explained by the relaxation of epigenetic control for genes whose monoallelic expression is not essential for brain function. However, this is an unlikely interpretation because Igf2 and Ig-DMR are critical for neurodevelopment. An alternative and more probable interpretation is that variation between monoallelic and biallelic expression of imprinted genes in the brain may be directly related to brain function and may reflect a response to changes in environmental cues, such as the maternal environment. In one study, Igf2 was shown to be abundantly expressed in the hippocampus and prefrontal cortex of adult rats and, in contrast to peripheral tissues, was maternally expressed (Ye et al. 2015). Given the established role of imprinted genes in maternal behavior and neurobiological processes (Lefebvre et al. 1998; Li et al. 1999), it is reasonable to conjecture that increased plasticity of imprinting in the brain may contribute to complex and heritable behavioral phenotypes, including exercise behavior. Additional insights are further discussed below.

Sexual dimorphism in parental imprinting in the brain

Earlier studies show that when provided access to a running wheel, females run significantly more, for more extended periods, and at higher average speeds than males (Swallow et al. 1998). This was also seen when offspring of HR females ran significantly longer distances, spent more time running, and ran at higher maximum speeds than the offspring of HR males (Kelly et al. 2010). As a result, females were also smaller and had lower percentage of body fat (Swallow et al. 1999) and generally regulated energy balance more precisely than males (Swallow et al. 1999). In the current study, sex further modified DMR methylation levels for Rasgrf1 in the cortex (Table 3) and Snrpn, Igf2, and Impact in the hippocampus (Table 4), all of which are paternally imprinted. These findings were expected, given the numerous sex-specific effects observed in the previous cross-fostering study, with cross-fostering increasing brain mass of female, but not male, HR mice (Cadney et al. 2021b).

Several reasons could account for sex differences in DNA methylation. First, the timing and duration of methylation marks are different in the male versus the female germlines (Bourc'his and Proudhon 2008). At the time of fertilization, paternal imprinting control regions are methylated in spermatozoa versus unmethylated in the oocyte. In contrast, maternal imprinting control regions are methylated in the oocyte but not in spermatozoa. Second, the developmental stage in which gametic methylation imprints are acquired is also sex-specific, with de novo methylation being initiated earlier in the male germline compared to the female germline. Methylation is acquired very early during spermatogenesis, in precursors of the self-renewing spermatogonial stem cells. Because of this, paternally methylated CpGs endure more replication cycles in spermatogonial stem cells and proliferating spermatogonia and persist several weeks to several years before the production of mature spermatozoa (Eichenlaub-Ritter et al. 2007). On the contrary, maternal methylation patterns are established just before ovulation. The longer duration of paternal methylation patterns and the higher number of replication cycles during gametogenesis may favor the erasure of paternal ICRs, while the shorter longevity of methylation patterns in the female germline may be better conserved. Third, while the total number of paternally and maternally imprinted genes is approximately even, most of the primary methylation marks acquired in the germline are of maternal origin. More than 15 imprinted clusters are dependent on ICRs harboring maternal germline methylation, while only three imprinted loci are controlled by paternal germinal marks: H19/Igf2, Gtl2/Dlk1, and A19/Rasgrf1 (Reik and Walter 2001). The rest of the paternally imprinted genes are controlled through the production in cis of a non-coding RNA, silencing itself on the maternal allele by DNA methylation inherited from the oocyte. Lastly, the maternal genome is required for the development of the embryo, while the paternal genome promotes the development of extraembryonic structures such as the placenta (Barton et al. 1984; McGrath and Solter 1984). However, it is difficult to address the specific influence of maternal and paternal germline imprints on the neurodevelopment of the offspring. Regardless, it is reasonable to deduce that sex-specific observations such as those found in our study may be informative when designing future studies to evaluate links between parent-of-origin genes and sexually dimorphic effects on wheel-running behavior.

Limitations and concluding remarks

The current study provides new insights into potential epigenetic mechanisms underlying the developmental origins of exercise and physical activity. However, some limitations will need to be addressed in future studies. First, the cross-fostering experimental design did not include non-fostering controls such that C offspring were reared by their biological C dam, and HR offspring were reared by their biological HR dam. This limitation was intentional, as the inclusion of such controls was not necessarily relevant to the experimental aims of the previous study. Additionally, resource constraints and early pandemic restrictions on in-person research were also concerns. Consequently, the current experimental design was limited to 20 litters. Second, the DNA in this study was extracted from hippocampal and cortical tissue made up of a heterogeneous collection of cells. DNA methylation can vary by cell type, so this must be carefully considered in future studies (Davies et al. 2005; Kravitz and Gregg 2019). Cell-type-specific imprinting effects are especially intriguing because they could serve as potential markers for novel subpopulations of brain cells controlling a particular neural circuit or behavior pattern. Third, we did not assess the ability of DNA methylation status to modulate gene expression; therefore, gene expression studies are needed to assess gene activity. Lastly, bisulfite sequencing does not differentiate between 5-methylcytosine and 5-hydroxymethylcytosine. This is critical because 5-hydroxymethylcytosine is enriched in the brain and is regulated during development (Kinney et al. 2011; Sherwani and Khan 2015). In the future, it will be essential to determine whether alterations in DNA methylation associated with physical activity and other early-life programming factors are mediated by enzymes that add not only 5-methylcytosines, but also 5-hydroxymethylcytosines.

Although we have observed line and sex-specific changes in DNA methylation patterns of several imprinted genes, it is prudent to question the relevance and potential functional consequences of small percentage changes in DNA methylation like those measured in this study (Breton et al. 2017). Indeed, several studies have found that small changes in CpG methylation are correlated with differential gene expression (Murphy et al. 2012b; Kile et al. 2013; Argos et al. 2015; Maccani et al. 2015; Montrose et al. 2017). Therefore, it is possible that even small shifts in DNA methylation that produce large effect sizes can directly impact the transcription and expression of these genes.

Because epigenetic marks can “drift” over time, it will also be critical to assess offspring later in life (Kochmanski et al. 2017) to determine whether DNA methylation changes are maintained, possibly amplify with age, and/or transmitted across generations. The expression of imprinted genes can also be regulated by additional epigenetic mechanisms, such as histone modifications, alternative promoter usage, and post-transcriptional processes (Jouvenot et al. 1999; Landers et al. 2004; Kim et al. 2018). Investigating the contribution of additional epigenetic mechanisms in the regulation of wheel-running behavior is also warranted.

In conclusion, using a cross-fostering paradigm to investigate non-genomic contributions to increased physical activity, we show that a line of mice selectively bred for high voluntary wheel-running behavior has stable alterations in DNA methylation levels of imprinted genes with known functions for fetal growth, development, and metabolism. These differential methylation patterns were also shaped by maternal upbringing and sex, supporting the hypothesis that genomic imprinting in the brain can contribute to complex and heritable behavioral phenotypes, and is further influenced by developmental programming factors during early life.