Examining age-dependent DNA methylation patterns and gene expression in the male and female mouse hippocampus.

DNA methylation is a well-characterized epigenetic modification involved in numerous molecular and cellular functions. Methylation patterns have also been associated with aging mechanisms. However, how DNA methylation patterns change within key brain regions involved in memory formation in an age- and sex-specific manner remains unclear. Here, we performed reduced representation bisulfite sequencing (RRBS) from mouse dorsal hippocampus - which is necessary for the formation and consolidation of specific types of memories - in young and aging mice of both sexes. Overall, our findings demonstrate that methylation levels within the dorsal hippocampus are divergent between sexes during aging in genomic features correlating to mRNA functionality, transcription factor binding sites, and gene regulatory elements. These results define age-related changes in the methylome across genomic features and build a foundation for investigating potential target genes regulated by DNA methylation in an age- and sex-specific manner.

Introduction

DNA methylation at the 5-position carbon of cytosine is a well-characterized heritable epigenetic mark that exists widely in mammals (Feng et al., 2010). Within the genome, methylation has been shown to be a critical regulator of biological function, and plays important roles in development, regulation of gene expression, suppression of repetitive element transcription and transposition (Taylor et al., 2019). Altered or aberrant DNA methylation patterns are hallmarks associated with a multitude of diseases, with implications in cancer, epigenome-wide association studies (EWAS), and cognitive variability in aging through dysregulation of methylation at specific genomic features (Mahmood and Rabbani, 2019; Zaghlool et al., 2020; Zhang et al., 2019; Maierhofer et al., 2017). Epigenetic age-prediction using methylated CpG sites identified from different tissues in mice, such as liver or brain, have shown that epigenetic age correlates are largely tissue-specific and do not share the same DNA methylation patterns across the genome, with varying correlation rates for age estimation depending on the origin of tissue sampled (Petkovich et al., 2017; Thompson et al., 2018; Wang et al., 2017; Field et al., 2018). As memory related markers are vital for understanding cognition in aging, and as epigenetic mechanisms such as histone modifications (Barrett el al., 2011; Kim and Kaang, 2017) are increasingly being shown as key modulators of memory formation, examining DNA methylation marks using tissue from a brain region directly involved in memory formation, such as the hippocampus, can provide more clear insight into understanding age-related methylation signatures.

Existing studies within mouse hippocampus investigated how patterns of DNA methylation found in the hippocampus compare to other tissues (i.e. blood, cortex), and how these patterns may be influenced by early-life environmental factors (Coninx et al., 2020; Harris et al., 2020; Zhang et al., 2018). However, the extent to which hippocampal DNA methylation can be used as a marker for cognitive variability in memory throughout the lifespan remains unknown. The study by Zhang et al., 2018 compared the DNA methylation signatures of the dorsal and ventral hippocampus as epigenetic correlates to experience-induced plasticity from environmental enrichment in young adult mice, but not aging mice. Both transcriptional and methylation differences were observed between the 2 subregions of the hippocampus despite similar rearing of the mice, demonstrating extreme specificity of gene expression and methylation across different hippocampal subregions. Furthermore, a study examining methylation at young and aging time points focused on whole hippocampal tissue rather than any specific subregion (Masser et al., 2017). Importantly, some site-specific changes in methylation were observed across age as well as lifelong sex differences, but limited regional annotation was performed to investigate whether methylation changes in aging are more prevalent or recurring within certain genomic features. To account for the lack of specificity in methylation dynamics across tissues and the genomic landscape, we focused our analyses specifically within the dorsal subregion of the hippocampus because it is necessary for the formation and consolidation of specific types of memories, and its function declines with age.

In the present study, we set out to achieve three goals: (1) examine mouse region-specific hippocampus methylation patterns across global genomic features; (2) characterize the regulatory context of sex-specific differentially methylated patterns; and (3) identify putative gene targets whose expression is impacted by age- and sex-associated changes in methylation. We performed reduced representation bisulfite sequencing (RRBS) from 2- and 20-month-old mouse dorsal hippocampus. In addition to thoroughly characterizing baseline methylation patterns that are age-and sex-specific, we also compare differentially methylated regions to previously published gene expression profiles from the mouse dorsal hippocampus. Together, these analyses reveal wide-spread age-associated hypomethylation, increased age-dependent methylation variability, sex-specific divergence in differentially methylated regions at specific genomic features, and putative regulatory mechanisms for age-associated differences in gene expression.

Materials and methods

Mice

Young adult (2 months old) male and female C57BL/6J mice were obtained from Jackson Laboratory and aged (20 months old) male and female C57BL/6 mice were obtained from the National Institute on Aging Aged Rodent colony. Mice were given ad libitum access to food and water and were group-housed while maintained on a 12h light/dark cycle for a minimum of one week for environment acclimation. Mice were kept group-housed in their homecage for at least 1 week before the time of tissue collection and were not subjected to any external manipulations in order to preserve baseline genomic patterns. All tissue collection experiments occurred during the light cycle and were approved by the Institutional Animal Care and Use Committee at the University of California Irvine.

Tissue collection

Mice were anesthetized through intraperitoneal (IP) injection of sodium pentobarbital (50mg/mL) and perfused transcardially with ice-cold PBS, pH 7.4, using a peristaltic perfusion pump (Thermo Fisher Scientific, Waltham, MA). After whole brain extraction, the dorsal hippocampus was bilaterally dissected and flash frozen immediately on dry ice. The ventral hippocampus, prefrontal cortex, hypothalamus, and blood from the right atrium of the heart were collected and stored in an identical manner. All perfusions and tissue collection were performed within the same 3-hour window to avoid circadian effects.

Sample processing

Bilateral dorsal hippocampal tissue was homogenized with a genomic lysis buffer for two minutes using a TissueLyser kit (Qiagen, Germantown, MD) and subsequently processed for DNA isolation according to the manufacturer’s instructions (Quick-DNA Miniprep Kit, Zymo Research, Irvine, CA). Quality of the samples was recorded using the Nanodrop 1000, where all samples held 260/280 ratios above 1.9. Genomic DNA at a concentration of 4.5ug/sample was digested overnight at 37 °C with the methylation-insensitive restriction endonuclease Msp1 (20U/uL, New England Biolabs, Ipswich, MA). Following digestion, samples were bead purified (AMPure XP, Beckman Coulter, Brea, CA) and size selected for 100–250bp using a double-sided size selection protocol. Size selected samples were washed with 85% ethanol and checked for DNA quality. Samples were then bisulfite converted using the EZ DNA-Methylation Gold kit also from Zymo Research.

Library construction

Library construction was completed using the Accel-NGS Methyl-Seq DNA Library kit (Swift Biosciences, Ann Arbor, MI). Quantities of bisulfite-converted DNA ranging between 50ng and 100ng were used as input. Multiple purification and size selection steps were performed after extension and ligation using ratios taken from the manufacturer’s protocol to retain fragments with a minimum size of 100bp. Size selection was followed by the addition of unique dual indices (Swift Biosciences) to each sample. Samples underwent five cycles of amplification and resulting PCR products were purified and size selected for optimal fragment length. Completed library concentration was recorded from a dsDNA High Sensitivity Assay on Qubit and run on a 2200 TapeStation (Agilent, Santa Clara, CA) for fragmented library size visualization and peak concentration.

Reduved representation bisulfite sequencing analysis

Library samples were sequenced on patterned flow cells (No-vaSeq 6000, Illumina, San Diego, CA). To analyze the quality of the sequencing data and call methylation levels, the Bismark (Krueger and Andrews, 2011) workflow pipeline (Ewels et al., 2020) was applied with GRCm38 (mm10) as reference genome. All samples used had PHRED quality scores over 34, indicating a base call accuracy of 99.9%. 11bp was set as adapter trimming length for both 5’ end and 3’end, read 1 and read 2, to achieve maximum mapping efficiency (up to ~70%). On average, sequencing produced 68 million reads per sample with 99.92% conversion rates and 58% mapping efficiency in alignment. After methylation levels were called, 40 million mapped reads and 2.2 million 5x CpGs were obtained on average. CpG sites captured at less than 5x were discarded from further analysis. 1 young female and 1 old male mouse were excluded from quantitative analysis due to insufficient read counts. SI Table 1 shows a summary of RRBS sample properties.

Differentially methylated region (DMR) Iientification

To identify the DNA methylation differences between the young and old mice, Differentially Methylated Cytosines (DMCs) were called by the R package MethylKit (Akalin et al., 2012). The output coverage files of Bismark were used as input for DMC calling. The CpGs in the input files were filtered for a minimum 5x read coverage depth for DMC calling. 3 tests in MethylKit (Fisher’s Exact test, logistic regression, and beta-binomial model-based methods) were used for DMC calling and to assess both across-group difference in methylation and heterogeneity within groups. The same cutoffs were used for all three DMC calling methods: q-value < 0.05, ≥ 30% absolute methylation percentage difference, with CpGs captured in a minimum of 5 samples for each group. DMCs outputted from the 3 methods were then merged and extended into 100 to 1,000 base pair (bp) sized DMRs based on methods previously described (Ziller et al., 2013). Young versus aged sample DMRs were generated using methylation data merged from both sexes according to the pipeline described above, while sex specific DMRs were obtained using the same pipeline but with CpGs captured in a minimum of 3 samples for each group (according to both age and sex). The overlapped DMRs across sex were defined as the DMRs with at least one base pair overlap across sex specific DMRs.

Clustering and enrichment analysis across genomic regions

The overall similarity between samples was characterized based only on their methylation profiles using clusterSamples() in MethylKit with the ‘correlation’ as distance metric, and ‘single’ as the agglomeration methods. To study the distribution and enrichment of 5x CpGs and DMRs across genomic regions, we retrieved genomic annotations from UCSC Table Browser. The annotation for genes, promoters, exons, introns, 3’UTRs were extracted from UCSC RefSeq track, whereas the LINEs, SINEs, LTRs were obtained from RepeatMasker track. The promoters were defined as regions 2000bp upstream and 200bp downstream of transcription start sites (TSS). CGIs, Enhancers, and Enhancer-Promoter Interactions were downloaded from CpG Island track and ENC+EPD Enhc-Gene track within whole brain tissue, respectively.

In calculating fold enrichment, only CpGs with 5x coverage were used. The DMR fold enrichment within a specific genomic feature, e.g. CGI, is defined as the ratio of the number of CpGs in DMRs falling within a CGI over the number of CpGs captured in RRBS that fall within a CGI (see formula below).

In order to evaluate the significance of DMR fold enrichment, a comparable background was generated. This was performed by calculating a separate dataset of fold enrichment values per genomic feature using random downsampling of CpGs (sampling dataset), where the number of CpGs in each sampling dataset was equal to the number of CpGs in DMRs of that respective feature. Sampled fold enrichment was defined as the ratio between the percentage of sampled CpGs in a feature across all sampled CpGs captured and the percentage of CpGs in a feature across all CpGs captured (see formula below).

The sampling process was performed 5,000 times per feature to calculate the feature-specific variability of sampled fold enrichment values. The resulting 5,000 different fold enrichment values in each feature were assumed as a Gaussian background distribution. p-values of feature-specific DMR fold enrichment were calculated based on z-score using the background fold enrichments obtained from the sampling dataset (Figs. 2B and C; Figs. 3C and D). A fold enrichment over 1 indicates that a hypermethylated or hypomethylated region is enriched for that respective feature relative to its representation in our data.

Fig. 2.

(A) Total number of differentially methylated regions (DMRs) from both males and females that are either hypermethylated or hypomethylated with age (≥30% absolute methylation difference, q < 0.05, 5 samples minimum per group). (B) Distribution of hypermethylated (dark green) or hypomethylated (light green) DMRs across various genomic features. Fold enrichment was normalized to the ratio of DMRs in that feature context, where a fold change over 1.0 (grey dotted line) is an overrepresentation of either hypermethylated or hypomethylated DMRs within that region. Significant enrichment by Z-test of hyper- or hypomethylated DMRs within that respective feature is indicated by the p-value cutoffs: * p < 0.05; ** p < 0.01, *** p < 0.001. (C) Fold enrichment for hypermethylated and hypomethylated DMRs showing the ratio between the percentage of histone modification peaks in DMRs vs. the percentage of peaks in all 5x CpGs. Fold enrichment greater than 1.0 indicates that the histone marks are overrepresented in DMRs, while less than 1.0 are underrepresented. Significant enrichment by Z-test of hyper- or hypomethylated DMRs within that respective feature is indicated by the p-value cutoffs: * p < 0.05; ** p < 0.01, *** p < 0.001. (D) Gene Ontology (GO) terms plotted against their associated gene distance to regulatory domain. Peaks used to generate GO term list are at 25kb and 77kb (shaded regions). (E) GO term list generated through GREAT analysis of hypermethylated regions only for the first, proximal peak and second, distal peak. Scale p = 0.1 to p = 0.0001. (F and G) Top motifs associated with hypermethylated (F) and hypomethylated (G) regions from merged male and female samples (HOMER software).

Fig. 3.

(A) Hypermethylated and hypomethylated DMRs that are shared between male and female samples. (≥30% absolute methylation difference, q < 0.05, 3 samples minimum per group). (B) Mean methylation difference (young to old) in male and female groups within different genomic features. (C) Across-sex comparison of DMR fold change enrichment for specific genomic features. Significant enrichment by Z-test of hyper- or hypomethylated DMRs within that respective feature is indicated by the p-value cutoffs: * p < 0.05; ** p < 0.01, *** p < 0.001. (D) Enrichment for histone modification marks in hypermethylated DMRs of male and female samples (left panel) and hypomethylated DMRs of male and female samples (right panel). Significant enrichment by Z-test of hyper- or hypomethylated DMRs within that respective feature is indicated by the p-value cutoffs: * p < 0.05; ** p < 0.01, *** p < 0.001. (E and F) Enriched motifs in hyper- and hypomethylated DMRs of male and female samples. Scale p = 0.1 to p = 0.001. (G) Number of GO terms plotted against their max distance to TSS. Peaks used for peak convergence are shown in males at 90kb, and females at 25kb and 80kb (shaded regions).

The histone modification peaks used for fold enrichment analysis in (Fig. 2C, Fig. 3D) was downloaded from the ENCODE/LICR ChIP-seq data (Roberston et al., 2007). The fold enrichment of DMRs in histone modification peaks is calculated based on the following formula, e.g. for H3K4me3:

MeDIP qPCR

A separate cohort of young adult (2 months old) male and female C57BL/6J mice were obtained from Jackson Laboratory and aged (20 months old) male and female C57BL/6 mice were obtained from the National Institute on Aging Aged Rodent colony. Mice were fed and maintained in the same manner as the sequencing cohort. The mice were sacrificed and tissue collected in the same manner as the sequencing cohort. DNA was isolated from dorsal hippocampus, ensuring 260/280 quality was above 1.9, and 5ug was used for chromatin shearing. Samples were processed with the MagMeDIP kit according to manufacturer instructions (Diagenode, Denville, NJ). Purified MeDIP DNA was used for quantitative PCR analysis using the Roche 480 Lightcycler II (Roche, Basel, Switzerland) and SYBR green. Inputs were normalized to the dilution factor used (10%). Percentage input was calculated with the % input method and enrichment was found using the equation AE^(adjusted input Ct - IP Ct) where AE is amplification efficiency. Primer sequences use to amplify DMR genomic regions were designed using Primer3 software and were as follows: Mvp L 5’-CGAGGACATTAACAGGGTGC-3’; Mvp R 5’-CCCGGGAGTAGATGTCTGTC-3’; Bmf L 5’-TGCTGAGATGGATGCAGACA-3’; Bmf R 5’-AGGCTTCAAGTCTGACTCGG-3’.

Statistical analyses

Welch’s t-tests were performed to determine statistical differences in the mean average methylation levels between young and aging mice groups as well as male and female groups, separately (Fig. 1D, left panel). This test does not assume equal population variance (Welch, 1947). Moreover, F-tests were applied for testing the variance differences between groups (Fig. 1D, right panel). These tests were applied to average methylation levels from samples, both globally (Fig. 1D) and within specific genomic features (Fig. 3B).

Fig. 1.

(A) Schematic of sample collection and processing. Created in Biorender software. (B) Number of CpGs mapped at ≥ 5x coverage depth captured in young and old male (YM and OM) and female (YF and OF) mice. Columns are split between CpGs captured that are specific to one age group versus shared (M and F Common) across age groups. (C) Average methylation level of ≥ 5x CpGs for each sample for male (top) and female (bottom) mice. (D) Mean (left) and variation (right) of global DNA methylation levels across samples grouped by age and sex. Welch’s t-test and F-test were performed for assessing global statistical differences in mean and variance of methylation levels, respectively, indicated by the p-value cutoffs: * p < 0.05; ** p < 0.01, *** p < 0.001 (E) Hierarchical clustering analysis of methylation levels across each individual sample within sex. Height represents the distance between samples using correlation as metric for sample clustering in methylKit.

To test whether the divergence in mean methylation changes across sex is statistically significant, Welch’s t-tests were performed across sex for CGIs and LINEs (Fig. 3B). In testing, we assumed distributions of mean methylation changes from young to aging between males and females were unequal. Therefore, the mean methylation changes used as input for t-tests were calculated by pairing all within-sex combinations of young samples and aging samples. To test for significant fold enrichment within genomic regions and histone modification marks, p-values were calculated using z-scores through a Z-test assuming a normal distribution obtained from the sampling dataset described above in Methods section 2.7.

Functional analysis

The Genomic Regions Enrichment of Annotations Tool (GREAT) v4.0.4 (McLean et al., 2010) was used to reveal the functional significance of the DMRs found in this study. DMRs were set as test regions and the whole genome mm10 was set as the background region in analysis. The gene regulatory domain in GREAT was set for 5kb upstream and 1kb downstream as proximal and 200kb as distal, with the max regulatory domain in Fig. 2D and Fig. 3G defined as the proximal plus distal domain. To detect the enrichment of known binding motifs in DMRs, HOMER (Hypergeometric Optimization of Motif EnRichment) v4.11 was used with the repeat sequence masked, 200bp fragment size for motif finding and genome mm10 as background (Heinz et al., 2010). Enrichment for histone modification marks was performed using ChIP-seq (Roberston et al., 2007) datasets obtained from ENCODE/LICR.

Gene expression analysis and integration with DNA methylation

To link the DNA methylation changes with its transcriptional consequences, we used RNA-Seq data from young and aging mice from our published work (Kwapis et al., 2018). Differential expression analysis was applied using Cyber T as previously described (Kayala and Baldi, 2012) to identify up- or down-regulated genes in the aged group. The paired t-test in Cyber T was used for young and aging group comparison with FPKM values as input. The log fold change for DEGs is defined as log2(Old Mean FPKM / Young Mean FPKM). Stanford’s Genomic Regions Enrichment of Annotations Tool (GREAT) software (McLean et al., 2010) (Version 4.0.4) was used with default parameters (basal plus extension/proximal 5Kb upstream, 1Kb downstream, plus distal up to 1000Kb) to find mm10 UCSC genes associated to input DMR files. From there, gene lists were intersected with DEGs identified between young and old mice. Overlapping genes were then submitted to STRING (Szklarczyk et al., 2019) (v11.0b) to identify protein-protein interaction (PPI) networks, with a provided PPI enrichment p-value for the generated network. Genes making up the most connected network were submitted to the PANTHER Classification System’s web interface (Mi et al., 2019) (v15.0) in order to determine enriched GO terms (FDR adjusted p-value ≤0.05).

Results

Aging mice show global hypomethylation and increased variability in DNA methylation across animals

Young and aging mouse hippocampus was dissected and prepared for RRBS (Fig. 1A). Specifically, we analyzed dorsal hippocampal tissue from a total of 40 mice, divided equally between young and aging male and female mice, with one old male and 1 young female excluded from analysis due to insufficient quality (10 young male YM, 9 old male OM, 9 young female YF, 10 old female OF). In total, there were approximately 7.8 million (M) CpGs captured in both young males and females, each with at least 5x read-depth. The total number of ≥ 5x CpGs in the aged groups was lower, with 4.8M in aging males and 6M in aging females (Fig. 1B). By overlapping the CpGs across the two age groups, we found 4.1M CpGs were shared across young and aging males. In contrast, a higher number of shared CpGs, around 4.8M, were found in females. The non-shared CpGs, meaning CpGs at sites in the genome that were captured as unique to 1 group, comprised 46% and 37% of the total in males and females, respectively, for the young group, and 13% of CpGs in males and 19% of CpGs in females labeled as unique to the aging group.

We calculated the average methylation level of CpGs captured across samples (Fig. 1C) and found the aging group, both in males and females, exhibited losses in global methylation overall and increased methylation level variability from sample to sample. Aging males showed significant (Welch’s t-test, p = 0.006, Fig. 1D, left panel) loss of average methylation levels compared to their young counterparts, whereas the loss in average methylation levels was non-significant in females (p = 0.056, Fig. 1D, left panel). Given the efficiency at which RRBS captures CpGs within the genome per number of reads sequenced, we were able to perform methylation analysis using a larger number of animals than done previously (Masser et al., 2017), while still capturing over 4M CpGs commonly shared across both age groups (albeit at lower sequencing depth). This allowed us to examine variability in methylation global patterns between individual animals. We observed a notable increase in the average variability in global methylation across aging samples compared to young samples in both males and females (Fig. 1D, right panel; F-test, p < 0.001 and p < 0.01, respectively).

To make an initial assessment of methylation divergence at single-CpG resolution across age, we performed a hierarchical clustering analysis on CpGs based on methylation ratio. This analysis confirmed that global methylation patterns in the dorsal hippocampus were indeed divergent across the lifespan and that aging animals generally showed more similarity in genome-wide methylation patterns compared to other aging animals rather than young animals (Fig. 1E). Interestingly, while correlation distance mostly ordered the samples according to age, many young and aging animals were closer in distance than some animals from their respective age groups. For example, old males 9 and 4 (OM04 and OM09) and one old female (OF09) did not collectively show divergent methylation patterns compared to young animals.

Hypomethylated regions in aged mice are associated with regulatory factor binding motifs related to plasticity

We next identified differentially methylated regions (DMRs) to investigate which features of the genome, as well as specific regulatory factor binding motifs, display altered methylation patterns in aging. We first calculated the total number of DMRs found between age groups, irrespective of sex. Regions with significantly higher methylation levels (≥ 0.3 methylation difference, q < 0.05, 5 sample minimum per group) in aging mice when compared to the same region within young mice were termed hypermethylated, whereas regions with significantly lower methylation levels in aging mice compared to young were termed hypomethylated (Fig. 2A). Over half (around 56%) of the total DMRs fell within the hypermethylated category, meaning regions with higher amounts of methylation in aging mice than compared to the same regions within young mice. Next, fold enrichment was normalized by the ratio of DMRs within a genomic feature and statistical comparisons were calculated through Z-tests. We found enhancers (ENCODE 3 UCSD/Ren) displayed the largest degree of fold enrichment in aging (Fig. 2B), with both hyper- and hypomethylated DMRs showing significant changes in these regions (p < 0.001). CpG islands (CGIs), areas that are densely populated with CpGs, contained a significant underrepresentation of hypomethylated regions in aging mice (p < 0.001). Hypomethylated DMRs were also highly underrepresented within promoter regions (p < 0.001). Of the genic regions analyzed, exons showed the greatest divergence in fold enrichment. Notably, DMRs were largely underrepresented within LINE and SINE repetitive regions. Histone modifications (Fig. 2C) showed only modest changes in DMR representation in aging, although H2K27me3 marks, typically associated with transcriptional repression, showed significant fold-enrichment in hypermethylated DMRs (p < 0.05).

Next, DMRs were compiled and run through Genomic Regions Enrichment Annotations Tool (GREAT) analysis (McLean et al., 2010) to identify genes associated with regions dynamically methylated during aging, and their subsequent gene ontology (GO) enrichment (Fig. 2D). The number of GO terms from hypermethylated DMRs were plotted against their associated gene(s) max distance from regulatory domain, as defined in GREAT from 0kb to 100kb, where multiple peaks identified as local maximums at 25kb and 77kb were observed. We compared the terms in the proximal first peak to terms found in the distal second peak (Fig. 2E). Only a small subset of the GO terms overlapped between the two peaks, namely covering receptor protein signaling pathways. Terms found only within the first peak included more neurogenesis and differentiation functions, as well as a handful of processes encompassing negative regulation of RNA functionality and DNA transcription, whereas the second peak encompassed cell morphogenesis in addition to adult behavior. Notably, hypomethylated regions did not generate any GO terms, similar to what was previously reported in a DMR study in whole hippocampal tissue (Harris et al., 2020). Overall, we observed an overrepresentation of GO terms associated with hypermethylated regions and a lack of annotations for hypomethylated regions, of which we describe potential implications in the discussion.

To further investigate whether methylation patterns can be indicative of changes in gene regulatory elements with age in the dorsal hippocampus, we performed motif enrichment analysis (HOMER) (Heinz et al., 2010) from both hypermethylated and hypomethylated regions to determine if the binding of specific regulatory factors could be impacted by changes in DNA methylation. The majority of motifs for hypermethylated DMRs include core homeobox development transcription factors and had similar sequence structure, as seen with the Lhx family (p < 0.01, Fig. 2F). Interestingly, motifs mapped to hypomethylated DMRs represent a neural activity related signature, with the majority composed of the AP-1 transcription factor complex family (JUN, FOS) and also displaying high sequence homology (Fig. 2G). Overall, hypermethylated and hypomethylated DMRs from our data reveal a new list of enriched transcription factor binding sites specific to the dorsal hippocampus that were not found in (Harris et al., 2020) and a stark contrast in the types of regulatory factors that bind to genomic regions that gain versus lose methylation with age.

Sex-specific divergence in genomic feature enrichment at differentially methylated regions

Abundance of hypermethylated regions in females in both humans and mice has been shown to be linked preferentially to autosomal regions, at least within liver tissue (Zhuang et al., 2020; McCarthy et al., 2014; Grimm et al., 2019). As age- and sex-specific difference were previously reported in (Masser et al., 2017), and given that sex differences in cognition and memory behavioral tasks as a function of age are prevalent (Jahanshad and Thompson, 2017), we investigated whether there are sexually divergent DNA methylation patterns in the aging dorsal hippocampus in a variety of genomic landscape features across the autosomal chromosomes. To achieve this, a new set of DMRs were established based on differentially methylated patterns found across age in females and males separately (≥ 0.3 methylation difference, q < 0.05, 3 sample minimum per group), although CpGs were still required to be captured in both males and females for direct comparison of sex-specific DMRs. Female samples had 805 hypermethylated DMRs and 653 hypomethylated DMRs, whereas male samples had 432 hypermethylated DMRs and 495 hypomethylated DMRs (Fig. 3A). Among the DMRs identified, only 14 hypermethylated and 22 hypomethylated regions were shared across sex. Interestingly, the number of hypermethylated DMRs in females was nearly 2-fold higher than in males, with the numbers of hypomethylated DMRs between the sexes slightly more similar. Comparison of mean methylation levels within genomic landscape features (Fig. 3B) showed that aging males exhibited significant hypomethylation in almost every feature (F-test, p < 0.01) except in CGIs and long interspersed nuclear elements (LINEs). The 3 largest mean changes were found in exons, introns and 3’ untranslated regions (3’UTRs). Conversely, aging females display less pronounced decreases in mean methylation, with 3’UTRs, exons, introns and short interspersed nuclear elements (SINEs) showing significant hypomethylation (F-test, p < 0.05). We then compared the differences between CGIs and LINEs, the only 2 contexts that slightly increased in methylation with age in males, to those same features in females and found LINEs were significantly higher in methylation in males rather than females (p < 0.001).

By focusing on hyper- and hypomethylated DMR enrichment within genomic features (Fig. 3C), we observed some notable differences between males and females. The largest sex-specific divergence in DMR enrichment was observed within 3’UTRs. Specifically, female fold enrichment was measured at 2.3 (left panel), while in male enrichment was measured at 0.7, illustrating a significant overrepresentation of hypermethylated DMRs in 3’UTRs within aging females only (p < 0.001). The opposite trend was observed regarding hypomethylated DMRs, with males showing a 2.1 fold enrichment while in female tissue it was 0.56, illustrating a significant overrepresentation of hypomethylated DMRs in 3’UTRs within aging males (p < 0.001). Of note, hypomethylated CGIs enrichment was significantly depleted in both males and females (p < 0.001), mirroring the absence of sex-combined hypomethylated DMRs within CGIs. However, female tissues did show significant overrepresentation in hypermethylated DMRs within CGIs (p < 0.001), while males remained significantly underrepresented (p < 0.05). Female-specific hypermethylated DMRs were also uniquely overrepresented within exons. In contrast, hypomethylated DMRs were significantly overrepresented in LINE elements within females (p < 0.01), while all other DMR contexts were underrepresented in LINE and SINE elements in males and females. Notably, hypermethylated and hypomethylated DMRs showed significant overrepresentation in enhancer regions across both sexes. All statistical comparisons were calculated through Z-tests.

We next performed an analysis of sex-specific DMR enrichment with respect to histone modifications using the same method described above in Fig. 2C. Here, female tissues again showed dominance in hypermethylated DMR overrepresentation across several histone mark contexts (e.g., H3K4me3, H3K27ac, H3K27me3), although both male and female hypermethylated DMRs were overrepresented in regions marked by H3K4me1, a mark for active enhancers (Fig. 3D). In contrast, hypomethylated DMRs showed only one significant change in representation across histone marks, with males showing significant overrepresentation in H3K4me1 (Z-test, p < 0.01). Furthermore, we also observed some sex-specific changes in methylation near genes encoding epigenetic regulators (SI Table 2).

Next, we sought to investigate whether there are potential transcriptional mechanisms impacted by the age- and sex-specific methylated patterns identified. To address this, binding motif analysis using HOMER was performed for male- and female-specific DMRs separately (Fig. 3E and Fig. 3F). In general, female samples showed greater enrichment in binding motifs in both hyper- and hypomethylated DMRs compared to males. Hypermethylated DMRs in both sexes were particularly enriched for a number of homeobox motifs, similar to sex-combined hypermethylated DMRs. While no motifs were found in male hypomethylated DMRs, the majority of female hypomethylated DMRs were related to the basic Helix-Loop-Helix(bHLH)/basic leucine zipper(bZIP) families and a neuronal-specific activity signature rather than general prenatal/postnatal development. Additionally, GREAT analysis was performed in males and females with input from hypermethylated regions, as hypomethylated regions again did not associate with any GO terms (Fig. 3G). Three peaks as local maximums were observed: two in females at 25kb and 80kb and one in males at 90kb. The GO terms from these peaks were compiled and listed and in a heatmap (SI Fig. 1C). DMRs near the 25kb peak, which were largely associated with female-specific hypermethylated DMRs, were enriched in terms associated with regulation of cellular signaling processes, RNA synthesis and general transcriptional functionality. By contrast, GO terms from the 80kb peak in females and the handful from the 90kb peak in males were more widely associated with development, ranging from neuron differentiation to ear, lung, sensory organ, nephron, genitalia and skeleton development, reinforcing the overarching developmental signature observed across our analyses.

RNA sequencing gene expression analysis with respect to differentially methylated regions

Given the prevalence of gene regulatory elements at regions undergoing dynamic changes in methylation levels throughout aging in our results, we next focused our analysis on determining correlations between our DMRs and differentially expressed genes (DEGs) derived from an RNA-Seq dataset generated from dorsal hippocampal tissue of young (3 months) and aging (18 months) homecage C57BL/6J mice (Kwapis et al., 2018). Among the genes captured in young and aging mice of both sexes, a total of 300 were significantly increased in expression in aged mice (paired t-test, p < 0.05), while 108 genes were significantly decreased in expression across the same age comparison (Fig. 4A; see SI Table 3 for list of DEGs). We next separated DEGs by sex and association with a DMR from our methylation data and found the total number of female-specific DEGs to be similar to that of males, however the number of DEGs associated to DMRs in females was higher than in males (Fig. 4B; SI Fig. 2A) and thus we focused our remaining analysis in (Fig. 4) on the 44 DEGs that overlapped with a DMR in females only. We created a protein-protein interaction (PPI) network to examine functional protein associations. The genes that comprised the strongest network for females were enriched with GO terms (FDR adjusted p-value ≤ 0.05) for learning, amyloidbeta, and inflammation/immune system categories (Fig. 4C), while males presented a smaller network with no GO terms (SI Fig. 2B). We then categorized female DMR associated DEGs by fold change and gene distance to any DMR, which resulted in a dense cloud of DEGs concentrated around 105 bp to 107 bp in distance from the nearest associated DMR (Fig. 4D). These genes exhibited variation in expression, showing both up- and downregulation covering a range of log2Fold Change between ± 5 or greater in both males (SI Fig. 2C) and females with a minimum fold change threshold of 1. Interestingly, assigning the DEGs to an enhancer or promoter context within a DMR (either hyper- or hypomethylated) resulted in a greater number of DEGs mapped to female DMRs than males (Fig. 4E and F; SI Fig. 2D and E), with a handful of female-specific genes in these regions listed as the same in our PPI network. More specifically, 10 DEGs mapped to female hypermethylated DMRs and 5 mapped to male hypermethylated DMRs (SI Fig. 2D). In contrast, mapped DEGs to hypomethylated DMRs were more even across the sexes, with 5 total in females (two in enhancers and three in promoters) and 6 total in males (four in enhancers and 2 in promoters). While many had smaller fold change values, these associations revealed a number of genes that may allow for meaningful insight into age-related changes in cognitive function in the dorsal hippocampus, a few of which are highlighted in the discussion. To validate a DMR within a promoter region, we performed MeDIP qPCR in a separate cohort of young and aging female mice of the same ages for Mvp, an upregulated gene in aging animals from the RNAseq data that also showed hypomethylation with age in our dataset at a location most proximal to its TSS (relative to all other DMR-associated DEGs) (Fig. 4D). Primers were designed to flank the sequence of the associated DMR. We found a trend of less methylation for this region in aging mice, though not significant (Fig 4G). Additionally, we chose a non-DEG associated DMR region in males for MeDIP qPCR that was identified as a hypomethylated DMR in the 3’UTR region of Bmf, a regulator of BCL2, as we found that 3’UTR regions showed highly divergent patterns in methylation during aging. Indeed, we found a significant decrease in 5mC enrichment in aging males for this region (SI Fig. 2F, p < 0.0235; two-tailed t-test, n =10).

Fig. 4.

(A) Differentially Expressed Genes (DEGs) with log2 fold change and p-value, where log2Fold change = log2(mean old FPKM / mean young FPKM). (B) Female specific DEGs that associate with a DMR (middle) compared to total DEGs (left) and total genes associated with DMRs (right). (C) Protein-protein interaction network (PPI) generated from STRING for female specific DMR associated DEGs. (D) Log2Fold Change of DEGs plotted against their distance (from TSS) to DMRs with a fold change threshold of 1.0. Colored dots above 1.0 or below −1.0 are categorized as DEGs. Genes that overlap as a target of either an enhancer or promoter within a DMR are colored as blue or purple, respectively. (E) Barplot of log2Fold change in expression of genes overlapped with enhancers and (F) promoters. (G) MeDIP enrichment for DMR within promoter of Mvp, p < 0.12; unpaired parametric 2-tailed t-test; n = 6.

Discussion

Sex-specific divergence of methylation patterns and gene expression across the genomic landscape

In this study, we found aging females show a different methylation signature overall than aging males. Though aging males have a significant drop in average global methylation levels and a greater increase in average global methylation variability from sample-to-sample (Fig. 1D and 1C), aging females show a greater amount of regional significant change, with a tendency for hypermethylation at genomic features and histone modifications (Fig. 3C and 3D). Of the total DMRs identified in autosomes in our data, females also had a greater proportion of DMRs than males in both hyper- and hypomethylated regions (Fig. 3A) despite similar coverage (Fig. 1B), leading to more DMR association with transcription factor binding motifs and GO terms (Fig. 3E, 3F, and 3G). Furthermore, our results mirror data that show sex-specific methylation differences in autosomes, particularly with males exhibiting higher global levels but females exhibiting higher methylation levels in DMRs (McCarthy et al., 2014; El-Maarri et al., 2007; Mamrut et al., 2015), and show that despite the aforementioned studies which used young to middle-age human tissue, these sex-differences can be found within the dorsal hippocampus of aging mice specifically and can be stable across the lifespan, as reported in (Masser et al., 2017).

With regard to sex-specific differences in enrichment for hyper-and hypomethylated DMRs across multiple genomic features, genomic 3’UTR regions were found to have the strongest, sex-dependent variation of methylation levels in aging, displaying extreme hypermethylation in females but hypomethylation in males to a similar degree (Fig. 3C). 3’UTRs are associated with translation, localization, and general functionality of mRNA transcripts (Bae and Miura, 2020; Wang et al., 2019), and have long been recognized as important regulators of gene expression. Changes in methylation dynamics at 3’UTRs have been largely understudied, with the majority of existing literature reporting that alteration of methylation levels at 3’UTR sequences does influence gene expression to some extent (McGuire et al., 2019), however an in-depth analysis of the role of methylation in this context, particularly in the aging hippocampus, has yet to be explained. The observation of sex-specific differences in 3’UTR methylation levels reported here creates an intriguing opportunity for future studies to focus on both sexually divergent and genomic landscape-specific properties of methylation that have molecular and behavioral implications.

Female-specific DEGs associated with a DMR also showed genes involved in hippocampal-dependent memory. For example, Fcgr2b (log2FC 4.361, BH p < 0.04) within our PPI network (Fig. 4C) and promoter hypermethylated DMRs (Fig. 4E) has been linked to expression-dependent modulation of long-term potentiation and amyloid beta neurotoxicity (Kam et al., 2013), highlighting the sex divergence seen in Alzheimer’s disease pathology (Ferretti et al., 2018). Interestingly, a gene found to be associated with enhancer hypermethylated DMRs in aging males, Gpr176 (SI Fig. 2D, log2FC −0.312, Benjamini-Hochberg (BH) p < 0.04), is a g-protein coupled orphan receptor involved in circadian rhythm patterns and negative regulation of cAMP activity (Doi et al., 2016). This is not the first appearance of gene expression variance in aging linked to circadian rhythm genes, which was reported previously in (Kwapis et al., 2018), where experience-induced expression of Per1, another major circadian clock gene, was restricted by the histone deacetylase HDAC3 in aging mice but rescued upon hippocampal deletion of HDAC3 in HDAC3flox/flox mice, along with age-related impairments in memory. A majority of the rest of the DEGs that are associated with DMRs in Fig. 4E and F and SI Fig. 2D and E are markers of inflammation and immune response, providing more functional pathways to investigate as potential target areas of age-related influence on methylation or other epigenetic mechanism dynamics.

Methylation dynamics influence on transcription factor binding

Motif enrichment analysis revealed the AP-1 complex (including the proteins JUN and cFOS), a major regulator of gene expression, as the most commonly associated transcription factor binding sites within hypomethylated regions in aged mice when compared to young animals in our results (Fig. 2G). In humans and mice, DNA methylation is known to influence the binding of transcription factors to their recognition motifs, either positively or negatively (Héberlé and Bardet, 2019; Stadler et al., 2011). Within the few studies examining AP-1 family binding affinity, 2 studies identified significant associations of transcription factor binding (i.e. BATF) with methylated CpGs at the AP-1 binding sequence. However, none explain potential effects of methylation at transcription factor binding sites with regard to aging dorsal hippocampal tissue (Qiao et al., 2016; Ji et al., 2019; Zuo et al., 2017; Ng et al., 2013), suggesting this is a key area of needed investigation.

Studies have shown that AP-1 regulates expression of immediate early genes, a gene class that is activated in response to neural activity and are associated with the stabilization of synaptic plasticity and long-term memory consolidation (Guzowski et al., 1999; Tischmeyer and Grimm, 1999). However, it remains unclear how AP-1 underlies mechanisms of long-term memory and long-term potentiation (LTP), especially with regard to DNA binding and interaction with epigenetic mechanisms. Thus, methylation dynamics at transcription factor binding sites (TFBSs) may serve as an important mediator of binding affinity and subsequent gene expression, particularly within plasticity associated transcription factors. For example, CREB-dependent regulation of c-Fos has been well described in long-term memory formation and in several brain regions including the hippocampus (Katche et al., 2010). Besides a key role in memory formation, FOS and JUN are bZIP proteins that are known to interact with other bZIP proteins in a cell-type and stimulus-specific manner in order to regulate expression of downstream targets via their promoters and enhancer elements. Our results suggest that investigating the mechanism by which hypomethylation occurs in the aging dorsal hippocampus at AP-1 controlled genes may reveal important insight into whether this DNA methylation change drives age-dependent cognitive dysfunction, or in contrast, is a possible compensatory mechanism that opens chromatin at these key sites to facilitate memory formation in the aging brain (Biddie et al., 2011). As proof of principle, we and others have demonstrated that opening chromatin can facilitate gene expression in the aging brain and allow subthreshold learning events to be encoded into long-term memory (McQuown et al., 2011; Kwapis et al., 2019; López et al., 2020).

Correlating methylation as an epigenetic mark for age-related cognitive impairment

Both human and rat memory studies have identified age-impaired (AI) individuals, where their cognition is significantly negatively affected with age, and age-unimpaired (AU) individuals, where there is efficient cognitive function similar to that of a young adult (Stark et al., 2013; Kosik et al., 2012; Castellano et al., 2012). The latter study showed that chronological aging regulated the epigenome, namely in basal and experience-dependent histone acetylation levels in rats that were characterized as being AI or AU in a spatial memory task dependent on the hippocampus, especially on the activity of the dorsal hippocampus. Numerous studies have used methylation marks to develop regression models to infer the epigenetic age of an individual but also to identify potential biomarkers and phenotypes associated with changes to the methylome (Petkovich et al., 2017; Thompson et al., 2018; Coninx et al., 2020), and there is a growing amount of evidence that epigenetic mechanisms are involved in age-related memory impairments (Fischer et al., 2007; Peleg et al., 2010; Oliveira et al., 2012). Thus, exploring variation in biological age between individuals, despite similar chronological age, is an intriguing question that may allow for exploration of epigenetic differences that can be investigated for correlational, or perhaps causal, relationships with cognitive function. E.g., our analyses showed an increase in animal-to-animal variation in global methylation levels in aging of both sexes but more in aging males - the extent of which was not observed in the younger groups. A fold enrichment analysis of CpGs that fall within DMRs (in various genomic contexts) across both sexes suggests that this difference is not due to DMR capture bias (SI Fig. 3B) and that variable global methylation levels by sex cannot be explained by sex-specific bias in DMR locations, which is discussed more below. Variability may, however, be indicative of the site-specific, rather than global, changes to methylation that may occur within one animal and differ from the next, although further analysis will have to be conducted in order to further develop our findings in methylation variance with age.

Additionally, though average methylation levels showed a global trend of hypomethylation with age (Fig. 1D), there were a greater number of hypermethylated DMRs captured, between aging animals and young, relative to hypomethylated DMRs (Fig. 2A). These data support prior studies that show hypermethylation arises with age in a genomic location-specific fashion, particularly in CpG islands, which correlate with age-related expression changes and pathological phenotypes such as those found in cancer (Xiao et al., 2016; Rakyan et al., 2010). Regions showing hypomethylation across age are numerous across the genome but are of much smaller abundance when only considering CpG sites that show statistically significant (based on q-value) losses in methylation. Conversely, sites that increase in methylation, though fewer, hold greater significance on a per-CpG level and subsequently drive hypermethylated DMRs. Both site-specific and global variances in methylation levels among samples have been described in both single- and multi-tissue studies. The correlations in methylation changes with age observed in our study with those found in other studies could be indicative of chronological and biological age patterns that arise in methylation specifically within the dorsal hippocampus, potentially creating new opportunities to measure biological age within brain regions important for memory formation and consolidation.

Previous results from Masser et al., 2017 utilized a bisulfite oligonucleotide-capture sequencing strategy that provided an extraordinarily high sequencing depth for targeted CpGs (20–40x) throughout the genome to measure changes in methylation across age within the hippocampus. This bait enrichment strategy enabled (by design) a high-resolution methylation analysis focusing specifically on CpG Islands (CGIs) and flanking regions of up to 4kb in both up and downstream directions from CGIs. Within this genomic context, no significant differences in methylation were observed across age (from 3 to 24 months). This suggests that “global” trends in methylation observed within a given tissue may be somewhat dependent on bisulfite sequencing CpG capture strategy (e.g., given that CGIs are typically hypomethylated and underrepresented within somatic tissues). Despite these differences, both studies (ours and Masser et al., 2017) observe site-specific changes in methylation across age as well as notable lifelong sex differences.

Our study expands on the findings of prior works that also performed methylation analysis in the hippocampus with single-base resolution by including specific genomic-region annotations for promoter and enhancer contexts, the latter of which can include distal enhancers at non-coding intergenic regions that bind transcription factors to act as a cis-regulatory element and are crucial for tissue-specific gene expression regulation (Tuvikene et al., 2021; Stadhouders et al., 2012; Ko et al., 2017). We also showed both hyper- and hypomethylation in H4K3me1 (Fig. 3D), a histone modification linked to primed and active enhancers and transcription, suggesting significant changes in DNA methylation can occur in regions associated with distal regulatory regions (Robertson et al., 2008).

It is important to note that while the RRBS performed in this study provides benefits in terms of lower cost and sample input for read depth, this method results in incomplete coverage of the genome, and has the possibility of capturing 5hmC. Though different methods of methylation analysis can create different target result biases, many of the studies aimed at creating epigenetic clocks with methylation values utilized RRBS and a regression analysis to do so (Petkovich et al., 2017; Thompson et al., 2018; Wang et al., 2017). In this study, we aimed to create a foundation for the addition of other age groups to identify an ‘epigenetic clock’ that is specific to the dorsal hippocampus, and behavioral groups to determine whether DNA methylation variability or aspects of a clock may predict whether an animal exhibits age-impaired or age-unimpaired cognitive function.

Conclusion

Previous studies have demonstrated that DNA methylation is involved in regulating gene expression required for memory formation (Penner et al., 2010; Penner et al., 2011; Moore et al., 2013; Miller and Sweatt, 2007) and that alterations in DNA methylation as well as other epigenetic regulatory mechanisms may underlie age-dependent cognitive dysfunction (Lubin et al., 2008; Lardenoije et al., 2015; Jarome et al., 2015). However, how DNA methylation patterns across the genome change in a brain region central for learning and memory in an age- and sex-specific manner remains unclear. In this study, we examined this question by selecting to focus on the dorsal hippocampus, which is a brain region that is pivotally involved in memory formation, becomes impaired with age, and exhibits sex-specific memory functions.

In general, we found that aging mice exhibited greater variability in average methylation levels. We also found several differences in the total number of DMRs as a function of age, across several genomic features and particularly at hypermethylated regions. Examination of what types of genes these DMRs may be affecting revealed that hypermethylated DMRs associated with homeobox-regulated genes, whereas hypomethylated DMRs associated with AP-1-regulated genes, suggesting there may be a shift towards open chromatin, which is intriguing as AP-1 is composed of FOS and cJUN–2 known transcriptional regulators essential for memory formation (Tischmeyer and Grimm, 1999).

The sex-specific results we identified in this study emphasized the divergence of methylation patterns seen throughout the genome. While both sexes exhibited an overall decrease in average methylation totals in aging, the general retention of methylation at more features in females than in males carried over to comparison of DMRs, where female hypermethylated DMRs were significantly overrepresented in more contexts, including 3’UTRs, CGIs and exons, than males. Sex differences were also observed in locations of specific regulatory factor motif binding, where male hyper-and hypomethylated DMRs were largely lacking in motif associations and females presented distinct changes in transcription factor functionality categories from hyper- to hypomethylated DMRs. Our study aimed to characterize the DNA methylation patterns in the dorsal hippocampus in young and aging mice of both sexes, and overall, we observed both previously identified aspects and novel age- and sex-specific findings. It will also be important to determine the sex-specific aspects of DNA methylation variability and clocks, especially in the context of neurodegeneration in which women exhibit higher incidence and more rapid cognitive decline (Jahanshad and Thompson, 2017) Ferretti et al., 2018; Conrin et al., 2018; Ritchie et al., 2018). Our sex-specific analyses provide new insight on regional- and genomic feature-specific evidence of contrasting methylation patterns in the dorsal hippocampus.