DNA methylation-based estimator of telomere length

Training and validation data from 3 cohorts

In stage 1 of our project, we evaluated data from n=3,334 individuals for whom both LTL and Illumina methylation array data were available. These were from three different studies: Framingham Heart Study offspring cohort (FHS, N=878), Women’s Health Initiative (WHI, N=818) and Jackson Heart Study cohort (JHS, N=1638, Table 1 and Supplementary Note 1). The same laboratory measured LTL by Southern blotting of the terminal restriction fragments [25]. DNA methylation levels were measured in different labs using the Illumina Infinium methylation array platform.

An overview of the data sets is found in Table 1. These US cohorts were comprised of two ethnic groups: 41% of European ancestry and 59% of African Ancestry. The age of the individuals ranged from 22 to 93 years. The training set used for constructing DNAmTL was comprised of N=2,256 individuals from the WHI and JHS cohorts for whom LTL and DNAm data were assessed from the same blood sample (collected at the same time). Although fewer than 20% of individuals in the training set were of European ancestry, our test data demonstrated that the resulting DNAmTL estimator applied equally well to individuals of European ancestry. We used two test data sets. The first test data set involved N=1,078 individuals comprised of N=100 from the WHI, N=100 from JHS, and N=878 from the FHS cohorts. The second test data set was collected in stage 2 of our analysis: it involved N=9,815 DNA methylation samples from additional cohorts (Bogalusa, Twins UK, Lothian Birth cohorts, InCHIANTI) to evaluate correlations between LTL and numerous age-related conditions and lifestyle factors. We also evaluated DNAmTL in publicly-available data from adipose tissue (N=648 from the Twins UK study [27,28]), liver (N=85) [28,29], and monocytes (n=1264 from the Multi-Ethnic Study of Atherosclerosis) [30]. Finally, we tested DNAmTL in in vitro studies to ascertain its applicability to cultured cells and to probe the nature of DNAmTL’s association with TL. Additional details of these studies can be found in Supplementary Notes 1 and 2.

DNAmTL versus measured TL in blood and adipose tissue

We restricted the analysis to CpGs that are present on both the Illumina Infinium 450K array and the Illumina EPIC methylation array (Methods). Using the training data (n=2,256), we regressed measured LTL (mean TRFs) on blood CpG methylations using an elastic net regression model [31]. This resulted in the automatic selection of 140 CpGs whose methylation levels best-predicted LTL (Supplementary data 1). The linear regression model allows a direct prediction of TL based on DNA methylation levels. The predicted TL value, also referred to as DNAmTL, possesses the same units (kilobase) as that of mean TRF. The correlation coefficient between DNAmTL and LTL in the training data was r=0.63, which was overly optimistic, as subsequent independent validation with test data sets produced lower correlations of r>0.40 (Figure 1). Further, using 12 large validation data sets, we found that the correlation between DNAmTL and LTL ranged from r=0.38 to r=0.5 (last column of Table 1) with the exception of the Lothian Birth cohorts (where it was close to zero). The correlations between DNAmTL and LTL were not confounded by age, as was evident from the high correlations between age-adjusted DNAmTL and age-adjusted LTL (e.g. r=0.34 in FHS and r=0.43 in Bogalusa Herat Study, Supplementary Figures 1A and 2A and E). A stratified analysis showed that the correlation between DNAmTL and LTL were neither confounded by sex (Supplementary Figures 12B and C) nor ethnicity (Supplementary Figure 2D and E). The DNAmTL biomarker was robust against potential effects of pre-processing steps in the DNAm data analysis as can be seen from the diverse normalization methods used by the different cohorts (Table 1). The DNAmTL measurement was also robust across time as can be seen with the FHS cohort where the blood samples for the LTL measurement (FHS exam 6) were collected 9.3 years earlier than those used for the DNAm measurement (FHS exam 8). This time lag biased the correlation toward the null hypothesis, i.e. generated a correlation (r=0.44, Figure 1B) that was overly conservative. A separate analysis of non-blood tissue revealed a higher correlation of r=0.65 between DNAmTL and TRF-based in adipose tissue samples from the Twins UK study (Supplementary Figure 3B).

Figure 1. Measured LTL versus DNAmTL in training and test datasets. Scatter plots of DNA methylation-based telomere length (DNAmTL, x-axis) versus observed LTL measured by terminal restriction fragmentation (y-axis). DNAmTL and LTL are in units of kilobase (kb). (A) Training data. (B) Test data from the Framingham Heart Study. (C) Test data from the Women's Health Initiative (BA23 sub-study). (D) Test data from the Jackson Heart Study. Each panel reports a Pearson correlation coefficient and correlation test p-value. Table 1 reports analogous results for additional cohorts (Bogalusa, Twins UK, etc).

DNAmTL across different blood cell types

To test whether DNAmTL differs across blood types, we used sorted blood cells and peripheral blood mononuclear cells (PBMCs) from 6 men (aged between 27 and 32 years old, Methods). We observed a statistically significant difference in median DNAmTL values (p=0.0033, Supplementary Figure 4) even though they were roughly comparable: CD8+ T cells (median=8.25), CD4+ T cells (median=7.64), B cells (median DNAmTL=7.43), PBMCs (median DNAmTL=7.55).

DNAmTL versus qPCR TL in the Twins UK study

We obtained leukocyte DNAm data from 792 participants (all women) from the Twins UK study whose LTL was measured by Southern blotting (N=346) and/or quantitative polymerase chain reaction (qPCR, N=779) (Supplementary Note 1). The correlation between DNAmTL and qPCR- based LTL (r=0.39, Supplementary Figure 5A) was similar to that of LTL measured by Southern blotting (r=0.40, Supplementary Figure 5B).

DNAmTL correlates more strongly with age than TL

Although DNAmTL was developed based on LTL, it displayed substantially stronger negative correlations with age at the time of blood draw (r ~ -0.80 to -0.62) than did measured LTL, based on our test datasets (r ~ -0.40 to -0.30, Figure 2). Multivariate regression models in the test data show that LTL shortened by 0.022 kilobases per year (p=2.3E-27) after adjusting for sex, race/ethnicity and other confounders (Table 2). Analogous multivariate regression models showed that DNAmTL reduced by 0.018 kilobases per year, but this was associated with a far more significant p value (p=6.0E-125) than that of measured LTL (p=2.3E-27). Although the DNAm-based biomarkers were derived from profiles of adults (22-93 years old), the resulting DNAmTL algorithm was equally applicable to profiles from children; even to those who were younger than 13 years of age, where a strong negative correlation of r=-0.81 was observed between DNAmTL of blood and age (Supplementary Figure 6A). Such expected negative correlations with age were also seen with DNAmTL in adipose tissue (r=-0.41. Supplementary Figure 3A), liver (r=-0.71, Supplementary Figure 7A), and in (sorted) monocytes (r=-0.60, Supplementary Figure 8A).

Figure 2. Chronological age versus measured LTL and DNAmTL. Chronological age versus measured LTL (panels A, C, E, in units of kilobase [kb]) and DNAmTL (panels B, D, F, in units of kb). (A, B) Test data from the FHS. (C, D) Test data from the WHI (N=100), (E, F) Test data from the JHS (N=100). Each panel reports a Pearson correlation coefficient and correlation test p-value.

Effect of sex and ethnicity

Because age would confound any potential relationship between DNAmTL and age-related traits such as health, it would be useful to derive an age-adjusted estimate of DNAmTL (referred to as DNAmTLadjAge). We therefore regressed DNAmTL on age and the resulting raw residual was defined as DNAmTLadjAge. A negative value of DNAmTLadjAge would indicate DNAmTL that is shorter than expected based on age, while a positive value would indicate the opposite. We noted that DNAmTLadjAge is heritable (heritability $h^2$=0.46, p=4.5E-11) according to a pedigree-based polygenic model analysis in the FHS cohort (N>2000, Methods).

Women tend to exhibit longer LTL than men of the same age [32]. Similarly, our multivariate regression models revealed that age-adjusted LTL and age-adjusted DNAmTL were indeed longer in females than in males. The p-values for age-adjusted LTL measured in this study (p=2.15E-4, N=1078) and that of a previous study [32] (p=5E-3, N ~ 730) were far less significant than those for age-adjusted DNAmTL (p=1.14E-15, Table 2). Roughly, women showed longer telomere length than men by 0.1 kilobase, given by DNAm or TRF measures (Table 2). We also found longer age-adjusted DNAmTL in female compared to male liver samples (P=0.017, Supplementary Figure 7C). With regards to ethnicity, age-adjusted DNAmTL of PBMCs (Table 2) and monocytes (Supplementary Figure 8B) revealed that US population of African ancestry have a longer LTL than those of European ancestry, consistent with previous observations made with TRF-based measured TL. Once again, the association of age-adjusted DNAmTL (p=1.6E-33, N ~ 1200) with ethnicity was stronger than those seen with age-adjusted LTL measured by TRF (p=1E-4) or quantitative polymerase chain reaction (p=1E-3, N ~ 2450) [33]. Overall, these results demonstrate that DNAmTL exhibits substantially more significant associations with age, sex and ethnicity than measured LTL.

DNAmTL is often superior to measured LTL in predicting mortality and health outcomes

Next, we compared the performance of DNAmTLadjAge with age-adjusted LTL in predicting time-to-death or time-to-heart disease in the training and test datasets (N=3,334) for which both measures were available (Figure 3). We found that longer DNAmTLadjAge was significantly associated with a lower hazard ratio (HR) for time-to-death, all-cause mortality (HR=0.31 and P=6.7E-9), time-to-coronary heart disease (CHD, HR=0.55 and p=9.5E-3), and time-to-congestive heart failure (CHF, HR=0.32 and p=9.7E-4, Figure 3). In women, later age at menopause was associated with significantly higher values of DNAmTLadjAge (p=0.025). Furthermore, physical activity was also positively associated with DNAmTLadjAge (p=0.013).

Figure 3. Comparing measured LTL with DNAmTL with respect to age-related conditions. Meta-analysis forest plots for relating age-related conditions (rows) to age-adjusted LTL (left panels) and age-adjusted DNAmTL (right panels). Panels in the first row (A, B) presents meta-analysis forest plots for Cox regression models of time-to-death. Meta-analysis of Cox regression models for (C, D) time-to-coronary heart disease (CHD) and (E, F) time-to-congestive heart failure (CHF). Rows in the forest plot correspond to training and test datasets (used for developing DNAmTL) stratified by race/ethnicity Each row presents the summary statistic at a (stratified) study da.taset and reports sample size (N), number of events, P value, hazard ratio and a 95% confidence interval resulting from a Cox regression model. (G, H) Meta-analysis for the association with age at menopause. (I, J) Meta-analysis for the association with self-reported physical activity status (yes/no). (G-J) Each row (study data set) presents the summary statistic, P value, beta coefficient and a 95% confidence interval resulting from a linear (mixed) regression model. In general, an insignificant Cochran Q test p-value (denoted by Het. P) is desirable because it suggests that results do not differ significantly across the strata. However, an insignificant Q test p-value could also reflect lack of statistical power.

By comparison, the results from age-adjusted LTLs (LTLadjAge) were far less significant in predicting lifespan (HR=0.81 and p=4.7E-3 compared to HR=0.31 and p=6.7E-9 for DNAmTL, Figure 3A, B) and were not significantly associated with time-to-CHD, time-to-CHF, age at menopause and physical activity (p>0.3, Figure 3).

We want to emphasize that our comparison between DNAmTLadjAge and LTLadjAge involved the same set of individuals for whom both measures were available, i.e. each association test used the same sample size and distribution in age, sex, and ethnicity. These results show that DNAmTLadjAge outperforms LTLadjAge when it comes to predicting important health-related conditions. However, our comparative analysis was subject to a limitation: the measures of LTLadjAge and DNAmTLadjAge in the FHS cohort corresponded to two different blood samples collected at different time points. We addressed this limitation in two ways. First, we repeated the analysis by omitting the FHS data. In the resulting test data (n=100 samples from the WHI and n=100 samples from the JHS), DNAmTLadjAge continued to outperform LTLadjAge (Supplementary Table 1). Second, we compared DNAmTLadjAge with LTL in a host of additional cohorts (Bogalusa Heart Study, Twins UK, Lothian Birth Cohorts) as detailed below.

Evaluating DNAmTL in large scale validation data

In the second phase of validation, we sought to test these associations with even larger, independent data sets. In total, we analyzed N=9,875 Illumina methylation arrays from blood samples of N=9,345 individuals from 9 cohorts across 7 studies: FHS, WHI BA23, WHI EMPC, JHS, InChianti, Lothian Birth Cohorts of 1921 and 1936 (LBC), UK Twins, and Bogalusa Heart Study (BHS, Table 1, Methods, and Supplementary Note 1). Of the samples, 4,039 individuals were available with measured LTL measurements based on Southern blot or quantitative polymerase chain reaction (qPCR). The data set was comprised of three different ethnic groups: European (77%), African (15%), and Hispanic (8%) ancestries. All but one cohort (BHS) were available for mortality analysis (N=9,044 methylation arrays on 8,514 individuals) with sufficient follow-up period. The mean chronological age at the time of the blood draw was 65.6 years and the mean follow-up time (for all-cause mortality) was 11.8 years (Supplementary Table 2). Once again, DNAmTL was negatively correlated with chronological age in all cohorts with sufficient variation in chronological age ($-0.83 \le r \le -0.42$), in which we excluded the Lothian Birth Cohort studies comprising individuals with similar ages (Supplementary Figure 9).

Further analyses of these data confirmed that higher values of DNAmTLadjAge were indeed associated with longer lifespan (Figure 4). Each kilobase increase of DNAmTLadjAge was associated with a hazard ratio of 0.37 for mortality (p=2.5E-20, Figure 4A), similar to what we observed in the training and test dataset (HR=0.31, Figure 3A). Higher values of DNAmTLadjAge were also associated with longer time-to-CHD (HR=0.51 and p=6.6E-5) and longer time-to-CHF (HR=0.27 and p=3.6E-6, Figure 4 D and G), mirroring yet again the results obtained from the training and test data sets.

Figure 4. Meta-analysis forest plots for predicting time-to-death due to all-cause mortality and time-to-cardiovascular disease in independent validation data. Meta-analysis forest plot for combining Cox regression hazard ratios for time-to-death, time-to-coronary heart disease (CHD), and time-to-congestive heart failure (CHF), based on age-adjusted DNAmTL (DNAmTLadjAge). The sample sizes for the analysis were up to 9,044 methylation arrays (8,541 individuals) across 8 cohorts. Left panels, middle panels, and right panels report meta-analysis results for (1) simple Cox regression models, (2) multivariate Cox models adjusted for blood cell counts, and (3) multivariate Cox model adjusted for traditional risk factors, respectively. Each row reports the hazard ratio associated with DNAmTLadjAge. (1) The simple Cox models (left panels) were adjusted for chronological age, sex and adjusted for intra-pedigree correlation and batch effects as needed. (2) The models in the middle panels involved additional covariates: imputed blood cell counts based on DNA methylation data. (3) The models in the right panels different from those of (1) by additional demographic characteristics, psychosocial behavior, and clinical covariates (Methods). Each panel reports a meta-analysis forest plot for combining hazard ratios associated with time to event. Each row presents the summary statistic at a (stratified) study dataset and reports sample size (N), number of events, hazard ratio and a 95% confidence interval resulting from a Cox regression model. In general, an insignificant Cochran Q test p-value (denoted by Het. P) is desirable.

Two types of multivariate Cox regression models demonstrated that these associations remained significant even after adjusting for (1) blood cell counts (Figure 4 B, E, H), and (2) classical risk factors (Figure 4 C, F, I) including body mass index, educational level, alcohol intake, smoking pack-years, prior history of diabetes, prior history of cancer, and hypertension status.

We went further and evaluated DNAmTLadjAge in different strata including age (younger/older than 65 years), prevalent clinical conditions at baseline and found that DNAmTLadjAge remained a significant predictor of time-to-death in each of these strata (Supplementary Table 3), e.g. HR=0.26 for individuals aged < 65 years and HR=0.41 for older individuals aged ≥ 65 years. DNAmTLadjAge also remained a significant predictor of time-to-CHF in most strata (Supplementary Table 4) and of time-to-CHD in specific strata such older age, normal BMI, or higher education attainment (Supplementary Table 5).

Our analyses revealed that higher DNAmTLadjAge values were associated with measures of physical fitness/functioning (P=7.6E-3), and disease-free status (p=0.019), while prior history of cancer was associated with lower DNAmTLadjAge values (p=0.053, Supplementary Figure 10). Interestingly, we also found a nominally significant association (p=0.026) between age-at-menopause and DNAmTLadjAge. We found that one year later age-at-menopause was associated 0.001 kilobase longer LTL (Supplementary Figure 10A). Our cross-sectional analyses, however, do not allow determination of cause-and-effect relationships, but we note that age-at-menopause is also associated with epigenetic aging [34].

Life-style factors and clinical biomarkers

To assess the effect of life-style factors and diet on DNAmTLadjAge in blood, we meta-analyzed large data sets from the FHS and WHI cohort (N up to 6,977, Methods) including their associations with clinical measurements. Age-adjusted DNAmTL was positively correlated with plasma-based estimates of mean carotenoid levels (robust correlation r=0.10, unadjusted p=1.5E-5), beta-Cryptoxanthin (r=0.10 and p=3.8E-6) and high-density lipoprotein (HDL, r=0.04 and p=1.3E-3) (Figure 5, Supplementary Figures 11 and 12). The positive correlation between DNAmTLadjAge and self-reported measures of carotene intake (p=4.0E-3, N=766 from the Lothian Birth Cohort from the year 1936) was consistent with these findings. Positive associations with DNAmTLadjAge were also observed for self-reported measures of fruit (p=3.1E-5), vegetable (p=1.2E-3), dairy (p=7.1E-3), fish (P=0.018), and carbohydrate consumption (p=0.02). Positive correlations were also evident between DNAmTLadjAge and socio-economic factors such as level of educational attainment (p=3.3E-8) and income (p=3.1E-5). These associations held for each sex separately in the FHS cohort (Supplementary Figure 11).

Figure 5. Cross sectional associations between age-adjusted DNAmTL versus lifestyle/dietary variables. Association analysis between age-adjusted DNAmTL (DNAmTLadjAge) and 43 variables including 15 self-reported diet, 9 dietary biomarkers, 14 variables related to metabolic traits and central adiposity, and 5 life style factors, based on the meta-analysis across the FHS WHI, LBC 1921 and LBC 1936 cohort. Robust correlation coefficients (biweight midcorrelation) analysis were performed on the FHS and WHI cohort while generalized linear regression analysis adjusted for sex was performed on the LBC 1921 and 1936 cohort, respectively. For each variable, we display number of datasets, number of total subjects, the robust correlation results from the meta-analysis, FHS, and the WHI cohort and the Z statistics for the LBC respectively. The meta-analysis was based on Stouffer’s method for the majority of the variables or fixed effect models. The 2-color scale (blue to red) color-codes bicor correlation coefficients in the range [-1, 1] or Z statistics. The green color scale (light to dark) applied to unadjusted P values. Cell entry "--" denotes not available. The correlation analysis results stratified by sex using the FHS cohort are listed in Supplementary Figure 11 and stratified by ethnic group using the WHI cohort are listed in Supplementary Figure 12, respectively.

There were also features that correlated negatively with DNAmTL. Smoking was strongly associated with lower DNAmTL values in leukocytes (p=2.3E-16) and in adipose tissue (P=0.036, Supplementary Figure 3D). C-reactive protein (p=4.3E-9), triglyceride levels (p=3.2E-4) and insulin levels (p=3.5E-4) were negatively correlated with DNAmTLadjAge in both FHS and WHI cohorts. There were also negative correlations of DNAmTLadjAge (in leukocytes) with waist-to-hip ratio (p=2.4E-3), body-mass index (BMI, p=3.6E-3) and physical exercise (p=0.02).

We caution the reader that our p-values are not adjusted for multiple comparisons.

DNAmTL of leukocytes exhibits stronger association with smoking than does measured LTL

Next, we used 4,039 subjects from our seven validation cohorts to interrogate the impact of smoking on telomere shortening in leukocytes (Methods). A detailed smoking history (smoking pack-years) was known for roughly half of these individuals (N=2216). The smoking variable was based on pack-years when available otherwise based on never versus ever smoking. We adjusted the smoking variable for potential confounders (age, sex and ethnicity) of the relationship with LTL. Our large-scale meta-analysis showed that DNAmTL greatly outperformed measured LTL (Stouffer’s meta p=1.2E-17 versus meta p=0.029) with regards to association with smoking (Table 3). According to fixed effect models analysis (Methods), we found DNAmTL was shortened by 0.022±0.0002 kilobases per smoking pack-year (P=5.9E-07) while TL was not significantly lengthened by 1.3E-04±5.0E-04 kilobases per smoking pack-year (P=0.79). In an analogous analysis, we found that smokers had significantly shorter DNAmTL by 0.063±0.009 kilobases associated with a very robust P value (3.2E-13). The smokers also exhibited shorter TL (0.01±0.003 kilobases) but with far less significant P value (7.5E-04). Focusing on individual cohorts corroborated these findings, e.g. the association of smoking with DNAmTL (p=1.70E-11) was significantly greater than the association of smoking with LTL (p=6.77E-2, Table 3) in the BHS cohort. We also carried out a sensitivity analysis that compared age-adjusted LTL between current smokers and former/never smokers in the BHS and subgroups defined by sex and ethnicity. Our sensitivity analysis confirmed that smoking was indeed significantly associated with shorter age-adjusted DNAmTL in each subgroup (1.95E-8 ≤ p ≤ 4.39E-3, Supplementary Figure 13). Conversely, no significant associations were observed between LTL and smoking in these subgroups (0.12 ≤ p ≤ 0.87, Supplementary Figure 14).

Omega 3 intake is associated with longer DNAmTL

We were interested in investigating the relationship between DNAmTL and omega-3 polyunsaturated fatty acid (PUFA) supplementation in a large observational cohort study. We found omega-3 supplement intake to correlate positively with age-adjusted DNAmTL (bicor r=0.088 and p=4.4E-5), and this correlation was far more significant than that exhibited by age-adjusted LTL (bicor r=0.085 and p=0.016, Supplementary Figure 15). For DNAmTLadjAge, the effect of omega 3 supplementation was more pronounced in males (r=0.08, p=0.012) than in females (r=0.047, p=0.11). A multivariate linear mixed effects model analysis resulted in suggestive evidence (p=0.09) that omega-3 intake is associated with longer DNAmTLadjAge even after adjusting for sex, BMI, educational levels, and smoking pack year.

DNAmTL relates to imputed blood cell composition

LTL is known to correlate with the abundance of naïve CD8+ T cells and other cell types [8]. Similarly, we found DNAmTL to be significantly correlated with several quantitative measures of leukocytes that were imputed using DNAm data (Methods) such as naïve CD8+ T cells ($r=0.42$, p=2.2E-151) and exhausted CD8+ T cells ($r=-0.36$, p=3.0E-102; Supplementary Figure 16 and Supplementary Data 2). A multivariate regression model revealed that 25% of the variation of age-adjusted DNAmTL and 4.6% of the variation in age-adjusted LTL could be attributed to imputed leukocyte cell composition in the FHS test data. Overall, DNAmTL exhibited substantially stronger correlations with imputed leukocyte cell composition than did LTL (Supplementary Figure 16 and Supplementary Data 2).

DNAmTL CpGs tend to be located near telomeres and enriched with cis-mQTL

To test if the 140 CpGs of DNAmTL model were enriched at regions near telomeres, we performed hypergeometric tests based on the regions threshold up to 3 Mb (Methods). Our analysis showed that the 140 CpGs had moderate enrichment (hypergeometric p=8.1E-3, Supplementary Table 6) within sub-telomeric regions. Sensitivity analysis based on different thresholds validated the result (p=5.9E-3 ≤ P ≤ 0.039, Supplementary Table 6). Furthermore, we tested whether the 140 CpG sites overlapped with 52,916 cis methylation quantitative trait loci (cis-mQTL [35], Methods). Strikingly, 51 out of the 140 CpGs markers were known cis-mQTL (hypergeometric p=2.6E-15, Supplementary Data 1) which is consistent with the high heritability observed for DNAmTL (additive h²=0.46). Several of the 51 SNPs (from cis-mQTL) were implicated in complex traits according to GWAS data bases (Methods), e.g. rs2147904 in 1p34.2 (P ≥ 1.0E-14), rs305082 in 16q24.1 (P ≥8.0E-131) and rs945631 in 1p22.1 (P=2.0E-8) are known to be associated with granulocyte composition in blood. Further, rs945631 is also implicated in circulating phospholipid trans fatty acids at a suggestive P value=3.0E-06 (Supplementary Data 3). Bipolar disorder (P ≥ 2.0E-09) was identified when we inspected the SNPs in linkage disequilibrium with the 51 markers (Methods and Supplementary Data 3).

DNAmTL and LTL exhibits different patterns of SNP associations

Alternatively, we briefly examined SNP association of DNAmTL across 10 loci implicated in measured LTL with genome-wide significance [36–38] (Methods). We previously found that SNPs in the TERT locus to be associated with epigenetic aging rates [19]. Inspecting the overlap regions between the 10 loci and the 140 CpGs (in DNAmTL, Supplementary Data 1), cg00580497 is located nearby the right arm of the TERT locus (within 500kb) and cg02282640 is located nearby the right arm of another highlighted gene MPHOSPH6 [38]. However, using the FHS cohort, we did not find these SNPs or those at other loci to exhibit any significant associations (P<0.05) with DNAmLTLadjAge (Supplementary Table 7).

Functional annotation of CpGs implicated in DNAmTL

We analyzed the genomic locations of the 140 CpGs underlying the DNAmTL using the GREAT software tool which assigns potential biological meaning to a set of genomic locations (here CpGs) by analyzing the annotations of nearby genes. Ten gene sets were identified as significant at a stringent Bonferroni corrected significance level (p<0.05) including “Cadherin, N-terminal” that mediate cell adhesion, polarization and migration (nominal p=1.0E-7) and SODD/TNFR1 signaling pathway (nominal p=5.2E-5, Supplementary Data 4). When focusing only on the subset of 72 CpGs that have negative coefficients in the DNAmTL model (increasingly de-methylated with telomere shortening), we found 10 statistically significant gene sets at a Bonferroni corrected p<0.05 including calcium-dependent adhesion and “Cadherin, N-terminal” (p=2.8E-10). When focusing only on the 68 CpGs with positive coefficients in the DNAmTL model we identified 4 gene sets at a Bonferroni corrected p<0.05 including the vascular endothelial growth factor (VEGF) signaling pathway (p=7.4E-5). The expression of these gene sets remains to be validated, and for moment their identities do not immediately proffer obvious clues as to how they may be linked with telomere attrition, making their involvement all the more intriguing.

DNAmTL is often inferior to epigenetic clocks

Numerous lines of evidence suggest that telomere attrition and epigenetic aging are distinct cellular features that are well-associated with the aging process. As such both TL and epigenetic clocks are potential estimators of biological age [9]. To ascertain how related they are to each other, we calculated, pairwise correlations between DNAmTLadjAge and four “age-independent” measures of epigenetic age acceleration derived from (1) the pan-tissue epigenetic clock by Horvath (2013) [12] (2) the blood-based clock by Hannum et al. (2013) [11], (3) the DNAm PhenoAge estimator by Levine et al. (2018) [15], and (4) the DNAm GrimAge estimator by Lu et al. (2019) [17]. Using N=2356 samples from the FHS, we found that DNAmTLadjAge exhibited moderate negative correlations with the four epigenetic age acceleration measures ($-0.44 \le r \le -0.20$, Supplementary Figure 17). We then compared the performance between DNAmTL and epigenetic clocks in predicting health outcomes. Using the validation data (N=6,850), we found that DNAmTLadjAge rivals measures of age acceleration based on Hannum’s [11] or Horvath’s clock [12]. but was clearly inferior to age-adjusted DNAm PhenoAge [15] or DNAm GrimAge [17]. In particular, we demonstrated that AgeAccelGrim (based on DNAm GrimAge) greatly outperformed DNAmTLadjAge with regards to predicting time-to-death (Cox P=5.8E-71 for AgeAccelGrim versus p=4.1E-15 for DNAmTLadjAge), time-to-CHD (p=2.5E-23 for AgeAccelGrim versus p=6.6E-5), and time-to-CHF (p=1.8E-35 versus p=3.6E-6, Supplementary Figure 18).

DNAmTL applied to in vitro cultured cells

All the above analyses have shown that associations of several human traits and health outcomes with DNAmTL were much stronger than with LTL. What then might be the biological meaning of DNAmTL? We examined this question by monitoring DNAmTL in cultured somatic cells without and with expressed telomerase activity. Primary keratinocytes isolated from healthy human skin were grown with serial passaging upon confluence, where cell numbers were obtained, and their population doublings determined. DNA methylation profiles of cells after different population doublings were measured using Illumina EPIC array. DNAmTL was correlated with population doubling in both telomerase- negative and telomerase-positive cells, which displayed no telomere shortening (Figure 6 and Figure 7). Collectively, these findings suggest that DNAmTL records cell replication rather than TL. In the following, we provide more details. Primary keratinocytes isolated from healthy human skin were grown with serial passaging upon confluence, where cell numbers were obtained, and their population doublings determined. DNA methylation profiles of cells after different population doublings were measured using Illumina EPIC array. As would be expected, DNAmTL of keratinocytes from five heathy donors reduced in function of cumulative population doubling (Figure 6A). Accordingly, while the DNAmTL values of neonatal donors 1 to 4 were mutually comparable, that of donor 5, which are keratinocytes from a 65 years-old donor was markedly smaller, which is in keeping with established telomere biology. However, when neonatal fibroblasts and adult coronary endothelial cells (EC) were compared, the DNAmTL values of adult EC were greater than those of neonatal fibroblasts (Figure 6B). Possible tissue-specific influence on telomere length was ruled out by TRF southern blot analyses which revealed that telomeres of neonatal fibroblasts were indeed longer than those of adult endothelial cells (Figure 6C). This challenges the notion that DNAmTL is a surrogate for telomere length in non-blood tissue. If it were, then DNAmTL of primary cells that are transduced with hTERT should be augmented since this enzyme replicates and stabilizes telomeres. The results in Figure 7A show that this was not the case, as DNAmTL values of primary human fibroblasts continued to decrease even though they became immortalized, and TRF Southern blotting (Figure 7B) confirmed that their telomeres were indeed significantly increased and stabilized after the introduction of hTERT. Collectively, these studies show that DNAmTL is not the same as telomere length (especially when it comes to in vitro studies), hence indicating that it is instead, a DNA methylation-based surrogate for biological outcomes that are linked to telomere length in blood samples from adults.

Figure 6. Application of DNAmTL on in vitro keratinocytes, neonatal fibroblasts and adult coronary artery endothelial cells. Panel (A) depicts decreasing DNAmTL values of keratinocytes from five heathy donors as a function of cumulative population doubling (y-axis, in units in kilobase). Panel (B) show that DNAmTL values of neonatal fibroblasts are smaller than those of adult coronary artery endothelial cells (EC). Both cell types exhibit decreasing DNAmTL in function of cumulative population doubling (x-axis). Panel (C) is the average telomere length measurement of neonatal fibroblasts and adult endothelial cells, which revealed that the telomeres of the former are longer than those of the latter.

Figure 7. Application of DNAmTL on hTERT-transduced primary human fibroblasts. Panel (A) DNAmTL of primary neonatal fibroblasts without (Donor A and B) or with hTERT (Donor A hTERT and Donor B hTERT) transduction demonstrated linear and continued decrease in value regardless of hTERT status. This contrasts with average telomere length measurement by TRF southern blotting (Panel B), which revealed substantial increase in telomere length of primary neonatal fibroblasts that were transduced with hTERT-expression vector.

Epidemiological cohorts

To establish DNAm-based telomere length in leukocytes we used N=2256 individuals from the WHI BA23 and JHS cohort in the training process and N=1078 individuals from the FHS cohort, WHI [54,55], JHS [56] cohort in the test process, as listed in Table 1. More details of the study cohorts are described in Supplementary Note 1.

Our validation analyses involved N=9,785 Illumina Infinium 450k or EPIC 850k arrays measuring blood methylation levels in N=9,345 individuals from seven independent cohorts across nine studies: the FHS dataset (N=2,356), WHI BA23 (N=1,389), WHI EMPC (N=1972), JHS (N=209), InChianti (N=924 from 1 to 2 longitudinal measures on 484 individuals), Lothian Birth Cohort of 1921 (N=404) and 1946 (N=906), Twins UK (N=794) and Bogalusa Heart Study (N=831, Table 1, Supplementary Table 1, and Supplementary Note 1). All statistical analyses were adjusted for the correlation structure due to pedigree effects or repeated measurements as described below.

LTL measurements in training and test datasets

The same lab generated the LTL data across the 3 cohorts [25]. LTL was measured using Southern blots of the terminal restriction fragment length. After extraction, DNA was inspected for integrity, digested, resolved by gel electrophoresis, transferred to a membrane, hybridized with labeled probes and exposed to X-ray film using chemiluminescence, as previously described in [25]. The inter-assay coefficient of variation for blinded pair sets was 2.0% for the WHI, 1.4% for the JHS and 2.4% for the FHS [25].

Average telomere length measurement of in vitro cultured cells

DNA from primary human fibroblasts were extracted according to the protocol provided by Zymo Research (USA) using the mini-prep kit (Cat No: D4004). 1.5 micrograms of DNA were digested with Hinf I and RsaI prior to being resolved through a 0.8% agarose gel. DNA was transferred to Hybond N+ nylon membrane by Southern transfer, after which it was baked for 20 minutes 120^oC. The subsequent steps of this process are as described in the protocol provided by TeloTAGGG (Cat No:12209136001, Sigma Aldrich, USA). Chemiluminescence signal was captured using a Kodak Gel Documentation apparatus and quantified using Kodak quantification software. Average telomere length was ascertained using the formula described within the protocol provided.

Transduction of primary human fibroblasts with hTERT

Primary human fibroblasts were isolated from neonatal foreskin and transduced with hTERT using recombinant retroviruses according to methods previously described [49].

Estimation of surrogate DNAm based telomere length in leukocytes

We developed an estimate of LTL based only on DNA methylation levels. The estimate was established using the elastic net regression model implemented in the R package glmnet [51]. The elastic net regression model corresponds to a choice of 0.5 for the alpha parameter in the glmnet function. Ten-fold cross validation was performed in the (WHI and JHS) training data to specify the underlying tuning parameter λ. The final model was based on lambda.1se, i.e., the λ value that led to the minimum cross validated error within one standard error.

DNAmTL applied to different blood cell types

To evaluate how DNAmTL (and epigenetic clocks such as DNAmAge [12], and DNAmAgeSkinBlood [14]) differ across different blood cell types, we analyzed DNA methylation data from sorted blood cells from 6 men aged from 27 and 32 years old. DNA methylation profiles were generated from peripheral blood mononuclear cells (PBMC) and single cell types: CD4+T, CD8+T, and B cells were measured using the Illumina Infinium 450k platform.

DNA methylation data

The DNA methylation profiling was based on the Illumina Infinium HumanMethylation450K BeadChip in the FHS and WHI cohort and was based on the Illumina Infinium EPIC 850K BeadChip in the JHS cohort. To ensure future use with EPIC arrays, we focused on the subset of 450,161 CpGs that were present on both platforms. We kept the original normalization methods to ensure consistency with previous publications. The WHI BA23 were normalized using the background correction method implemented in the software GenomeStudio. WHI EMPC were normalized based on BMIQ [57] for beta-mixture quantile normalization. The JHS data were normalized using the "noob" normalization method implemented in the minfi R package [58,59].

Statistical models used in validation analysis

Our validation analysis involved several regression models. Cox regression for various censored outcomes such as time to death (all-cause mortality), time-to-CHD, time-to-CHF, and time to any cancer. Multivariate linear regression for our DNAm based measures (independent variable) associated with and number of age-related conditions (dependent variable) and physical function score, respectively. Linear regression was used to relate age at menopause (independent variable) with DNAmTL. Logistic regression analysis for binary outcomes allowing us to estimate odds ratios for the binary variables such as cancer status, hypertension, type 2 diabetes, and disease-free status. The multimorbidity index was defined as the number of age-related conditions including arthritis, cataract, cancer, CHD, CHF, emphysema, glaucoma, lipid condition, osteoporosis, type 2 diabetes (see Supplementary Note 1). In our validation analysis, we used the age-adjusted variable, DNAmTLadjAge, which is not correlated with chronological age. All regression models included the following covariates: chronological age, sex, and batch effect as needed. To avoid bias due to familial correlations from pedigrees in the FHS cohort or intra-subject correlations resulting from repeated measurements, we used the following techniques. For censored time variables, we used robust standard errors, the Huber sandwich estimator, implemented in the R function "coxph". We used linear mixed models with a random intercept term, implemented in the R function lme. We used generalized estimation equation models (GEE), implemented in the gee function, for our logistic regression models. Additional covariates related to demographic characteristics, psychosocial behaviors and clinical covariates were adjusted in multivariate Cox models analysis. Those additional covariates include BMI, educational attainment (category), alcohol consumption (gram/day), self-reported smoking pack-years, three medical covariates: status of (any) cancer, hypertension and type 2 diabetes at baseline. BMI was categorized into 3 groups: 18.5 -24.9 (normal), b) 25 to 29.9 (overweight), and c) >=30 (obese). The categories associated with educational attainment were a) less than high school, b) high school degree, c) some college, and d) college degree and above. Smoking pack-years and educational variables were not available in the JHS cohort. Smoking status (never, former and current) was used in the analysis of the JHS cohort. Our stratified analysis was conducted in strata defined by age (<65 versus ≥ 65 years), BMI, education, prior history of hypertension, type 2 diabetes or cancer. All models used in the stratified analysis adjusted for age, sex, and (possibly) batch effect.

Meta-analysis

We mainly used fixed effects meta-analysis models weighted by inverse variance to combine the results across validation study sets into a single estimate. Toward this end, we used the metafor R function. Alternatively, we used Stouffer’s meta-analysis method (weighted by the square root of sample size) to combine results for variables whose scale/definition differed across study sets, e.g. multimorbidity (number of age-related conditions), disease-free status and physical function scores.

Cox models that include blood cell counts

We also fit multivariate Cox regression models that adjusted for imputed blood cell counts in addition to chronological age, batch, and pedigree structure, for predicting time-to-death and time-to-CHD. The blood cell counts were imputed based on DNA methylation levels as described elsewhere [60,61]. To avoid multi-collinearities between blood cell counts, we only included the following seven blood cell counts into the multivariate model: naïve CD8+T, exhausted cytotoxic CD8+ T cells, plasma blasts, CD4+T, natural killer cells, monocytes and granulocytes.

Heritability analysis

In general, DNAm based biomarkers are highly heritable [19,62,63]. To evaluate whether DNAmTLadjAge is heritable as well, we estimated the narrow sense heritability ${\mathit{h}}^2$ using the polygenic models defined in SOLAR and its R interface solarius [64]. Heritability is defined as the total proportion of phenotypic variance attributable to genetic variation in the polygenic model. The quantitative trait DNAmTL was adjusted for both age and sex. The robust polygenic model (with the option of a t-distribution) was used to estimate heritability. The heritability estimate correspondents to the variance component associated with the kinship coefficient. We used all individuals from the FHS cohort for whom DNA methylation data were available (irrespective of the availability of the observed LTL measure).

SNP associations of DNAmTL

We performed 14 SNP associations across 10 distinct susceptibility loci associated with LTL reported from three large-scale studies: (I) meta-analysis association of LTL in chromosome 5 TERT only (N=53,724) [36], (II) a genome-wide meta-analysis of LTL (N=37,684) [37], or (III) a genome-wide meta-analysis of LTL (N=26,089) [38]. SNP associations were performed for DNAmTLadjAge and LTLadjAge respectively, using the individuals of the FHS cohort available for both measures (N=811). The association analysis was based on linear mixed models whose random covariance matrix was determined by kinship coefficients of the pedigree structure, adjusted for sex and three principal components as fixed effects. We conducted the association tests using the function "relmatLmer" from the R library "lme4qt".

LTL measures versus blood cell composition

The imputed blood cell abundance measures were related to TRF based and DNAm based LTL measures, using the training datasets from the WHI BA23 and the JHS cohort and the test dataset from the FHS cohort, involving n=3,134 individuals. The following imputed blood cell counts were analyzed: B cell, naïve CD4+ T, CD4+ T, naïve CD8+ T, CD8+ T, exhausted cytotoxic CD8+ T cells (defined as CD8 positive CD28 negative CD45R negative), plasma blasts, natural killer cells, monocytes, and granulocytes. The blood cell composition imputation of the naive T cells, exhausted T cells, and plasma blasts was based on the Horvath method [65]. The remaining cell types were imputed using the Houseman method [61]. More details are described in Supplementary Methods. To avoid confounded by age, we used the age-adjusted DNAmTL (DNAmTLadjAge) variable for analysis. The correlation results were combined across studies via the same fixed effects meta-analysis model.

GREAT analysis

We applied the GREAT analysis software tool to three sets of CpGs: (1) all the 140 CpGs underlying DNAmTL model, (2) the 72 CpGs with negative coefficients in the model, and (3) the other 68 CpGs with positive coefficients in the model. CpGs in non-coding regions typically lack annotation with respect to biological functions. GREAT assigns biological meaning to a set of non-coding genomic regions (implicated by the CpGs) by analyzing the annotations of the nearby genes. Toward this end, the GREAT software performs both a binomial test (over genomic regions) and a hypergeometric test over genes when using a whole genome background. We performed the enrichment based on default settings (Proximal: 5.0 kb upstream, 1.0 kb downstream, plus Distal: up to 1,000 kb) for gene sets associated with GO terms, MSigDB, PANTHER, KEGG and InterPro pathway. To avoid large numbers of multiple comparisons, we restricted the analysis to the gene sets with between 5 and 3,000 genes. We report nominal P values and two adjustments for multiple comparisons: Bonferroni correction and the Benjamini-Hochberg false discovery rate.

Diet and lifestyle factors

We performed a robust correlation analysis (biweight midcorrelation, bicor) or generalized linear regression analysis between DNAmTLadjAge and 43 variables including 15 self-reported dietary variables, 9 dietary biomarkers, 14 variables related to metabolic related traits and central adiposity, and 5 life style factors, using the FHS, WHI, LBC 1921 and/or LBC 1936 cohort (N up to 6977), as listed in the first two columns in Figure 5. In the FHS cohort, we conducted the robust correlation analysis in males and females, respectively. Next we combined the results via fixed effect models weighted by inverse variance. In the FHS cohort, we used linear mixed effects models to account for pedigree structure. In the WHI cohort, we conducted the robust correlation analysis in each ethnic group separately. Next we combined the results via fixed effects meta-analysis models. In the LBC cohorts, we performed generalized linear regression analysis adjusted for sex for the 1921 and 1936 follow-up, respectively. Most of the final results were combined based on the Stouffer’s method weighted by the square root of sample size across the study cohorts. A few variables that were only available in both FHS and WHI cohort and presented in the same scaling of effect sizes were combined based on the fixed effect models weighted by inverse variance. In the FHS cohort, the diet and the clinical variables were based on the datasets archived in dbGAP pht002350.v4.p10 and pht000742.v5.p10, respectively. Both datasets were collected at exam 8, aligned with the blood drawn for DNA methylation profiles. In the WHI cohort, blood biomarkers were measured from fasting plasma collected at baseline. Food groups and nutrients are inclusive, including all types and all preparation methods, e.g. folic acid includes synthetic and natural, and dairy includes cheese and all types of milk. The individual variables are explained in [66].

Interrogating smoking impact on telomere length in blood

We used the 4,039 subjects across seven studies from our validation cohorts to compare the smoking association with (1) DNAm TL and with (2) measured LTL. Multivariate linear regression analysis of telomere length on smoking was performed on each study, adjusted for chronological age and adjusted sex, ethnicity and pedigree effects as needed. Smoking information was based on pack years when available (N=2216) otherwise based on a binary variable (never smokers versus ever smokers). The Stouffer’s method weighted by square root of sample size was performed to combine the results across all studies. Fixed effect models weighted by inverse variance were applied to the studies with pack years information and the other studies only with the binary smoking variable, respectively.

Enrichment analysis of CpGs in DNAmTL near telomeres

We performed hypergeometric analysis for evaluating the overlap between the 140 CpGs comprising our DNAm TL model and the CpGs nearby telomere regions from 3 mega base (Mb) at each chromosome tail. In the hypergeometric test, the margin for the total number of CpGs across whole genome was based on 453093 CpGs present in both 450k and Epic arrays. We also performed sensitivity analysis based on a variety of regions thresholded at 2, 4, 6 and 8 Mb respectively.

Enrichment analysis of CpGs in DNAmTL annotated with DNA methylation quantitative trait locus (mQTL), conducted in McRae et al. [35]. The hypergeometric test was used to evaluate the overlap of the 140 CpGs with 52915 cis-mQTL. The mQTL CpGs were identified using a stringent criterion described in the original reference: a cis-mQTL was identified either in the Brisbane Systems Genetics Study or in the Lothian Birth Cohorts at P < 1.0E-11 and was replicated in the other cohort at P < 1.0E-6.

The resulting 51 SNPs underlying cis-mQTL loci were cross referenced to GWAS Catalog tool (version v1.02 database, see URL). We not only report SNPs that met genome-wide (P < 5.0E-8) significance levels but also suggestive significance (P< 5.0E-6). Similarly, we characterized the 51 SNP using the HaploReg (version 4.1) tool [67] (see URL) which allows one to visualize SNPs in linkage disequilibrium (LD r² ≥ 0.8 based on European ancestry).

In vitro cultured cell procedure

URLs

GWASCatalog, https://www.ebi.ac.uk/gwas/search?query=28346442

HaploReg, http://www.broadinstitute.org/mammals/haploreg/haploreg.php