
==== Front
eBioMedicine
EBioMedicine
eBioMedicine
2352-3964
Elsevier

S2352-3964(24)00315-3
10.1016/j.ebiom.2024.105279
105279
Articles
NMR metabolomics-guided DNA methylation mortality predictors
Bizzarri Daniele abc
Reinders Marcel J.T. bc
Kuiper Lieke de
Beekman Marian a
Deelen Joris afg
van Meurs Joyce B.J. dh
van Dongen Jenny ijk
Pool René ik
Boomsma Dorret I. ijk
Ghanbari Mohsen l
Franke Lude m
BBMRI ConsortiumGeleijnse J.M.
Boersma E.
van Spil W.E.
van Greevenbroek M.M.J.
Stehouwer C.D.A.
van der Kallen C.J.H.
Arts I.C.W.
Rutters F.
Beulens J.W.J.
Muilwijk M.
Elders P.J.M.
't Hart L.M.
Ghanbari M.
Ikram M.A.
Netea M.G.
Kloppenburg M.
Ramos Y.F.M.
Bomer N.
Meulenbelt I.
Stronks K.
Snijder M.B.
Zwinderman A.H.
Heijmans B.T.
Lumey L.H.
Wijmenga C.
Fu J.
Zhernakova A.
Deelen J.
Mooijaart S.P.
Beekman M.
Slagboom P.E.
Onderwater G.L.J.
van den Maagdenberg A.M.J.M.
Terwindt G.M.
Thesing C.
Bot M.
Penninx B.W.J.H.
Trompet S.
Jukema J.W.
Sattar N.
van der Horst I.C.C.
van der Harst P.
So-Osman C.
van Hilten J.A.
Nelissen R.G.H.H.
Höfer I.E.
Asselbergs F.W.
Scheltens P.
Teunissen C.E.
van der Flier W.M.
van Dongen J.
Pool R.
Willemsen A.H.M.
Boomsma D.I.
Suchiman H.E.D.
Barkey Wolf J.J.H.
Beekman M.
Cats D.
Mei H.
Slofstra M.
Swertz M.
Reinders M.J.T.
van den Akker E.B.
Boomsma D.I.
Ikram M.A.
Slagboom P.E.
no
Slagboom Pieternella E. af
van den Akker Erik B. e.b.van_den_akker@lumc.nl
abc∗
a Molecular Epidemiology, Department of Biomedical Data Sciences, Leiden University Medical Center, Leiden, the Netherlands
b Leiden Computational Biology Center, Department of Biomedical Data Sciences, Leiden University Medical Center, Leiden, the Netherlands
c Delft Bioinformatics Lab, TU Delft, Delft, the Netherlands
d Department of Internal Medicine, Erasmus MC, Rotterdam, the Netherlands
e Center for Nutrition, Prevention and Health Services, National Institute for Public Health and Environment (RIVM), Bilthoven, the Netherlands
f Max Planck Institute for the Biology of Ageing, Cologne, Germany
g Cologne Excellence Cluster on Cellular Stress Responses in Aging Associated Diseases, University of Cologne, Cologne, Germany
h Department of Orthopaedics & Sports, Erasmus Medical Center, Rotterdam, the Netherlands
i Department of Biological Psychology, Vrije Universiteit Amsterdam, Amsterdam, the Netherlands
j Amsterdam Reproduction and Development (AR&D) Research Institute, Amsterdam, the Netherlands
k Amsterdam Public Health Research Institute, Amsterdam, the Netherlands
l Department of Epidemiology, Erasmus MC, Rotterdam, the Netherlands
m Department of Genetics, University Medical Center Groningen, Groningen, the Netherlands
∗ Corresponding author. Leiden Computational Biology Center, Leiden University Medical Center, Einthovenweg 20, 2333, ZC, Leiden, the Netherlands. e.b.van_den_akker@lumc.nl
n BIOS Consortium, see Consortium Banner Supplementary Material.

o BBMRI-NL Consortium, see Consortium Banner Supplementary Material.

17 8 2024
9 2024
17 8 2024
107 10527914 3 2024
25 7 2024
29 7 2024
© 2024 The Authors
2024
https://creativecommons.org/licenses/by/4.0/ This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).
Summary

Background

1H-NMR metabolomics and DNA methylation in blood are widely known biomarkers predicting age-related physiological decline and mortality yet exert mutually independent mortality and frailty signals.

Methods

Leveraging multi-omics data in four Dutch population studies (N = 5238, ∼40% of which male) we investigated whether the mortality signal captured by 1H-NMR metabolomics could guide the construction of DNA methylation-based mortality predictors.

Findings

We trained DNA methylation-based surrogates for 64 metabolomic analytes and found that analytes marking inflammation, fluid balance, or HDL/VLDL metabolism could be accurately reconstructed using DNA-methylation assays. Interestingly, a previously reported multi-analyte score indicating mortality risk (MetaboHealth) could also be accurately reconstructed. Sixteen of our derived surrogates, including the MetaboHealth surrogate, showed significant associations with mortality, independent of relevant covariates.

Interpretation

The addition of our metabolic analyte-derived surrogates to the well-established epigenetic clock GrimAge demonstrates that our surrogates potentially represent valuable mortality signal.

Funding

BBMRI-NL, X-omics, VOILA, 10.13039/501100022216 Medical Delta , 10.13039/501100003246 NWO , ERC.

Keywords

DNA methylation predictors
NMR metabolomics
Ageing biomarkers
Epigenetic clock
Metabolic risk score
Epidemiology
==== Body
pmc Research in context

Evidence before this study

The development of a comprehensive “biomarker of ageing” capable of identifying frailty and functional decline in the global ageing population would represent a significant leap forward in geroscience and healthcare prevention. Omics clocks, sophisticated prediction models condensing extensive molecular data into manageable multi-biomarker scores, hold immense promise in elucidating the biological trajectory of ageing. Focusing our attention on second generation of clocks, which draw predictive value from longitudinal endpoints, 1H-NMR metabolomics and DNA methylation are emerging as pivotal omics-layers, with prominent scores such as MetaboHealth, PhenoAge, and GrimAge. DNA methylation emerges as a remarkably precise estimator of chronological age, with additional potential to encode for risk factors like as smoking, consequently also predicting mortality and disease susceptibility. Nonetheless, also 1H-NMR metabolomics exhibits robust predictive capabilities in tracking various risk factors and endpoints, particularly when concerning cardiovascular mortality. Recent investigations underscore the distinct and complementary nature of mortality-oriented scores derived from both -omics layers. They appear to independently contribute to the prediction of frailty and mortality among older adults, highlighting the significance of integrating two omics layers.

Added value of this study

In this study, we investigate the potential of DNA methylation profiles to serve as epigenetic predictors for 64 metabolomics features and one metabolomics mortality score (MetaboHealth). Leveraging a dataset comprising of four large Dutch epidemiological cohorts, and designing rigorous calibration and feature selection techniques, we achieved a successful quantification of metabolomics markers related to inflammation, fluid balance, HDL-, and VLDL-related markers, and MetaboHealth. These markers exhibited distinct signals independent of previously trained methylation-base surrogates, both quantitatively and as mortality indicators, as validated in an independent test set. Consequently, a mortality model incorporating nine of our DNAm-metabolomics features demonstrated superior predictive capability compared to current state-of-the-art models.

Implications of all the available evidence

While previous research has successfully harnessed models to leverage clinical health variables, cell counts, and proteomics markers from methylation profiles, our study unveil the added value of integrating them with faithful metabolomics surrogates. These additional DNAm-based features not only expand the scope of epigenetic clocks to encompass cardiometabolic signal, but also demonstrate complementary signal to GrimAge and its surrogates in predicting mortality. Moreover, through a meticulous evaluation of the epigenetic content within our DNAm-based metabolomics models, we observe an optimised selection of highly informative CpG sites. These sites resonate with metabolic and developmental processes, cell differentiation dynamics, and the intricate landscape of cardiometabolic diseases, thus enriching our understanding of ageing-related pathways and facilitating more precise predictive modelling.

Introduction

A common goal in geroscience is to identify mechanisms that drive ageing and design interventions that might slow down or even reverse the rate of ageing.1 For this purpose, it is essential to have indicators not only quantifying ageing, but simultaneously marking the trajectory of overall health decline.2 While calendar age is a core risk factor for almost any common disease, it has many limitations for capturing the variability in health-span. Crucially, calendar age does not capture the effects of an individual's lifestyle, nor incorporates readouts of functional decline. Instead, faithful markers of biological age would allow to quantify the vulnerability to acute and chronic diseases irrespective of an individual's calendar age, and to develop and monitor effective healthy lifestyle advice and anti-ageing interventions. The earliest approaches to construct such markers of biological age relied on clinical measures of physiological capacity.3 Later molecular and -omics approaches gained popularity, initially including markers such as leukocyte telomere length,4 followed by multi-marker algorithms based on high-throughput platforms, such as DNA methylation,5 transcriptomics,6 metabolomics,7 and proteomics.8 Importantly, these algorithms were trained to estimate cross-sectional chronological age. Of these omics approaches, particularly DNA methylation-based algorithms exhibited remarkably high accuracies in predicting calendar age,5 and were named ‘DNA methylation clocks’. Nonetheless, while interesting by itself, this observation highlighted a fundamental limitation in this design of markers of biological age. Since nearly-perfect age predictors would arrive to similar observations as chronological age, they would lose their characteristics as age-independent health status indicators.9

Concomitantly, a second generation of -omics markers was introduced, which instead were trained to predict the mortality risk. Prominent examples of these mortality-trained multivariate markers include the DNA methylation-based PhenoAge10 and GrimAge11 and the 1H-NMR metabolomics-based MetaboHealth.12 These predictors were trained quite differently. The wide availability of the Nightingale Health 1H-NMR metabolomics in large prospective population studies, in combination with its relatively narrow though informative content (∼250 analytes), allowed for a more classic and direct approach. Deelen et al. trained MetaboHealth as a linear combination of 14 metabolic features, showing a strong predictive value, not only for mortality risk, but also for other outcomes, including pneumonia,13 and frailty.14 Conversely, the DNA methylation platform by Illumina contains hundreds of thousands of features, and thus requires additional guidance to robustly capture the mortality signal. Hence, the PhenoAge and GrimAge were trained using the so-called two-stage approaches, in which more widely-available markers associated with mortality were leveraged to help extract the mortality signal.10,11 PhenoAge achieved this by first training an all-cause mortality predictor based on clinical measures (e.g., glucose, C-reactive-protein), which was then re-estimated using DNA methylation. Similarly, DNAm-GrimAge is composed by a combination of DNA methylation-based surrogates for molecular or phenotypic markers known to associate with mortality. Interestingly, both two-step training strategies yielded DNA methylation-based scores that can associate not only with mortality, but also with a wide diversity of disease outcomes. These developments indicate that mortality-trained predictors for biological age can be trained using different omics platforms, and moreover, that DNA-methylation might serve as a platform to integrate these signals captured by different data sources. This latter concept was recently further substantiated by the work of Gadd et al., who systematically developed DNA-methylation-based predictors for 109 plasma proteins, called EpiScores. Their findings demonstrated significant associations with various incident morbidities over a span of 14-years.15 Moreover, these surrogates were afterwards employed, in an attempt to refine the GrimAge score in a clock known as bAge.16

In a recent study we demonstrated that mortality-based predictors such as MetaboHealth and GrimAge are instrumental in predicting frailty in studies of middle-aged and elderly individuals.14 Importantly, we also showed that these scores confer mutually independent information for predicting five frailty indicators and mortality (specific details in Kuiper et al.14). Viewing these developments in the field, we thus pose the question to what extent the mortality signal captured by 1H-NMR metabolomics could be transferred and integrated with the mortality signals captured by the DNA-methylation platform. For this purpose, we will evaluate both strategies for training two-stage DNA-methylation based mortality predictors. On one hand, we will train a DNA methylation-based predictor re-estimating directly MetaboHealth, akin the strategy of PhenoAge. On the other hand, we will train DNA-methylation surrogates for single metabolomics features, and combine these in an overall score, akin GrimAge. Moreover, we will evaluate to what extent DNA-methylation surrogates features from different origins capture mutually independent signals, also with respect to predicting mortality risk.

Methods

Dataset description

Cohorts

This study was performed using DNA methylation data (DNAm, Illumina 450 k array) and 1H-NMR metabolomics (Nightingale Health, platform version 2020) from 4 Dutch cohorts: LifeLines-Deep (LIFELINES), Leiden Longevity Study (LLS-PAROFFS), Netherlands Twin Register (NTR) and Rotterdam Study (RS-I), all part of the BIOS consortium.17,18 For the current study, the BIOS multi-omics compendium was further extended with 1145 samples from the NTR for which the entire process of array measurement to quality control and normalization was done together with the other BIOS-NTR samples,19 and 904 samples from the Rotterdam Study (RS-II).14 A thorough description of all cohorts and their ethics statement are provided in the Supplementary Materials. The datasets were realised by the Dutch part of the Biobanking and BioMolecular Resources and Research Infrastructure (BBMRI-NL). The final dataset contained data for 5238 individuals.

Metabolomics data

The metabolomics data were generated by the BMBRI-NL Metabolomics Consortium for the cohorts LIFELINES, LLS-PAROFFS, NTR and RS-I. The current study employs a total of 4334 EDTA plasma samples quantified with the high-throughput proton Nuclear Magnetic Resonance (1H-NMR) platform made available by Nightingale Health Ltd., Helsinki, Finland (platform version 2020). This technique can quantify over 250 metabolic features, including also ratios and derived features.20,21

DNA methylation data

DNA methylation data were generated by the subsection of BBMRI-NL named Biobank-based Integrative Omics Study (BIOS) Consortium for the cohorts LIFELINES, LLS-PAROFFS, NTR, RS-I, and RS-II (total of 5238 samples). The DNAm was assessed from whole blood samples with an Illumina iScan BeadChip according to the manufacturer's protocol: the Illumina HumanMethylation450 BeadChip (450 k array). For compatibility with the following versions of the Illumina array, we only considered CpG sites which are available in the Illumina HumanMethylation450 BeadChip and the MethylationEPIC BeadChip. We analysed the DNAm β values, which range from 0 to 1, to indicate the proportion of methylated sites at a specific CpG in a sample.

Mortality data

We evaluated the associations of the DNAm-based features with all-cause mortality in a subsample of the Rotterdam Study (RS-I and RS-II) comprising a total of 1544 samples, of which a subset of 640 (RS-I, 104 deaths) had also Illumina 450 k and Nightingale Health metabolomics. The information on the vital status of the participants in RS was last updated on the 20th of October 2022. The dataset comprehends 1544 samples, 285 of which are deceased. All the DNAm-based features were z-scaled within the RS.

Pre-processing

Quality control of the metabolomics dataset

To ensure the quality of our data, we applied standardised quality control processes, which have been described in previous publications (summarised in Figure S1).7,22 First, we limited our analyses to a subset of 65 features (out of 250), previously selected to be a mutually independent subset.7,12,22 This selection includes fatty acids, routine lipid concentrations, lipoprotein subclasses and low molecular weight metabolites. A complete list of the variables can be found in the Supplementary Materials Table S3. In addition, pyruvate was excluded due to its high missingness in NTR (80%). Despite a small percentage of values under detection limit for acetoacetate (8% in NTR), and an even smaller percentage of outliers in glucose and xl_hdl_c (less than 0.15%), we decided to retain all other variables (Figure S1C–E). Samples with more than 1 outlier (2 from Lifelines and 1 from RS-I) were further removed. We then used nipals (from the package pcaMethods) to impute the 584 missing values, which accounted for 0.211% of the remaining values. The final dataset included 4334 samples and 64 metabolic measures.

Quality control of the DNA methylation dataset

The quality control and normalization of the DNA methylation (DNAm) was performed using a workflow developed by the BIOS Consortium for each cohort and thoroughly described in DNAmArray (https://molepi.github.io/DNAmArray_workflow/). In brief, sample-level QC was performed with the R package MethylAid.23 Probes were set to missing based on the number of available beads (≤ 2), intensity equal to zero, or the detection p value (p < 0.01). Probes with more than 5% missing were excluded from all samples. The remaining missingness was imputed using impute.knn from the R package impute.24 Functional normalization was then applied as implemented in minfi. Finally, we removed an ulterior set of ∼60,000 underperforming probes as suggested by Zhou et al.25

Calibration of 1H-NMR-metabolomics

To minimise any bias that may arise from batch effects among the four cohorts included in our study, we performed a cross-cohort calibration. We followed the assumption that similar phenotypic characteristics result in similar metabolomics profiles.26 We used sex, age, and BMI as matching characteristics, given their well-known association with the metabolomic features in the Nightingale Health Platform.7,22,26, 27, 28 We considered LIFELINES as our reference cohort, as it spanned a broad range of age, and BMI and had an equal representation amongst sexes. To further minimise the impact of sex on our results, we selected the subset of samples used for cross-cohort matching to have equal numbers of men and women. Following this strategy, we identified the following subsets of participants used for matching: 73 men and 73 women in LLS-PAROFFS; 140 men and 140 women in NTR; 37 men and 37 women in RS (Figure S2).

Based on these matching samples between cohorts, we calculated the shift in mean and standard deviation for each metabolic feature required to transform the distribution of values observed in a cohort to match the distribution in the reference cohort. We then applied this transformation to all samples of each cohort (see Supplementary Materials). The final dataset was log-transformed and standard normalised (zero mean and unit standard deviation) across all samples to obtain normally distributed concentration with comparable ranges across all metabolic features.

T-distributed stochastic neighbour embedding (tSNE, R package Rtsne) was used to visually inspect the effect of this calibration, comparing the sample similarities before and after calibration. Moreover, K-nearest neighbour batch effect test (kBET, R package kBET) was applied to the matching samples of each biobank before and after calibration to quantitatively evaluate the mixing of the samples.29 Finally, we used principal variance Component Analysis (PVCA, R package pvca), to determine if the calibration disrupted the sources of variability of the dataset.30

Application of previously trained multivariate models

MetaboHealth

The MetaboHealth model is a mortality predictor based on Nightingale Health metabolomics concentration.12 We applied this model both on the uncalibrated and calibrated version of the 1H-NMR metabolomics dataset using the R-package MiMIR.31

Epigenetic clocks

We projected the Horvath, Hannum, DNAm PhenoAge in our data using the R package methylclock and the DNAm GrimAge clocks using Python scripts provided by Lu et al.11,32 About 1000 out of 30,000 CpG sites required for calculating these biological age scores were missing or removed during QC and were imputed using the “datMiniAnnotation3_GOLD.csv” file,33 as outlined by the authors of the original papers.11 The accuracy of the epigenetic clocks projected in BIOS is indeed preserved (Figure S13).

EpiScores

We projected the EpiScores and bAge using the code available from the work of Bernabeu et al.16

Estimation and evaluation of the epigenetic-based metabolic features

We derived prediction models for the 64 metabolomics features and the MetaboHealth score using blood methylation data. NTR, and LIFELINES, respectively the largest cohort and the cohort used as a reference for calibration, were employed for model development using the internal loop of a nested 5-Fold-Cross-Validation (5-Fold CV; Supplementary Material SM2). We employed ElasticNET regression from the R package glmnet to train the models (Supplementary Material SM3).

Other studies show the benefit of pre-selecting the features before using ElasticNET regression.16,34 For this reason, we performed Epigenome Wide Association Studies (EWASes) to identify CpG sites showing linear association with each feature separately in NTR and LIFELINES (metabolic feature ∼ CpG site) and pre-selected the CpG probes with a consistent association sign (positive or negative in both cohorts) and nominal significance (p value <0.05) to enrich for potentially predictive CpGs, while minimizing the risk of prematurely excluding valuable signal. Subsequent training of the ElasticNET regression models determined the final set of CpG sites included in each model (Supplementary Material SM2). An overview of the number of CpGs selected during the EWAS pre-selection and ElasticNET regression phases can be found in Supplementary Material Table S4. CpG selection and optimization of the penalization parameter λ were performed in the internal loop of the nested 5-Fold-Cross-Validation. The mixing parameter alpha was fixed at 0.5, based on previous work.11,22 Model performances in NTR and LIFELINES are reported using the outer loop of the nested 5-Fold-Cross-Validation and should give an unbiased impression of the performances of models created with our training procedure. The final models were obtained by training the ElasticNET models with the optimised parameters on the whole of NTR and LIFELINES. Resulting models where then applied to the held-out LLS-PAROFFS and RS datasets. Finally, reported measures for model performance are the Pearson correlation (R) and the root mean square error (RMSE) of the predicted DNAm metabolic features with their measured concentrations.

CpG sites characterization

To gain more insight into the biological phenomena that characterise our DNAm metabolomics models we evaluated their fully data-driven selection of 22,145 CpG sites.

EWAS enrichment analysis

We utilised the MRC-IEU EWAS Catalog35 and the EWAS Atlas36 to assess the previously known phenotypic annotation (traits) of the CpG sites selected by our models. Both the EWAS Catalog and EWAS Atlas are online databases that compile results from Epigenome Wide Association Study results. By merging these two resources, both downloaded on March 13th, 2023, we aimed to gather a comprehensive list of previous EWASs. We gathered only the associations conducted with Illumina 450 k and accompanied with a PMID. We then unified redundant or possible synonymous traits between the two catalogues (EWAS Atlas and EWAS Catalog) and selected only the trait-to-CpG associations significant after Bonferroni correction (using the number of CpGs in the Illumina 450 k ∼480,000). This process yielded 742,635 CpG-trait associations.

Next, we employed Fisher's exact test to assess the enrichments of each of the CpG sites selected by our DNAm models. To account for multiple testing, we used the Benjamini-Hochberg correction.

Annotation of the genomic position of the CpG sites

We used the R package annotatr to annotate the genomic features. CpG sites from the 450 k array were annotated using CpG Island (CGI) centric categories. The annotations we utilised are as follows: CGI (annotated in the R package AnnotationHub), shores (2 Kb upstream or downstream the CGI), shelves (2 kb flanking the CpG shores), interCGI (the rest of the CpGs). As for the genic annotations we considered regions 1–5 Kb upstream of the transcription starting site (TSS), promoters (<1 Kb upstream of the TSS), 5′UTR, 3′UTR, exons, introns, boundaries between introns and exons, and intergenic regions. Additionally, we report the annotations of active enhancers determined by Anderson et al.37

For the enrichment analyses of the annotations described above we calculated odds ratios (OR) of the CpGs included in each model compared to the rest of the 450 k array. Statistical significance was evaluated using the Fisher's exact test. The significance of the associations was established with an FDR <0.05.

Gene ontology enrichment analyses

To gain further insights into the genetic context of the set of CpG sites selected, we investigated the genes in cis, considering a maximum distance of 100 kB of distance.

Next, we utilised the genes associated with each CpG selections from our models to perform a functional enrichment using Gene Ontology. We employed GOfuncR package to explore the Biological Processes and Molecular Functions. The significance of the associations was established with an FDR <0.05. This analysis resulted in 2365 significant associations between CpG sites and genes.

Associations with mortality in the Rotterdam Study (RS-I and RS-II)

Univariate mortality associations

We used Cox Proportional hazard to univariately associate our 65 DNAm metabolomics features, the pre-trained DNAm clocks (e.g., PhenoAge, GrimAge), and 109 protein EpiScores with mortality (see Supplementary Materials). All models were corrected for age at blood sampling and sex. Additionally, we evaluated the association with mortality of our DNAm metabolomics features when correcting for sex, age and GrimAge. All p values were corrected using Benjamini Hochberg and considered significant if the FDR <0.05. We used the R-package survival to calculate the Cox regressions.

Multivariate mortality models

We then combined the DNAm features with sex and age in 4 different stepwise Cox regression models (see Supplementary Materials). The first base model included our DNAm metabolomics features. The second and third model added to the first model respectively DNAm-GrimAge and the DNAm-based components of the GrimAge model. Finally, the fourth model is based on a combination of our DNAm metabolomics features, the DNAm-based components of the GrimAge and the DNAm-based protein EpiScores.

To select the interesting DNAm surrogate, we used a stepwise (backward/forward) procedure for each Cox regression model. For each of the above-described selections, we started from a model containing the full set of variables and we removed or added an unselected metabolic surrogate at each round based on the improvement on the model calculated from the C-index, taking also into account the significance of the p value of each variable included in the model.

To compare the performances of the Cox regression models we used the R package survcomp within the Rotterdam Study.38 We compared the C-indices of the newly developed models with baseline (GrimAge) using a Student t-test as described in Haibe Kans et al.39 Moreover, we plotted the ROC curves at 5 and 10 years of all the models.

Role of funders

BBMRI-NL and BIOS contributed to the generation of the metabolomics data and the data sharing and computational resource infrastructure. The funding sources had no role in the design of this study, and did not have any role during its execution, analyses, interpretation of the data, or decision to submit results.

Ethics

The complete ethical statements for each cohort are available in “Supplementary Materials: BIOS Consortium, Ethics Statements”.

Results

Cross-cohort calibration of 1H-NMR metabolomics data

To derive DNA methylation-based models predicting metabolic features we analysed data gathered by partners of the BIOS consortium,17,18 totalling 4334 individuals for whom both DNA methylation (Illumina 450 k) and 1H-NMR Metabolomics (Nightingale Health Plc) data have been assayed. The resulting dataset had contributions of four independently collected population studies: LIFELINES-DEEP (LIFELINES), Leiden Longevity Study Partners-Offspring (LLS-PAROFFS), Rotterdam Study (RS), and Netherlands Twin Register (NTR), each with their own inclusion criteria, as reflected by differences in subject characteristics that range from the younger and leaner population of NTR (mean age = 37.57 years and mean BMI = 24.32 cm/kg2) to the older and heavier population of RS (mean age = 67.15 years and mean BMI = 27.71 cm/kg2) (Fig. 1, Table S1). A reduced dimensionality projection using a t-distributed neighbour embedding (tSNE) suggested that the interindividual variance in metabolomics data could not only be attributed to interindividual phenotypic variability but was also capturing some systematic differences between studies (Figure S3a–c). Following Makinen et al., we implemented a calibration technique suitable for cross-cohort harmonization, which starts with the assumption that individuals with similar phenotypic characteristics should on average exhibit similar metabolomics profiles.26 For this purpose, we identified pairs of samples across cohorts with matching age, sex, and BMI, and used LIFELINES as a common reference to calibrate the other studies (Figure S2, more details in methods). A t-SNE projection of the calibrated data revealed a substantial reduction of the systematic differences between studies, as also quantified by k-BET (k-nearest neighbour Batch Effect Test) (Fig. 2a and b, Figure S3). Principal Variance Component Analyses (PVCAs) further confirmed this observation indicating that the variation attributable to study differences was attenuated, while maintaining the variation attributed to relevant biologically variability (Fig. 2c).Fig. 1 Study and methods overview. a) Study overview. (i) We employed 4334 samples, from 4 cohort of the BIOS Consortium, DNAm methylation and metabolomics to train and test our surrogates (orange border). (ii) Coupled with 1544 samples from the Rotterdam Study (black border) to evaluate their associations with mortality. (iii) We applied a calibration to harmonise all the metabolomics data, using Lifelines as a reference dataset. (iv) We then train ElasticNET models, on LIFELINES and NTR. Using the DNA methylation data, we predict two types of outcomes: 1) the pre-trained metabolomics mortality predictor (MetaboHealth), and 2) the 64 metabolic features. (v) The DNAm models are evaluated using 1) the hold-out valuation sets (LLS and RS) and 2) a 5-Fold Cross Validation on the training sets (NTR and LIFELINES). (vi) Finally, we use the DNAm models to generate surrogate metabolomics features in the RS dataset (1544 samples) and 1) evaluate their univariate associations to mortality (while correcting for age, sex), and 2) trained a complete Cox regression combining our DNAm metabolomics features and the pre-trained DNAm surrogates. b) Data availability, usage, and main phenotypic characteristics of the individuals in each cohort.

Fig. 2 Harmonization of the metabolomics data. a) Distribution of glucose in NTR and LIFELINES before (upper figure) and after (lower figure) calibration. b) tSNE of the metabolomics dataset after calibration and coloured by the four biobanks (LIFELINES, LLS, RS and NTR). c) Principal Variance Component Analysis (PVCA) before (red) and after (green) calibration, estimating the variance explained in the dataset by available clinical variables (e.g., sex, age, BMI, diabetes). d) Bar-plots showing the differences in men and women in the calibrated MetaboHealth in the four cohorts. e) Observed mean values of age, BMI, eGFR, hsCRP and pressure, and f) alcohol consumption, current smoking, and diabetes ordered following the calibrated MetaboHealth in different percentiles over the entire BIOS population. On top of each panel the Spearman correlation ρ.

The construction of the MetaboHealth score as published by Deelen et al.12 does not include a cross-cohort calibration, but instead standardises the individual metabolic features per study prior to computation of the score (Figure S4). While this does make the MetaboHealth score more comparable across cohorts, and satisfactory for most meta-analysis purposes, it does also negate any real biological differences that may exist between studies. Conversely, when computing the MetaboHealth score on the calibrated data, i.e., after removing unwanted study differences and supplying all data as one dataset, it produced scores with interpretable differences and consistent trends across cohorts (Fig. 2d–f). For instance, consistent with our expectation, the calibrated MetaboHealth scores now tend to be higher among the studies with the older individuals RS and LLS-PAROFFS (Figure S4a). In addition, it showed a more pronounced age-associated increase in men than in women, consistently over all cohorts (Fig. 2d). Lastly, higher calibrated MetaboHealth percentiles correlated with increasing age, BMI, high-sensitive CRP, and increasing prevalence of diabetes and alcohol usage (Fig. 2e and f).

DNA methylation-based predictors recapitulate metabolic markers previously associated with mortality

Our first objective was to determine if DNA methylation could simulate the MetaboHealth score, our metabolomics-based mortality predictor (Fig. 1). To enforce the selection of consistent signal in different cohorts, we implement an output-specific pre-selection of consistent DNA methylation sites in the studies reserved for model development, NTR and LIFELINES (methods, Supplementary Material SM2). This Epigenome Wide Association Study (EWAS) yielded 17,705 CpG sites showing a consistent univariate association with MetaboHealth, both in direction of association and nominal significance (p value < 0.05). Pre-selected sites were then used as input for the ElasticNET regression model predicting the MetaboHealth values (Supplementary Material SM2, Table S4)). The resulting model, indicated as “DNAm-MetaboHealth” comprised ∼1000 sites and showed good accuracy in the 5-Fold Cross Validation test sets (5-FCV) (median R∼0.52, RMSE∼0.43), which was slightly lower, but stable, in the replication sets (LLS-PAROFF: R∼0.34, RMSE∼0.38; RS: R∼0.33, RMSE∼0.5) (Fig. 3).Fig. 3 DNAm metabolites accuracies. Circular heatmap representing the accuracies of the DNAm-based models for 64 1H-NMR metabolic features by Nightingale Health and MetaboHealth. The outer ring shows the correlation between measured and DNAm-based metabolomics features, while correlation between DNAm surrogates with age and sex are shown in the middle and inner ring, respectively. Mean CV states for mean cross validation results in the cohorts LL and NTR together in the training (train) and test (test) sets, while results in the left-out set are indicated with RS and LLS. Moreover, the metabolomics features are annotated for their metabolomics group type (e.g., amino acids, fatty acids etc.) and if they were or were not included in MetaboHealth. Finally, we indicated with asterisks the tertiles of mean correlations with the measured metabolites over the test sets.

In parallel, we built distinct predictors for 64 metabolic features from Nightingale Health Plc, following the same training design as for DNAm MetaboHealth (Fig. 1a and Supplementary Material SM2). The resulting DNAm-based surrogates for the metabolomic features showed diverse mean accuracies over the different test sets (5-FCV, LLS-PAROFFS and RS), with 23 models being accurate (mean R across test sets >0.35), 20 mildly accurate (0.2 >mean R across test sets ≤0.35), and 21 low accuracy models (mean R across test sets <0.2) (Fig. 3 and Figure S5). In the latter group we find 5 out of 8 amino acids, several LDL-related variables, all the ketone bodies and all the glycolysis related markers. The middle group is enriched with IDL related markers, 6 out of 14 fatty acids, and 2 out of 3 glycolysis related metabolites. The accurate group of DNAm-metabolomic features included, 8 out of 10 HDL-related markers, 4 out of 8 VLDL-related molecules, glycoprotein acetyls, creatinine, 3 out of 8 amino acids, and several fluid balance markers (e.g., MUFA%, Omega6%). Higher accuracies are often accompanied by a higher correlation with age (e.g., DNAm Leucine and Isoleucine) or sex (e.g., DNAm Creatinine) (Fig. 3, inner circles), similar to what is observed for the GrimAge DNAm based components.11 Notably, only 7 of the most accurate surrogate markers were part of the 14 original metabolomic features composing the MetaboHealth score. Nevertheless, for 18 of the 23 most accurate surrogate markers, it was previously shown that the respective metabolic features significantly associated with mortality.12

DNAm metabolomics surrogates confer a unique and relevant signal

To foster the concept that DNA-methylation measurements might potentially serve as a platform to integrate biomarker signals captured from various data sources, we conducted two types of experiments. First, we ensured that the signals conveyed by our DNAm surrogates of metabolomic features, constitute mutually independent markers, and not a multitude of highly similar signals (Figure S6a). Then, comparing the correlations between our surrogates with the original metabolomic features (Figure S6a, upper triangle), we observed a structure remarkably congruent with the correlation structure observed between the original metabolites (Figure S6a, lower triangle), albeit overall at slightly lower magnitude. This indicates that, apart from the correlation structure between the original markers, no systematic high inter correlations are observed, which thus suggests that the information in DNA methylation measurements are sufficiently rich to reconstitute many closely related biomarker signals, without introducing artificial interdependency.

Secondly, we explored to what extent our DNAm surrogates of metabolomic features report signal compared to previously constructed DNA-methylation estimates, after regressing out age. Overall, modest correlations (mean|R|∼0.13, min|R|∼0.01, max|R|∼0.5) are observed between our DNAm-based metabolomics surrogates and DNAm-based multivariate clocks (Horvath, Hannum, PhenoAge, GrimAge, and bAge) (Fig. 4). Furthermore, correlations with other pre-trained DNAm-based surrogate molecular markers (GrimAge components and the 109 protein EpiScores) are generally modest, with a few notable exceptions. Particularly, the GrimAge surrogates DNAm-leptin and DNAm-adm, and 4 Episcores (2771.35 [Gene: IGFBP1], 4929.55 [Gene: SHBG], 3505.6 [Gene: LTα], CD6), present a relatively high positive correlation with HDL related surrogate markers and a relatively high negative correlation with the amino acids (DNAm-Leucine, DNAm-Isoleucine, and DNAm-Valine). We observe the inverse pattern for the GrimAge surrogate DNAm-PAI-1, and 4 other EpiScores (4930.21 [Gene: STC1], 2516.57 [Gene: CCL21], 3343.1 [Gene: ACY1], 3470.1 [Gene: SELE]) which also show a high correlation with VLDL surrogate markers. Notably, these correlations might suggest a link between the metabolome and protein markers related to immune signalling (LTα, CD6, CCL21, SELE), energy balance and metabolism-related hormones (IGFBP1, SHGB, ADM, leptin, and STC1), and atherosclerosis/thrombosis inhibitor (PAI1). Nonetheless, most of the markers exhibit limited correlations with their predecessors (Fig. 4), implying the presence of previously unexplored information in DNA methylation patterns.Fig. 4 Correlations with pre-trained DNAm scores after regressing out the ageing signal. Correlations between our DNAm metabolic features and previously trained clocks (Hannum, Horvath, PhenoAge, bAge, and GrimAge), the DNAm surrogates included in GrimAge and the 109 DNAm-based surrogates for proteins (EpiScores) by Gall et al. These correlations were calculated on the age regressed DNA methylation-based features (indicated as “AgeAccel”).

Metabolomic surrogates improve the mortality predictions of GrimAge

Next, we evaluated the DNAm-metabolomics features for their predictive value for all-cause mortality in the Rotterdam Study (RS-I and RS-II). We emphasise that this whole cohort was solely employed for evaluation purposes through the analyses. For this purpose, we utilised a total of 1544 samples from this cohort (mean age at baseline of 64 years, 251 deceased, and a median follow-up of 11 years; Fig. 1), by incorporating an additional 863 samples with available Illumina 450 k and mortality information, but not 1H-NMR metabolomics. In accordance with previous studies,11,14,15 we assess univariate Cox proportional hazard models (in years of follow up) adjusted for relevant covariates, specifically sex, together with age, BMI and cell counts levels at blood sampling (Fig. 5A and Figure S7A).

First, we evaluated our DNAm-MetaboHealth predictor, which showed a significant association with all-cause mortality (HR = 1.29, p = 5.57 × 10−04 [cox regression]) (Fig. 5a), in line with the original metabolomics-based MetaboHealth score (Figure S8C). Next, we evaluated the individual DNAm metabolomics features and observed significant associations with mortality for 15 out of the 64 surrogate metabolites, of which 6 (out of 14) metabolomics features were included in the original MetaboHealth score. In addition, our estimated effects were overall consistent with those found by the study performed by Deelen et al. which employed a considerably larger dataset of 44,168 individuals (Supplemental S7c), with our most significant findings being amongst their strongest effects. We observed an increased risk for higher estimates of 8 DNAm-based features, with the strongest being DNAm-Glucose (HR = 1.22, p = 7.29 × 10−03 [cox regression]), DNAm-Glycoprotein Acetyls (HR = 1.35, p = 1.74 × 10−05 [cox regression]). Conversely, we observed protective effects for 7 features, such as DNAm-PUFA% (HR = 0.81, p = 1.63 × 10−03 [cox regression]), DNAm-Histidine (HR = 0.82, p = 1.67 × 10−03 [cox regression]), DNAm-Valine (HR = 0.83, p = 1.6 × 10−02 [cox regression]), and DNAm-Albumin (HR = 0.83, p = 9.5 × 10−03 [cox regression]). In addition, 5 nominal significant additional metabolites showed discordant mortality associations between sexes. Explicitly, for males we observed mortality associations with DNAm-Glutamine, whereas DNAm-Tyrosine, DNAm-Leucine, DNAm-Total Fatty acids, and DNAm-PUFA associated with mortality in women (Figure S7B).Fig. 5 Associations with time to death. a) Significant univariate associations of the DNAm metabolomics features with time to all-cause mortality in RS (N = 1542 with 285 reported deaths). The associations are grouped based on the metabolomics groups coloured by the significant associations or the metabolites with mortality in Deelen et al. The asterisks (∗) separates nominal significant DNAm metabolomics features from the FDR significant ones [cox regression]. b) Stepwise Cox regression predicting of time to all-cause mortality optimised in RS, composed combining age, 3 DNAm surrogates included in GrimAge, and 9 DNAm metabolic models and 12 protein EpiScores. Finally, c) presents the ROC curves and the accuracies (AUC) at 5-years mortality for our newly developed clocks as compared to previously trained scores. Similarly, d) reports the Hazard Ratios for each of the Age accelerated multivariate scores, regressing also for sex, BMI and cell counts. In this last figure the scores as split for containing our DNAm metabolomics features (red) or not (blue).

Almost all the pre-trained DNAm clocks that we considered (Hannum, Horvath, PhenoAge, GrimAge) and some of their intermediate surrogates exert mortality associations in RS (Figure S9a–c). Next, we attempted to refine the current standard for biological age estimation, specifically GrimAge (CI = 0.79, p = 4.6 × 10−77 [cox regression]), by training a multivariate all-cause mortality predictor including our DNAm metabolomics features. As a first exploration, we trained a Cox regression model with age at blood sampling, sex, DNAm-GrimAge and DNAm-MetaboHealth, which only showed minor, but significant, improvements in the C-index (CI = 0.8, p = 4.6 × 10−77 [cox regression]) (Figure S10b). As a second exploration, we performed a stepwise (backward/forward) Cox regression model to identify a minimal set of features including age, sex, our 64 DNAm-metabolomic features and DNAm-GrimAge (CI = 0.81, p = 1.7 × 10−83 [cox regression]) (Figure S10d). Nonetheless, the best performing model was obtained when including age, sex, 3 out of 8 GrimAge components (predicting Leptin and ADM and TIMP_1), 12 of the 109 protein EpiScores, and 9 out of 64 DNAm metabolites (CI = 0.82, p = 1 × 10−85 [cox regression]) (Fig. 5b). The selected DNAm metabolic features included DNAm-Tyrosine, DNAm-S-VLDL-L, DNAm-S-LDL-L, DNAm-M-LDL-L, DNAm-APOB, DNAm-LA, DNAm-omega3 and DNAm-omega6, and DNAm-MUFA. Analyses of Schoenfeld residuals highlighted no violation of the proportional hazard assumption (Figure S14). In any case, all the newly introduced scores exhibited a significantly improved C-index (Figure S10e)39 and a higher AUC at 5 and 10 years compared to the GrimAge (Fig. 5c and d). Overall, this indicates that DNAm-surrogates from different origin, phenotypic, proteomic, or metabolomic, might confer mutually independent information for mortality prediction.

DNAm metabolomics models introduce relevant CpG selections

After establishing the value of our DNAm metabolomics features in predicting mortality, we explored the nature of the signal included in our models by investigating the CpG sites picked by the ElasticNET regression, which can shrink contributions of unnecessary features to zero. Predictors selected a median of ∼750 CpG-sites, with a minimum of 234 CpG-sites for DNAm-Acetoacetate and a maximum of 1569 for DNAm-ApoA1 (Fig. 6B, and Table S4). A total of 22,145 probes were included in at least 1 model. Comparison of the genomic positions of the selected CpG-sites with the rest of the 450 k array highlighted an underrepresentation of probes positioned in CpG Islands, and a preferential selection for CpG shelves and shores, known to be more dynamic areas (Fig. 6a).40,41 Noteworthy is the higher tendency to select CpGs co-locating with enhancers, cis-acting short regions of DNA that control the temporal and cell-specific activation of gene expression (Fig. 6a).37Fig. 6 CpG selections of the ElasticNET models. a) Log2 Odds ratio indicating the enriched in annotations of the CpG by our ElasticNET models. b) The central heatmap reports the log 10 p values [Fisher test] of the enrichments CpG sites selected by our models (rows) and the 50 most significant traits in the EWAS Catalog and Atlas enriched (rows). Bottom: the median coefficients in each DNAm model, and the number of CpGs per model. Right: the median coefficients given by our DNAm model to the overlapping CpGs with each trait. c) The nine most used probes (rows) over the 65 ElasticNET models (columns), coloured by metabolic groups. Top: The models were ordered by the mean accuracy over the test sets (CV, LLS, and RS). Right: The number of models which include each CpG and their nearest genes. d) Manhattan plot-like figure indicating the Variable importance of the single CpG probes in the DNAm metabolic models.

Functional enrichment analyses using the most proximal genes to the selected CpG-sites highlighted pathways associated to “developmental processes”, “cell differentiation”, and “regulation of metabolic processes” from Biological Processes in Gene Ontology (Figure S11c). Concomitantly, enrichment analyses of phenotypic annotations in the EWAS Catalog and EWAS Atlas (Fig. 6b), indicated that the CpG-sites are known to be largely related to peripheral tissue differentiation,42 foetal brain development43 and gestational age.44 Nonetheless, the CpG sites with the highest median coefficients across all our models were the ones annotated for metabolite-related traits, such as “Triglycerides”, and “Fasting Glucose” (Fig. 6b, Figure S11a and b). In total, 203 traits exhibited a significant enrichment for the CpG selections made by our models. Notably we find also highly significant associations with “Ageing”, and “all-cause mortality”, indicating that we do identify CpGs related to age-related processes.

Despite their interesting overarching signal, the DNAm-metabolic models show little overlap with each other in their CpG selections, with the majority showing overlaps well below 15%, apart for a few exceptions of highly correlated metabolites (e.g., 83% between DNAm-Total_cholines and DNAm-phosphoglycerides; Figures S6 and S12). Nonetheless, a handful of CpG probes were chosen in more than 30 models with largely consistent coefficient signs (Fig. 6c). Interestingly, while some of these 9 features have a higher importance weight on the DNAm metabolomics models (e.g., cg00574958, or cg06500161), others only exert a more minor influence (e.g., cg14938561, cg00461022). The 9 CpG sites with higher importance weight don't favour one specific metabolic group but seems to be relevant to many metabolic markers (Fig. 6d). Not surprisingly, also the nearest genes to these 9 probes are noteworthy. For instance, TXNIP, which includes cg19693031 (chosen in 43 DNAm-metabolomics models), was previously associated to hyperglycaemia and insulin resistance, and ABCG1, nearby cg06500161 (in 42 DNAm-metabolomics models) was associated to plasma lipid levels and stroke (Fig. 6d).

Discussion

A comprehensive quantification of biological ageing, as a way to assess the overall, holistic health status and disease susceptibility of individuals,14 would constitute a major advance for healthcare and preventive research. A diversity of molecular markers has been proposed as indicators of biological age relating to health- and lifespan. Here we integrated well-established DNA methylation-based and 1H-NMR metabolomics resources for biological age prediction with mortality as a primary endpoint. To our knowledge, the potential synergistic effects arising from combining these two molecular sources remained thus far largely unexplored, and we believe that a collection of models predicting metabolomics features may be relevant within the rapidly growing repertoire of DNA methylation-based estimates.10,11,15,45 A structured training and evaluation design aided us to demonstrate the robustness of our features. We highlighted the distinct signal expressed by our models and their feature selection. Finally, we explored the use of our DNAm-based surrogates of metabolomics features in combination with previously trained DNAm-based surrogates (e.g., GrimAge constituents) suggesting that these confer complementary information.

We applied ElasticNET regression, a widely utilised algorithm to train epigenetic clocks, to data of four large population cohorts to derive DNAm-based surrogates for a previously derived multi-analyte score indicating mortality (MetaboHealth), and for 64 individual metabolomics features. The direct estimation of metabolomics-based mortality by constructing a DNAm surrogate for the MetaboHealth score showed promising results (mean R in test-sets = 0.397). Moreover, we were able to construct DNAm surrogates for many, but not all, metabolomics features with good replication accuracies (mean R in test-sets >0.35), including health markers for HDL and VLDL metabolism, inflammation, and fluid balance. Less accurate were the DNAm surrogates for amino-acids, ketone bodies, glycolysis, and LDL-related markers (mean R in test-sets <0.2). Nevertheless, considering the limited number of available markers and the low accuracy thresholds previously used for DNAm scores (R >0.1 in test sets),15,45 we continued evaluating all 65 models. This decision was further corroborated by a previous report by Stevenson et al. who suggested that their DNAm surrogate for CRP was a more reliable indication of chronic inflammation than its measured counterpart, even when considering the modest correlation between CRP and its surrogate.46 Overall, our DNAm metabolomic features conveyed a signal coherent with the quantified metabolomics variables and independent from most of the previously reported DNA methylation-based clocks and molecular surrogates.

Great emphasis was given to the harmonization of metabolomic data collected across different cohorts, prior to training our DNAm-based models for individual metabolites or the MetaboHealth score. Non-biological variability that may originate from inter cohort differences in sample collection, storage, or handling could confound model training. Typically, this challenge in epidemiology is addressed by applying a z-scaling per cohort prior to conducting a meta-analysis, which in effect discards all differences, both technical and biological, between cohorts. In other words, while allowing to draw conclusions on the similarities in associations with endpoints between cohorts, this strategy does not allow for a direct comparison of the underlying molecular profiles between cohorts. To address this issue we applied a calibration technique, which we developed adapting methodologies previously applied in longitudinal studies.26 This calibration technique showed its merit in harmonizing the metabolomics profiles, while preserving the natural biological heterogeneity within and between the different study populations. Importantly, this approach allowed for an evaluation of the MetaboHealth score across cohorts, showing consistent age and sex specific trends per study, and global predictive power for established clinical variables, such as hsCRP and diabetes.

Previous studies have shown advantages of pre-selecting CpGs when training ElasticNET regression models.16,47, 48, 49 Following this example, we implemented a pre-selection of CpG sites showing a high variability and consistent association with the outcome of interest during the training phase of our 5-Fold Cross Validation procedure. During this selection process, we avoided imposing stringent criteria (it was based on nominal p value) and refrained from correcting for covariates to prevent an excessive loss of potentially meaningful signal. Approximately 22,000 CpG sites were included in at least one DNAm-based models. Enrichment analyses showed that the selected CpGs are more likely to be enhancers in CpG shelves and shores and are in the proximity of genes enriched for regulation of metabolic and developmental processes, or cell differentiation. This finding resonates with a longstanding hypothesis, that the ageing methylome reflects processes underlying intricate cellular and molecular changes linked with development and differentiation.50 Furthermore, CpG sites selected for our surrogates were also previously associated to age (e.g., Ageing, all-cause mortality), inflammatory (C-reactive proteins), or metabolically related traits (e.g., triglycerides and metabolic syndrome). Strikingly, we found a highly recurrent selection of 9 CpGs in at least 30 distinct DNAm surrogate models, suggesting that these CpGs form a fundamental link between the blood metabolome and DNA methylome. All these loci have been previously found associated with metabolic traits and processes,51 and most of these 9 CpGs and their nearest genes are considered powerful classifiers for diabetes stratification.52, 53, 54 Remarkably, 3 of these 9 CpG probes showed significant univariate association with mortality within the Rotterdam Study (Figure S8D). This reassures over the valuable cardiometabolic content latent in our DNAm models.

Besides, our main intent was to evaluate the possibility to extrapolate the mortality signal from the metabolome to DNA methylation. To do so, we tested which of our surrogates might be indicative of all-cause mortality in a subset of the Rotterdam Study (1544 persons, 285 deaths). Notably, we observed a successful detection, albeit partial, of the mortality signal exerted by the metabolomics platform. We could successfully derive a DNAm-based version of MetaboHealth, which significantly associates with all-cause mortality, although it showed a lower hazard ratio than the original score.14 This might in part be explained by the fact that only 6 of the 14 DNAm surrogates for the metabolites constituting the MetaboHealth showed associations with all-cause mortality. Overall, we observed significant associations with mortality for 15 out of 64 DNAm-based metabolites. The detected effects are consistent with the results previously reported by Deelen et al. in a large study using the original metabolomic features measured in 44.168 individuals. This consistency further underpins that DNAm surrogates for metabolomic features could potentially be leveraged as epigenetic markers of biological ageing.

To further explore this concept, we trained a multivariate model for all-cause mortality, that was allowed to select from all available DNAm surrogates using a stepwise forward/backward regression. This final model included 9 DNAm metabolomic features together with the competing covariates age, 3 GrimAge components and 12 plasma protein EpiScores. The resulting model combining DNAm surrogates from different origin showed a significantly improved mortality prediction (C-index = 0.82) compared to the GrimAge score (C-index = 0.79) (Fig. 5 and Figure S10). Our composite scores showed a substantial refinement of the AUC at 5 and 10 compared to the original GrimAge (Figure S10G and H). Overall, this suggests that a broader collection of DNAm-surrogates of independent origin, such as proteomics, phenotypes, and now also metabolomics, might confer a more comprehensive indication on epigenetic-based biological ageing.

An important limitation of the current study for leveraging mortality signals is its limited sample size, which is modest when compared to the large dataset that Deelen et al. employed to evaluate the mortality associations of the metabolomics features and to build a multi-analyte predictor for mortality. Despite the limited power, we found significant associations with mortality for the DNAm surrogates of the multi-analyte score MetaboHealth and 15 individual metabolic features, which were consistent with those observed by Deelen et al. A second limitation consists in the inclusion of only Dutch population cohorts, which does not ensure a correct replication in other populations. A third limitation is the usage of a single endpoint, mortality, for evaluating the potential applications of our DNAm surrogates as marker for biological age. We acknowledge that ageing and its associated decline in overall health is a complex multi-factorial process, that is only partially captured by mortality risk. Previous work reported the merits of the 1H-NMR metabolomics in estimating several different types of endpoints,7,13,55,56 or even end-of-life related-phenotypes such as frailty,14 leading us to speculate that our DNAm surrogates for metabolomic features might also be instrumental for capturing these ageing endophenotypes.

In conclusion, we have demonstrated that metabolite markers previously associated with mortality could be leveraged to help extract the mortality signal captured by the DNA methylation platforms. Moreover, we showed that our DNAm surrogates capture mortality signal that is independent of the mortality signal captured by previous DNAm scores, such as GrimAge or its separate DNAm surrogate constituents. Overall, this does suggest that even more mortality signal could be extracted given the availability of proper mortality-associated biomarkers.

Contributors

EBvdA, DB, MJTR and PES conceived and wrote the manuscript. DB performed the analyses. EBvdA and MJTR verified and supervised the analyses. PES, MB, JBJvM, JvD, DIB, RP, MG, LF were involved in data acquisition of the cohort data. All authors (DB, MJTR, LMK, MB, JD, JBJvM, JvD, RP, DIB, MG, LF, PES, EBvdA) discussed the results read, and approved the final version of the manuscript. Data used in this manuscript was mostly funded by the BBMRI-NL Metabolomics Consortium.

Data sharing statement

BBMRI-nl and BIOS-nl data are available upon request at https://www.bbmri.nl/services/samples-images-data. All DNAm metabolomics scores can be obtained with a script at: https://github.com/DanieleBizzarri/DNAm_metabolomics_scores.

Declaration of interests

Authors declare no competing interests.

Appendix A Supplementary data

Figure S1

Preprocessing of the metabolomics dataset. Percentages of a) missing values. c) zeros, and e) outliers, in each 64 metabolomics features in the 4 BIOS cohorts (VUNTR, RS, LLS_PAROFFS and LIFELINES). Bar-plots representing the number of b) missing values, d) zeros, and f) outliers in the samples divided per cohort. Finally, g) indicates the Z-score distributions of the 1168 imputed missing values in the dataset.

Figure S2

Distributions of the matches used to perform the calibration on the metabolomics dataset. The number of matching samples between a) b) and c) are the age distributions of the matching samples. d), e) and f) are the BMI distributions of the matching samples.

Figure S3

Comparisons of the metabolomics datasets before and after the calibration. tSNE plots coloured by biobanks, sex, and age before (respectively a, b, c) and after (respectively h, i, l) the calibration. kBET shows an improved mixing in the matching samples between LIFELINES and VUNTR (d before and h and after), LLS LIFELINES and LLS_PARTOFFS (e before and n after), and RS and LIFELINES (f before and o after). PVCA before g) and after p) calibration.

Figure S4

Calibrated and uncalibrated MetaboHealth. Box-plot comparing the a) calibrated and b) uncalibrated MetaboHealth values in each biobank, c) Bar-plots showing the differences in men and women in the uncalibrated MetaboHealth in the 4 cohorts. d) Observed mean values of age, BMI, eGFR, hsCRP and pressure and e) alcohol consumption, current smoking, and diabetes ordered following the uncalibrated MetaboHealth percentile over the entire BIOS population. On top the Spearman correlations (ρ) and its p value f) Correlation chart comparing the calibrated and uncalibrated MetaboHealth, with age, sex, BMI and biobanks. The upper triangle part of the figure indicates the mutual correlations and their relative p value.

Figure S5

DNAm metabolomics features accuracies divided in tertiles. Correlation of the DNAm metabolomics features with their quantified counterpart over the test sets (5Fold Cross-Validation test sets, LLS, and RS) divided in tertiles. Explicitly a) shows features with mean R >0.35, b) 0.2>mean R<0.35, c) mean R<0.2.

Figure S6

Metabolites intercorrelations and DNAm metabolites intercorrelations. Intercorrelations of the DNAm metabolomics features. A) Clustered intercorrelations between the DNAm metabolic features in the upper triangle and intercorrelations between the measured metabolomics features in the lower triangle. On the right there is also an heatmap indicating the number of CpGs used per model.

Figure S7

Univariate mortality associations in RS. a) Complete Univariate associations of the DNAm metabolomics features with time to all-cause mortality in RS (N = 1542 with 285 reported deaths). The associations are grouped based on the metabolomics groups and coloured by the significant associations or the metabolites with mortality in Deelen et al. b) Univariate associations of the DNAm metabolomics features split for sex. c) Comparisons of the univariate DNAm metabolomics associations to mortality to the associations of the metabolomics features in Deelen et al. d) Univariate mortality associations of the most used CpG sites in the DNAm metabolomics features. On the right side of each forest plot there are the p values for each univariate association [cox regression].

Figure S8

Mortality associations of the measured metabolomics in the 664 samples (99 deceased) of the Rotterdam Study with metabolomics and mortality data. a) Univariates associations with mortality of the metabolomics features divided in metabolomics groups. b) Comparison of the HRs of the metabolomics features in the RS with what previously reported by Deelen et al. c) Comparison of the univariate mortality associations with MetaboHealth and DNAm_MetaboHealth. On the right side of each forest plot there are the p values for each univariate association [cox regression].

Figure S9

Univariate mortality associations in RS of the pre-trained scores. a) Univariate associations of each of the evaluated DNAm-based clocks (GrimAge, PhenoAge, Hannum and Horvath). b) The mortality associations with the DNAm surrogates included in GrimAge. c) The mortality association for each of the 109 plasma protein EpiScores. On the right side of each plot there are the p values for each univariate association [cox regression].

Figure S10

Multivariate mortality models built in the Rotterdam Study. a) Cox regression using GrimAge. b) Cox regression including age, sex, DNAm MetaboHealth score and GrimAge (DNAm_MetaboHealth+DNAm_GrimAge). c) Stepwise cox regression built with the DNAm metabolomics features and the GrimAge score (DNAm_metabolites+ DNAm_GrimAge). d) Stepwise cox regression built with the DNAm metabolomics features and the DNAm proteins surrogates included in GrimAge (DNAm_metabolites+ DNAm_proteins). e) p values evaluating the significance of the improvement in the C-indices of the newly developed multivariate mortality models compared with GrimAge. f) ROC curves and the accuracies (AUC) at 10-years mortality for our newly developed clocks as compared to previously trained scores.

Figure S11

Enrichment analysis of the CpGs selected by DNAm models. a) Log10 p value [fisher test] of all the significant associations in the enrichment analysis over the selected CpG in the EWAS Catalog and Atlas. b) Median coefficients of the CpGs in the DNAm metabolomics models for the overlapping CpGs with the first 50 significantly enriched traits in the EWAS Catalog and Atlas. c) Log10 p value [fisher test] of the first 50 Gene Ontology traits enriched with the nearest genes to the selected CpG sites.

Figure S12

Percentage overlap of the CpG sites selected by each DNAm model.

Figure S13

Correlation of the projected epigenetic clocks in BIOS. a) Horvath clock, b) Hannum, c) PhenoAge, d) GrimAge, and e) bAge.

Figure S14

Schoenfield residuals for the model comprising DNAm metabolomics, GrimAge surrogates, and protein Episcores. Each figure represents the Schoenfield residuals for each of the features selected within our model, with the p value [Schoenfield] on the top, demonstrating the absence of patterns with time.

Supplementary Materials

BBMRI-consortium-author

Acknowledgements

This work was performed within the BBMRI Metabolomics Consortium funded by: BBMRI-NL (financed by NWO 184.021.007 and 184.033.111), X-omics (NWO 184.034.019), VOILA (ZonMW 457001001) and 10.13039/501100022216 Medical Delta (METABODELTA: Metabolomics for clinical advances in the Medical Delta). EBvdA is funded by a personal grant of the 10.13039/501100003246 Dutch Research Council (NWO; VENI: 09150161810095). Acknowledgements for all contributing studies can be found in the Supplementary Material-BIOS Consortium. Additional NTR samples were funded by the 10.13039/501100000781 European Research Council (ERC-230374) project Genetics of Mental Illness (DIB).

Appendix A Supplementary data related to this article can be found at https://doi.org/10.1016/j.ebiom.2024.105279.
==== Refs
References

1 López-Otín C. Blasco M.A. Partridge L. Serrano M. Kroemer G. Hallmarks of aging: an expanding universe Cell 186 2023 10.1016/j.cell.2022.11.001
2 Partridge L. Deelen J. Slagboom P.E. Facing up to the global challenges of ageing Nature 561 2018 45 56 30185958
3 Comfort A. Test-Battery to measure ageing-rate in man Lancet 294 1969 1411 1415
4 Blackburn E.H. Greider C.W. Szostak J.W. Telomeres and telomerase: the path from maize, Tetrahymena and yeast to human cancer and aging Nat Med 12 2006 1133 1138 17024208
5 Horvath S. DNA methylation age of human tissues and cell types Genome Biol 14 2013 R115 24138928
6 Peters M.J. Joehanes R. Pilling L.C. The transcriptional landscape of age in human peripheral blood Nat Commun 6 2015 8570 26490707
7 van den Akker E.B. Trompet S. Barkey Wolf J.J.H. Metabolic age based on the BBMRI-NL 1H-nmr metabolomics repository as biomarker of age-related disease Circ Genom Precis Med 13 5 2020 541 547 10.1161/CIRCGEN.119.002610 33079603
8 Menni C. Kiddle S.J. Mangino M. Circulating proteomic signatures of chronological age J Gerontol A Biol Sci Med Sci 70 2015 809 816 25123647
9 Zhang Q. Vallerga C.L. Walker R.M. Improved precision of epigenetic clock estimates across tissues and its implication for biological ageing Genome Med 11 2019 54 31443728
10 Levine M.E. Lu A.T. Quach A. An epigenetic biomarker of aging for lifespan and healthspan Aging (Albany NY) 10 2018 573 591 29676998
11 Lu A.T. Quach A. Wilson J.G. DNA methylation GrimAge strongly predicts lifespan and healthspan Aging (Albany NY) 11 2019 303 327 30669119
12 Deelen J. Kettunen J. Fischer K. A metabolic profile of all-cause mortality risk identified in an observational study of 44,168 individuals Nat Commun 10 2019 1 8 30602773
13 Nightingale Health UK Biobank InitiativeJulkunen H. Cichońska A. Slagboom P.E. Würtz P. Metabolic biomarker profiling for identification of susceptibility to severe pneumonia and COVID-19 in the general population Elife 10 2021 e63033
14 Kuiper L.M. Polinder-Bos H.A. Bizzarri D. Epigenetic and metabolomic biomarkers for biological age: a comparative analysis of mortality and frailty risk J Gerontol A Biol Sci Med Sci 78 2023 1753 1762 37303208
15 Gadd D.A. Hillary R.F. McCartney D.L. Epigenetic scores for the circulating proteome as tools for disease prediction Elife 11 2022 e71802
16 Bernabeu E. McCartney D.L. Gadd D.A. Refining epigenetic prediction of chronological and biological age Genome Med 15 2023 12 36855161
17 Zhernakova D.V. Deelen P. Vermaat M. Identification of context-dependent expression quantitative trait loci in whole blood Nat Genet 49 2017 139 145 27918533
18 Bonder M.J. Luijk R. Zhernakova D.V. Disease variants alter transcription factor levels and methylation of their binding sites Nat Genet 49 2017 131 138 27918535
19 van Dongen J. Gordon S.D. McRae A.F. Identical twins carry a persistent epigenetic signature of early genome programming Nat Commun 12 2021 5618 34584077
20 Soininen P. Kangas A.J. Würtz P. Suna T. Ala-Korpela M. Quantitative serum nuclear magnetic resonance metabolomics in cardiovascular epidemiology and genetics Circ Cardiovasc Genet 8 2015 192 206 25691689
21 Würtz P. Kangas A.J. Soininen P. Lawlor D.A. Davey Smith G. Ala-Korpela M. Quantitative serum nuclear magnetic resonance metabolomics in large-scale epidemiology: a primer on -omic technologies Am J Epidemiol 186 2017 1084 1096 29106475
22 Bizzarri D. Reinders M.J.T. Beekman M. Slagboom P.E. Bbmri-Nl null van den Akker E.B. 1H-NMR metabolomics-based surrogates to impute common clinical risk factors and endpoints eBioMedicine 75 2022 103764
23 van Iterson M. Tobi E.W. Slieker R.C. MethylAid: visual and interactive quality control of large Illumina 450k datasets Bioinformatics 30 2014 3435 3437 25147358
24 Hastie T. Tibshirani R. Narasimhan B. Chu G. impute: impute: imputation for microarray data 2023 10.18129/B9.bioc.impute
25 Zhou W. Laird P.W. Shen H. Comprehensive characterization, annotation and innovative use of Infinium DNA methylation BeadChip probes Nucleic Acids Res 45 2017 e22
26 Mäkinen V.-P. Karsikas M. Kettunen J. Longitudinal profiling of metabolic ageing trends in two population cohorts of young adults Int J Epidemiol 51 2022 1970 1983 35441226
27 Telle-Hansen V.H. Christensen J.J. Formo G.A. Holven K.B. Ulven S.M. A comprehensive metabolic profiling of the metabolically healthy obesity phenotype Lipids Health Dis 19 2020 90 32386512
28 Ala-Korpela M. Lehtimäki T. Kähönen M. Cross-sectionally calculated metabolic aging does not relate to longitudinal metabolic changes-support for stratified aging models J Clin Endocrinol Metab 108 2023 2099 2104 36658689
29 Büttner M. Miao Z. Wolf F.A. Teichmann S.A. Theis F.J. A test metric for assessing single-cell RNA-seq batch correction Nat Methods 16 2019 43 49 30573817
30 Li J. Bushel P.R. Chu T.-M. Wolfinger R.D. Principal variance components analysis: estimating batch effects in microarray gene expression data Batch effects and noise in microarray experiments 2009 John Wiley & Sons, Ltd 141 154
31 Bizzarri D. Reinders M.J.T. Beekman M. Slagboom P.E. van den Akker E.B. MiMIR: R-shiny application to infer risk factors and endpoints from Nightingale Health's 1H-NMR metabolomics data Bioinformatics 38 2022 3847 3849 35695757
32 Pelegí-Sisó D. de Prado P. Ronkainen J. Bustamante M. González J.R. methylclock: a bioconductor package to estimate DNA methylation age Bioinformatics 37 2021 1759 1760 32960939
33 Mukherjee S. Stamatis D. Bertsch J. Genomes OnLine Database (GOLD) v.8: overview and updates Nucleic Acids Res 49 2021 D723 D733 33152092
34 Higgins-Chen A.T. Thrush K.L. Wang Y. A computational solution for bolstering reliability of epigenetic clocks: implications for clinical trials and longitudinal tracking Nat Aging 2 2022 644 661 36277076
35 Battram T. Yousefi P. Crawford G. The EWAS catalog: a database of epigenome-wide association studies Wellcome Open Res 7 2022 41 35592546
36 Xiong Z. Yang F. Li M. EWAS open platform: integrated data, knowledge and toolkit for epigenome-wide association study Nucleic Acids Res 50 2022 D1004 D1009 34718752
37 Andersson R. Gebhard C. Miguel-Escalada I. An atlas of active enhancers across human cell types and tissues Nature 507 2014 455 461 24670763
38 Schröder M.S. Culhane A.C. Quackenbush J. Haibe-Kains B. survcomp: an R/Bioconductor package for performance assessment and comparison of survival models Bioinformatics 27 2011 3206 3208 21903630
39 Haibe-Kains B. Desmedt C. Sotiriou C. Bontempi G. A comparative study of survival models for breast cancer prognostication based on microarray data: does a single gene beat them all? Bioinformatics 24 2008 2200 2208 18635567
40 Ziller M.J. Gu H. Müller F. Charting a dynamic DNA methylation landscape of the human genome Nature 500 2013 477 481 23925113
41 Irizarry R.A. Ladd-Acosta C. Wen B. The human colon cancer methylome shows similar hypo- and hypermethylation at conserved tissue-specific CpG island shores Nat Genet 41 2009 178 186 19151715
42 Islam S.A. Goodman S.J. MacIsaac J.L. Integration of DNA methylation patterns and genetic variation in human pediatric tissues help inform EWAS design and interpretation Epigenetics Chromatin 12 2019 1 30602389
43 Spiers H. Hannon E. Schalkwyk L.C. Methylomic trajectories across human fetal brain development Genome Res 25 2015 338 352 25650246
44 Bohlin J. Håberg S.E. Magnus P. Prediction of gestational age based on genome-wide differentially methylated regions Genome Biol 17 2016 207 27717397
45 Chen Q. Dwaraka V.B. Carreras-Gallo N. OMICmAge: an integrative multi-omics approach to quantify biological age with electronic medical records bioRxiv 2023 10.1101/2023.10.16.562114
46 Stevenson A.J. McCartney D.L. Hillary R.F. Characterisation of an inflammation-related epigenetic score and its association with cognitive ability Clin Epigenetics 12 2020 113 32718350
47 Choi H. Joe S. Nam H. Development of tissue-specific age predictors using DNA methylation data Genes (Basel) 10 2019 888 31690030
48 Bergersen L.C. Ahmed I. Frigessi A. Glad I.K. Richardson S. Preselection in lasso-type analysis for ultra-high dimensional genomic exploration Frigessi A. Bühlmann P. Glad I.K. Langaas M. Richardson S. Vannucci M. Statistical analysis for high-dimensional data 2016 Springer International Publishing Cham 37 66
49 Croiseau P. Legarra A. Guillaume F. Fine tuning genomic evaluations in dairy cattle through SNP pre-selection with the Elastic-Net algorithm Genet Res 93 2011 409 417
50 Seale K. Horvath S. Teschendorff A. Eynon N. Voisin S. Making sense of the ageing methylome Nat Rev Genet 23 2022 585 605 35501397
51 Gomez-Alonso M.D.C. Kretschmer A. Wilson R. DNA methylation and lipid metabolism: an EWAS of 226 metabolic measures Clin Epigenetics 13 2021 7 33413638
52 Soriano-Tárraga C. Jiménez-Conde J. Giralt-Steinhauer E. Epigenome-wide association study identifies TXNIP gene associated with type 2 diabetes mellitus and sustained hyperglycemia Hum Mol Genet 25 2016 609 619 26643952
53 Krause C. Sievert H. Geißler C. Critical evaluation of the DNA-methylation markers ABCG1 and SREBF1 for Type 2 diabetes stratification Epigenomics 11 2019 885 897 31169416
54 Lai C.-Q. Parnell L.D. Smith C.E. Carbohydrate and fat intake associated with risk of metabolic diseases through epigenetics of CPT1A Am J Clin Nutr 112 2020 1200 1211 32930325
55 Buergel T. Steinfeldt J. Ruyoga G. Metabolomic profiles predict individual multidisease outcomes Nat Med 28 2022 2309 2320 36138150
56 Ahola-Olli A.V. Mustelin L. Kalimeri M. Circulating metabolites and the risk of type 2 diabetes: a prospective study of 11,896 young adults from four finnish cohorts Diabetologia 62 2019 2298 2309 31584131
