Introduction
Optimal treatment of
isocitrate dehydrogenase (
IDH)-wildtype glioblastoma is particularly challenging due to its infiltrative nature and aggressive behavior. Current standard treatment includes maximal safe resection and adjuvant combined radiochemotherapy for newly diagnosed glioblastoma [
18]. Despite this multimodality treatment, long-term local tumor control is achieved only in rare cases, and the vast majority of patients’ relapse. To date, there is no standardized treatment regimen for recurrent glioblastoma, and it is unclear which patients benefit from local or systemic therapies at this time [
44]. Ringel and colleagues reported the survival benefit of resection of recurrent glioblastoma and established recurrent surgery as an option for second-line therapy when this can be done safely [
39]. Additional studies confirmed the prolonged survival after recurrent surgery [
3,
50,
62]. Nevertheless, a major challenge in the search for a more effective surgical and therapeutic regimen is the extensive inter- and intratumoral heterogeneity, which is considered one of the main factors for treatment failure in glioblastoma [
36,
37].
Genome-wide DNA methylation profiling has recently emerged as a tool enabling accurate molecular classification of central nervous system (CNS) tumors and has the potential to further stratify patients according to their molecular pathological characteristics [
6,
7]. It allows the subdivision of glioblastoma into different methylation subclasses such as the most abundant
receptor tyrosine kinase I (
RTK I),
receptor tyrosine kinase II (
RTK II), and mesenchymal (MES) subclasses [
49]. Recently, methylation-based classification of glioblastoma has become increasingly important for predicting therapeutic response and aids in clinical decision-making for distinct subclasses [
11‐
13,
38,
60]. One issue that needs to be addressed when generally classifying glioblastoma by subgroups is heterogeneity within the tumor itself [
23,
34]. In addition, the extent of heterogeneity appears to influence patient prognosis [
23,
34]. Focusing on DNA methylation subclasses of glioblastoma, spatial heterogeneity of these subclasses has been analyzed in newly diagnosed glioblastoma, with studies reporting varying degrees of heterogeneity [
53,
59].
Drawing on these findings, an important consideration is whether temporal changes in DNA methylation subclasses occur during glioblastoma progression, what factors contribute to a potential subclass transition, and the extent to which such a transition influences patient outcome. To explore this, we conducted global DNA methylation profiling in 47 patients, comparing tissue and serum samples from the initial and recurrent surgeries.
Methods
Study population
Matched tissue samples from 47 patients diagnosed with
IDH-wildtype glioblastoma, who underwent initial and recurrent surgeries at University Medical Center Hamburg-Eppendorf and University Medical Center Freiburg (both Germany), were analyzed. Diagnosis was based on the current WHO classification [
29]. The extent of resection (EOR) was stratified into gross total resection (GTR), near-GTR, and partial resection (PR). A GTR was defined as a complete removal of contrast-enhancing parts, a near-GTR as a removal of more than 90% of the contrast-enhancing parts, whereas a resection of lower than 90% was defined as PR/biopsy. The EOR of contrast-enhancing parts was evaluated by magnetic resonance imaging (MRI) performed within to 48 h after index surgery. Overall survival (OS) was calculated from diagnosis until death or last follow-up, and progression-free survival (PFS) from diagnosis until progression according to Response Assessment in Neuro-Oncology (RANO) criteria based on local assessment [
58]. The progression-to-progression survival (PPS) was calculated from recurrence surgery until next progression, and progression-to-overall survival (POS) from recurrence surgery until death or last follow-up.
DNA methylation profiling
DNA was extracted from tumor tissue and bulk plasma and analyzed for genome-wide DNA methylation patterns using the Illumina EPIC (850 k) array. The tumor tissue of interest for performing DNA methylation was chosen by a board-certified neuropathologist of the Department of Neuropathology, University Medical Center Hamburg-Eppendorf, Germany. Processing of DNA methylation data was performed with custom approaches [
7].
Inclusion criteria
Methylation profiling results from first and recurrent surgery were submitted to the molecular neuropathology (MNP) methylation classifier v12.8 hosted by the German Cancer Research Center (DKFZ) [
6]. Patients were included if the calibrated score for methylation class family glioblastoma,
IDH-wildtype was > 0.84 at time of diagnosis [
7]. In addition, patients with a score below 0.84 but above 0.7 with a combined gain of chromosome 7 and loss of chromosome 10 or amplification of epidermal growth factor receptor (
EGFR) were included in accordance with cIMPACT criteria [
5]. Furthermore, a class member score of ≥ 0.5 for one of the glioblastoma subclasses was required. In addition, the following clinical criteria were defined: age above 18 years at time of diagnosis, local tumor progression at first recurrence, and availability of data for OS and PFS.
t-SNE embeddings
To compute t-SNE embeddings, first, 50 principal components (PCs) were computed for 25,000 of the most variable CpG sites. Subsequently, the class TSNE() from the Python package scikit-learn (v1.2.2) was fitted on the PCs with the perplexity argument set to 10, and was used to transform the PCs into t-SNE embeddings. The t-SNE results were visualized using the scatter function from the Python package Matplotlib (v3.8.2).
Differentially methylated probes and gene set enrichment
Differentially methylated probes in newly diagnosed glioblastoma with subclass transition were computed using the function dmpfinder from the R package minfi (v1.40.0). The resulting data frame was filtered with a p value cutoff of 0.01 and an absolute beta-value difference of 0.1. This led to the identification of 1962 differentially methylated probes (1415 hypermethylated probes and 547 hypomethylated probes). Gene annotations for these CpGs were extracted from the Illumina EPIC Manifest file (v1.0 B5). Only the following gene region feature categories were retained: TSS200, TSS1500, 1stExon, and 5’UTR. The resulting 493 hypermethylated and 113 hypomethylated genes were then used for gene set enrichment analysis based on GO biological processes (2023) with the R Package clusterProfiler (v4.2.0).
Copy number alterations
Genome-wide DNA methylation profiling was further used to analyze chromosomal copy number alterations and to provide information regarding gene amplifications, gains and losses as routinely done by the tumor classifier of the DKFZ. Genes and chromosomal regions were manually evaluated for differences in copy numbers from baseline and compared with other indicator genes on the array. For the assessment of relevant deviations, we followed the recommendations described by Capper et al. [
7].
Cell state composition analysis
To infer the abundance of cell type and cell state in the samples, we applied the Silverbush et al. deconvolution method [
47] to each sample that underwent bulk DNA methylation analysis using EPIC V2.0 arrays. The deconvolution method is a reference free method that uses a hierarchical matrix factorization approach inferring both cell types and the cell states therein. The method was trained on the DKFZ GBM cohort and tested on the TCGA GBM cohort and was shown to be able to infer the abundance of cell types in the microenvironment (immune cells, glia, and neurons) as well as malignant cell states (malignant stem-like cells component and two differentiated cells components). The stem-like component provides a deconvoluted estimation of the fraction of stem-like cells in a glioblastoma sample, specifically those exhibiting a transcriptional cell state of OPC-like or NPC-like [
32]. Both Differentiated 1 and Differentiated 2 represent the estimated fraction of differentiated cells in a glioblastoma sample, particularly those with a transcriptional cell state of MES or AC-like [
32]. We applied the method as described in Silverbush et al. using the cell type and cell state encoding provided in the manuscript and via the engine provided in EpiDISH [
63] package, with Robust Partial Correlations (RPC) method and maximum iterations of 2000.
Cell type deconvolution
Processing of methylation arrays
All Idats corresponding to methylation array data of tumor tissue and patients serum were processed similarly using the minfi package in R (version 1.40.0). The data was processed using the preprocessIllumina function. Only probes with detection p values < 0.01 were kept for further analysis. Also, probes with < 3 beads in at least 5% of samples, as well as all non-CpG probes, SNP-related probes, and probes located on X and Y chromosomes were discarded. The CpG intensities were converted into beta values representing total methylation levels (between 0 and 1).
Cell type deconvolution
Cell type deconvolution was applied to methylation arrays of tumor tissue and serum. Non-negative least square (NNLS) linear regression was used in deconvolving the beta values of methylation arrays into cell type components [
30,
41,
51]. As a reference, a publicly available signature was obtained from Moss et al. (2018) consisting of gene expressions for 25 cell type components (Monocytes_EPIC, B-cells_EPIC, CD4T-cells_EPIC, NK-cells_EPIC, CD8T-cells_EPIC, Neutrophils_EPIC, Erythrocyte_progenitors, Adipocytes, Cortical_neurons, Hepatocytes, Lung_cells, Pancreatic_beta_cells, Pancreatic_acinar_cells, Pancreatic_duct_cells, Vascular_endothelial_cells, Colon_epithelial_cells, Left_atrium, Bladder, Breast, Head_and_neck_larynx, Kidney, Prostate, Thyroid, Upper_GI, Uterus_cervix) and 6,105 unique CpGs[
30].
Proteomics
Proteomic processing of human glioblastoma samples
Formalin-fixed paraffin embedded (FFPE) samples of tumors were obtained from tissue archives from the neuropathology unit in Hamburg. Tumor samples were fixed, dehydrated, embedded in paraffin, and sectioned at 10 μm for microdissection using standard laboratory protocols. For paraffin removal FFPE tissue sections were incubated in 0.5 mL n-heptane at room temperature for 30 min, using a ThermoMixer (ThermoMixer
® 5436, Eppendorf). Samples were centrifuged at 14.000 g for 5 min and the supernatant was discarded. Samples were reconditioned with 70% ethanol and centrifuged at 14.000 g for 5 min. The supernatant was discarded. The procedure was repeated twice. Pellets were dissolved in 150 µL 1% w/v sodium deoxycholate (SDC) in 0.1 M triethylammonium bicarbonate buffer (TEAB) and incubated for 1 h at 95 °C for reverse formalin fixation. Samples were sonicated for 5 s at an energy of 25% to destroy interfering DNA. A bicinchoninic acid (BCA) assay was performed (Pierce
™ BCA Protein Assay Kit, Thermo Scientific) to determine the protein concentration, following the manufacturer’s instructions. Tryptic digestion was performed for 20 µg protein, using the Single-pot, solid-phase-enhanced sample preparation (SP3) protocol, as described by Hughes et al. [
20]. Eluted Peptides were dried in a Savant SpeedVac Vacumconcentrator (Thermo Fisher Scientific, Waltham, USA) and stored at -20° until further use. Directly prior to measurement dried peptides were resolved in 0.1% FA to a final concentration of 1 μg/μl. In total 1 μg was subjected to mass spectrometric analysis.
Liquid chromatography–tandem mass spectrometer parameters
Liquid chromatography–tandem mass spectrometer (LC–MS/MS) measurements were performed on a quadrupole-ion-trap-orbitrap mass spectrometer (MS, QExactive, Thermo Fisher Scientific, Waltham, MA, USA) coupled to a nano-UPLC (Dionex Ultimate 3000 UPLC system, Thermo Fisher Scientific, Waltham, MA, USA). Tryptic peptides were injected to the LC system via an autosampler, purified and desalted by using a reversed phase trapping column (Acclaim PepMap 100 C18 trap; 100 μm × 2 cm, 100 A pore size, 5 μm particle size; Thermo Fisher Scientific, Waltham, MA, USA), and thereafter separated with a reversed phase column (Acclaim PepMap 100 C18; 75 μm × 25 cm, 100 A pore size, 2 μm particle size, Thermo Fisher Scientific, Waltham, MA, USA). Trapping was performed for 5 min at a flow rate of 5 µL/min with 98% solvent A (0.1% FA) and 2% solvent B (0.1% FA in ACN). Separation and elution of peptides were achieved by a linear gradient from 2 to 30% solvent B in 65 min at a flow rate of 0.3 µL/min. Eluting peptides were ionized by using a nano-electrospray ionization source (nano-ESI) with a spray voltage of 1800 V, transferred into the MS, and analyzed in data dependent acquisition (DDA) mode. For each MS1 scan, ions were accumulated for a maximum of 240 ms or until a charge density of 1 × 16 ions (AGC target) were reached. Fourier-transformation-based mass analysis of the data from the orbitrap mass analyzer was performed by covering a mass range of 400–1200 m/z with a resolution of 70,000 at m/z = 200. Peptides with charge states between 2 + –5 + above an intensity threshold of 5,000 were isolated within a 2.0 m/z isolation window in top-speed mode for 3 s from each precursor scan and fragmented with a normalized collision energy of 25%, using higher energy collisional dissociation (HCD). MS2 scanning was performed, using an orbitrap mass analyzer, with a starting mass of 100 m/z at an orbitrap resolution of 17,500 at m/z = 200 and accumulated for 50 ms or to an AGC target of 1 × 105. Already fragmented peptides were excluded for 20 s.
Raw data processing
Proteomics samples were measured with liquid chromatography tandem mass spectrometry (LC–MS/MS) systems and processed with Proteome Discoverer 3.0 and searched against a reviewed FASTA database (UniProtKB: Swiss-Prot, Homo sapiens, February 2022, 20,300 entries). To cope with protein injection amount differences, the protein abundances were normalized at the peptide level. Perseus 2.0.3 was used to obtain log2 transformed intensities. The imputation was performed using the Random Forest imputation algorithm (Hyperparameters: 1000 Trees and 10 repetitions) in R, version 4.3.
Differential protein abundance
The protein abundances of the nine samples were corrected for the batch effect using ComBat [
24]. Protein level normalization was performed by centering protein abundances per file around the median value. For cases, measured multiple times, the mean value of abundances across all measurements was calculated. Values are log2 normalized (Peptide level), RF-imputed, batch-effect corrected normalized (Protein level) protein abundance (Mean value per case, across all measurements). For differential abundance analysis,
p values and log2 fold changes for each quantified protein (total 4716) were estimated using the empirical Bayes statistical test as implemented in the limma R package [
35].
Weighted correlation network analysis (WGCNA)
The WGCNA package in R (version 1.70.3) was used to identify protein co-expression modules [
28]. The minimum module size was set to 15 and a merging threshold of 0.40 was defined. Based on the assessment of scale-free topology, soft-power of 18 was selected. To construct modules, first we corrected for any technical batch effect using empirical Bayes-moderated adjustment using empiricalBayesLM function of WGCNA. Modules were assessed based on their correlation with traits (No transition and transition) and their levels of significance (associated with two-tailed Student’s
t test). The significant modules (
p < 0.05) were used for further analysis. All genesets within a module were used for overrepresentation analysis using clusterProfiler package in R (Version 4.2.0) [
61]. Further, to identify cell type enrichment within each module, enrichr API in python was used (Package maayanlab_bioinformatics, version 0.5.4 with PanglaoDB library available within the package) [
40]. To assess the module scores on single cells, Scanpy’s score_genes function was used to calculate module scores using core glioblastoma single-cell atlas [
27].
Detection of soluble factors in patient serum
Plasma from glioblastoma patients was taken 1 h prior primary surgery and then isolated by double spin centrifugation of whole blood. Samples were aliquoted and stored at − 80 °C before use. Soluble factors were detected using the LEGENDplex Neuroinflammation Panel 1 (Biolegend, San Diego, CA, USA) according to the supplier’s protocol were simultaneously determined using a multiplex bead-based immunoassay (LEGENDPlexTM Human Neuroinflammation Panel 1, Cat. No. 740795, Biolegend, USA). Data were acquired using the BD LSR Fortessa and Beckman Coulter Cytoflex LX flow cytometer and analyzed with the BioLegend LEGENDplex software.
DNA tumor purity
Tumor-purity was calculated using the RF_purify Package in R [
22]. This package uses the “absolute” method which measures the frequency of somatic mutations within the tumor sample and relates this to the entire DNA quantity [
8].
3D volumetric segmentation
We analyzed T1-weighted as well as T2-weighted FLAIR (Fluid attenuated inversion recovery) from MRI using the program BRAINLAB. To measure tumor volume, the tumor region of interest was delineated with the tool “Smart Brush” enabling a multiplanar 3D reconstruction. Volume of contrast enhancement and FLAIR hyperintensity was assessed in cm3.
Ethics statement
This study was approved by the medical ethics committee of the Hamburg chamber of physicians (PV4904). Informed written consent was obtained from all patients.
Statistical analysis
Gaussian distribution was confirmed by the Shapiro–Wilk normality test. For parametric data, unpaired two-tailed Student’s t test or one-way ANOVA with Tukey’s post hoc tests to examine pairwise differences were used as indicated. Survival curves were visualized as results from the Kaplan–Meier method applying two-tailed log rank analyses for analyzing statistical significance. In general, a p value less than 0.05 was considered statistically significant for all experiments. Statistical analyses were performed using SPSS Inc. (Version 29, Chicago, IL, USA). Data illustrations were performed using GraphPad Prism 10. Alluvial plots were graphed with R.
Discussion
A factor contributing to the aggressive behavior of glioblastoma is the spatiotemporal heterogeneity, which poses a major challenge in finding an optimal therapeutic regimen and targeting tumor cells [
37,
52,
57]. Although the spatial heterogeneity of DNA methylation subclasses in newly diagnosed glioblastoma has been described previously [
53,
59], their robustness and its clinical significance during disease progression remains of great interest. Our study provides the following important findings: (1) Temporal changes in DNA methylation subclasses were observed in 28.2% of tumors, with the majority of transitions occurring towards the mesenchymal subclass. (2) Transition to the dominant DNA methylation subclass was more frequent after incomplete tumor resection; however, there was no association with adjuvant treatment modalities or the time between initial and recurrent surgery. (3) Glioblastomas with subsequent mesenchymal transition demonstrated a higher stem cell-like state but a decreased immune cell state, as well as upregulated metabolic and catabolic processes at the time of diagnosis. (4) Despite a decreased immune cell state in tissue, significantly higher circulating immune signatures, as well as the cytokines IL-6 and IL-18, were observed in serum at the time of diagnosis in patients with a mesenchymal transition. (5) Conversely, matched tumor tissue at recurrence showed an increased immune cell state but decreased stem cell-like state after a mesenchymal transition, revealing an opposite cell composition than at the initial diagnosis. (6) The survival outcome of glioblastoma patients with a subclass transition was comparable to that of patients without such a transition.
Tumor profiling based on DNA methylation facilitates more accurate differentiation of brain tumor subgroups and allows subclassification of glioblastoma, with
RTK I,
RTK II, and MES being the most common subclasses [
6,
49]. For these DNA methylation subclasses, two studies investigated the spatial heterogeneity of malignant gliomas [
53,
59]. Wenger and colleagues analyzed 38 biopsies from 12 patients with newly diagnosed glioblastoma and reported intratumoral heterogeneity of the dominant DNA methylation subclass in 5 patients, demonstrating the existence of different DNA methylation subclasses within one tumor [
59]. A more recent study by Verburg et al. failed to confirm this degree of heterogeneity and reported more stable DNA methylation subclasses across tumors when analyzing 133 biopsies from 16 patients, 7 of whom diagnosed with glioblastoma [
53]. While these studies investigated spatial differences of DNA methylation subclasses, we focused on possible temporal heterogeneity and found that the dominant DNA methylation subclass changed in 28.2% of 39 tumors between first and recurrent surgery, with the transition to the mesenchymal subclass being most likely. These findings indicate that DNA methylation subclasses exhibit greater stability compared to transcriptional subtypes, since large-scale studies investigating the transcriptional glioblastoma subtypes reported about a transition of the dominant subtype in about 50.0% of the patients [
26,
52]. Varn and colleagues observed a most frequent switch to the transcriptional mesenchymal subtype which is in accordance with our findings [
52]. In the past years, the mesenchymal transition was investigated extensively and various drivers, such as radiation, and immune cell interactions were identified [
14,
17,
32,
43]. When investigating potential surgery- or treatment-related factors contributing to a DNA methylation subclass transition, we identified incomplete removal of the contrast-enhancing tumor. In contrast, treatment characteristics and time between surgeries as well as genetic alterations had no impact. To obtain a more comprehensive understanding of the epigenetic profile, we conducted an analysis of CpG sites in newly diagnosed glioblastoma. Through this analysis, we identified specific methylation signatures that were differentially abundant and associated with metabolic and catabolic processes in glioblastoma cases exhibiting a subclass transition. While the link between DNA methylation subclass heterogeneity and increased tumor metabolism remains unclear, changes in metabolic profiles have been identified as a hallmark of high-grade gliomas [
16] and have demonstrated prognostic impact in glioblastoma patients recently [
42]. An underlying mechanism could be the interaction with stem cells, as demonstrated previously [
45], which is also related to another finding of our study. Integrating cell state composition analysis showed a markedly increased stem cell-like state in tumors with mesenchymal transition at time of diagnosis, highlighting the importance of stem cells in transition processes and progress [
15,
52,
57]. This result is consistent with our finding that increased residual contrast enhancement after first surgery favors subclass transition since peripheral, contrast-enhanced tumor areas are considered as one niche of stem cells [
33]. However, at time of recurrence, there was a comparatively reduced stem cell percentage, which agrees with the results of Varn et al. [
52].
The interplay between stem cells and immune cells in the tumor environment has been illustrated several times [
15,
31,
46]. Here, we found a decreased immune state in newly diagnosed glioblastoma with mesenchymal transition, in opposition to the stem cell state, which is inverted in recurrent tumor tissue after mesenchymal transition. Interestingly, signatures of increased circulating immune cells in the serum of patients with a mesenchymal transition were already observed at the time of diagnosis, so that here a stem cell-associated immune evasion preceding mesenchymal transition during the course of the disease can be hypothesized and could be subject of future studies. The relevance of the change in immune cell composition observed here in tumors with a mesenchymal transition is consistent with two recent studies [
10,
56]. Furthermore, our study identified possible soluble factors which are differentially concentrated in patients serum at time of diagnosis. The cytokines IL-6 and IL-18 were highly elevated in serum of mesenchymal-transitioning tumors. This might be explained with an interplay with stem cells given the existing literature [
2,
9,
21,
55], and for both factors it appears reasonable to consider these as biomarkers in further studies with larger cohorts.
Since recurrent glioblastoma and especially the mesenchymal transcriptional subtype are considered particularly aggressive, we asked whether a DNA methylation subtype transition is relevant to patient survival [
14]. Previously, studies demonstrated worse PFS and OS in patients with a transition to the mesenchymal transcriptional subtype at recurrence [
26,
57]. This might be explained in conjunction to our study with the increased stem cell-like state at time of diagnosis in these transitioning tumors, as already demonstrated in recurrent gliomas [
52].
However, in our study, survival of patients with and without a subclass transition was comparable. The lack of prognostic relevance of methylation subclass transition is most plausibly consistent with previous studies that reported comparable survival between
RTK I,
RTK II, and MES tumors at the time of initial diagnosis [
11,
13,
25,
60]. Although methylation subclass does not appear to be a general prognostic marker, its clinical importance in predicting the probability of glioblastoma-associated seizure [
12,
38] and the benefit of surgical resection [
13] has been previously highlighted. This underscores the need for additional rapid [
1,
19] and intraoperative [
54] methods to detect the methylation subclass.
Publisher's Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.