DNA methylation-based age acceleration observed in IDH wild-type glioblastoma is associated with better outcome—including in elderly patients

Elderly patients represent a growing proportion of individuals with glioblastoma, who however, are often excluded from clinical trials owing to poor expected prognosis. We aimed at identifying age-related molecular differences that would justify and guide distinct treatment decisions in elderly glioblastoma patients. The combined DNA methylome (450 k) of four IDH wild-type glioblastoma datasets, comprising two clinical trial cohorts, was interrogated for differences based on the patients’ age, DNA methylation (DNAm) age acceleration (DNAm age “Horvath-clock” minus patient age), DNA methylation-based tumor classification (Heidelberg), entropy, and functional methylation of DNA damage response (DDR) genes. Age dependent methylation included 19 CpGs (p-value ≤ 0.1, Bonferroni corrected), comprising a CpG located in the ELOVL2 gene that is part of a 13-gene forensic age predictor. Most of the age related CpGs (n = 16) were also associated with age acceleration that itself was associated with a large number of CpGs (n = 50,551). Over 70% age acceleration-associated CpGs (n = 36,348) overlapped with those associated with the DNA methylation based tumor classification (n = 170,759). Gene set enrichment analysis identified associated pathways, providing insights into the biology of DNAm age acceleration and respective commonalities with glioblastoma classification. Functional methylation of several DDR genes, defined as correlation of methylation with gene expression (r ≤ −0.3), was associated with age acceleration (n = 8), tumor classification (n = 12), or both (n = 4), the latter including MGMT. DNAm age acceleration was significantly associated with better outcome in both clinical trial cohorts, whereof one comprised only elderly patients. Multivariate analysis included treatment (RT, RT/TMZ→TMZ; TMZ, RT), MGMT promoter methylation status, and interaction with treatment. In conclusion, DNA methylation features of age acceleration are an integrative part of the methylation-based tumor classification (RTK I, RTK II, MES), while patient age seems hardly reflected in the glioblastoma DNA methylome. We found no molecular evidence justifying other treatments in elderly patients, not owing to frailty or co-morbidities. Supplementary Information The online version contains supplementary material available at 10.1186/s40478-022-01344-5.


Introduction
Patients over the age of 70 represent an increasing fraction of individuals with glioblastoma (GBM) [40], yet elderly patients are often excluded from clinical trials due to the estimated poor prognosis and a short overall Page 2 of 18 Bady et al. Acta Neuropathologica Communications (2022) 10:39 survival (6-9 months). Usually, a shorter duration radiation therapy schedule in association with temozolomide (TMZ) is favored for elderly patients aiming to shorten the treatment duration in view of the overall short life expectancy. Hypofractionated radiotherapy (15 × 2.66 Gy) has been shown clinically equivalent to standard 30 × 2 Gy fractionation [43], in combination with TMZ chemotherapy this regimen has been superior to hypofractionated radiotherapy (RT) alone [42]. In tumors with a methylated promoter of the O6-methylguanine methyl transferase (MGMT) gene, exclusive or initial TMZ chemotherapy alone may be an alternative to radiotherapy. Median overall survival remains unsatisfactory with treatment regimens. The choice of the treatment regimen depends on the patient's biological rather than chronological age and frailty, and the tumor MGMT promoter methylation status [59]. Yet, a fraction of elderly patients have a more favorable outcome with disease control beyond 18-24 months. Little is known about the molecular make-up of GBM in older patients (7th and 8th decade of life) and whether it differs from a "general" GBM IDH wild-type (IDHwt) population.
We have reported previously from the Nordic trial for elderly GBM patients that the MGMT promoter methylation status is not associated with patient age [33]. The predictive value of the MGMT status for benefit from TMZ is supported by the aforementioned trials. Expectedly, mutations in the isocitrate dehydrogenase gene 1 (IDH1), a hallmark of gliomas in young adults, are rare in the elderly [33,60]. Moreover, no obvious GBM-specific genetic alterations have been associated with the age of the patients, once removing cases with IDH mutant (IDHmut) tumors that are recognized and classified as a distinct disease entity since 2016 [31]. Hence, systematic molecular characterization is required, in order to justify distinct treatment modalities not owing to frailty of the patients, or suggesting new, more adapted (targeted) approaches that can be tested in clinical trials in this elderly patient population. The DNA methylome of cells is known to contain age related information reflected in age associated epigenetic changes. These are thought to arise from innate biological mechanisms, such as deterioration of epigenetic maintenance resulting in changes of DNA methylation over time that allowed the construction of accurate predictors of chronologic age. Such predictors are known as epigenetic clocks or DNA methylation clocks [12,16,21,22,62]. Furthermore, certain biological features or environmental factors have been associated with an acceleration of the DNA methylation age (DNAm age) in the blood, certain tissues, and cancer. This difference between chronologic age of the patient and estimated DNAm age is called DNAm age acceleration, and has been found associated with obesity, and smoking, but also Down syndrome or cancer, as reviewed by Horvath and Raj [23]. It has been proposed that biological age acceleration may hold the potential for disease specific biomarkers or frailty measures for individuals [5,23,34]. In addition, the tumor DNA methylome comprises information on cell of origin, plus tumor development related alterations that allow accurate tumor classification and identification of new tumor entities [9,37]. Previous studies have reported on associations of tumor DNAm age with chronologic age of the patients in adult glioma and GBM [30,65]. However, these studies included glioma of variable grades and subtypes, including IDHmut glioma that constitute an epigenetically and clinically distinct entity, associated with a glioma CpG island methylator phenotype (G-CIMP +) [10,39,51].
In this translational research study we aimed at identifying age-related molecular differences that may unravel novel therapeutic options and guide rational management of elderly GBM patients. To this end we set out to investigate the DNA methylome of adult IDHwt GBM for age related differences, including the patients' age, DNAm age acceleration of the tumor, DNA methylationbased tumor classification (Heidelberg), and CpG methylation entropy. Further, we paid special attention to the potential implication of DNA damage response (DDR) genes, given the current standard of care, treating GBM patients with DNA damaging therapies. We stablished the methylome for two clinical trial cohorts treating newly diagnosed GBM patients [33,48,49], whereof one recruited elderly patients only [33]. The study included an additional two external GBM data sets [8,9], and retained only tumors meeting the criteria of GBM grade 4 of the recent update of the WHO CNS5 classification 2021 [32].

Patient samples
Two sets of newly diagnosed GBM samples were obtained from patients treated in clinical trials according to pre-specified clinical criteria: from the Nordic trial (n = 116), and from the trial conducted by the European Organisation for Research and Treamtent of Cancer (EORTC) and the National Cancer Institute of Canada (NCIC) (EORTC 26,981/NCIC CE.3), pooled with the samples from the Lausanne Pilot clinical trial (LN-Pilot; Biobank of the Brain and Spine Tumor Center, BB_031_BBLBGT, of the Centre Hospitalier Universitaire Vaudois, CHUV, Lausanne, Switzerland) (n = 219). The constituted clinical trial cohorts were restricted to patients for whom enough frozen or paraffin embedded tissue was available, thus excluding cases with diagnostic biopsies only. The cohorts overlap largely with those for which the MGMT promoter methylation status was reported in the original trials [19,20,33]. The Nordic phase III trial recruited patients 60 years or older with a WHO performance score (PS) [58] of 2 or less, allowing PS 3, if it was due to neurological deficits. Patients were randomized to one of two regimens of RT (60 Gy in 30 fractions, or 34 G in 10 fractions) or standard dose TMZ (200 mg/m 2 , days 1-5 every 28 days) [33] (trial registration number ISRCTN81470623). The LN-Pilot trial and the EORTC 26,981/NCIC CE.3 trial (NCT00006353) recruited patients between the ages of 18 and 70, and a WHO PS of 2 or less [48,50]. Patients in the uncontrolled phase II LN-Pilot trial received the TMZ/RT→TMZ standard regimen (the current standard of care), and patients in the pivotal trial EORTC 26,981/NCIC CE.3 were randomized to RT (60 Gy in 30 fractions) or TMZ/ RT→TMZ. Patients signed informed consent for translational research according to institutional and international guidelines and regulations.

DNA methylation analysis
For genome-wide DNA methylation analysis DNA was isolated from macro-dissected tumor tissue (frozen samples: DNeasy Blood & Tissue Kit, Qiagen; formalin fixed paraffin embedded (FFPE) samples: EX-WAX ™ Paraffinembedded DNA Extraction Kit, S4530; Merck KGAa) and quantified (Quant-iT ™ PicoGreen ® dsDNA Assay Kit, #P7589, Life Technologies). DNA samples were analyzed on the Human Methylation 450 K BeadChip (Illumina, San Diego CA, USA) at the Genomics platform of the University of Geneva. FFPE-derived tumor DNA samples were processed after passing a PCR-based quality control (Infinium HD FFPE QC Assay Protocol). DNA samples were subjected to bisulfite treatment (EZ DNA Methylation-Gold ™ Kit, Zymo Research) as previously described, and FFPE samples were analyzed in separate batches after pretreatment with the restauration kit as recommended (llumina) [3,26].

Data availability of own datasets
The datasets from the trial cohorts are available at the Gene Expression Omnibus database (GEO) (http:// www. ncbi. nlm. nih. gov/ geo/) under the accession numbers GSE195684 for the Nordic trial samples, and the samples from the LN-Pilot trial and the EORTC 26,981/NCIC CE.3 trial are available at GSE195640, or GSE60274. The latter comprises data from a subset of GBM samples of the EORTC-NCIC & Pilot trials and 5 non-tumoral brain tissue samples that we have previously published [26].
Methylation data from an additional 5 non-tumoral brain tissues used, is available under GSE104293 [3].

Preprocessing DNA methylation data
The CpG probes with detection p-values > 0.01, located on the sex chromosomes, or in SNPs were removed from each dataset. The functional normalization [13] for Illumina 450 k arrays includes noob (normal-exponential out-of-band) background correction, dye-correction (chemistry I vs II) and RUV-2 step (removing unwanted variation) based on control probes. This normalization was performed by the function preprocess-Funnorm from the R package minfi. DNA methylation was summarized by M-values [11]. The ComBat procedure [24] based on common CpG probes was used to aggregate the four datasets to limit experimental variation and batch effects across the four datasets.

Copy number variation
Copy Number status for each marker was assessed using the combined intensities for methylated and unmethylated signals and circular binary segmentation to detect copy number aberration (CNA) events as previously described [2]. The homozygous deletion status (HD) of CDKN2A was determined using copy number probe means (CpGs) located in the CDKN2A gene and applying a mixture model [55].

Sample purity
The tumor purity (HMpurity) of each sample was estimated as previously proposed [3] using the GBM TCGA datasets.

MGMT promoter methylation status
The DNA methylation status of the MGMT promoter and the MGMT score (logit-transformed probability) were determined using the MGMT-STP27 regression model based on HM-450 k data [2,4]. In brief, the

Molecular subtypes
The G-CIMP status was determined by unsupervised clustering (Ward's algorithm with Euclidean distance) as previously reported [3] and served as approximation for the IDH mutation status, as this information was not available for all samples in any of the datasets. All G-CIMP + samples were removed. Molecular subtypes were obtained using the classification procedure of central nervous system tumors based on the analysis of DNA methylation patterns [9] (version v11b4, www. molec ularn europ athol ogy. org). For the main analyses of this study, we considered only GBM classified as mesenchymal (GBM_MES), RTK I (GBM_RTK_I) or RTK II (GBM_RTK_II).

DNAage and age acceleration
DNA methylation-based estimates of age (DNAm age) were calculated using ElasticNet regression [66] using 353 CpG sites selected by the Horvath clock [21,23]. The DNA methylation data were calibrated before the computation of DNAm age as recommended [21]. The metric called Age Acceleration (Accel) was obtained by subtraction of chronological age (age) from DNA methylation age (DNAm age). For subsequent analysis, 306 clock probes were conserved after filtering and aggregation of the four datasets (detection p-values > 0.01).

DNA methylation entropy (HME)
Estimation of the DNA methylation entropy (HME) is given by the normalized Shannon entropy [46] adapted for methylation fraction (p) given by Beta values [16]. For the i th methylation marker, two states were possible: unmethylated (1 − pi) or methylated (pi) and the maximal entropy is given by log (2). The DNA methylation Entropy (HME) for N methylation markers is computed as follow: The HME metrics were defined for all CpGs (global HM-entropy) and for the 12 strata constituted by the Island regions (CpG islands, shores, shelves or open sea) and promoter location status (promoter or not in promoter). The HME table was described by variation partitioning [7] for age, age acceleration and GBM classification. The results are illustrated in a Venn diagram containing the variation fractions for the three supplementary variables.

Expression
Gene expression from RNA sequencing (RNA-seq) data (Level 3), from the TCGA GBM dataset selected for this study, was quantified for the transcript models using RSEM [29] and normalized within samples to a fixed upper quartile for TCGA. Further details are available at the DCC data portal of TCGA. Gene-level data were restricted to genes expressed in at least 70% of samples. The complete dataset was normalized by the VOOM procedure [28].

Detection of DMP
The associations of CpG-probes (CpGs) with age or age acceleration (Accel) were investigated by comparison of the model including the variable of interest (e.g. CpG ~ age or CpG ~ Accel) and the null model (e.g. CpG ~ 1) based on linear models with F-test. The Bonferroni procedure was used to account for multiple testing comparisons. A differentially methylated position (DMP) was defined as a candidate for which the q-value was less than 0.1.

Detection of functional methylation
The correlation of methylation with the expression level for each CpG-probe located in a gene promoter within 1500 nucleotides up-or down-stream of the transcription start site (TSS) was estimated by Spearman correlation test. A methylated position was defined as functional, when the q-value was less than 0.05 and the correlation coefficient was inferior to − 0.3 (negative effect on gene expression). The gene locations were based on the Homo sapiens (human) genome assembly GRCh37 (hg19) available at the Genome Browser of the University of California, Santa Cruz (UCSC build hg19, http:// genome. ucsc. edu) [25].

Cox regression model and other tests
For the continuous variables Wilcoxon test (t) or Kruskal and Wallis test (a) were used to test the differences between two or more groups. The independence between qualitative variables and groups was tested with Pearson's Chi-squared with Yates' continuity correction. Survival univariate and multivariate models were computed by Cox proportional hazards regression model [54]. The association tests with GBM classification, study origins and the interaction between these both variables were performed by model comparison using Wald's test with sandwich estimation of the covariance matrix and F-statistic. The covariance of the models is estimated by sandwich methods with the type version HC3 [64] to compensate heteroscedasticity. Principal component analyses (PCA) and the permutation multivariate analyses of variance (ADONIS) [1] using Euclidean distances were used to investigate the association between additional variables (e.g. age, age acceleration or GBM classification) and DNA methylation data. Analyses and graphical representations were performed using R-4.1.2 and the R package rms and survival [17,54].

Patient characteristics
Age-dependent DNA methylation features of GBM were investigated in samples from patients enrolled in clinical trials for newly diagnosed GBM. They comprised samples from the Nordic trial (n = 116), treating elderly patients, and the EORTC 26,981/NCIC CE/3 and the Lausanne Pilot (LN-Pilot) trials (n = 219). The DNA methylome was established on the 450 k platform. The baseline description of the full cohorts is presented in Additional file 1: Table S1. The DNA methylation-based tumor classification [9] and its association with age is visualized in Additional file 1: Fig. S1. Of note, the patient cohort constituted of the EORTC/NCIC & LN-Pilot trials exhibited more diversity in glioma subtypes than the Nordic cohort. Detailed annotation of all samples is provided in Additional file 2. Two additional external GBM datasets were included in the analyses. The age ranges of the patients in the four full GBM datasets, after removing all IDHmut /G-CIMP + samples, were as follows: Nordic (n = 115) 60 to 83 years (median 70, SD 4.8 years), EORTC/ NCIC & LN-Pilot (n = 195) 26 to 70 years (median 55, SD 9.39 years), TCGA (n = 113) 21 to 85 years (median 62, SD 11.9 years), and DKFZ (n = 235) 18 to 86 years (median 59, SD 13.3 years). The tumors of the four datasets were classified according to the Heidelberg DNA methylation-based classifier [9]. Most GBM were classified as MES, RTK I or RTK II (88% for TCGA, 91% for the other 3 datasets), with few samples belonging to GBM-MID or GBM-G34 (no GBM-G34 in Nordic, Additional file 1: Fig. S1) or others that are now considered distinct tumor entities in the updated WHO classification 2021, and few non-GBM classifications. This study was subsequently restricted to IDHwt GBM subtypes MES, and RTK I or II, for which the baseline characteristics of the patients are summarized in Table 1 for each dataset.
The distribution of the tumors in the three GBM subclasses was significantly different between the four datasets (p = 0.014, chi-squared). There were less MES GBM comprised in the DKFZ dataset as compared to the others (22% less), likely due differences in patient selection of the study. There were no differences in the proportion of female and male patients, or the frequency of MGMT promoter methylation between the datasets.
After initial filtering, 361,745 CpGs were retained for subsequent analyses. After batch correction no effect of the four datasets was observed on the DNA methylation based organization of the samples as illustrated in a PCA (R 2 = 0.002, p = 1.00, Fig. 1a-b). In contrast, the global organization of DNA methylation was significantly affected by GBM methylation subclasses (R 2 = 0.083, p = 0.01, Fig. 1c). Similarly, no differences were observed for sample purity between the data sets ( Table 1), but between the GBM subclasses (Wald's test, p < 0.001, Table 2; Fig. 1e).
The lowest purity was associated with the MES GBM subtype (Fig. 1e), in line with a more pronounced fraction of tumor infiltrating cells that has been associated with this subtype [57]. The Nordic trial, recruiting only elderly patients (Table 1), introduced a significant difference of age among the four datasets (Wald's Test, p < 0.001, Table 2). However, no age related association with the three GBM subtypes was observed (Wald's Test, p = 0.805, Table 2, Fig. 1d). Finally, no association was observed between the GBM subtype and the WHO performance score at study entry (PS, scale 0 to 4, with higher values indicating greater disability [58]) (Cochran-Mantel-Henszel chi-squared test with stratification by study origin, p = 0.720). These data was available for the patients treated in the EORTC 26,981/NCIC CE/3, the LN-Pilot study and the Nordic trial, respectively (Additional file 1: Table S1).

Age related differential methylation (DMP)
First, we analyzed the DNA methylation data for associations with the patients' age. Age dependent methylation identified 19 CpGs (p-value ≤ 0.1, after Bonferroni correction; Additional file 1: Table S2). Of these, ELOVL Fatty Acid Elongase 2 (ELOVL2) methylation at cg16867657 has previously been published as a biomarker for chronologic age (r = 0.92) [14] and is part of forensic age predictors [38,63]. Similarly, methylation levels of Tripartite Motif Containing 59 (TRIM59) have been associated with chronologic age [34,63]. Several of the probes were associated with cancer relevant genes.
Three CpGs met our criteria of functional methylation that we defined as negative correlation of methylation with expression (≤ − 0.3 and p-value adjusted for multiple testing ≤ 0.1) of the corresponding, annotated gene.

Table 1 Baseline description and parameters of GBM datasets
Description of datasets used in the analyses that were restricted to IDHwt GBM WHO grade 4, 2021, and corresponding parameters determined in this study. Significant differences between the datasets were observed for patient age and DNA methylation age of GBM (DNAm age), and GBM subtypes.
Age Acc, DNAm age acceleration; GBM subgroup, methylation-based classification into mesenchymal (MES), RTK I and RTK II; Global HM entropy, human methylation based entropy of all genomic regions; HME Prom CpG, human methylation based entropy in promoter CpG islands; HM purity, human methylation based determination of sample purity; MGMT promoter methylation status, unmethylated U, methylated M, classified by MGMT-STP27 procedure; MGMT score, calculated with MGMT-STP27; SD, standard deviation.   (Fig. 2). The small age range of this older population, may explain the latter ( Table 1). The age acceleration was significantly different to zero (t test, p < 0.001) with an averaged acceleration of 36.81 years and a standard deviation of 23.99 years. No differences in DNAm age acceleration were observed between male and female patients (p = 0.774, ANOVA, stratified by study). It is of note that the association between the numerous copy number variations (CNVs), characteristic for GBM, and the interrogated clock probes and DNAm age was weak. Only 7% of the variation of the DNAm age (p < 0.001) was explained, based on the regression of the four first axes from PCA of the clock probes, hence, excluding a strong confounding influence by CNVs.

Nordic (%) TCGA (%) P-Value
A large number of CpGs associated with age acceleration were identified (DMP age accel, n = 50,551). Most of the age related CpGs (16 of 19) were also associated with age acceleration. Although by construction, DNAm age acceleration is not associated with the patients' age (DNAm age acceleration = chronologic age of the patient minus DNAm predicted age of the tumor, see methods section). Interestingly, over 70% age acceleration associated CpGs (n = 36,348) overlapped with those associated with tumor classification (n = 170,759) into the three methylation-based GBM subtypes (MES, RTK I, RTK II) (Fig. 3a). Hence, it was not surprising that age acceleration was significantly associated with the GBM subtype (Wald's test, p < 0.001, Table 2; Fig. 1f ). The tumors classified as GBM RTK II, trended to exhibit higher age acceleration than the two other GBM methylation subclasses. However, the sample purity constituted a weak confounding factor, with a Spearman correlation coefficient between DNAm age acceleration and purity of r = 0.290 (p < 0.001), illustrated in Additional file 1: Fig. S2.
The overlap of the CpGs fulfilling all three criteria, age, DNAm age acceleration, and methylation-based subclassification, comprised 4 CpGs that also included the probe in the ELOVL2 promoter (Additional file 1: Table S2). This suggests that DNA methylation features of DNAm age acceleration are an integrative part of methylationbased tumor classification, while age seems only slightly reflected in the tumor DNA methylome. The breakdown of the selected CpGs by genome regions and function (e.g. promoter, gene body, enhancer, etc.); the three variables, age, DNAm age acceleration and classification, and affiliation with the Horvath clock, is detailed in Fig. 3c. It is of note that the observed significant association of 129 clock probes with the GBM classification was not dependent on their location on the 22 autosomes (Fisher's test with p-value estimated by Monté-Carlo simulation, p = 0.558).

Pathway analysis of DNAm age acceleration and tumor classification
Next, we were interested in the pathways to which the genes belonged that were associated with DNAm age acceleration or classification. Performing signature analysis using MiSigDB (molecular signatures database; gene set enrichment analysis [GSEA], adjusted p ≤ 0.1) and the CpGs associated with DNAm age acceleration yielded 1220 pathways. The top pathways were dominated by gene-sets characterizing epigenetic properties of neural precursor cells, or gene sets for midbrain neurotypes, and other progenitor cells, and axon development (Fig. 3d) invoking developmental features and cell of origin. The tumor classification associated CpGs yielded 23 pathways, of which most (n = 20) overlapped with those from DNAm age acceleration, in line with the large overlap of CpGs between the two criteria.
The common pathways were also dominated by gene sets characterizing epigenetic features of neural progenitors, included sets for transmembrane transporters and channels, and cancer gene sets (Module 11 and 64) (Fig. 3d). The few CpGs associated with age were linked with gene sets of downstream targets of STAT5B (Signal Transducer And Activator Of Transcription 5) and RSRFC4 (alias of MEF2A, Myocyte-Specific Enhancer Factor 2A), with 3 of the 4 genes overlapping in the two sets, the latter was also associated with DNAm age acceleration.
Subsequently we looked only at the functional CpGs (Fig. 3c), as they may shed light on associated biological mechanisms, and determined the associated pathways (GSEA, adjusted p ≤ 0.1). The "functional" pathways associated with DNAm age accel (n = 167) mostly (119, 71%) overlapped with those associated with classification (n = 294) (Additional file 1: Suppl Fig. S2). The overlapping "functional" pathways were dominated by gene-signatures characterizing epigenetic properties of neural precursor cells, or gene sets for midbrain neurotypes, and other progenitor cells [27,35,36]. These comprised gene sets associated with specific chromatin marks (histone H3 dimethylation mark at K4, H3K4me2, open chromatin, and trimethylation mark at K27, H3K27me3, repressive mark) associated with neural progenitor cells (NPC) [35,36], and distinct gene signatures derived from single cell sequencing of human embryonic midbrain cells [27]. In addition, signatures of immune cells, and some cancer related signatures were comprised, the latter including the Verhaak expression signature for mesenchymal GBM (Additional file 1: Fig. S2). Interestingly, but not surprising, among the top 30 functional pathways associated with the methylation-based classification comprised three signatures of the expression based GBM classifier defined by Verhaak [56] (signatures for mesenchymal, classical, and proneural GBM). These observations support the overall consistency and biological relevance of the findings. Fig. 3 Associations of DNA methylome with age, DNAm age acceleration, and GBM subclassification. Methylome-wide association studies of the samples were performed for age (green), DNAm age acceleration (pink) and GBM classification (blue) and visualized in a Venn-diagram, for the interaction between the three sets of candidate CpG-probes a and the associated genes b. The genomic location of the selected CpGs, their association with age, DNAm age acceleration, GBM classification, the CpG functionality, and association with the Horvath clock are summarized in c.
Gene set enrichment analysis (GSEA), using the MSigDB database, was established for candidate CpG-probes related to age, DNAm age acceleration and GBM classification. The list of gene sets d corresponds to the pathways significantly enriched for age and GBM classification. The number (NB) of genes per gene set is indicated and the p-value is represented by the color code. The gene lists overlapping with those associated with functional methylation (Additional file 1: Fig. S3) are highlighted in bold Human methylation entropy (HME) and age Subsequently, we investigated the genome-wide variation of DNA methylation that may affect regulatory functions and genomic/epigenomic stability of the tumors. For this purpose, we used the measure of the Human Methylation Entropy (HME) that quantifies the methylation complexity for a given CpG or a given genomic region [16]. Low heterogeneity/high similarity corresponds to low entropy  scores (range 0 to 1; e.g. 100% methylation, 0 entropy; 50% methylation, organized in fully methylated and fully unmethylated alleles, 0.25 entropy; 50% methylation, organized randomly, 0.65 entropy, for more details see [45]). The entropy at 12 distinct genomic regions (promoter, shore, etc.) tested, showed characteristic features by region, with strikingly lower entropy levels within CpG promoter islands (Fig. 4). This may suggest that CpG promoter islands present more homogenous DNA methylation patterns, e.g. methylated or unmethylated states, likely due to their direct regulatory function in gene expression. HME at different genomic regions (Additional file 1: Table S3) and global HME were associated with GBM subtype (Wald's Test, p < 0.001, Table 2, Fig. 4e). Indeed, HME was highest in RTK I for most genomic regions. In contrast, no difference between tumor subtypes was detected for HME in the promoter CpGs Islands (Wald's Test, p = 0.917, Table 2, Fig. 4f ). No differences were observed between datasets at the distinct regions and over all regions combined. Finally, we performed variation partitioning of HM-entropy metrics to evaluate the contribution of age, age acceleration and GBM classification. This revealed some weak association of HM-entropy with GBM classification, explaining 13% (R 2 = 0.131) of the variance, while age had very little impact (R 2 = 0.005), and the contribution of DNAm age acceleration was also small (R 2 = 0.050).

DNA damage response (DDR) and DNAm age acceleration
Since all GBM patients, including the elderly, are treated with genotoxic therapy we had a closer look at the involvement of DNA damage response (DDR) genes in DNAm age acceleration. Evaluating only CpG probes associated with the promoter region of DDR genes as input (list as defined by Pearl et al. [41], 3947 CpGs, associated with 403 genes that might reveal age dependent treatment responses, did not yield any age related candidate CpGs.  Table S4). DMPs associated with tumor classification consisted of 272 CpGs in 138 DDR genes whereof 22 in 12 genes were functional (Additional file 1: Table S5). Four genes with functional CpGs were associated with both DNAm age acceleration and tumor classification and comprised CDKN2A, MGMT, the gene for the Helicase like Transcription Factor (HTLF), and SMC1B, as annotated in Fig. 5. CDKN2A is a prominent tumor suppressor gene, frequently affected by homozygous deletions (HD) in IDHmut GBM [8], Therefore we tested the methylation of the CpG-probe cg07562918 (M-value) located in the CDKN2A gene that was significantly associated with age acceleration and tumor classification, for association with the HD-status. Methylation of this CpG was significantly associated with HD of the CDKN2A gene (Wilcoxon test p < 0.001). However, no significant direct effect of the HD status of the CDKN2A gene was identified on DNAm age acceleration (Wilcoxon test, p = 0.163; Additional file 1: Fig. S4).

Age acceleration and outcome
Finally, we integrated DNAm age acceleration into the multivariable model for outcome for the two clinical trial cohorts, and included the following parameters: treatment (TMZ/RT→TMZ for EORTC/NCIC & LN-Pilot; TMZ for Nordic), MGMT status (MGMT methylation), and the interaction between treatment and MGMT methylation (predictive factor). DNAm age acceleration was significant in both cohorts, EORTC/NCIC & LN-Pilot p = 0.0138 and Nordic p = 0.00161 (Table 3). The predictive value of MGMT methylation (interaction between treatment and MGMT methylation) was confirmed in the EORTC/NCIC & LN-Pilot cohort (p < 0.0001). The interaction term in the Nordic cohort was not significant (p = 0.09). However, the number of MGMT methylated patients in the TMZ-arm was very small (n = 13), suggesting lack of power for this test.

Discussion
In this study, we investigated the DNA methylome of IDHwt GBM for age related differences that could hint at predictive or prognostic factors or indicate particular vulnerabilities for treatments in elderly patients. However, we found no strong direct associations, although we identified the methylation probe in ELOVL2 that has reportedly the highest association with age as a single marker [5,14] and TRIM59 methylation, another robust marker for aging [34,63]. ELOVL2 is a GBM-relevant gene [15,44]. However, methylation of this probe was not functional according to our analysis. TRIM59 was among the genes associated with functional methylation. It encodes a protein with ubiquitin-transferase activity, and it has been associated with various regulatory processes and maybe involved in innate immune regulation.
We then used the metric of DNAm age and age acceleration proposed by the Horvath DNA methylation clock that has been trained on multiple tissues using the HM-27 k array (does not comprise CpG probes associated/annotated with ELOVL2), and has been reported to be highly accurate to predict chronologic age (r = 0.92) [14]. This methylation clock is considered independent of tissue type and mitotic potential [21,23]. Various epigenetic clocks have been developed with potential as biomarkers. The aim is not only to determine accurate chronologic age, but to develop biological clocks for specific purposes, such as tissue specific sensors of disease, risk of disease, including cancer risk, or external stress, such as smoking history, or all-cause mortality (Grim-Age) (reviewed in [5,12,23,34]).
DNAm age acceleration observed for GBM averaged at almost 40 years, with a wide variability. No association was found with the patients' sex.
The acceleration associated gene sets, identified using the MiSigDB, comprised several distinct cell type-specific signatures obtained by single cell sequencing of cells derived from embryonic ventral midbrain [27]. Other signatures reflected epigenetic features at developmentally regulated genes with high CpG density promoters that are characterized by bivalent histone marks (both active Heatmap of DNA methylation of functional CpGs in DDR genes significantly associated with DNAm age acceleration. The clustering was obtained by Ward's algorithm using Euclidean distance based on centered and scaled data. The samples and the methylation of the functional CpGs associated with DDR genes were classified into 3 clusters (consensus k-means clustering for 200 repetitions) and 5 clusters (consensus k-means clustering for 200 repetitions), respectively. The samples are annotated with the patients' age (age in years), DNAm age acceleration (Accel, in years), MGMT methylation score (MGMTscore), Purity index based on DNA methylation (HMpurity), study origin (study) and GBM classification (Class). The probes are annotated for statistical significance associated with age (Age), DNAm age acceleration (Accel) and GBM classification (Class); yes, red; no, blue. The Spearman's correlation of methylation with expression (measure of functionality) is given for each CpG, Corr (E/M). The clustering is dominated by the high correlation of DNA methylation among the functional CpGs of individual genes H3K4me3 and repressive H3K27me3 marks) that are poised but silent in ES (embryonic stem) cells and active in NPC, such as the genes encoding the Oligodendrocyte Transcription Factor 1 and 2 (Olig1 and Olig 2). Changes in the pattern of histone marks have a strong influence on DNA methylation. The signatures identified here, associated with DNAm age acceleration reflected a change, with loss of the active mark H3K4me3 and retention of H3K27me2 and/or the repressive mark H3K27me3 that has been associated with partial increase of DNA methylation [35]. Hypermethylation of such regulatory genes that are active under differentiation are found preferentially aberrantly methylated and thereby silenced in cancer. Such genes are reported to be regulated by thritorax-group and/or polycomb-repressive complex 2 (PRC2) proteins. Furthermore, CpGs whose methylation increases with age at specific locations, have been found overrepresented near polycomb target genes and bivalent chromatin domains. These genes are of importance for stemness and cell differentiation and have been found frequently methylated in cancer [21,53]. Beside these insights into epigenetic mechanisms associated with age acceleration and the neurotypes specific signatures [27,35,36], suggestive of developmental features and cell of origin, the analyses also yielded signatures for immune cells, cell death, and cancer. The biological relevance is further supported by the fact that a subset of these signatures overlapped with those associated with functional methylation, implicating that the observed increased methylation was negatively associated with expression of the corresponding genes, affecting related pathways ( Fig. 3d; Additional file 1: Fig. S2). The investigation of mechanisms and biological consequences reflected in methylation aging clocks is an active field of research that has been extensively reviewed elsewhere [5,12,47]. When associating the observed DNAm age acceleration of the GBM with patient outcome, we observed a significant effect, when analyzing the cohorts treated in clinical trials. Increased age acceleration in these IDHwt GBM was associated with better outcome, raising the question of the biological meaning of DNAm age acceleration in the context of these tumors. Therefore, we consider tumor related DNAm age acceleration as a measure for epigenetic distance associated with tumor development and GBM subtype. This is in accordance with the observation that DNAm age acceleration associated Table 3 Multivariable Cox regression models including age acceleration, GBM classification, global HM-entropy (HME) and HME at promoter CpGs Cox models, adjusted for Treatment (TRT), methylation status of MGMT promoter (MGMT), and the interaction of these two variables (TRT X MGMT) HME, human methylation entropy; HR, hazard ratio; MGMTmeth, methylated MGMT; Pr(( >|z|), p-value for z-statistics; P-value <0. 5 probes overlap largely with those associated with methylation-based GBM classification. It has been suggested that deterioration of epigenetic maintenance contributes to age related changes [23], and more recently, it has been proposed that DNA breakinduced epigenetic drift may contribute to aging [18]. Given the gross structural changes observed in GBM this is an interesting hypothesis. However, the numerous CNVs, characteristic of GBM, only explained a minor part (7%) of DNAm age of the tumors. Along the same lines, while a subset of clock probes was significantly associated with the methylation-based GBM classification, this was not dependent on their chromosomal location. Hence, this measure seems to comprehensively capture GBM related epigenetic changes that are associated in part with the GBM subclassification. The DNA methylation-based GBM classification is constituted of multidimensional information reflected by DNA methylation that are contributed by biological features such as purity, age acceleration, HM-entropy (variability of the DNA methylation), and other factors. The purity, reflecting infiltrating immune cells, measured by means of a DNA methylation variation, contributes to the differentiation of MES GBM samples, while DNAm age acceleration facilitates the differentiation of the RTK II GBM samples from those classified as RTK I and MES, and the HME discerns RTK I from the two others. Interestingly, among the significant pathways associated functional methylation and tumor classification, we identified previously reported expression-based GBM classification signatures [56]. This underlines the strong functional implication of DNA methylation on the expression phenotypes of the tumors and the coherence of the results.
The HM entropy metrics (12 distinct regions, and global entropy) exhibited genome-location dependent variation explained in part by the GBM subtypes (13%), but was basically not affected by age (0.5%). As an exception, HM entropy at promoter associated locations was lowest, and was not different across the GBM subtypes, indicating little permissiveness for variation in accordance with direct functional implications of methylated versus unmethylated status of promoters on gene expression.
Finally, we were interested in the DDR genes with functional methylation associated with DNAm age acceleration and classification that may yield mechanistic insights in the context of the genotoxic treatments that are part of the standard of care. Enhanced POLE4 methylation (8 functional probes) was associated with DNAm age acceleration (Additional file 1: Figure S3). Deletion of Pole-4 that encodes an accessory subunit of the DNA polymerase epsilon complex, has been reported to have no strong effects on its own in worms, however, in absence of the gene encoding the regulator of telomere elongation helicase 1 (rtel-1) apparently led to synthetic lethality due to impaired homologous recombination (HR) [6]. Similarly, six functional probes of MGMT were associated with DNAm age acceleration, whereof three were also associated with GBM subclassification. Functional probes in CDKN2A, HLTF, and SMCB1 were also associated with both. Other functional probes were only associated with tumor classification, e.g. the gene encoding the Fanconi anemia complementation group M protein (FANCM) or the gene encoding the AlkB Homolog 1, Histone H2A Dioxygenase (ALKBH1). The former is involved in homology directed DNA repair, and the latter takes part in the repair of DNA alkylation damage and has recently been associated with the regulation of the level of N 6 -methyladenine (N 6 -mA) DNA modifications, implicated in epigenetic regulation of gene expression relevant in GBM [61]. While MGMT methylation is a known predictive factor for responsiveness to the alkylating agent TMZ in GBM, it remains to be explored, whether any of the other identified probes and their associated genes and pathways, indicate potentially actionable vulnerabilities.
In conclusion, our efforts to identify epigenetic differences in GBM of elderly patients revealed that once all high grade gliomas not classified as GBM IDHwt WHO grade 4 are removed, the distribution and spectrum of the GBM subtypes seems to be comparable across adult age, at least from a DNA methylation point of view (Fig. 1). DNAm differences of the tumors quantified as age acceleration or HME overlap with features of tumor classification, while age is hardly reflected. Interestingly, the epigenetic distance measured as DNAm age acceleration was associated with better outcome in the cohorts treated in clinical trials also for elderly GBM patients. However, our analyses yielded no molecular evidence for age related differences that would advocate different treatment modalities in elderly GBM patients.