New targeted approaches for epigenetic age predictions

Abstract

Background

Age-associated DNA methylation changes provide a promising biomarker for the aging process. While genome-wide DNA methylation profiles enable robust age-predictors by integration of many age-associated CG dinucleotides (CpGs), there are various alternative approaches for targeted measurements at specific CpGs that better support standardized and cost-effective high-throughput analysis.

Results

In this study, we utilized 4647 Illumina BeadChip profiles of blood to select CpG sites that facilitate reliable age-predictions based on pyrosequencing. We demonstrate that the precision of DNA methylation measurements can be further increased with droplet digital PCR (ddPCR). In comparison, bisulfite barcoded amplicon sequencing (BBA-seq) gave slightly lower correlation between chronological age and DNA methylation at individual CpGs, while the age-predictions were overall relatively accurate. Furthermore, BBA-seq data revealed that the correlation of methylation levels with age at neighboring CpG sites follows a bell-shaped curve, often associated with a CTCF binding site. We demonstrate that within individual BBA-seq reads the DNA methylation at neighboring CpGs is not coherently modified, but reveals a stochastic pattern. Based on this, we have developed a new approach for epigenetic age predictions based on the binary sequel of methylated and non-methylated sites in individual reads, which reflects heterogeneity in epigenetic aging within a sample.

Conclusion

Targeted DNA methylation analysis at few age-associated CpGs by pyrosequencing, BBA-seq, and particularly ddPCR enables high precision of epigenetic age-predictions. Furthermore, we demonstrate that the stochastic evolution of age-associated DNA methylation patterns in BBA-seq data enables epigenetic clocks for individual DNA strands.

Introduction

During aging, DNA methylation (DNAm) is continuously lost or gained at specific CG nucleotides (CpG sites) of our genome. Conversely, the DNAm levels at multiple CpG sites can be combined to estimate age, and these models are often referred to as “epigenetic clocks” [1]. Epigenetic clocks raise hopes as a biomarker in forensic medicine, to determine the donor age of an unknown specimen or of a person with an allegedly unknown age [2]. On the other hand, accelerated epigenetic aging has been shown to be associated with shorter life expectancy [3,4,5,6,7], and it is liable to be affected by environmental exposure, gender, specific mutations, and diseases [5, 8,9,10]. Therefore, epigenetic clocks seem to reflect aspects of biological age, which opens perspectives as a surrogate for intervention studies. It is even conceivable that therapeutic regimen in future medicine will rather be stratified by epigenetic age than chronological age. To translate epigenetic biomarkers into an approved medical test, it is advantageous to select a manageable set of informative genomic regions, which can be targeted by DNA methylation assays that are sufficiently fast, cheap, robust, and widely available for clinical diagnostics [11, 12].

Initially, models for epigenetic age predictions were based on Illumina BeadChip data [13, 14]. This microarray platform enables cost-effective and relatively precise DNAm measurements at single-base resolution. In contrast, whole genome bisulfite sequencing (WGBS) or reduced representation bisulfite sequencing (RRBS) does not always cover the same CpG sites and a limited number of reads may entail lower precision of DNAm levels [15, 16]. Another advantage of the Illumina BeadChip technology is that a multitude of publicly available datasets can easily be integrated into the analysis. Various different epigenetic age-predictors have been described that consider up to several hundreds of CpGs [17,18,19]. The most commonly used age predictor for multiple different tissues is based on 353 CpGs to facilitate epigenetic age predictions with a median error of 3.6 years [19].

For targeted analysis of specific age-associated CpGs, various studies described epigenetic age-predictors based on bisulfite pyrosequencing [18, 20] or by the Sequenom’s EpiTYPER assay [21]. Single Base Primer Extension Assay (SNaPshot) was also used by many laboratories for epigenetic age prediction, but the accuracy is apparently lower [22, 23]. Recently, droplet digital PCR (ddPCR) was reported to enable precise DNAm measurements [24, 25], and hence it might facilitate epigenetic age predictions without PCR bias. Barcoded bisulfite amplicon sequencing (BBA-seq), which is based on next generation sequencing, enables multiplexed analysis of PCR amplicons [26, 27]. The strength of BBA-seq lies within the parallelization of multiple DNA sequences on one lane with relatively long amplicons (up to 500 bases), and with very high coverage (usually > 1000 fold) [28, 29]. So far, only few groups described ddPCR [30] and BBA-seq [31, 32] for epigenetic age predictions and a direct comparison of these methods is still elusive.

It is largely unclear, how age-associated DNAm is regulated and if it is functionally relevant, per se. Transcription factors (TFs) or long non-coding RNAs (lncRNAs) might target epigenetic writers, such as DNA methyltransferases (DNMTs) or ten-eleven translocation family enzymes (TETs), to specific sites in the genome [33]. This process may also involve alternative splicing of DNMTs [34]. CCCTC-binding factor (CTCF), as an insulator TFs that is involved in chromatin architecture, has also been shown to be methylation sensitive [35]. If age-associated DNAm was directly regulated by epigenetic writers, it would be anticipated that the DNAm pattern of neighboring CpGs on the same DNA strand is coherently modified. Alternatively, it has been proposed that age-associated DNAm is evoked by “epigenetic drift,” which may occur due to stochastic accumulation of errors, e.g., in copying DNAm patterns during cell replication [36,37,38]. In this case, DNAm on individual DNA strands might rather follow stochastic patterns.

In this study, we have further optimized and directly compared epigenetic age predictors based on pyrosequencing, ddPCR, and BBA-seq of specific age-associated regions. Furthermore, our data indicate that the correlation of age-associated DNAm with chronological age peaks close to CTCF binding sites. Age-associated DNAm is not coherently modified on individual DNA strands and this enabled alternative single-read age-predictors that reflect heterogeneity in epigenetic aging within a specimen.

Results

Selection of age-associated CpGs for blood

To select age-associated CpGs for further comparison, we used a training set of 973 DNAm profiles of healthy human blood samples (1 to 101 years old), which are derived from seven different studies and based on the 450 k Illumina BeadChip platform (Additional file 1: Figure S1A). For further analysis, we excluded CpGs that were associated with a single-nucleotide polymorphism (SNP) or reported to be cross-reactive [39, 40]. To reduce the impact of leukocyte composition [41], we excluded CpGs with high variation in DNAm across hematopoietic subsets (R > 0.02 in six different cellular subsets; Additional file 1: Figure S1B) [42]. Furthermore, we excluded CpGs on sex chromosomes, CpGs that are significantly affected by smoking, and those which were no more comprised on the new EPIC Illumina BeadChip microarray. Thus, 313,728 CpGs were used for further selection of age-associated CpGs (Additional file 1: Figure S1C) [39, 40, 43, 44]. To identify candidates that are best suited for targeted DNAm analysis, we selected CpGs by linear correlation with chronological age using Pearson correlation > 0.5 or < − 0.5. Age-associated DNAm might also follow a non-linear monotonic function and therefore we alternatively used Spearman’s rank correlation > 0.5 or < − 0.5. Since age-associated DNAm was shown to rather follow a logarithmic pattern, particularly in children [19, 45], we also selected CpGs that correlated with logarithmic age across all samples (Pearson’s correlation > 0.5 or < − 0.5). It needs to be taken into account that there is notoriously variation between different Illumina BeadChip datasets due to modifications in sample preparation, handling, and measurement in different studies (Additional file 1: Figure S1D-G). Therefore, we have only considered CpGs that met the corresponding criterion (R > 0.5 or < − 0.5) in each of the seven individual studies of the training set. Overall, 65 CpGs passed at least one of these three filter criteria—and despite the different filter criteria, the overlap was remarkably high—with 18 CpGs passing all three thresholds (Fig. 1a). Subsequently, we trained a multivariable model based on the 65 CpGs using functions of transformed chronological age with logarithmic dependence for pediatric and linear dependence for adult donors, as previously described by Horvath et al. [19]. This model provided a high correlation with chronological age in the training set (R² = 0.95; median error = 3.0 years; Fig. 1b) and in an independent validation set of 3674 blood samples of five different studies (R² = 0.82; median error = 3.3 years; Fig. 1c). Moreover, the accuracy of our 65 CpG model even outperformed the accuracy of age-predictions using the epigenetic aging signatures with 353 CpGs [19] or 71 CpGs [17], often referred to as Horvath and Hannum clocks, respectively. It even outperformed the more recent skin and blood clock of Horvath et al. [46]; however, the later performed slightly better if we only considered datasets that provided DNAm values for all 391 relevant CpGs (Additional file 1: Figure S2 and S3) [17, 19, 46].

Fig. 1

Selection of age-associated CpGs and targeted analysis with pyrosequencing. a Illumina BeadChip profiles of 973 blood samples of 7 studies (all 450 k) were used to select age-associated CpGs by Pearson’s correlation (blue), Spearman’s rank correlation (red), or Pearson’s correlation after logarithmic transformation of age (green). Sixty-five CpGs passed at least one of these thresholds (correlation coefficient: R > 0.5 or R < − 0.5) and the Venn diagram depicts a very high overlap. b, c A multivariable linear model for these 65 age-related CpGs revealed high correlation with chronological age in the training set (b; n = 973), and in an independent validation set (c; n = 3674; color code corresponds to different studies as indicated in Supplementary Fig. S3A). d, e Epigenetic age prediction based on pyrosequencing of 6 CpGs (d) and 9 CpGs (e) in blood samples of the training set (n = 40; blue) and an independent validation set (n = 40; red)

Epigenetic age predictor based on pyrosequencing of specific CpGs

To further narrow down the best-suited CpGs for targeted DNAm measurement, we selected among the 65 age-related CpGs those with the best linear correlation and highest slope in Illumina BeadChip training sets (Additional file 1: Figure S4A). The selected CpGs were associated with the genes of Elongation Of Very Long Chain Fatty Acids Protein 2 (ELOVL2), which is well known to reveal high age-association [21], Coiled-Coil Domain-Containing Protein 102B (CCDC102B), Four And A Half LIM Domains Protein 2 (FHL2), Immunoglobulin Superfamily Member 11 (IGSF11), Collagen Type I Alpha 1 Chain (COL1A1), and MEIS1 Antisense RNA 3 (MEIS1-AS3). When we analyzed DNAm at these CpGs in 40 blood samples by pyrosequencing, we observed high correlation with chronological age, which was overall higher than for the independent validation set of 450 k Illumina BeadChip profiles (Additional file 1: Figure S4B). A multivariable linear regression model based on pyrosequencing of 40 samples revealed a good correlation with chronological age in an independent validation set of 40 blood samples (R² = 0.86; Fig. 1d). The correlation was even increased, when we integrated three additional CpGs, which were selected in our previous work (associated with the genes ITGA2B, ASPA, and PDE4C [18]), into the model (R² = 0.89; Fig. 1e). However, measurement of the independent validation set, which was measured several months later, revealed a small systematic offset at several CpGs (COL1A1, FHL2, and IGSF11; Additional file 1: Figure S4B), and therefore the median error was considerably larger in the validation set (median error = 6.8 years with 6 CpG model and 5.2 years with 9 CpG model). We could not identify the influencing parameters that evoked this small offset, but pyrosequencing appears to be sensitive to small technical variations due to possible PCR bias of methylated versus non-methylated sequences.

Analysis of age-associated DNAm with droplet digital PCR

Droplet digital PCR (ddPCR) is based on randomly distributing a single sample into small droplets, which are then processed for PCR amplification separately. We anticipated that this technology reduces PCR bias, because sequences are amplified with different efficiency and therefore affect the ratio in conventional PCR. For each amplicon two probes were designed that either detect the methylated or the non-methylated sequence and DNAm level was subsequently determined by Poisson distribution (Fig. 2a). Reliable ddPCR assays could be established for the sequences of CCDC102B, COL1A1, MEIS1-AS3, FHL2, PDE4C, ASPA, and IGSF11, while this was hampered for ITGA2B and ELOVL2 due to neighboring CpG sites. The ddPCR results revealed a clear correlation with chronological age at all tested CpGs, which was often even slightly higher than for pyrosequencing (Table 1; Fig. 2b–h). We then generated a multivariable model based on ddPCR measurement of these seven CpGs. This provided relatively reliable age-predictions in an independent validation set of 40 blood samples (R² = 0.89; Fig. 2i). The median error of these ddPCR age-predictions was only 2.9 years, which was significantly lower than for pyrosequencing (P = 0.00005), indicating that age-predictions in independent sets of samples might be more accurate with ddPCR.

Fig. 2

Age-associated DNA methylation measurements with droplet digital PCR. a Two-dimensional amplitude analysis of duplex ddPCR (blue: positive droplets for methylated CCDC102B; green: positive droplets for non-methylated CCDC102B; orange: double-positive droplets; black: negative droplets). b–h DNAm measurements by ddPCR in the training set (n = 40; blue) and independent validation set (n = 40; red). i Epigenetic age prediction based on ddPCR measurements of 7 CpGs in blood samples of the training and validation set

Table 1 Correlation of DNA methylation and chronological age with different targeted approaches

Epigenetic age-predictions with bisulfite barcoded amplicon sequencing

Subsequently, we analyzed DNAm of the nine age-associated regions with deep sequencing technology (Illumina MiSeq). Bisulfite barcoded amplicon sequencing (BBA-seq) was initially performed on a training set of 38 blood samples. Overall, DNAm levels of BBA-seq and pyrosequencing revealed good correlation (R² = 0.92; Additional file 1: Figure S5), indicating that DNAm measurements were feasible for all amplicons. Furthermore, the correlation of DNAm levels in BBA-seq with chronological age was in a similar range as observed for pyrosequencing or ddPCR (Table 1; Additional file 1: Figure S6). A multivariable linear regression model was then established for the nine CpGs with the highest correlation with chronological age per amplicon. This approach provided high accuracy for age-predictions in the training set (R² = 0.95; median error = 2.8 years). Unexpectedly, the accuracy of epigenetic age-predictions was even higher in the independent validation set of 39 blood samples (R² = 0.87; median error = 2.4 years; Fig. 3a). Alternatively, we used all CpGs within the nine amplicons for Lasso (Fig. 3b) and elastic net regression models (Fig. 3c) with 10-fold cross-validation. Both machine learning approaches performed better than the linear regression model on the training set, but the precision was lower on the validation set, which might be attributed to small technical offsets between different BBA-seq runs or overfitting due to the relatively small sample size.

Fig. 3

Epigenetic age predictions using BBA-seq. a–c Nine amplicons with age-associated CpGs were analyzed by bisulfite barcoded amplicon sequencing (BBA-seq) in a training set of 38 blood samples (blue) and an independent validation set of 39 samples (red). Age predictions were based on a multivariable linear model of 9 CpGs within 9 amplicons (a), Lasso regression model of 17 CpGs within 8 amplicons (b), or elastic net regression model of 26 CpGs within 8 amplicons (c). d, e In analogy, genomic DNA was isolated from buccal swabs and analyzed by BBA-seq (training set: n = 46, blue; independent validation set: n = 49, red) with a multivariable model of 9 CpGs (d), by Lasso regression model of 27 CpGs within 7 amplicons (e), and by elastic net regression model of 26 CpGs within 7 amplicons (f)

We have then analyzed if the method is also applicable for buccal swab samples, which are a widely used specimen in legal medicine. Due to the different cellular composition, the multivariable model for the 9 CpGs was retrained to provide a good correlation in a training set of 46 buccal swab samples (R² = 0.93; median error = 4.0 years), albeit the correlation remained lower in the independent validation set of 49 buccal swab samples (R² = 0.75; median error = 6.9 years; Fig. 3d). The accuracy could be slightly increased by Lasso (Fig. 3e) and elastic net algorithms (Fig. 3f) that were generated based on all CpGs of the amplicons, but it remained lower than for blood samples. This might be due to the heterogeneous composition of buccal epithelial cells and leukocytes in buccal swab samples [47].

Age-associated DNA methylation changes peak close to CTCF binding sites

One of the advantages of BBA-seq is that amplicons are longer than in pyrosequencing and include measurement of more neighboring CpGs. We have compared the correlation of neighboring CpGs with chronological age, particularly in the amplicons of ELOVL2, PDE4C, and FHL2, which harbor 36, 26, and 18 CpGs, respectively. Plotting of correlation coefficients against the genomic locations revealed curvy distributions (Fig. 4a). Similar distributions, but less distinct, were also observed for the BBA-seq data of buccal swabs (Additional file 1: Figure S7). Notably, age-associated DNAm changes at these three age-associated regions peaked in vicinity of binding sites of CCCTC-binding factor (CTCF), which is involved in organization of chromatin structure. Furthermore, chromatin immune precipitation (ChIP)-seq data of human embryonic stem cells (hESC; GSM822297), K562 (GSM822311), and A549 cell lines (GSM822289) indicated that CTCF binds close to the peak of age-associated DNAm changes (Fig. 4b). When we analyzed enrichment of our previously identified 65 age-associated CpGs within ChIP-Seq read peaks of CTCF, they were significantly enriched in each of the three ChIP-seq experiments (Fig. 4c). These data may suggest that regulation of age-associated DNAm changes is related to CTCF binding and/or the three-dimensional chromatin conformation [48].

Fig. 4

Age-associated DNAm changes peak close to CTCF binding sites. a Pearson’s correlation of age with DNAm levels of CpGs within the amplicons of ELOVL2, PDE4C, and FHL2 are plotted for the blood samples of the training set (n = 38, blue) and validation set (n = 39, red). X-axis represents the position of CpGs within the amplicons. b Enrichment of CTCF binding at the position of these amplicons (gray shaded region) was then analyzed in chromatin immune precipitation (ChIP) sequencing data of hESC (GSM822297), K562 (GSM822311), and A549 cells (GSM822289). Peak heights were automatically trimmed by IGV tool (indicated in brackets). The positions of predicted CTCF binding motives are also presented. c Boxplot of normalized read counts of CTCF ChIP-seq data of A549, hESC, and K562 cell lines either at the 65 age-associated CpGs or at 1000 randomly chosen CpGs from 450 k BeadChip array. Read counts from ChIP-seq data were normalized by quantile normalization and analyzed within a window of 500 base pairs (p value was estimated by Mann-Whitney rank test)

Analysis of age-associated DNAm patterns within individual BBA-seq reads

Age-associated DNAm peaks at specific sites in the genome, but it remained unclear whether neighboring CpGs are coherently modified or not. If DNAm was regulated by a targeted protein complex, it would be expected that neighboring CpGs are conjointly modified on individual strands. To address this question, we analyzed the DNAm pattern of individual reads in BBA-seq, which reflect the binary code of methylated and non-methylated cytosines in individual DNA strands. In fact, the DNAm patterns at the age-associated regions were very heterogeneous (Fig. 5a). When we simulated random DNAm patterns based on the likeliness of DNAm at specific CpGs at a given age, the patterns were similar to the experimental results (Fig. 5b). The correlation in DNAm at individual CpGs was overall very low (Fig. 5c). Thus, age-associated DNAm at neighboring CpGs is apparently not coherently modified and rather seems to be acquired in a stochastic manner.

Fig. 5

Analysis of age-associated DNAm patterns within individual BBA-seq reads. a Heat map to exemplarily depicts frequencies of DNAm patterns within the 36 neighboring CpGs of the ELOVL2 amplicon in BBA-seq data of a young (21 years old) and old sample (72 years old). b For comparison, heatmaps are presented based on random simulation of DNAm patterns under the assumption that DNAm at neighboring CpGs occurs entirely independent (simulations correspond to 21 and 72 year old donors). c Pearson’s correlation of DNAm levels between neighboring CpG sites within ELOVL2 amplicon (BBA-seq data of training set). d For each BBA-seq read of ELOVL2 training set, we estimated the epigenetic age based on the binary sequel of methylated and non-methylated CpGs. The plot depicts relative read count of every donor in the training set that were classified for predicted ages between 0 and 200 years (relative read count normalized by read count per sample). e–g The mean age-predictions based on individual BBA-seq reads were determined for each sample and then plotted against the chronological age of the samples of the training (blue, n = 38) and validation set (red, n = 39). This analysis was performed independently for the amplicons of ELOVL2 (e), PDE4C (f), and FHL2 (g)

We then reasoned that DNAm patterns within individual reads might also be used for age-predictions. To this end, we developed a mathematical model based on the BBA-seq data by assigning each DNAm pattern the most likely corresponding age (between 0 and 200 years; Fig. 5d). As anticipated for younger donors, the model revealed a higher number of young read predictions, whereas older donors had more reads that were predicted older. Notably, the mean of strand-specific age-predictions correlated well with the chronological age of the donors in the training and validation sets (Fig. 5e–g). In fact, in the validation set, the correlation of the mean of single read predictions with chronological age for ELOVL2, PDE4C, and FHL2 was even better than the correlation of the mean DNA methylation levels at individual CpGs in amplicons of the BBA-seq data (Table 1). This supports the notion that single-read age-predictors are applicable and that epigenetic clocks tick independently within cells of the same sample.

Discussion

In the advent of new technologies for DNAm measurements, there is a continuous need to revisit, optimize, and validate epigenetic clocks. In this study, we used the linear correlation of age-associated DNAm with chronological age as a proxy for the precision of DNAm measurements. While all methods tested revealed good or even very good correlations with chronological age, they all have their advantages and limitations for the application in epigenetic clocks.

So far, most epigenetic clocks for human samples were derived from Illumina BeadChip datasets [15]. It provides unparalleled opportunity to measure DNAm at single-base resolution, albeit not all CpGs of the genome are represented on these platforms (for instance about 1.7% of the human CpGs are covered by the 450 k BeadChip). Our analysis clearly demonstrated that DNAm levels in 450 k Illumina BeadChip data are highly correlated with chronological age at various specific CpGs. For development of our targeted assays, we focused particularly on blood samples and filtered for CpGs with low variation between leukocyte subsets to reduce the impact of the cellular composition. Our predictor with 65 age-associated CpGs provides similar or even better accuracy than other commonly used predictors [17,18,19]. Furthermore, our signature is also applicable for the current EPIC BeadChip version. While normalization regimen, integration of more CpGs, and machine learning algorithms can further enhance age-predictions, the goal of our selection was to identify suitable CpGs for the subsequent targeted approaches. Focusing on a smaller number of CpGs in targeted assays is a tradeoff between the applicability and accuracy [15]. Targeted analysis is usually faster, more cost-effective, and better applicable for laboratories that do not have immediate access to Illumina BeadChip technology. Furthermore, targeted assays are less dependent on availability of specific BeadChip versions, which is important given the long timeframes needed to develop a standardized procedure towards certification according to in vitro diagnostics regulations (IVDR).

Despite the remarkable linear correlation of age-associated DNAm with age, there is evidence for rather logarithmic association, particularly in childhood [19, 45, 49]. We selected candidate CpGs either for linear correlation (Pearson’s correlation), continuous non-linear association (Spearman’s rank-order), or by linear correlation with the logarithm of age. To our surprise, there was a very high overlap of selected CpGs with the three filter criteria, indicating that age-associated DNAm changes are generally accelerated in early life and are then rather linearly acquired in adulthood and the elderly. In fact, age-predictions for pediatric cohorts were clearly improved by age-transformation with a logarithmic adjustment, as elegantly described by Horvath [19], and this approach should therefore also be considered for targeted epigenetic age predictors if applied to pediatric cohorts.

Bisulfite pyrosequencing is currently the most popular targeted method for epigenetic age predictions in forensics [50]. This method is relatively simple and it has been shown to have a very high precision in DNAm measurement [11]. The accuracy of our epigenetic age-predictions with pyrosequencing was in a similar range as for Illumina BeadChip models. Nevertheless, the conventional PCR reaction before pyrosequencing can evoke amplification bias for methylated or non-methylated strands [50].

Droplet digital PCR might reduce this technical PCR amplification bias, because the individual droplets are either scored as positive or negative independent of the PCR efficacy. So far, only few studies used ddPCR for DNAm analysis [51, 52], and only one recent study used it for measurement of age-associated DNAm changes in a pediatric cohort [30]. In this study, we performed a direct comparison of epigenetic age-predictions of ddPCR and pyrosequencing, and ddPCR indeed provided significantly lower median error of age-predictions in the validation set. However, fluorescent probes could not be designed for all sequences, particularly if many neighboring CpGs were located in the target region. Overall, DNA methylation measurements with ddPCR appear to be more precise than pyrosequencing and this should be further validated by interlaboratory comparison.

Bisulfite barcoded amplicon sequencing makes use of the technical advances in massive parallel sequencing. Recently, Naue and coworkers described a similar age predictor using massive parallel sequencing [32]. Our comparative approach substantiates the notion that BBA-seq is a powerful method for epigenetic clocks, particularly if multiple samples need to be analyzed in parallel. The correlation of DNAm measurements with chronological age at individual CpGs was similar in BBA-seq data as compared to pyrosequencing or ddPCR. In analogy to pyrosequencing, we observed a systematic off-set between the BBA-seq results of different sequencing runs, which might be attributed to PCR bias. On the other hand, the long BBA-seq amplicons provided insight into DNAm patterns on the same DNA strand. We demonstrate that the correlation of DNAm at neighboring CpGs with chronological age follows a bell-shaped curve. Notably, this correlation often peaked close to CTCF binding sites, which resembles one of the best characterized architectural proteins for the 3D chromatin conformation [48]. This is in line with previous studies indicating that age-associated DNAm changes are enriched at CTCF binding sites [53, 54]. Furthermore, it has been suggested that age-associated hypomethylation, but not hypermethylation, is related with CTCF binding sites across various human tissues [55, 56]—albeit we observed this particularly at the hypermethylated regions in ELOVL2, PDE4C, and FHL2. Another genome scale study supported the notion that CTCF enrichment occurred both in age-associated hyper- and hypomethylated regions [54]. Occupancy with CTCF is tissue-specific and it is influenced by DNAm [57, 58]. In turn, CTCF was also capable to change DNAm status, for instance, by reducing the DNAm level at the vicinity of its binding positions [59]. Thus, the relevance of CTCF-binding for age-associated DNAm changes needs to be further explored.

Furthermore, our analysis of DNAm patterns in individual reads of BBA-seq demonstrated that neighboring CpGs are modified rather independently. We have recently described similar findings for DNAm changes during long-term culture of cells in vitro [29, 60]. While the stochastic changes at neighboring CpGs challenge the view of directed regulation of age-associated DNAm, they may support the notion that this process is evoked by “epigenetic drift,” possibly caused by changes in chromatin conformation. The finding of stochastic DNAm changes provided also the basis for our mathematical approach of epigenetic age predictions for individual BBA-seq reads. It is yet unclear if epigenetic aging is accelerated synchronously at different genomic regions within the same cell. Single Cell DNAm analysis has been described using WGBS [61], but the coverage at individual CpGs (including age-associated CpGs) is notoriously very low, and thus analysis of epigenetic clocks is hampered in such datasets. On the other hand, there is evidence that age-associated DNAm is coherently modified in cancer, which possibly reflects the epigenetic state of the tumor-initiating cell [62]. Thus, the epigenetic age predictions based on individual BBA-seq reads seems to reflect heterogeneity of epigenetic aging within a sample.

Conclusions

Our comparative approach demonstrates that targeted analysis of age-associated DNAm via pyrosequencing, ddPCR, and BBA-seq enables similar correlation with chronological age as described for larger signatures on Illumina BeadChip profiles. In the independent validation set, the highest accuracy of epigenetic age-predictions was reached with ddPCR, which might be due to lower PCR bias. Thus, this relatively new method should be considered for further round robin tests to validate the precision and accuracy for epigenetic age-predictions in clinical and forensic application. Furthermore, our analysis of BBA-seq data demonstrated that age-associated DNAm peaks close to specific sites in the genome, particularly at CTCF binding sites. The finding that DNA methylation at neighboring CpGs is not coherently modified opened the perspective for our alternative single read predictors. The mean of strand-specific age-predictions correlated very well with chronological age. This method will gain additional power with longer reads, as provided by nanopore sequencing, and it will shed additional light into the heterogeneity of cellular aging within a sample.

Methods

Sample collection

Peripheral blood samples of healthy donors (n = 80) were obtained from the Department of Transfusion Medicine at the University Hospital of RWTH Aachen. Buccal swab samples were collected with Mastaswab MD555 (Mast Group ltd., Reinfeld, Germany) at the Institute for Legal Medicine of the Heinrich Heine University in Düsseldorf, Germany (n = 95). This study has been performed according to the guidelines approved by the local ethics committees of RWTH Aachen University (EK 041/15) and Heinrich Heine University of Düsseldorf (Permit number 4939).

Selection of age-associated CpG sites

In total, 4647 DNAm profiles of human blood samples of 12 different studies (all analyzed on the HumanMethylation450 [450 k] BeadChip platform; no samples with known malignancies) were retrieved from Gene Expression Omnibus Database (GEO; Additional file 1: Table S1) [17, 45, 63,64,65,66,67,68,69,70,71]. To further select age-associated CpGs for targeted analysis, we excluded (i) 105,454 SNP-associated (provided by the updated ‘MASK_general’ subset) [40] and 29,233 cross-reactive CpGs [39]; (ii) 26,426 CpGs with high variation across the different hematopoietic subsets in (GSE35069; variances > 0.02) [42]; (iii) 2901 CpGs that are significantly affected by smoking [43, 44]; (iv) 11,648 CpGs of sex chromosomes; and (v) 31,396 CpGs that were not included in the new Illumina BeadChip EPIC platform to facilitate better comparison with future data. The remaining 313,728 CpGs were then filtered for their correlation with chronological age. It has been demonstrated that most age-associated DNA methylation changes can be accurately modeled as a function of logarithmic age, particularly in children [45]. Therefore, we selected CpGs by their correlation with the logarithm of chronological age (cutoffs for all comparisons: R < − 0.5 or > 0.5). Alternatively, we used Pearson’s correlation or Spearman’s correlation with chronological age (also R < − 0.5 or > 0.5), since those approaches were used in other studies before and we anticipated that selection for linear and non-linear DNA methylation changes would provide complementary subsets of CpGs. Due to variation in beta-values between different datasets the corresponding criterion (Pearson’s log age, Pearson’s age, or Spearman’s age; R < − 0.5 or > 0.5) needed to be fulfilled in each of the seven different datasets of the training set, individually.

Before establishing the multivariable model for epigenetic age predictions based on the 65 age-associated CpGs, we used a transformed age instead of chronological age, as described before [19]. In brief, adult.age was assigned as 20 years, F function was used for age transformation as follows:

$$ F\ \left(\mathrm{age}\right)=\log\ \left(\mathrm{age}+1\right)-\log\ \left(\mathrm{adult}.\mathrm{age}+1\right),\mathrm{if}\ \mathrm{age}\le \mathrm{adult}.\mathrm{age}. $$

$$ F\ \left(\mathrm{age}\right)=\left(\mathrm{age}-\mathrm{adult}.\mathrm{age}\right)/\left(\mathrm{adult}.\mathrm{age}+1\right),\mathrm{if}\ \mathrm{age}>\mathrm{adult}.\mathrm{age}. $$

Subsequently, we established multivariable regression model (F (transformed age)) based on transformed age on the training sets from Illumina BeadChip (Additional file 1: Table S2) [19]. After estimation of transformed age, the results need to be inversely transformed for epigenetic age-predictions:

$$ \mathrm{DNAm}\ \mathrm{age}=\mathrm{inverse}.F\ \left(\mathrm{transformed}\ \mathrm{age}\right)=\left(1+\mathrm{adult}.\mathrm{age}\right)\times \exp\ \left(\mathrm{transformed}\ \mathrm{age}\right)-1,\mathrm{if}\ \mathrm{transformed}\ \mathrm{age}<0 $$

$$ \mathrm{DNAm}\ \mathrm{age}=\mathrm{inverse}.F\ \left(\mathrm{transformed}\ \mathrm{age}\right)=\left(1+\mathrm{adult}.\mathrm{age}\right)\times \mathrm{transformed}\ \mathrm{age}+\mathrm{adult}.\mathrm{age},\mathrm{if}\ \mathrm{transformed}\ \mathrm{age}\ge 0 $$

To investigate batch effects of the different datasets of the training and validation microarrays, we performed simple mean imputation and principal component analysis with the python package scikit-learn. We used all available CpGs of the datasets within the training (385,587 CpGs) and validation (462,889 CpGs) groups (only 1000 samples of GSE55763 were included for the validation datasets) or used the 65 CpGs of our age predictor.

Isolation of genomic DNA and bisulfite conversion

Genomic DNA was isolated from 50 μl blood with the QIAamp DNA Mini Kit (Qiagen, Hilden, Germany), or from buccal swab with NucleoSpin Tissue Kit (Macherey-Nagel, Düren, Germany). DNA was quantified with a Nanodrop 2000 Spectrophotometer (Thermo Scientific, Wilmington, USA) and 200 ng genomic DNA was bisulfite converted with the EZ DNA Methylation Kit (Zymo Research, Irvine, USA).

Pyrosequencing

Bisulfite converted DNA was used for PCR amplification by PyroMark PCR Kit (Qiagen, Hilden, Germany). Twenty micrograms of PCR products was immobilized to 5 μl Streptavidin Sepharose High Performance Bead (GE Healthcare, Piscataway, NJ, USA), and subsequently annealed to 1 μl sequencing primer (5 μM) for 2 min at 80 °C. PCR and pyrosequencing primers (Metabion, Planegg-Martinsried, Germany) are provided in Additional file 1: Figure S8 and Additional file 1: Table S3. Pyrosequencing was performed on PyroMark Q96 ID System and analyzed with PyroMark Q CpG software (Qiagen). To estimate epigenetic age, we either used a multivariable model based on six CpGs (Additional file 1: Table S4), or a nine CpG model that also considered the three CpGs of our previous work [18] (Additional file 1: Table S5).

Droplet digital PCR (ddPCR)

Droplet digital PCR was performed with a QX200™ Droplet Digital™ PCR System (Bio-Rad, CA, USA). Primers and dual-labeled probes were designed by Primer3Plus software (Additional file 1: Table S6). The reaction mixture consisted of 10 μl of 2X ddPCR Supermix (no dUTP; Bio-Rad), 1 μM of the forward and reverse primers, 250 nM of the probes targeting the methylated and unmethylated DNA sequences, and 25 μg of bisulfite converted DNA in a final volume of 20 μl. Together with 70 μl of droplet generation oil, it was then subjected into a DG8 disposable droplet generation cartridge (Bio-Rad). The water-in-oil droplets were produced by QX200 Droplet Generator (Bio-Rad). Forty microliters of the generated droplets were transferred to the ddPCR 96-well plate (Bio-Rad). After heat sealing with the PX1 PCR Plate Sealer (Bio-Rad), the plate was placed in the C1000 Touch Thermal Cycler (Bio-Rad) for PCR run. The thermal cycling conditions were 95 °C for 10 min, followed by 40 cycles of 94 °C for 30 s and 1 min (2.5 °C/s ramp rate) at 58 °C (FHL2 and PDE4C), 54 °C (COL1A1) or 53 °C (CCDC102B, MEIS1-AS3, ASPA, and IGSF11) with a 10 min step at 98 °C for enzyme deactivation and a final hold at 4 °C. Subsequently, the plates were read on the QX200 droplet reader (Bio-Rad) and data were analyzed by QuantaSoft 1.7.4 software (Bio-Rad). The percentage methylation of each reaction was calculated by Poisson statistics according to the fraction of positive droplets for methylated and non-methylated probes. The multivariable model for ddPCR is provided in Additional file 1: Table S7.

Bisulfite barcoded amplicon sequencing (BBA-seq)

Target sequences with candidate CpG sites were amplified by PyroMark PCR kit (Qiagen) (Additional file 1: Figure S9). The forward and reverse primers contain handle sequences for the subsequent barcoding step (Additional file 1: Tables S8). PCR conditions were 95 °C for 15 min; 45 cycles of 94 °C for 30 s, 58 °C for 30s, 72 °C for 30s; and then final elongation 72 °C for 10 min. The amplicons of each donor (e.g., of 9 different CpGs) were pooled at equal concentrations, quantified with Qubit (Invitrogen), and cleaned up with paramagnetic beads from Agencourt AMPure PCR Purification system (Beckman Coulter). Four microliters of PCR products were subsequently added to 21 μl PyroMark Master Mix (Qiagen) containing 10 pmol of barcoded primers (adapted from NEXTflexTM 16S V1-V3 Amplicon Seq Kit, Bioo Scientific, Austin, USA) for a second PCR (95 °C for 15 min; 16 cycles of 95 °C for 30 s, 60 °C for 30s, 72 °C for 30s; final elongation 72 °C for 10 min). PCR products were again quantified by Qubit Kit (Invitrogen), combined in equimolar ratios, and cleaned by Select-a-Size DNA Clean & Concentrator Kit (Zymo Research, USA). Ten picometers DNA library was diluted with 15% PhiX spike-in control and eventually subjected to 250 bp pair-end sequencing on a MiSeq lane (Illunima, CA, USA) using Miseq reagent V2 Nano kit (Illumina). Two samples of the training set and one sample of the validation set were not successfully amplified and therefore not considered for further analysis.

To estimate DNAm levels for each CpG based on BBA-seq data, we used the Bismark tool [72]. The average number of reads per sample and genomic region was approximately 25,000 and only sequences that occurred at least 10 times were further considered. Multivariable models for epigenetic age predictions based on those CpGs that revealed highest correlation with chronological age per amplicon are provided in Additional file 1: Table S9 and S10 (for blood and buccal swabs, respectively). Alternative age prediction models were generated by machine learning as described before [73]. In brief, we applied a penalized lasso and elastic net regression model from glmnet R package on the training set from BBA-seq data. The best-fitted candidate CpGs and model was chosen by 10-fold cross-validation on the training set (Additional file 1: Tables S11 - S14).

Association with CTCF binding sites

Chromatin immune precipitation sequencing data (ChIP-seq) for CTCF in A549 cell lines (GSM822289), H1 human embryonic stem cells (GSM822297), and K562 cell lines (GSM822311) were analyzed. Each CpG site was extended for 250 bp in both directions and quantile normalized read counts of the ChIP-Seq experiments were compared for the 65 age-associated CpG sites in comparison to 1000 randomly chosen CpGs from the 450 k array. Enrichment was estimated by Mann-Whitney rank test. The CTCF binding motifs around the target 65 age-associated CpG sites were predicted by RGT-motif analysis (www.regulatory-genomics.org/motif-analysis) with default parameters.

Simulation of stochastic DNAm patterns

We simulated randomly generated DNAm patterns under the assumption that methylation at neighboring CpGs occurred independently. The probability that the CpG site i of gene X is methylated in the generated patterns is given by a linear function f_X,S(i), which was based on correlation of chronological age versus DNAm at i in the training set. For the simulations we used the function random from the python 3 library random.

Epigenetic age predictions for individual BBA-seq reads

The following algorithm was developed to estimate epigenetic age based on the binary sequel of methylated and non-methylated CpGs within individual reads of BBA-seq data: Let X be a gene with n_X CpG sites. Using the training set, we estimated the probability for each CpG to be methylated at a given age. Using linear regression, we approximated the methylation frequency of a site i of gene X as a function of age a. This yields a set of n_X linear functions F_X,i (1 ≤ i ≤ n_X), where F_X,i(a) approximates the methylation frequency of site i of gene X at age a. Since the functions F_X,i are linear, they approach infinity or minus infinity if a approaches infinity or minus infinity. Since methylation frequencies assume only values between zero and one, we defined a set of n_X functions p_X,i (1 ≤ i ≤ n_X) as follows: p_X,i(a) = F_X,i(a) if 0 ≤ F_X,i(a) ≤ 1; if F_X,i(a) < 0, then p_X,i(a) = 0; if F_X,i(a) > 1, then p_X,i(a) = 1. We interpreted p_X,i(a) as the probability that site i of gene X is methylated in a donor of age a. For a given methylation pattern P, we could then calculate the probability Pr(P,a) that pattern P comes from a donor of age a. For this we assumed that methylation of different sites occurs independent from each other. Formally, we obtain Pr(P,a) = q₁∙ … ∙q_nx, where q_i = p_X,i(a) if site i is methylated in pattern P and q_i = 1 − p_X,i(a) otherwise. To each detected pattern, we assigned the age a_P that maximizes Pr(P,a) (if the maximum is not unique we use the average of the ages that maximize Pr(P,a)). The estimated age of the donor corresponded to the average of the a_P, where we included a given value a_P multiple times if the pattern P is detected multiple times. To practically determine the a_P, we calculated Pr(P,a) for a between 0 and 200 years in steps of 1 year.