Blood samples from DBA/2J and C57BL/6 mice were collected at the University of Ulm by submandibular bleeding (100–200 μL) of living mice or postmortem from the vena cava or from the heart, for both training and validation sets. All mice were accommodated under pathogen-free conditions. The experiments were approved by the Regierungspräsidium Tübingen, Germany.

Genomic DNA was isolated using a Qiagen QIAamp DNA Mini and Blood Mini Kit (Qiagen, Hilden, Germany) and measured with a NanoDrop2000 spectrophotometer (Thermo Scientific, Wilmington, United States). 500 ng per sample of isolated DNA was bisulfite converted and hybridized with the Infinium Mouse Methylation BeadChip (Illumina Inc., San Diego, CA, United States), from 12 C57BL/6 mice for the training set (all female), additional 12 C57BL/6 mice for the validation set (all female) and 12 DBA/2J mice (6 male, 6 female). The platform interrogates more than 285,000 methylation sites per sample at single-nucleotide resolution. The association to CpG islands, shelve, and shore regions was taken from the Illumina BeadChip annotation. Raw data from IDAT files was read and processed in R with the ENmix package, and the beta values were normalized using the ENmixD method (ENmix background correction and RELIC dye-bias normalization) (Xu et al., 2016). Probes with detection p-values > 0.01 or more than 10% of NA values were filtered out. One sample with more than 10% missing values was not considered for further analysis. Furthermore, CpGs on X and Y chromosomes were not considered for epigenetic clocks. Probe-type bias adjustment was performed with the Regression on Correlated Probes method (Niu et al., 2016). Initially, quantile normalization for ENmixD was included, which provided a similar selection of age-associated CpGs. However, since this normalization regimen masked the general shifts in age-associated DNAm, we pursued without quantile normalization. To identify significantly differentially methylated CpGs between B6 and DBA mice, the Wilcoxon rank-sum test with Bonferroni-correction (p < 0.01) for the total amount of CpGs in the Mouse Methylation BeadChip was used, after removing the probes in X and Y chromosomes, cross-reactive probes, and probes with single nucleotide polymorphism (SNP). The methylation profiles are available at the Gene Expression Omnibus (GEO accession number: GSE200527).

In this study, the initial selection of relevant CpGs was based on the slope of linear regressions of DNA methylation and chronological age to filter for CpGs that have large changes in absolute DNAm levels. Since the two mouse strains have different life expectancy, we normalized the data according to the sex and strain of the mice (Yuan et al., 2009). The CpGs were then filtered independently for B6 and DBA mice by the slope (greater than 30%). The slope of the correlation for each CpG was calculated: s l o p e C p G = ∑ ( x i − x ¯ ) ⋅ ( y i − y ¯ ) ∑ ( x i − x ¯ ) 2 Where x i is the normalized age of the sample i , x ¯ is the mean normalized age of the samples, y i is the normalized beta value of the sample i for the CpG, and y ¯ is the mean of the normalized beta values of the CpG. The overlap of the CpGs with slopes >30% in both mouse strains was 105 CpGs. Since the sample number was not sufficient to generate a multivariable model for this relatively large signature we used age-estimations for the individual CpGs with linear models using the normalized beta values (slope and intercept with age). The average of all these single-CpG estimations was considered as the final age prediction.

Homology alignments of human and murine genomes for the 105 CpGs were performed with the Ensembl genome browser with 121 nucleotide windows. A MATLAB script was used to find CpGs in the Illumina HumanMethylation450 BeadChip (450k) that were located either inside these homolog regions (categorized as “homolog”) or we determined the distance in base pairs (bp) to the center of the closest homolog region (categorized as “distance”). To further investigate age-associated changes at these CpGs the GSE40279 dataset was employed (Hannum et al., 2013). Data was analyzed with the R package GEOquery. Beta-values were used to determine the Pearson correlation with chronological age for all CpGs. The p-value of the correlations (α < 0.01) was estimated with a t-test for linear correlation with Bonferroni-correction for the total amount of CpGs on the 450k array.

Genomic DNA (about 500 ng) was bisulfite converted with the Zymo Research Group EZ DNA Methylation Kit (Zymo Research, Irvine, United States). Primers were designed with PyroMark Assay Design 2.0 software (Metabion, Planegg-Martinsried, Germany; Supplementary Table S1). One of the PCR primers was biotinylated. Age-associated regions were amplified using the PyroMark PCR Kit (Qiagen, Hilden, Germany) with manufactures instructions: 1 µL (approximately 25 ng) of the CT-converted DNA sample was used, with 12.5 µL of the Pyromark MasterMix 2X, 2.5 µL of Coral Load 10X, 0.4 µL of MgCl2, 1.3 µL of forward and reverse primers each and 6 µL of nuclease-free water (25 µL in total per tube). The PCR settings were: initial activation (95°C, 15 min), followed by 50 cycles of denaturation (95°C, 30 s) + annealing (56°C, 30 s) + extension (72°C, 30 s), and one cycle of final extension (72°C, 10 min). Pyrosequencing was then performed on the PyroMark Q48 Autoprep system using the PyroMark Q48 Advanced Reagent Kit according to the manufacturer instructions: 10 µL of amplified sample was sequenced with manual primer loading with 2 µL of sequencing primer. The results were analyzed using PyroMark Q48 Advanced software. The pyrosequencing reads covered neighboring CpGs (Aspa and region 1: 2 CpGs; Wnt3a, Prima1, and Kcns1; 3 CpGs; Tbc1d16 and Hsf4: 4 CpGs) and their position (pos) was numbered consecutively in the direction of sequencing. For the final selection of a multivariate regression model, the sklearn package in Python was used to identify the combination of 4 CpG sites from different amplicons that minimizes the MAD.

To investigate age-associated DNAm with the Mouse Methylation BeadChip, we used blood samples of 12 B6 mice (11–117 weeks old) and 12 DBA mice (6–109 weeks old) and we detected that the proportion of methylated CpGs decreases with age (Supplementary Figure S1). Pearson correlation between age and DNAm levels revealed that there was more age-associated hypomethylation than hypermethylation, and this age-associated demethylation was more pronounced in B6 than in DBA mice (Figure 1A). In B6 mice, age-associated hypermethylation was enriched in CpG islands, whereas hypomethylation was rather associated with shore, shelf, and open-sea regions. In contrast, in DBA mice there was rather age-associated hypomethylation at CpG islands, indicating that there are pronounced differences in age-associated DNAm patterns between these mouse strains.

With the aim of selecting candidate CpGs for epigenetic clocks with most prominent changes in DNAm levels in both strains, we filtered on the slope_CpG > 0.3 or < -0.3 in linear regression with age: 505 CpGs reached this threshold for B6 mice and 324 CpGs for DBA mice. 105 CpG were in the overlap of both strains (Figure 1B). As anticipated, these CpGs also revealed overall high Pearson correlation with age. When we directly compared the age-associated DNAm between B6 and DBA mice, it became evident that there are marked differences in age-associated DNAm. Several CpGs with positive age-associated DNAm in one strain revealed negative age-associated DNAm in the other strain, and vice versa (Figure 1C). Thus, there may even be antagonistic age-associated DNAm changes among mice from distinct strains.

We reasoned that a generally applicable epigenetic clock for mice should focus on the intersection of age-associated DNAm across different mouse strains. Our predictor was therefore trained for B6 mice but based on the 105 overlapping CpGs that also revealed marked age-associated DNAm changes in DBA mice (Figure 1D). The signature based on single-CpG linear predictions was trained on the above mentioned DNAm profiles of the 12 B6 mouse blood samples. The list of 105 CpGs with the slope and intercept for the predictions can be found in Supplementary Table S2. For validation, we analyzed an additional set of 12 B6 mouse blood samples (7–91 weeks old; one sample with >10% missing values was removed from the analysis). The age predictions correlated with chronological age (R ² = 0.60) with a mean average error (MAE) of 18.14 weeks and a mean absolute deviation (MAD) of 16.77 weeks. In analogy, the predictor was applied to the 12 DBA mouse blood samples that were initially considered for the selection of age-associated CpGs but not for training of the model. Despite the different background, the predictor gave relative precise estimates for the DBA cohort (R ² = 0.96; MAE = 7.11 weeks; MAD = 7.57 weeks).

The marked differences in age-related DNAm prompted us to directly compare the methylome of the two mouse strains. Therefore, we first performed a pairwise correlation of all CpGs on the Mouse Methylation BeadChip for each of the individual samples. The DNA methylation profiles of B6 and DBA correlated particularly within the strains in two clear clusters, independent of mouse age, indicating that there are marked differences in their methylome (Figure 2). Furthermore, 17,718 CpGs revealed significant differential methylation (adjusted p-value < 0.01) between B6 and DBA mice. Thus, the methylome of the two inbred mouse strains differs extensively.

In order to investigate if age-associated DNAm across B6 and DBA mice are also preserved in human, we utilized DNAm profiles of 656 human blood samples for comparison (GSE40279; 19–101 years old). When analyzing age-associated CpGs in relation to specific genomic regions we observed enrichment of hypermethylation in CpG islands (Supplementary Figure S2), as described before (Christensen et al., 2009). Homologous regions between human and mice were then identified with Ensembl. Out of the 105 overlapping CpGs in B6 and DBA, 83 were located in a region in the mouse genome with a homolog in the human genome. Only 14 CpGs in these homolog regions were also found in the HumanMethylation450 BeadChip, of which seven showed a significant age-methylation correlation in the human homolog region (Table 1; Supplementary Table S2). Furthermore, age-associated hyper- and hypomethylation was consistent in both mouse strains and humans for all of the seven selected CpGs. These genomic regions with inter-species conserved age-associated DNA methylation changes were located in the genes aspartoacylase (Aspa), WNT family member 3A (Wnt3a; cg21934230), heat shock transcription factor 4 (Hsf4) and a disintegrin and metalloproteinase with thrombospondin motifs 6 (Adamts6).

The homologous genomic regions of Aspa, Wnt3a and Hsf4 comprised several neighboring CpGs on the Mouse Methylation BeadChip and the HumanMethylation450 BeadChip. The highest correlated CpGs in these genes were at only few bases apart at the homologous sequences and they revealed similar age-associated hyper- or hypomethylation across species (Figure 3). Notably, the corresponding human region for ASPA was already selected as one of three CpGs in our previous aging signature (Weidner et al., 2014). Furthermore, Hsf4 was one of the three CpGs in our previous aging signature for mice (Han et al., 2018). In addition, two CpGs in the gene proline rich membrane anchor 1 (Prima1; cg31044702 and cg31044701) were amongst the top candidates in B6 and DBA (Supplementary Table S2) and this region was also in our previous three CpG signature (Han et al., 2018). Thus, the Mouse Methylation BeadChip identified very similar regions as previous selected from RRBS data (Petkovich et al., 2017), as well as conserved age-associated regions between mice and human.

Finally, we investigated if the top candidate CpGs that were selected with the Mouse Methylation BeadChip would also provide robust and reliable assays for targeted analysis with pyrosequencing. To this end, the list of 105 CpGs was ranked by mean R ² between B6 and DBA mice to then select “Region1” (cg42528232) that is not associated with a specific gene, Aspa (cg29748675), and TBC1 Domain Family Member 16 (Tbc1d16, cg30271979). Furthermore, the top CpGs with homologies in human were considered (Aspa, cg29748675; Wnt3a, cg29601161; and Hsf4, cg46095458), as well as Potassium Voltage-Gated Channel Subfamily S Member 1 (Kcns1), and Prima1 that were derived from our previous signature (Han et al., 2018). For all these seven regions we established a pyrosequencing assay that comprised a total of 21 CpGs. When we analyzed an independent set of 15 B6 mice (10–125 weeks; Supplementary Table S3) all CpGs revealed very high correlation with age (R ² > 0.8 for each CpG analyzed; Supplementary Figure S3).

To identify the best combination of these age-associated CpGs with regard to MAD of age-predictions, we analyzed every possible multivariate regression model with combinations of 4 CpG sites from different amplicons. The best multivariate model for age-prediction included the variables α (Aspa pos1, cg29748675), ß (Hsf4 pos3, 5 bp upstream cg46095458), γ (Wnt3a pos2, 9 positions upstream cg29601161) and δ (Prima1 pos1, cg31044702): P r e d i c t e d   a g e B 6 i n   w e e k s = 167.4533 + 1.2421 α + 0.9824 β − 1.3110 γ − 1.6088 δ

This model was validated with an independent set of B6 mice (age range 7–115 weeks; Supplementary Table S3) and the results revealed very high precision of epigenetic age predictions (R ² = 0.98, MAE = 5.03 weeks, MAD = 5.69 weeks; Figure 4A). Furthermore, the predictor gave high precision, without clear offset, in DBA mice (n = 12; R ² = 0.91, MAE = 10.7 weeks, MAD = 11.84 weeks; Figure 4B). The DNAm levels of the selected CpG sites showed similar association with age in both mouse strains with consistently very high correlations (R ² > 0.89) (Supplementary Figure S4). To determine the impact of the sex of mice, we performed a Mann-Whitney U test in the B6 training and validation sets and there was no significant impact of sex on epigenetic age-predictions (p-value = 0.33; Figure 4C).