Many chronological aging clocks can be found throughout the epigenome: Implications for quantifying biological aging

The magnitude of epigenetic clock site (ECS) methylation changes across the lifespan is significantly smaller than the largest age‐related changes within non‐clock sites, and smaller still than known differences in DNA methylation between tissues (Figure 1). Looking at age‐related sites by a variety of metrics, we see that age‐predictive sites by linear regression and non‐linear mutual information regression follow the same pattern of small magnitude changes found in clock sites. We then analyze the highest weighted locus by each metric (e.g., epigenetic clock site in red) (Figure 1b). Even the most “age‐related” locus, as identified by multiple methods, shows a small magnitude difference over the lifespan. This finding leads to many questions. How do such small age‐related changes occur consistently across individuals and tissues composed of multiple cell types? Are these small differences occurring in a restricted set of loci, or are there large numbers of these sites from which a variety of clocks can be constructed? Perhaps methylation values reach an asymptote and regress toward the young values?

A logical starting point in determining the number of possible ECSs is to recreate Horvath's original clock. We independently collected as much of the original data that was publicly available from NCBI Gene Expression Omnibus (GEO), obtaining ~75% of the training data. We selected a set of loci whose methylation state can predict age using the same preprocessing, including imputation, normalization, and elastic net modeling.

We replicated the results of the original Horvath model (Figure 2a, 2b). We also trained a model on 9699 samples of Illumina 450k methylation array data deposited in NCBI GEO (Figure 2C), which allowed us to use data from more sites (~450,000) than the 21,369 sites used in the Horvath report. We only included samples from experiments with over 100 samples, as we found that including smaller experiments quickly diminishes clock performance. Our full 450k model outperformed the original Horvath model, but this is expected as it had access to more features in both the training and test sets. While some data from Horvath's original paper were unavailable and/or lacked age data that were publicly available, we demonstrated roughly equivalent performance across the three models with mainly different sets of sites in each model (Figure 2d).

Because clock models trained on almost identical data can select different loci as the most predictive set and perform equivalently, we examined how many loci could be used as clock sites via an iterative “knockout” approach. After using elastic net regression to identify the most informative loci, these sites were then removed and new clock models trained on the remaining sites. We fit 5 new clock models after every removal and removed any locus used by any of the 5 clocks. By iteratively removing the most predictive sites every round, subsequent models should become increasingly hindered in their performance. Interestingly, rather than a rapid depletion of predicted performance, we saw a gradual linear decrease in prediction accuracy (as defined by Pearson correlation coefficient between predicted and actual age, Figure 2e) and rise in model error (Figure 2e) until about 20% of sites have been removed from the training data. These data illustrate the breadth of CpGs that are strong age predictors from the sites measured by the Illumina 450k microarray.

DNA methylation has known roles in specifying tissue identity (Koh & Rao, 2013; Macaluso & Giordano, 2004) via differentially suppressing and activating specific regions of the genome in a cell type‐specific manner. Our own data and previous reports (Horvath, 2013; Lowe et al., 2016; Thompson et al., 2018) on pan‐tissue epigenetic clocks demonstrate that chronological age‐predictive models work equally well across most tissues when trained on many tissues. However, clocks trained on blood alone have poor performance at predicting tissue‐specific impairments, such as cognitive decline(Starnawska et al., 2017).

This raises the question of whether pan‐tissue clocks either still capture enough tissue‐specific/relevant changes by inclusion of multiple tissues in the training set, or if pan‐tissue clocks are missing tissue‐relevant changes, and potentially sacrificing insight into tissue‐specific aging and disease.

To look for tissue‐specific age association, we used an ordinary least squares regression to model the effects of age, tissue, and their interactions. 39823 loci were significantly (q<0.05) associated with tissue, while 9587 were associated with age. Of the 6226 sites with significant main effects of both tissue and age, 3939 had a significant interaction effect, suggesting a tissue‐specific aging change. The rest of the age‐associated sites had tissue‐independent effects (as defined by non‐significant tissue:age interactions), suggesting common age‐related changes across tissues (Figure 3a, c). Looking closer at the loci with significant age main effects, the tissue‐specific aging interactions are weak (Figure 3b). This suggests that there is some common process across tissues that causes these changes regardless of whether it is a regulated program or a systemic form of epigenomic entropy. Dysregulations of epigenetic machinery have been implicated as the driver of the relationship between DNA methylation and aging(Bell et al., 2019), but our prior work in mouse found no changes in the direct methylation machinery with aging(Hadad et al., 2016).

Observing that most loci with a main effect of age were similar across tissues (except for saliva) was surprising in light of the strong relationship between methylation and tissue identity. This leads us to question whether elastic net would identify ideal pan‐tissue predictors as the best sites to predict age even in a single tissue, or if they would be too weak of a signal compared to tissue‐specific signals. When we trained clocks on data from only one tissue at a time, they perform markedly worse on predicting age in other tissues (Figure 3d and e). This indicates that the tissue‐specific signals are in fact strong enough to drown out the non‐tissue‐specific signals, which are more abundant. Notably, saliva sample ages were uniquely well‐predicted by all other tissue models, but poor predictors of other tissue ages. When comparing the coefficients of loci identified by elastic net for multiple tissues, we do see moderate associations (r > 0.5) for all tissues—indicating that some loci are being chosen that covary even by tissue‐specific models (Figure 3f). The loci chosen by these tissue‐specific clocks are largely uniformly distributed across genomic features, but are notably enriched near GTEX (Lonsdale et al., 2013) eQTLs. Perhaps the enrichment of these loci near regions where point mutations are sufficient to affect disease‐related gene expression is indicative of a potential biological link, but determining whether these changes are causal, compensatory, or silent in the aging process will require further studies in epigenome editing/manipulation.

We have made two key observations regarding methylation aging, namely 1) non‐linear site selection outperforms linear for training epigenetic clocks (Figure [Link], [Link], [Link]) and 2) age‐predictive loci have very small changes over the lifespan (Figure 1). This raised the question of how sites that would need to change by a fraction of a percent per year can consistently predict age. These small changes are observed even in a pan‐tissue context. The changes are much larger than we would expect from senescent cells, which are relatively rare and appear at different rates between tissues (Tuttle et al., 2020), and also are unlikely explained by some change in blood constituent cells (Chen et al., 2016) as they would vary wildly in their abundance with vascularity of the tissue, such as between whole blood and saliva. To answer this question, we performed an experiment to see if the most predictive loci from full age range models are identifiable in windows within the lifespan (Figure 4a). We trained three epigenetic clocks separately using data from three groups: young (25–50), middle (50–75), and aged (75–100) samples, then tested them on each age group. We found sites most predictive of aging in younger ages (25–50) were poor predictors later in life (75–100). Based on our observations of the linear and non‐linear relationships, we then looked at how clock sites from both Horvath and our own clocks change throughout the lifespan. We found that 99.5% of loci exhibit a distribution with at least one inflection point with aging that could be regressed into an overall linear trend (Figure 4d, Figure [Link], [Link], [Link]), alongside loci where locally fitted regression matched the sites’ linear fits (Figure 4e). While other reports (e.g., ‐ (Marioni et al., 2019)) have noted a non‐linear change with age, the parabolic distribution we observe appears to be novel. These changes may be attributable to changes in variance with aging (Slieker et al., 2016), but they describe an overall trend to increasing variance that would not necessarily lead to a parabolic distribution. Since these parabolic sites are also selected by epigenetic clocks on the basis of an extrapolated linear trend, further exploration of the trajectory of methylation aging may yield understanding of how clock predictions would respond to aging interventions and diseases. Similar to tissue‐specific clocks, clocks trained on discrete sections of the lifespan are largely uniformly distributed with an enrichment near eQTLs (Figure 4f).

To understand the biological relevance of different loci that are predictive of or related to aging, we performed genomic enrichment analyses on three sets of methylation loci. These sets were as follows: 1) 6666 loci we determined that change with age by linear regression, 2) 2624 that serve as predictive clock sites, and 3) 79921 sites related to aging as measured by mutual information. When exploring the 6666 age‐related sites identified by traditional statistical analysis (linear regression) we found strong enrichments for genomic regions with limited known biological function for DNA methylation such as intergenic regions and sites outside of CpG islands, with under‐representation in regions such as promoters and CpG islands where DNA methylation is understood to regulate gene expression/genomic accessibility (Figure 5). The results were similar but weaker for age‐related sites selected using mutual information. In contrast, epigenetic clock sites were previously reported as enriched in VISTA (Visel et al., 2006) enhancers and near GTEX (Lonsdale et al., 2013) eQTLs. We attempted to further interrogate the activity state of promoters near age‐predictive cytosines, but it is impossible to generate a meaningful multi‐tissue prediction of activity with current data. Performing the enrichment against a target tissue of interest did not yield any significant results (Figure [Link], [Link], [Link]). They were also enriched in open chromatin regions as defined by aggregated DNase hypersensitivity data. However, the regulatory element enrichment appears to be driven by enhancers as aggregate regulatory elements, because subsets of enhancers, super enhancers, and repressors show neither enrichment nor depletion. All methods of identifying age‐related methylation sites found they were depleted in CpG island bodies and gene promoters while being enriched in intergenic regions.

We then analyzed the clock sites’ relationship to DNA‐binding proteins by comparing the probe sequence used to detect methylation loci to DNA binding motifs using HOMER. We attempted to interrogate the enrichments with respect to the “targeted” cytosines but found no significant results after multiple testing (Supplement [Link], [Link], [Link]). The top 3 most significant sequence enrichments for known DNA‐binding protein motifs were all TEAD proteins, which are enhancer binding proteins that aid in the initiation of transcription and have been linked to cellular senescence (Xie et al., 2013) and age‐related disease (Tsika et al., 2010). These enhancer enrichments were previously reported (Bell et al., 2019), but in their analysis were under the threshold for significance after adjusting for the array background. Nonetheless, these provide a potential mechanism to further explore the biological impact of altered clock site methylation.

If epigenetic age acceleration is truly a biomarker of aging, we would expect samples from patients with age‐related diseases and conditions tend to have higher predicted ages than healthy samples. To this end, we annotated 1767 samples for a variety of conditions including multiple sclerosis, obesity, Alzheimer's disease, and smoking status (Figure 6). We compared measures of “age‐acceleration” in these groups using both our own chronological 450k clock and the published PhenoAge(Levine et al., 2018) “biological” clock. PhenoAge predicted most samples were age‐decelerated, with multiple sclerosis, obesity, and depression being significantly “accelerated” compared to controls (Figure 6a). Meanwhile, smoking status and HIV decreasing age acceleration seem contradictory to prior reports (Esteban‐Cantos et al., 2021; Levine et al., 2018). Since the PhenoAge clock is regressed against 10‐year mortality risk coerced into units of years, perhaps it is not surprising that even “healthy” control samples are predicted to be age‐accelerated. In contrast, our chronological 450k clock had much closer to zero average age acceleration. Notable exceptions are the significant increase in “age acceleration” of patients with NAFLD or NASH (Figure 6b). In both cases, these analyses have an additional uncertainty term in their age predictions that make interpreting age acceleration values difficult, with most samples falling within the range of model error. A potential explanation for age acceleration predictions going the opposite direction one might expect is found in looking at trajectory of aging methylation values over the lifespan. Namely, some show a non‐linear fit that is “hidden” to the linear elastic net regression (Figures S4 and S5). While not all epigenetic clock sites exhibit this non‐linear relationship, we noted in our primary feature selection that selecting sites with the best linear relationship reduced performance much more than selecting sites by non‐linear feature selection (Figure [Link], [Link], [Link]). This may explain why some age‐related diseases appear to be “decelerated” due to the same methylation values being present earlier and later in life.

Raw data were collected from NCBI’s Gene Expression Omnibus (GEO). For replicating Horvath's experiments, individual datasets were extracted and manually curated from metadata using geoquery (Davis & Meltzer, 2007). Of these, only 20 datasets were available and matched the given age distributions and sample numbers as reported in Horvath's original report (Horvath, 2013). These data make up dataset 1, used for the direct Horvath replication (Figure 2b, d) and the iterative trimming models (Figure 2e). Our second dataset consists of all publically available data from the Illumina 450k Human Methylation BeadChip Array that could be annotated by our label extraction program ALE (Giles et al., 2017), with accompanying metadata as presented in geometadb (Zhu et al., 2008). Sample labels for all samples’ sex, age, and tissue were extracted from text using our previously published tool ALE (Giles et al., 2017). We then exclude any samples with annotated ages under 25 for two reasons. First, we see aging as a process that begins post‐development, with human development ending around 25 years of age. Second, our annotation model is often incorrect about units, and as a result, our confidence in ages <25 is much lower than the rest of the lifespan due to common age labels with units of weeks/months falling in this range. Disease annotations were hand curated by reading their metadata and the metadata of the associated GSEs. A table of annotations can be found in Supplemental Data [Link], [Link], [Link].

The replication dataset for Horvath's clock model was preprocessed using code from Horvath 2013 to ensure identical normalization and imputation. We also used the same age transformation as reported in Horvath 2013(Horvath, 2013). The full 450k data were imputed in a similar pipeline, using KNN imputation in sets of 120,000 probes with subsequent normalization. For our linear models, batch effect correction was performed using ComBat (Leek et al., 2012) to control for experiment ID after removing samples which did not fall in the beta distribution.

Due to the large number of samples and probes included in the full 450k dataset, modeling could not be performed on the full dataset simultaneously. As such, we tested many feature selection pipelines in the smaller Horvath subset to determine their ability to preserve useful sites for age prediction (Figure,,). This led us to using mutual information as the primary feature selection method, from which we selected the top ~20% of sites and using those 79931 loci for all downstream analyses. [Link] [Link] [Link]

F‐regression and mutual information were performed using sklearn's(Fabian Pedregosa et al., 2011) implementation on one quarter of all loci at a time due to memory constraints. Top 10% delta mean and variance were computed by comparing the mean and variance differences for each locus between the young (25–35) and aged (65–100) sets of samples.

Elastic net regression was utilized to generate clock models, as in Horvath 2013. Horvath's replication was performed using R’s glmnet(Friedman et al., 2009), while the full 450k model was performed using sci‐kit learn's ElasticNetCV (Fabian Pedregosa et al., 2011). In both cases, 10‐fold cross‐validation modeling was used with the L1 vs L2 selection parameter at 0.5 (elastic net), and the selectivity parameter adjusted based on cross‐validation.

Age‐related locus analysis was performed using OLS regression in the model methylation ~age + tissue +age:tissue. After FDR multiple testing correction, 676 loci were left for downstream analyses. PhenoAge predictions were determined as described in Levine, et. al. (Levine et al., 2018).

Identifying gene, gene‐related, and CpG island‐related enrichments were performed by using bedtools to annotate sets of loci with features from UCSC. Feature enrichments were computed based on hypergeometric tests comparing the various groups with a 10% FDR correction for multiple testing.

Motif enrichments were performed using both HOMER (Heinz et al., 2010) and MEME (Machanick & Bailey, 2011). Fasta files were constructed with the 50 bp sequence of each probe on the Illumina 450k Array for all loci, and age‐related loci from mutual information, elastic net, and OLS regression models. These loci were then analyzed for known and de novo motifs, with enrichments calculated against the background of all sites on the array.

Identifying enrichments in genomic features were based on data from UCSC genome browser. CpG islands and genic regions were defined using UCSC’s annotated island bodies and genes, with promoters being 2kb upstream of the TSS and shore/shelves being defined as 2kb blocks up and downstream of the island body. Simultaneous comparison of all regulatory elements was performed using ORegAnno regulatory elements (Griffith et al., 2007). Vista enhancers were used to look specifically at known human enhancers (Visel et al., 2006). eQTLs were taken from GTEx (Lonsdale et al., 2013). Chromatin density was inferred using DNase hypersensitivity data from ENCODE DNase‐seq (Consortium & E.P., 2004).

The datasets analyzed in the current study are available in NCBI GEO. These datasets are all from GPL13534↗ and can be obtained from https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL13534↗.