Comparing 14 biological aging measures and their links to 174 diseases

Cox proportional hazards regression was run for each clock-disease pairing, first adjusting for just age and sex (Supplementary Data 2) followed by a fully-adjusted model that also controlled for body mass index, smoking, alcohol consumption, education, and socioeconomic deprivation (Supplementary Data 3). Interaction models were considered for both the basic and fully-adjusted models to test for differences by smoking status (ever/never smoked) or sex (Supplementary Data 2). Age, estimated cell proportions and relatedness were regressed out of each clock prior to the Cox regression analyses. The fully-adjusted covariate setup might directly capture information that is included in the clocks, e.g., a DNAm proxy for smoking is included in the construction of GrimAge, a predictor of mortality⁷. However, these covariates are easy to assess and commonly collected in clinical settings. Thus, in terms of disease prediction, it is valuable to note what added information epigenetic clocks might contribute. A total of 176 Bonferroni significant associations were found for 13 of the clocks across 57 diseases in the fully-adjusted models (P < 2.9 × 10⁻⁴). In addition to the Cox regressions, logistic regression 10-year classification was conducted. Differences in AUC were tested for between nested models (first adjusting for the full set of aforementioned covariates and then adding in an epigenetic clock).

There were between 12 (Lin clock) and 72 (both GrimAge v2 and DunedinPACE) statistically significant associations (P < 0.05/174) in the age- and sex-adjusted Cox regression models (Table 1 and Supplementary Data 2).

An R Shiny app displays the results from both the fully-adjusted Cox and logistic regression analyses (https://shiny.igc.ed.ac.uk/Epigenetic_Clock_and_Disease_Association/↗). Users can visualise findings for specific or multiple clocks and disease outcomes or diseases within a grouping, e.g, psychiatric or cardiovascular traits.

To help contextualise the magnitude and significance of the clock-disease associations, we first ran survival models for 10-year all-cause mortality (n_cases = 842; Supplementary Data 6)–one of the best-established health outcomes to be related to epigenetic age acceleration¹². The largest and most significant association was with GrimAge v2 (Hazard Ratio (HR) per SD of age acceleration = 1.54, 95% CI [1.46, 1.62], P = 7.1 × 10⁻⁶²). All clocks had significant associations at P < 0.05 apart from Horvath's skin and blood clock, Ying's AdaptAge and Lin's age predictor. Using logistic regression, an AUC of 0.851 was obtained for the model with covariates only. This increased by up to 0.014 (1.4%) upon the addition of GrimAge v2 (AUC = 0.865) as a covariate.

There were 27 unique disease outcomes where the clock-disease hazard ratio was Bonferroni-significant and exceeded the magnitude of the corresponding clock-mortality association (Supplementary Data 7). A diverse set of diseases were highlighted, however there was a strong focus on respiratory/smoking-related and liver-related outcomes, including primary lung cancer (HR_GrimAgev1 per SD = 1.56 [1.42, 1.72], P = 5.3 × 10⁻¹⁹) and cirrhosis (HR_GrimAgev2 = 1.86 [1.57, 2.21], P = 8.9 × 10⁻¹³). Also of note, were associations with diabetes (HR_DunedinPACE = 1.44 [1.33, 1.57], P = 9.6 × 10⁻¹⁹), Crohn's disease (HR_PhenoAge = 1.39 [1.19, 1.64], P = 4.7 × 10⁻⁵) and delirium (HR_Zhang10 = 1.44 [1.23, 1.68], P = 6.7 × 10⁻⁶).

In the fully adjusted sex- and smoking-stratified analyses and interaction models, there were only two significant findings from analyses where the number of cases per stratum exceeded 30 (Supplementary Data 2): the associations between GrimAge v1 and v2 and the risk of post-viral fatigue syndrome (HR per SD of GrimAge v1 in ever smokers = 1.51 compared to 0.81 in never smokers, P_interaction = 6.2 × 10⁻⁵).

Again, several clocks were linked to respiratory/smoking-related disease outcomes, including COPD, lung cancer and respiratory failure.

All components of Generation Scotland received ethical approval from the NHS Tayside Committee on Medical Research Ethics (REC Reference Number: 05/S1401/89). All participants provided broad and enduring written informed consent for biomedical research. Generation Scotland has also been granted Research Tissue Bank status by the East of Scotland Research Ethics Service (REC Reference Number: 15/0040/ES), providing generic ethical approval for a wide range of uses within medical research. This study was performed in accordance with the Helsinki declaration.

Self-reported biological sex was considered in this study–this was cross-referenced via genetic data. The distribution of sex by disease outcome is provided in Supplementary Data . Analyses were conducted on the full cohort to maximise statistical power; an interaction with sex and sex-stratified models were also considered. Participants did not receive compensation for participation. 1

Generation Scotland (GS) is a family-based research study comprised of 24,084 volunteers¹⁴. Volunteers joined the cohort between 2006 and 2011 with the majority of recruitment taking place via invitation from the individuals' general practitioners (GPs). Initially, individuals aged 35–65 years were recruited from across five geographic areas (Glasgow, Aberdeen, Dundee, Ayr and Arran). They were then encouraged to recruit family members living in Scotland to join the study. This resulted in a baseline cohort aged between 17 and 99 years. All individuals were invited to complete questionnaires and to attend a baseline clinic visit for further questionnaires and testing, resulting in a detailed collection of lifestyle, clinical, health, cognitive and sociodemographic data. The majority also agreed to provide biosamples, including blood, from which DNA has been extracted and profiled for genetic and epigenetic data. Most individuals also provided consent for researchers to access their medical records via data linkage to primary and secondary care codes.

Methylation data have been quantified for 18,869 GS volunteers at a total of 851,610 genomic (CpG) sites after quality control via the Illumina EPIC850k array¹¹. Following 10 study withdrawals, the analysis sample size was 18,859. Epigenetic age and age acceleration estimates were derived via the Biolearn platform. Biolearn is an open-source python library that enables easy and versatile analyses of biomarkers of aging data. It provides tools to easily load data from publicly available sources and contains reference implementations for common aging clocks such as the Horvath clock, DunedinPACE, and many others that can easily be run in only a few lines of code¹⁵.

Age acceleration for each clock/predictor was calculated as the residual from a linear mixed model where the clock/predictor was regressed on age and estimated cell proportions derived from the Biolearn platform using the Salas et al. approach¹⁶. To avoid collinearity (as the cell proportions sum to one), basophil proportion was not included as a covariate. A pedigree kinship matrix was included as a random effect to control for the family structure present in the cohort. The regression models were run using the lmekin function from the coxme R package (2.2-22)¹⁷. As there are multiple approaches to estimating white cell proportions from DNA methylation data, a sensitivity analysis was conducted where each clock was adjusted for proportions derived from the EpiDISH R package¹⁸ (version 2.22.0; reference panel set to "cent12CT.m", method = "RPC"; basophil proportion was excluded to avoid collinearity as the sum of the proportions is one). The correlation between these residuals and those used in the analyses ranged between 0.90 (Zhang clock) to 0.99 (YingDamAge), suggesting minor differences between the two approaches.

Using a list of consensus definitions based on primary and secondary care codes¹⁹ 308 disease outcomes were derived from the electronic health records data in Generation Scotland from January 1980 to April 2022. Despite having volunteer consent to access primary care data, these were only accessible for ~40% of the cohort due to consent constraints with individual GP surgeries (the data holders). These data were subsequently filtered to only consider the first diagnosis for each disease, which could be made in either a primary or secondary care setting. Those with a diagnosis prior to joining the study were excluded from downstream analyses for that specific disease outcome, in order to ensure investigation of incident as opposed to prevalent disease outcomes. Not all diseases were observed in the cohort. The latest date of linkage with primary and secondary healthy records was April 2022. Filtering to diseases with at least 30 cases over the first 10-years of follow-up resulted in 174 outcomes for the main analyses. A list of the 174 diseases observed and the number of incident events (events occurring after the blood draw) along with mean age at diagnosis are presented in Supplementary Data 1.

Mortality records were obtained via Community Health Index linkage to the National Records of Scotland. These data formed part of our censoring criteria where individuals either survived and remained unaffected by a disease or died during the follow up period to April 2022.

Demographic, lifestyle and health variables were considered as covariates in the regression models with each clock and disease outcome. These included: age, sex, body mass index (BMI, calculated as weight in kg divided by squared height in metres), smoking pack years (calculated by self-reported packs per day multiplied by number of years as a smoker), alcohol intake, which was assessed as units consumed over the last week in a self-report questionnaire, years of education (an ordinal variable with 11 self-report categories: (1) 0 years, (2) 1–4 years, (3) 5–9 years, (4) 10–11 years, (5) 12–13 years, (6) 14–15 years, (7) 16–17 years, (8) 18–19 years, (9) 20–21 years, (10) 22–23 years, (11) 24 or more years), and an area-based rank of socioeconomic deprivation (Scottish Index of Multiple Deprivation, SIMD). SIMD measures relative deprivation across 6976 data zones, calculated across seven domains: income, employment, education, health, access to services, crime and housing.

K-nearest neighbours' imputation–using k = 10 and a maximum missingness per person of 4/7 of variables– was used to generate complete data for covariates via the impute package in R (version 1.80.0)²⁰. There was complete data on age and sex with a range of missingness for the other variables (SIMD: n = 1152, BMI: n = 120, Education: n = 1012, Pack years: n = 395, Alcohol: n = 1734). Natural log transformations were applied to alcohol and pack years (adding a constant of 1 to all values to account for non-drinkers/smokers) and BMI. All covariates were then standardised to mean zero and unit variance prior to imputation.

For the primary analysis, Cox proportional hazards regression models were run using the survival package in R (version 3.8-3)²¹, adjusting for the aforementioned covariates, and were used to relate age acceleration for each epigenetic clock to 10-year onset for the list of 174 incident disease outcomes. The analyses were subset to the most frequently affected sex where the proportion of male or female cases was >90% for any given disease. A total of 174 models were run for each of the 14 clocks, giving 2436 models in total. P value filtering was conducted per individual clock-regression model, resulting in a Bonferroni threshold of P < 0.05/174 = 2.9 × 10⁻⁴. A further set of analyses were run where time to all-cause mortality was the outcome.

Additional analyses included basic-adjusted Cox models (controlling only for age and sex) as well as sex- and smoking (ever/never)-stratification and interaction testing for both the primary analysis (fully-adjusted Cox models) and the basic-adjusted Cox models. When stratifying by smoking or testing the interaction formally, the continuous pack years variable was excluded as a covariate.

The proportional hazards assumption was tested by examining the Schoenfeld residuals. A test of the residuals for both the overall model and the predictor variable of interest (age acceleration) was performed using the cox.zph function in the survival package. This was done for both the basic- and fully-adjusted Cox models as well as for the stratified analyses (male, female, ever smoker, never smoker).

A formal comparison of the log Hazard Ratios for the clocks across all 174 disease outcomes in the fully adjusted model (primary analyses) was then conducted using a linear mixed effects model (using the lme4 and lmerTest R packages, versions 1.1-35.4 and 3.1-3, respectively)^22,23. The logHR was the outcome while the clock was the predictor (GrimAge v1 was the reference level) and disease was included as a random effect on the intercept. The emm function from the emmeans R package (version 1.10.2)²⁴ was then used to obtain the marginal means for each clock, prior to formal pairwise comparisons with a Tukey adjustment to account for multiple comparisons.

In addition to the Cox models, logistic regression was used to generate area under the curve (AUC) estimates for 10-year disease classification. Null models contained the aforementioned covariates, while a full model also included an age acceleration measure. A formal comparison of the nested models was conducted using the roc.test function from the pROC R library (version 1.19.0.1)²⁵.

Two-sided statistical tests were considered for all analyses.

Further information on research design is available in the linked to this article. Nature Portfolio Reporting Summary