Epigenetic drift underlies epigenetic clock signals, but displays distinct responses to lifespan interventions, development, and cellular dedifferentiation

Introduction

Changes in DNA methylation with age, a component of “epigenetic aging”, are widely observed across the tree of life. Age-associated DNA methylation patterns manifest as two general phenomena; one leading to stereotypical shifts in mean methylation levels at individual cytosines that can be modeled to predict individual age with high accuracy [1], and the other leading to increased variability or “disorder” in DNA methylation states due to the erosion of the epigenetic landscape [2–4]. These phenomena are hypothesized to be linked as average methylation values of individual cytosines are reported to drift from hyper- or hypo-methylated (e.g., ≥80%, ≤20%) states to more intermediate levels (e.g., 20–80%) with age [5]. However, the extent to which age-associated changes to the DNA methylome reflect distinct or similar underlying processes remains unresolved.

Over the last decade, dozens of epigenetic clocks have been developed for a range of taxonomic groups including humans [6], rodents [7], fish [8], birds [9], and trees [10]. Epigenetic clocks are typically constructed as linear models that predict chronological age or age-related phenotypes using mean methylation levels from a relatively small number of individual cytosines. The rate of epigenetic aging, measured as the discrepancy between chronological age and epigenetic age estimates, is associated with environmental conditions [6], life history traits (e.g., age at first menarche [11] and menopause [12]), and has become a widely used indicator of biological age and attendant disease risk [1, 13, 14]. More recently, epigenetic clocks have been applied to understanding epigenetic rejuvenation events occurring either naturally during early embryonic development or as a consequence of cellular reprogramming. Whereas epigenetic age estimates of induced pluripotent stem cells (iPSCs) are typically reset to zero [6], transient treatments with Yamanaka factors that do not fully induce dedifferentiation also reduce epigenetic age estimates and have been recognized as a promising anti-aging therapeutic avenue [15]. Kerepesi, et al. have also reported a period of epigenetic rejuvenation occurring during early development in which epigenetic age estimates decrease after conception until reaching a “ground zero” state coinciding with gastrulation [16]. Yet, age estimates derived from epigenetic clocks may not fully capture other facets of epigenetic aging, and here, we integrate multiple measures of age-associated DNA methylation patterns to examine these phenomena more broadly.

The mechanistic underpinnings of epigenetic clock signals are still unclear, but with millions of CpG dinucleotides in the genome [17], and minimal overlap of individual CpGs included across different epigenetic clocks [18], the “ticking” of epigenetic clocks is suggested to be the product of a more general epigenetic maintenance system than can be reflected by the clock sites alone [1]. Commonly referred to as epigenetic “drift”, the failure of this maintenance system has many references in the recent literature [2, 5, 19–22]. Yet, despite an abundance of reports examining age-related epigenetic drift [5, 21, 23], a consensus definition is lacking, with studies often defining drift to mirror the analytical approach employed [24]. For example, “epiallele frequency” [25], “discordance” [26], “disorder” [3], “entropy” [27], and “heterogeneity” [3] have all been used to assess epigenetic drift and reflect different analytical approaches. Perhaps the most inclusive definition of epigenetic drift is a change in the status of DNA methylation over time [21, 23]. Yet, according to this definition, even programmed changes which guide developmental processes could be considered epigenetic drift, and it is likely more useful to define epigenetic drift as a stochastic, rather than a deterministic change in methylation states. One popular approach for assessing stochastic changes in methylation is using Shannon’s Entropy [28]. Originating in information theory, this metric measures the amount of uncertainty in an occurrence or event. However, when applied to DNA methylation, Shannon’s Entropy simply reflects average methylation values (whether genome wide or at a specific CpG) and is also likely influenced by heterogeneity among cells. Heterogeneity of epigenetic patterning within cells requires analyzing single cells or in the case of bisulfite sequencing experiments, can be inferred from linked CpGs occurring on individual reads [3, 29].

Herein, we apply novel read-based strategies to resolve age-associated epigenetic disorder across the mouse genome. By considering methylation states between individual CpGs and their immediate neighbors, we directly assess epigenetic disorder and investigate its relationship to epigenetic clock signals, embryonic development, lifespan interventions, and cellular reprogramming. Borrowing from the conceptual framework of Waddington’s epigenetic landscape, we hypothesize that low levels of epigenetic disorder characterize robust epigenomic states and that gains in disorder occurring with age lead to “erosion” of this landscape [2, 22, 30–32]. We find that approximately 30% of the genome is disproportionately affected by age-related epigenetic disorder. Loci which act as predictors in conventional epigenetic clocks based on mean methylation levels appear to be enriched in regions that both accumulate and lose disorder with age, suggesting a direct link between epigenetic disorder dynamics and clock signals. We subsequently develop epigenetic clocks based on our regional disorder (RD) metric and compare age predictions with those produced using conventional epigenetic clocks and those based on entropy. Upon exploring the influences of development, lifespan interventions, and cellular dedifferentiation, we identify similarities as well as clear divergence between epigenetic clock signals based on either mean DNA methylation or regional DNA methylation disorder. Contrary to predictions based on prior studies, we find that disorder increases during early development and global levels of disorder are unaffected after cellular rejuvenation. Collectively, our findings suggest that DNA methylation disorder dynamics are a key contributor to epigenetic clock signals, yet also highlight a fundamental decoupling of disorder dynamics from canonical epigenetic aging that is likely to inform the potential of lifespan intervention strategies.

Materials and Methods

Results

Disorder in DNA methylation patterns are strongly correlated with age on a regional and global scale. Of the 249,015 regions assessed, RD was significantly correlated with age in 76,353 regions (30.7%), with RD increasing with age in 70,094 genomic regions (91.8%; Figure 1B) and decreasing with age in 6,259 genomic regions (8.2%; Figure 1B). The average RD across all regions, or global disorder, increases with chronological age according to a quadratic relationship (R² = 0.51, p < 2.2e-16; Figure 1C). Consistent with increases and decreases in RD being driven by distinct processes, regions experiencing increases in RD (cor ≥0.25) display a quadratic relationship to age (R² = 0.74, p < 2.2e-16; Figure 1D), whereas regions experiencing decreasing RD (cor ≤0.25) display a linear relation with age (R² = 0.77, p < 2.2e-16; Figure 1E).

With the exception of the Y chromosome, every chromosome incurs significant age-related accumulation of RD (Figure 1F). Given that more than 30% of the genome experiences age-associated RD, only sites with a p-value ≤ 0.05 and a correlation coefficient greater ≥0.5 (n = 4149) or ≤−0.5 (n = 286) were considered as age-associated for enrichment tests, with all other regions considered background (n = 244,580). Regions accumulating RD with age were significantly enriched in genes (p < 2.2e-16; Figure 1G) and promoters (p < 2.2e-16; Figure 1H), depleted in enhancers (p < 2.2e-16; Figure 1I), enriched in PRC2 target genes (p < 2.2e-16; Figure 1J), and depleted in both transcription factor (p = 0.00015; Figure 1K) and CTCF (p = 1.35e-12; Figure 1L) binding sites. Although enrichment scores were less robust, regions losing RD with age were significantly enriched in genes (p = 8.725e-05; Figure 1G), promoters (p = 7.073e-11; Figure 1H), and PRC2 target genes (p = 0.020; Figure 1J) and were depleted in enhancers (p = 0.030; Figure 1I). The mean CpG density did not differ in age-associated regions when compared to background (Figure 1M). Genes which accumulate disorder with age (n = 1,635) were significantly enriched in 552 different biological processes (GO:BP) with the most significant terms relating to nervous system development and differentiation (Supplementary Table 1; Figure 1N). Genes losing disorder with age (n = 197) were enriched in 14 different biological processes, with the most significant terms relating to multicellular organismal development (Supplementary Table 2; Figure 1N).

We next examined the relationship between RD and regional averages of Shannon’s Entropy, a commonly used measure of epigenetic drift. Regional entropy (RE) is calculated directly from mean methylation values, and thus has a strong relationship to regional methylation (RM), even when regional values are averaged across all individuals (Figure 2A). However, in loci where RE reaches its maximum (RE = 1, mean methylation = 50%), RD spans from fully ordered to fully disordered (RD = 0–1). We found that the relationship between average RD and average RE is best explained by a quadratic relationship (R² = 0.83, p < 2.2e-16; Figure 2B) with increasing RE generally indicating increases in RD. While the relationship between global RD and global RE does not change with age (Supplementary Figure 2), the relationship becomes increasingly variable at greater values of RE (Supplementary Figure 3). Age dependent changes to RM and RD are linked as 78.7% (n = 38,926) of regions experiencing modest age-associated RD (cor ≥ |0.25|) also incur modest age-associated RM (cor ≥ |0.25|); however, the remaining 21.3% (n = 10,531) of age-associated changes in RD do not correspond with RM, and 34.1% (n = 20,152) of age-associated RM occur independently of changes in RD (Figure 2C).

We also aimed to understand how signals underlying epigenetic clocks relate to epigenetic disorder. Interestingly, there is a clear enrichment of Petkovich epigenetic clock loci in regions which increase and decrease in RD with age (Figure 3A), with the absolute correlation coefficient of RD and age being significantly higher in Petkovich epigenetic clock regions when compared to those not included in the clock (p < 2e-16; Figure 3B). However, 37 of the 90 total clock CpGs fall into the same 200 bp genomic region. To more thoroughly resolve the relationship between epigenetic clock signals and epigenetic disorder, we built a series of epigenetic clocks based on CpG, RM, RE, and RD states. Over the 153 LOO iterations for each clock type, there was no difference in absolute error across clocks, suggesting that each methylation metric is capable of predicting chronological age with equivalent accuracy (Figure 3C). Similarly, there was no difference in the mean absolute error produced by the representative clocks (Supplementary Figure 1).

To further compare the influence of methylation context on clock composition, we assessed the overlap of loci incorporated into each clock type as well as the frequency in which they were selected (referred to as robustness). Of the LOO iterations, the CpG clocks selected 312 different regions with an average robustness of 0.11 (Figure 3D), RM clocks selected 106 different regions with an average robustness of 0.05 (Figure 3E), RE clocks selected 330 different regions with an average robustness of 0.11 (Figure 3F), and RD clocks selected 483 different regions with an average robustness of 0.13 (Figure 3G). Interestingly, the mean absolute RD correlation coefficients for age were significantly higher for CpG, RM, RE and RD clock regions when compared to non-clock regions (CpG p < 2e-16, RM p < 2e-16, RE p < 2e-16, RD p < 2e-16; Figure 3H–3K). The majority of clock sites (86.5%) were specific to each clock type; however, seven regions were selected across all clock types. Pan-clock regions are all associated with genes (Map10, Nlrp5-ps, Rasef, Rnf220, Evx2, Gm21297, and Apba1), with five (71%) regions located within promoters, and all increase in disorder with age (Figure 3L). While all four datatypes produce low errors in age prediction (Figure 3C), the discordance of chronological age with the age prediction (or “delta epigenetic age”) is highly correlated with age adjusted global disorder in all datatypes (CpG R² = 0.09, p = 7.53e-05; RM R² = 0.07, p = 0.0007, RE R² = 0.12, p = 6.89e-06, RD R² = 0.18; p = 3.00e-08; Figure 3M).

We then tested the influence of common lifespan manipulations on epigenetic age estimates across different clock types. Caloric restriction led to a reduction in age-associated RD, but the effect varied across strains. Male C57BL/6 mice fed a calorie restricted diet had significantly younger epigenetic ages when compared to controls as determined by all clock types (CpG p = 0.00043, RM p = 2.16e-07, RE p = 4.62e-06, RD p = 5.75e-08; Figure 4A). However, mean age adjusted global disorder appeared unaffected (p = 0.62; Figure 4B). Conversely, male B6D2F1 mice fed a calorie restricted diet only had significantly younger epigenetic ages as determined by the RM epigenetic clock (CpG p = 0.46, RM p = 0.018, RE p = 0.073, RD p = 0.39; Figure 4C). However, there was a slight trend for calorie restricted individuals to have greater mean age adjusted global disorder when compared to mice on a standard diet (p = 0.097; Figure 4D). Genetic interventions which extend lifespan resulted in a general decrease in epigenetic age. Snell Dwarf mice had significantly younger epigenetic ages when compared to controls according to the RD and RM clocks, but not the CpG or RE clocks (CpG p = 0.078, RM p = 0.0059, RE p = 0.23, RD p = 0.011; Figure 4E). Snell Dwarf mice also showed reduced mean age adjusted global disorder compared to control mice (p = 0.027; Figure 4F). GHR knock out also resulted in significantly younger epigenetic ages according to the RM and RE clocks, but not the CpG or RD clocks (CpG p = 0.088, RM p = 0.0080, RE p = 0.031, RD p = 0.64; Figure 4G), and no difference in age adjusted regional disorder was observed between GHRKO and control mice (p = 0.39; Figure 4H).

We next investigated the impacts of cellular dedifferentiation on epigenetic disorder by comparing DNA methylomes of iPSC cells and their differentiated precursors. A significant reduction in epigenetic age predictions after dedifferentiation was observed across all clock types except for RE (CpG p = 7.11e-08, RM p = 3.41e-06, RE p = 0.20, RD p = 3.93e-08; Figure 5A). Interestingly, there was no difference in the global disorder between kidney or lung fibroblasts when compared to their respective iPSCs (p = 0.28; Figure 5B). Given that dedifferentiation led to a reduction in epigenetic age estimates but did not affect global disorder, we sought to identify those regions in which RD is affected by dedifferentiation. Upon comparing RD across all fibroblasts and iPSCs, 26,512 regions significantly increase in RD after differentiation and 19,419 regions significantly decrease in disorder after dedifferentiation, but these regions do not disproportionately acquire age-associated RD relative to background (Figure 5C). Interestingly, the influence of dedifferentiation on RD epigenetic clock estimates are driven by differences in RD at just several clock sites (Figure 5D), with four regions contributing 35.7% of the overall effect.

Consistent with a previous report identifying an epigenetic rejuvenation event occurring during early development [16], we observed a significant decrease in epigenetic age predictions occurring between embryonic days 4 and 6 using the Stubbs CpG methylation epigenetic clock [7] (CpG p = 1.32e-08; Figure 5E) and the RM clock (RM p = 0.0019; Figure 5F). Conversely, we observe a significant increase in epigenetic age predictions during early embryonic development when using the RD and RE clocks (RD p = 3.81e-05; RE p = 1.90e-08; Figure 5G, 5H). In addition, global disorder is strongly increased during development (Figure 5I). Upon comparing RD across all samples from embryonic days 0.5–3.5 (n = 24) and embryonic days 6.5–7.5 (n = 12), 31,687 regions significantly increase in RD during development and 368 regions significantly decrease in RD during development. These regions were not significantly enriched in regions with age-associated RD (Figure 5J). Similar to our findings examining the influence of de-differentiation, the predictions of the epigenetic clock appear to be driven strongly by differences in RD at just several clock sites, with one region having a total effect size of 2.51 months, contributing 14.6% of the difference in ages between groups (Figure 5K). Given that the clock has 86 of the regions included as predictors represented, we would expect each region to contribute just 1.16% to the overall effect size.

Discussion

Epigenetic drift is broadly hypothesized to be a primary contributor to epigenetic aging. However, drift is a multifaceted phenomenon encompassing both stochastic and deterministic processes and is unlikely to be fully captured by a single metric. In this study, we report an approach for spatially resolving genomic patterns of DNA methylation disorder, which is distinct from measures of both average methylation and entropy. Age-associated changes in regional disorder (RD) are found in approximately one third of the genome and generally reflect the accumulation of disordered methylation; however, the opposite pattern is observed in a subset of regions in which DNA methylation patterns become more ordered with age. Given that epigenetic drift or disorder is thought to be driven in part by stochastic processes, we hypothesize that the directionality of changes in RD represent different functional pathways. Yet, the specific biological mechanisms that mediate losses of RD with age remain unclear. Age-associated gains and losses of RD were disproportionately observed in coding regions, promoters, and regions harboring PRC2 target genes, and age-associated increases in RD were strongly enriched in developmental genes, especially those functioning in neural development. These gains in disorder with age support the deleteriome model of aging [37], wherein small deleterious errors accumulate in the epigenome without effect until later in life, when epigenomic stability is compromised [2]. We suggest that disorder accumulates across the genome until it reaches a critical threshold – this may explain why the majority of regions across the genome are characterized by relatively low disorder (RD <0.5; data not shown). The value of this hypothetical threshold and the factors which contribute to the accumulation of disorder have the potential to explain the rate of aging and possibly maximum lifespan across species [32].

Consistent with the hypothesis that disorder in DNA methylation patterns underlies signals in conventional CpG clocks (i.e., those based on mean CpG methylation levels), we find that loci comprising clocks constructed using RD, RM, RE, and CpG contexts are all enriched for regions in which disorder changes with age, and loci which are shared across clocks are involved biological processes such as cell cycle, tumor suppression, and development. Interestingly, Evx2 is shared across all clock contexts, and is frequently reported as being age-associated across a variety of vertebrate species [38–41]. However, we identify notable distinctions and minimal overlap across clock contexts. For example, while a subset of CpGs selected as predictors in an epigenetic clock were enriched in regions with age-associated disorder, many CpG clock sites also fell into regions lacking age-associated changes in disorder. Thus, while disorder underlies some components of traditional CpG epigenetic clocks, other components may be attributed to other processes like coordinated changes in methylation or cell type composition. By contrasting the effects of caloric restriction, genetic manipulations, cellular reprogramming, and development across different clock types, we further identify both similarities and clear distinctions according to DNA methylation context and genomic scale. For example, while traditional lifespan extending treatments, such as caloric restriction, broadly affect RD epigenetic clocks, there is no observable effects on global disorder. This may be because the magnitude of the effect of caloric restriction is too small to be seen when averaged across the entire genome. Similarly, while CpG and RM clocks demonstrate a “ground zero” occurring during mid-development [16], we see an increase in RE and RD clock predicted ages during the same period, suggesting that average methylation states may not fully reflect how the DNA methylome changes throughout development. Collectively, these findings demonstrate the connections between epigenetic drift and other aspects of epigenetic aging, while also highlighting a complexity that should be considered when assessing read outs from epigenetic clocks alone.

In mice, global disorder changes with age according to a quadratic function, with decreases in disorder occurring rapidly earlier in life prior to a steady increase with age. This pattern is consistent with previous findings of a quadratic relationship between global DNA methylation entropy and age in the naked-mole rat [42]. The initial high level of global disorder suggests that development, as well as aging, may be characterized by a disorganized epigenetic landscape – possibly due to a transitionary period between methylation states. Given the dynamic nature of the DNA methylome during development [6, 43, 44], it is likely that RD metrics, like other measures of DNA methylation that provide temporal snapshots, capture this transition as high disorder. While data from embryonic samples suggest that disorder increases during early development, the trajectory of global disorder throughout development, and whether it corresponds with previous findings of an epigenetic “ground zero” during development [16], will require a more complete developmental series to fully resolve.

Age estimates derived from epigenetic clocks are ideal for predicting chronological age (i.e., forensics, conservation and management applications [45, 46]) as well as identifying the consequences of accelerated epigenetic aging (i.e., biomarkers in biomedical approaches [1]). Yet, collapsing mean methylation levels into a single value presents challenges for understanding the drivers and biological pathways responsible for epigenetic aging. Given the push towards targeted, high-throughput approaches (e.g., bead-based assays) for acquiring data on age-associated methylation [38, 47], critical biological information is missed. While CpG level resolution has been integral in developing our understanding of epigenetic aging, clocks built using regionally averaged methylation perform with similar accuracy to those trained on individual CpGs. We further demonstrate that the effect size of individual clock sites varies widely, and thus, changes in methylation states of just one or several clock loci can be misinterpreted as wholesale changes in epigenetic age. This is important especially when age estimates are compared across studies and different datasets. For example, we report that a single RD clock region accounted for nearly 10% of the difference in age estimation between fibroblasts and iPSCs. While the age predictions generated corroborate previous findings [6, 33], the inclusion (or exclusion) of this single region vastly changes our interpretation of the effects of de-differentiation on epigenetic age. Thus, using epigenetic clocks of any kind gives us a narrow, and potentially easily skewed, understanding of epigenetic aging at the genomic scale.

Overall, this study provides robust empirical evidence that epigenetic drift, as measured by epigenetic disorder, accumulates with age in non-random places of the mouse genome. Our analyses suggest that epigenetic disorder underlies aspects of traditional epigenetic clocks and highlights critical gaps in our interpretation of epigenetic aging. Although more work needs to be done to better resolve the drivers of epigenetic disorder – we provide an empirical basis for testing assumptions about this emerging phenomenon.

Supplementary Materials

Supplementary Figures

Supplementary Tables