Pan‐Epigenetic Age Prediction in Mammals

To characterize epigenetic change during aging, we collected histone modification (ChIP‐seq) and DNA methylation (whole‐genome bisulfite sequencing; WGBS) profiles from nine previous studies (Bujold et al. 2016; ENCODE Project Consortium et al. 2020; Fernández et al. 2016; Hillje et al. 2022; Meer et al. 2018; Petkovich et al. 2017; Signal et al. 2024; Stubbs et al. 2017; Yang, Occean, et al. 2023). These 3491 profiles were drawn from 482 humans and 523 mice with each individual donor contributing profiles from one or more of seven epigenetic layers across up to twelve tissues (average 3.5 profiles per donor; Figure S1). Each histone profile reflected the proportion of cells bearing that particular histone modification at each genomic locus. For DNAm, each profile encoded the percent of cells in the measured tissue for which each CpG was methylated. To harmonize features across studies, all measured sites within the same gene body were pooled and their signal averaged to produce a single continuous value per epigenetic layer per gene (Figure 1a; Methods).

Correlating the epigenetic signal in each gene with age, we found that certain epigenetic layers exhibited more age‐related change than others. For H3K27me3, we observed a significant association of the occupancy of this epigenetic mark with age for 77.3% of human genes, whereas for H3K9me3, a significant age association was seen for only 0.6% of genes (Figure S2a,b). Nevertheless, there was a clear coordination between the age‐related changes among epigenetic layers, with repressive marks H3K27me3, H3K9me3, and DNA methylation aligned in their direction of change, and activating marks (H3K36me3, H3K4me1) exhibiting an opposing pattern (Figure 1b,c). Unexpectedly, while H3K27ac and H3K4me3 had been largely associated with the activation of transcription, in our aging study they showed a modest but significant trend in the same direction as the repressive marks (Figure 1c). Overall, we found that for 12 of the 21 possible pairs of layers, the direction of age‐related change between the two layers was associated (Figures 1d–i and S3a–d). For instance, DNA hypomethylation co‐occurred with a loss of H3K9me3, whereas the accrual of H3K9me3 co‐occurred with DNA hypermethylation (Figure 1d). In contrast, H3K4me1 and DNAm exhibited the opposite relationship, with an increased deposition of one associating with a loss of the other (Figure 1e).

We found that the genes exhibiting the strongest associations with age were largely consistent across epigenetic layers (Figure 2a). For example, the genes for which DNA methylation state was most correlated with age were highly likely to be the same genes for which H3K9me3 occupancy was correlated with age (34‐fold enrichment, p = 5.5 × 10⁻¹²⁰; Figure 2a). This same convergence among epigenetic layers was seen in mice (Figure 2b) involving largely the same genes as in humans (Figure 2c). In humans, the strongest coordination was seen for H3K4me3 and H3K27ac, while in mice H3K4me3 and H3K27me3 were most tightly coupled, although nearly all epigenetic layers showed significant convergence on a common set of genes in both species. At a false discovery rate (FDR) of 20%, a total of 143 genes exhibited age‐related changes across all seven epigenetic layers (Figure 2d).

Genes exhibiting the greatest epigenetic repression with age—those that accumulated repressive marks and lost activating marks—were highly enriched for developmental functions related to HOX gene expression and organ formation (Figure 2e,f; Methods). Such genes included SOX1, WNT1, and HOXB5. In contrast, genes with the greatest epigenetic activation with age were related to inflammatory pathways (Figure 2f), including SIPRA and CAMP.

Taken together, these analyses (Figures 1 and 2) indicated that aging is accompanied by a coordinated shift in histone modifications and DNA methylation, with both activating and repressive marks changing together at shared genomic loci. While the direction of change varied by epigenetic mark, the loci that underwent the largest modifications during aging were conserved across species, suggestive of an intertwined "pan‐epigenetic" process that unfolds over the lifespan.

Having identified core genes at which aging reflects coordinated changes across epigenetic layers, we next asked whether these various layers are equally predictive of age. To address this question, we trained a collection of "single‐layer" clocks, each designed to estimate age using genes as input features, with values of each gene drawn from one particular epigenetic layer (Figure 3a). Of the seven layers, H3K36me3 and H3K9me3 were omitted from clock‐related analyses due to sparse coverage in mice, yielding a total of five epigenetic clock models (Figure S1c). Using a 10‐fold cross‐validation procedure to assess these models (Methods), we found that the prediction accuracy varied by epigenetic layer (lowest ρ: H3K4me1 = 0.54; highest ρ: DNAm = 0.91, Figure S4a–c). Training single‐layer clocks for increasing numbers of donors revealed that some marks had a greater age‐prediction capacity than others, even when the number of training donors was matched (Figure S4d). As the training sample size increased, H3K27me3 and DNAm showed the most rapid improvements in model performance (scaling rate: DNAm = 0.112, H3K27me3 = 0.109; Figure S4d; Methods), while H3K4me1 showed the slowest (H3K4me1 scaling rate = 0.056). On the other hand, each of these single‐layer clocks exhibited similar performance whether applied to humans or mice (Figure S4a–c), suggesting that, like DNA methylation (Lu et al. 2023; Wang et al. 2020), histone modifications follow an evolutionarily conserved trajectory of change during aging.

We next compared the results of the single‐layer clocks to a "pan‐epigenetic" clock designed to predict donor age from any of the five epigenetic layers analyzed (Figure 3a; Methods). This model had an architecture identical to the single‐layer clocks but was trained using profiles from all epigenetic layers. This pan‐epigenetic clock accurately predicted age in both humans (Spearman ρ = 0.67; Figure 3b) and mice (Spearman ρ = 0.80, Figure 3c) from any epigenetic mark (Figure 3d). Notably, the performance of the pan‐epigenetic clock closely matched that of the separate single‐layer clocks in each species (Spearman ρ = 0.64, p = 0.048; Figure 3e) and tissue (Spearman ρ = 0.95, p = 2.0 × 10⁻⁶; Figure 3f) when assessed on held‐out donors in 10‐fold cross‐validation. When we left out all profiles from one epigenetic layer during training, the model largely retained the ability to estimate the age of held‐out donors using that particular layer (Figure 3g; Methods). H3K4me1 had the smallest decrease in performance when left out (−23%), while H3K27me3 had the largest (−62%). These results demonstrated that the aging signals contained within each epigenetic layer are similar enough to be effectively represented by a single model.

If there is a coordinated process driving epigenetic change across layers, individuals appearing older/younger than their chronological ages according to one layer would likewise appear older/younger according to other layers. To test this hypothesis, we applied the pan‐epigenetic clock to donors profiled for all five epigenetic modifications (n = 109 human donors). We found that the degree of over‐ or under‐prediction of chronological age was indeed synchronized across layers (Figure 4a–d; Methods). For instance, for individuals whose H3K27ac profile indicated that they were 1 year older than their chronological age, their H3K4me3 profile produced an overprediction of 0.52 ± 0.1 years (mean ± s.d., Figure 4a). This directional synchronization of the rate of epigenetic change between modifications was not limited to associations within histone marks, with DNAm patterns also yielding age predictions that were significantly associated with those of histone marks (Figure 4b,d).

Reduced methylation age is strongly promoted by chronic caloric restriction (CR) (Wang et al. 2017), which has long been known to extend lifespan in many species (Speakman and Mitchell 2011). We thus evaluated whether this intervention would likewise slow the progression of age‐related changes across the various layers of the epigenome. Applying the pan‐epigenetic clock to tissues from C57BL/6 mice fed a CR or control diet (Hillje et al. 2022; Petkovich et al. 2017), we found that every layer tended to produce younger age predictions in CR mice, but that these reductions reached statistical significance for only DNAm, H3K4me1, and H3K27ac (Figure 4e,f). Ages inferred from DNAm and H3K27ac profiles of CR mice were on average 4.6 and 0.9 months younger than their chronological ages, respectively. Thus, CR appears to exert a youth‐preserving influence on multiple layers of the epigenome.

ChIP–seq narrowPeak files were gathered from six public datasets—ENCODE, BLUEPRINT, CEEHRC, Signal et al. (2024), Yang, Occean, et al. (2023), and Hillje et al. (2022). The locations of autosomal gene bodies in humans and mice were retrieved from the Ensembl v109 catalog (Harrison et al. 2024). Then, for each histone profile separately, the mean narrowPeak score of all peaks overlapping a given gene was taken as the epigenetic value of that gene, while genes without any overlapping peaks in that sample were assigned a value of zero. Mouse gene identifiers were converted to their one‐to‐one human orthologs via pybiomart (Smedley et al. 2009), placing every profile in a common human‐gene coordinate system. Genes with no overlapping peaks in any sample were discarded. This procedure yielded gene‐level histone modification profiles aligned between species.

Whole‐genome bisulfite sequencing files were obtained from six public datasets—ENCODE, BLUEPRINT, CEEHRC, Stubbs et al. (2017), Meer et al. (2018), and Petkovich et al. (2017). CpG sites were intersected with autosomal gene bodies defined in Ensembl v109 via pybedtools (Dale et al. 2011). For each sample, the mean methylation fraction of all CpGs falling within a given gene was recorded. Genes lacking overlapping CpGs were assigned a value of zero. Mouse gene identifiers were mapped to their one‐to‐one human orthologs using pybiomart (Smedley et al. 2009), placing every methylome on the same human‐gene coordinate system.

Gene‐level histone mark and DNA methylation datasets were first scaled separately. Within each {dataset × epigenetic‐mark} stratum: epigenetic signal values were min–max normalized to the interval [0, 1] and, for DNA methylation datasets, missing values (< 0.5% of all values) were imputed as the mean across all other samples. To mitigate inter‐study batch effects, each dataset was inverse normal transformed (INT) (McCaw et al. 2020). The epigenetic signal values of each gene were ranked across samples; these ranks were converted to uniform percentiles, and the inverse standard normal cumulative distribution function was applied to map these percentiles to z‐scores. During the training of all machine learning models, INT and mean‐imputation were instead performed within each cross‐validation fold, to prevent information leakage across folds. In each cross‐validation fold, the INT and mean‐imputation parameters from the training set were applied to the held‐out test data before prediction. After normalization, all datasets, including both histone marks and DNA methylation, were combined and genes for which epigenetic values were present in all samples were retained (n = 17,602).

Within the batch‐corrected matrix, gene‐wise associations between chronological age and epigenetic values were quantified separately for every {species × epigenetic‐mark} stratum using Spearman correlation. Ranking the correlation values by their unsigned magnitude, pairwise overlaps among the 1000 genes with the strongest age association were computed for all histone‐mark combinations within each species (Figure 1d,e). The expected overlap under random expectation was taken to be:Expected overlap=Ntop×NtopG2where Ntop = 1000 genes and G = 17,602 (the total number of genes). Fold enrichment was calculated as the observed divided by expected overlap. Statistical significance was assessed with two‐sided binomial tests followed by Bonferroni correction. Concordance between human and mouse age‐associated genes was evaluated analogously (Figure 1f).

For each pair of epigenetic layers A and B: (1) the 1000 genes with the strongest positive (age‐increasing set) or negative (age‐decreasing set) age‐association (Spearman ρ) in layer A were selected; (2) the distribution of Spearman ρ's in layer B was extracted for the genes in the age‐increasing and age‐decreasing sets of layer A; (3) the median difference in Spearman ρ's between these sets (age‐increasing—age‐decreasing) was calculated; and (4) a two‐sided Mann–Whitney U test was used to assess if there existed statistical significance of medians between the distribution of Spearman ρ values in the age‐increasing versus age‐decreasing sets. These p values were Bonferroni‐corrected across all (7 choose 2=) 21 comparisons (α = 0.05). A significant median difference signified cross‐talk between epigenetic layers, with positive values indicating that genes gaining signal in layer A with age also tend to gain signal in layer B, whereas negative values denoted opposing trends (Figures 1g–i and S3a–d).

Genes showing the strongest age‐associated epigenetic repression were defined by taking the intersection of the top 1000 most positively age‐associated genes according to each repressive mark (H3K27me3, H3K9me3, DNAm) and the top 1000 most negatively age‐associated genes according to each activating mark (H3K27ac, H3K36me3, H3K4me1, H3K4me3) in each species (n = 778 genes). The converse was done to define the most epigenetically activated genes with age (n = 427 genes). Enrichments of these genes for Gene Ontology Biological Process and Reactome pathways were measured using the Enrichr API via the gseapy (Fang et al. 2023) package, testing against the "GO_Biological_Process_2023" and "Reactome_Pathways_2024" libraries, respectively. The universe for each test comprised all genes quantified in the age–correlation analysis. For each gene set (repressed or activated), terms containing fewer than two genes or with Benjamini–Hochberg–adjusted p ≥ 0.05 were excluded. The Combined Score from Enrichr (log p × z‐score) was used as the enrichment metric.

Developmental and aging time points were aligned between humans and mice by representing the age of each donor as a percentage of the maximum lifespan of its respective species, as described by Lu et al. (2023). Briefly, each donor's age (in years) was divided by the species‐specific maximum lifespan reported in GenAge (de Magalhães et al. 2024) (120 years for humans; 4 years for mice) to obtain a unitless relative age between 0 and 1. These values were then transformed with the negative log–log transformation:Scaledage=−log−logchronologicalagemaximum lifespan

The reverse transformation was used to present values in units of years in all figures.

Epigenetic profiles originating from cell lines, sex‐specific tissues (prostate, breast or other reproductive organs), placenta and any tissue represented by < 50 libraries were excluded from both model training and evaluation. Calorically restricted samples (n = 298 mice) were likewise withheld from model fitting but retained for subsequent inference. Due to their limited representation in mice (Figure S1c), H3K36me3 and H3K9me3 profiles were omitted from all clock construction steps, leaving five epigenetic layers (H3K4me3, H3K27ac, H3K27me3, H3K4me1, and DNA‐methylation). The resulting dataset comprised 2598 profiles from 1005 unique donors.

Five independent gradient boosted decision‐tree regressors (Chen and Guestrin 2016) were fit—one for each epigenetic layer (H3K4me3, H3K27ac, H3K27me3, H3K4me1 and DNA‐methylation) to predict the age of mouse and human samples. Models used the 17,602 genes retained after layer integration as predictors (see "Integration and Normalization of Epigenetic Layers"); no additional covariates were included. Training employed stratified ten‐fold cross‐validation in which the number of profiles from donors of each age quintile was balanced across folds, while ensuring that all profiles from the same donor resided in a single fold. Default XGBoost hyper‐parameters were used, and performance was computed over the ten held‐out folds. Within each training fold, we retained only features whose association with age, quantified by Spearman's ρ in human and mouse samples separately, was stronger than |ρ| > 0.2 in both species and had the same sign in humans and mice.

The pan‐epigenetic clock was likewise an XGBoost model, trained in an identical way to the single‐layer clocks. The same feature set, feature selection, and samples were used. However, instead of samples from each epigenetic layer being trained on and predicted by separate clocks, a single model was trained over all profiles in 10‐fold CV.

To determine whether the pan‐epigenetic clock leverages shared aging information across layers, we retrained the model five times, each time leaving all profiles of one epigenetic layer out of the training set and performing 5‐fold CV on the remaining samples. Then, in each CV fold, the performance of the model was evaluated on (1) all samples from the left‐out layer and (2) the held‐out samples from the layers used in training.

Using donors from the CEEHRC, ENCODE, and BLUEPRINT datasets who had been profiled for multiple epigenetic marks (n = 219 donors), single‐layer clocks were trained on training cohorts of varying size. For each mark, we repeatedly (ten times) sampled 20–200 donors in 20‐donor increments without replacement, trained a single‐layer clock on each training set exactly as described above, and recorded the age‐prediction accuracy (Spearman ρ). The scaling rate for a given epigenetic layer was taken as the slope of a zero‐intercept ordinary‐least‐squares regression of these Spearman ρ values on the log₁₀‐transformed number of training donors. The log transformation was used as the scaling rate was nonlinear.

All analyses were performed in Python 3.10. Data analysis was conducted using Pandas 1.5.3, SciPy 1.10.0, Pingouin 0.5.3, and Statsmodels 0.13.5. Data were visualized with Seaborn 0.12.1 and Matplotlib 3.7.1.