DNA methylation heterogeneity attributable to a complex tumor immune microenvironment prompts prognostic risk in glioma

ABSTRACT

Gliomas are malignant tumours of the human nervous system with different World Health Organization (WHO) classifications, glioblastoma (GBM) with higher grade and are more malignant than lower-grade glioma (LGG). To dissect how the DNA methylation heterogeneity in gliomas is influenced by the complex cellular composition of the tumour immune microenvironment, we first compared the DNA methylation profiles of purified human immune cells and bulk glioma tissue, stratifying three tumour immune microenvironmental subtypes for GBM and LGG samples from The Cancer Genome Atlas (TCGA). We found that more intermediate methylation sites were enriched in glioma tumour tissues, and used the Proportion of sites with Intermediate Methylation (PIM) to compare intertumoral DNA methylation heterogeneity. A larger PIM score reflected stronger DNA methylation heterogeneity. Enhanced DNA methylation heterogeneity was associated with stronger immune cell infiltration, better survival rates, and slower tumour progression in glioma patients. We then created a Cell-type-associated DNA Methylation Heterogeneity Contribution (CMHC) score to explore the impact of different immune cell types on heterogeneous CpG site (CpG^ct) in glioma tissues. We identified eight prognosis-related CpG^ct to construct a risk score: the Cell-type-associated DNA Methylation Heterogeneity Risk (CMHR) score. CMHR was positively correlated with cytotoxic T-lymphocyte infiltration (CTL), and showed better predictive performance for IDH status (AUC = 0.96) and glioma histological phenotype (AUC = 0.81). Furthermore, DNA methylation alterations of eight CpG^ct might be related to drug treatments of gliomas. In conclusion, we indicated that DNA methylation heterogeneity is associated with a complex tumour immune microenvironment, glioma phenotype, and patient’s prognosis.

KEYWORDS:

• Glioma
• DNA methylation heterogeneity
• tumor immune microenvironment
• prognosis
• T-lymphocyte infiltration
• clinical drug treatment

Introduction

Glioma is an extremely serious type of brain tumour, and glioblastoma (GBM) is the most aggressive form with higher cellular and molecular heterogeneity [1–5]. Lower-grade gliomas (LGG) have the potential to progress to higher-grade tumours over time [6,7]. GBM has a median overall survival of approximately 15 months [8–10]. LGG has a better prognosis than GBM, but treatment is also challenging [11,12]. The diversity of GBM may be influenced by multiple factors including genetic, epigenetic, and tumour immune microenvironment [13]. Tumour heterogeneity is associated with enhanced genomic proliferation, altered tumour cell growth morphology, and inconsistent gene mutation frequency across cells [14,15]. High intra- and inter-tumour heterogeneity contribute to drug resistance and even treatment failure in some GBM patients [16]. Differences in DNA methylation patterns between cells affect different cancer phenotypes, giving rise to widespread tumour heterogeneity in clinical and diagnosis [17,18].

Several studies have shown that 60–80% of CpG sites in the human genome are methylated [19,20]. Approximately 2% of the 26.9 million CpG sites show intermediate DNA methylation signals [21,22]. Methylation of a single CpG site in a cell type is either fully methylated or unmethylated [23]. Mixing of different cell types in tumour tissues may result in elevated intermediate methylation, reflecting DNA methylation heterogeneity in the complex tumour immune microenvironment [23]. Some quantification methods of DNA methylation heterogeneity, such as the Proportion of Discordant Reads (PDR) and epiallele, capture the methylation status of individual CpG sites in different cells, facilitating the analysis of DNA methylation heterogeneity within samples [17,21,24,25]. Association between DNA methylation heterogeneity and tumour clonal progression has been demonstrated in hepatocellular carcinoma and lower-grade glioblastoma [26,27]. However, deciphering the influences of different immune cell types mixed in tumour immune microenvironment on DNA methylation heterogeneity remains inadequate.

Our study revealed a bimodal distribution of DNA methylation profiles in single purified human immune cells, while we observed that intermediate methylation levels were elevated in the tumour tissue. We utilized the Proportion of sites with Intermediate Methylation (PIM) to quantify DNA methylation heterogeneity, which showed increased levels with immune cell infiltration. Furthermore, our analysis of GBM and LGG samples from The Cancer Genome Atlas (TCGA) indicated that lower PIM scores were associated with better survival outcomes and slower tumour progression in glioma patients. We identified cell-type-associated heterogeneous CpG site (CpG^ct) for B cell, CD4⁺ T cell, CD8⁺ T cell, granulocyte, monocyte, and neutrophil, and constructed a Cell-type-associated DNA Methylation Heterogeneity Contribution (CMHC) score to dissect the effect of different immune cell types on the DNA methylation level of a given CpG^ct. In GBM, LGG samples from TCGA, and other validation datasets [28,29], we selected the prognosis-related CpG^ct to construct the risk score, i.e., the Cell-type-associated DNA Methylation Heterogeneity Risk (CMHR) score, which was independent of the age, gender, tumour grade, O-6-methylguanine-DNA methyltransferase (MGMT) promoter status, and isocitrate dehydrogenase (IDH) status in glioma patients. CMHR predicted the IDH status and glioma histological phenotype better than T-lymphocyte cytotoxic score. Furthermore, we found that DNA methylation level alterations of prognosis-related CpG^ct may be associated with the drugs Temozolomide, Bevacizumab, and radiation therapy in glioma patients.

The results of this study focus on the link between DNA methylation heterogeneity and the intricate tumour immune microenvironment, suggesting that DNA methylation heterogeneity could prove to be a valuable tool for monitoring glioma progression, predicting treatment response, stratifying subtypes, and conducting phenotype studies. These findings have potential implications for the diagnosis and treatment of glioma patients and may lead to more personalized treatment options in the future.

Materials and methods

Data acquisition and preprocessing

We obtained The Cancer Genome Atlas (TCGA) level 3 gene expression data, standard β-value DNA methylation profiles (based on Illumina Infinium HumanMethylation450 BeadChip), and phenotype-survival data of glioblastoma (GBM) and lower-grade glioma (LGG) from the University of California Santa Cruz (UCSC) Xena browser (Genomic Data Commons [GDC] hub, available at https://gdc.xenahubs.net). Mutation annotation format (MAF) files of GBM and LGG were obtained via the TCGA data portal (https://www.cancer.gov/tcga/). Other validation cohorts of glioma patients were downloaded from Gene Expression Omnibus (GEO, accession numbers: GSE103659; GSE60274). We removed the genes that were not expressed in more than 70% of the samples in the gene expression data. The standard β methylation profiles were processed in the same way. The missing values of standard β profiles were supplemented using the ‘knnInputation’ function in the R package DMwR.

We download the processed β-values of purified human immune cells from GEO with identifiers: GSE103541, GSE118144, and GSE166844 (based on Illumina Infinium HumanMethylationEPIC Beadchip). These three EPIC datasets include six immune cell types, i.e., B cells, CD4⁺ T cells, CD8⁺ T cells, granulocytes, monocytes, and neutrophils. The six immune cell types are involved in immune-inflammatory response and anti-tumour immunity [30,31]. We used whole genome bisulphite sequencing (WGBS) and reduced representative bisulphite sequencing (RRBS) from GEO to show the distribution of β-value for single cell type in Figure 1c. In addition, we download three methylation arrays (GSE207426, GSE85087, and GSE195684) to explore the association between DNA methylation alterations of prognosis-related CpG^ct and drug therapies for glioma patients. We supplemented the number of samples and data types for each dataset mentioned above in Table S1.

Figure 1. The association between the PIM score and the tumor immune microenvironment. (a). The bubble plot shows different correlations in GBM (upper) and LGG (bottom) samples from TCGA. The immune score, stromal score, and tumor purity were assessed by the ESTIMATE algorithm [32]. (b). The charts show the correlation degree (see methods) between the PIM score corresponding to different intermediate methylation levels and cell marker expressions in various cell types (table S2). A higher correlation degree indicates a stronger association between marker gene expressions and intermediate methylation levels. The range of correlation degrees at different intermediate methylation levels is larger, and the color bar is deeper. (c). The tumor immune microenvironment consists of multiple cell types (left, created with biorender.com, https://biorender.com/). The distributions of DNA methylation profiles varied in single cell type and tumor bulk samples (right, single cell type shows 0% or 100% methylation status, intermediate methylation level was higher in tumor tissue). The lollipop plot shows glioma tumor tissue had the highest intermediate methylation level.

Classification of tumor immune microenvironmental subtypes

Firstly, we calculated the single-sample gene set enrichment score (ssGSEA) in 34 cell types for GBM samples (n = 168) and LGG tumours (n = 529) from TCGA by the R package gsva (gsva function with method = ‘ssgsea’) [33]. Of these, 33 cell types were chosen from Xcell (https://xcell.ucsf.edu/) [34], and we added ‘microglia’ after retrieving the literature [35]. The marker genes of 34 cell types were listed in Table S2. ssGSEA scores for the total number of 697 TCGA samples were shown in Table S3.

Secondly, we used the ssGSEA score matrix of GBM and LGG samples from TCGA, performing Non-negative Matrix Factorization (NMF) clustering using the R package NMF [36]. The parameter k was set ranged from 3 to 8, and we repeated the NMF analysis 50 times. The best k = 3 was selected when the coupling coefficient was maximum. Therefore, we grouped GBM and LGG patients from TCGA into three tumour immune microenvironmental subtypes (NMF subtypes, NMF-1, NMF-2, and NMF-3).

Evaluation of DNA methylation heterogeneity

The Proportion of sites with Intermediate Methylation (PIM) calculates the number of CpG sites with intermediate methylation, which has been used to measure DNA methylation heterogeneity [1,23]. We defined a correlation degree (see follow methods) to select the ranges of intermediate β-values, and the results were shown in Figures 1a,b. We then used the intermediate methylation levels ranges from 0.2 to 0.6. In addition, the β-value ranges from 0.2 to 0.6 have been used in Wenger et al. study to measure DNA methylation heterogeneity [1]. We thus calculated the proportion of CpG sites with β-values from 0.2 to 0.6 ($\mathit{Cp}{G}^{\mathit{inter}}$) across genome-wide as follows,

$\mathit{PIM}=\frac{\mathit{num}\left(\mathit{Cp}{\mathit{G}}^{\mathit{inter}}\right)}{\mathit{N}},$

where N represents the total number of genome-wide CpG sites for each patient. A higher PIM score means more intermediate methylation sites were enriched in tumour tissue, and reflects stronger DNA methylation heterogeneity. We used the PIM score to compare intertumoral DNA methylation heterogeneity and dissect the association with the clinical presentation of glioma patients.

Cell-type-associated DNA methylation heterogeneity contribution score

We constructed a Cell-type-associated DNA Methylation Heterogeneity Contribution (CHMC) score to dissect the effect of different immune cell types on the β-values of the cell-type-associated heterogeneous CpG sites (CpG^ct). The steps were displayed as follows:

Step 1. Construction of DNA methylation reference matrix

We used EPIC datasets of purified human immune cells from GEO (Table S1), and there were six immune cell types including B cell, CD4⁺ T cell, CD8⁺ T cell, granulocyte, monocyte, and neutrophil due to their immunity functions in tumours [30,31]. We performed the Wilcoxon signed-rank test (H0: the β-value of $\mathit{Cp}{G}_i$ in cell-type A is equal to that in cell-type B, where i represents the i^th CpG site), and the P-values were corrected by the Benjamini-Hochberg multiple tests (BH). β-values for each CpG site in a cell type were compared pairwise to every other cell type (there are 15 paired cell combinations for 6 immune cell types). CpG sites with adjusted p-values less than 0.05 were retained as cell-type-associated CpG sites (*).

To investigate the association between DNA methylation heterogeneity and the tumour immune microenvironment, we further selected the cell-type-associated heterogeneous CpG sites (CpG^ct) in glioma tumour tissues. Based on the cell-type-associated CpG sites (*) identified above, we filtered the CpG^ct by the following steps: (i) CpG sites (*) with mean β-values between 0.2 and 0.6 in tumour tissues; (ii) CpG sites (*) with the coefficients of variation (CV, calculated by the ratio between standard deviation and mean) [23] of β-values in the top 75%; (iii) CpG sites (*) with the proportion of intermediate methylated tumour samples ranked in the top 75%.

We then created the methylation reference matrix (M). The rows of M are K (1, 2, … k) CpG^ct, and the columns of M are six immune cell types (C = 6, each cell type has a DNA methylation profile $m_c$). M could be expressed as (1) or (2), where $m_{\mathit{ij}}$ denotes the mean β-value of the i^th CpG^ct within j^th cell type in EPIC datasets of purified human immune cells (Table S1), which was expressed as (3). Where T means the number of samples of the same cell type in EPIC datasets, ${\beta }_{\mathit{it}}$ represents the β-value of i^th CpG^ct in the t^th cell sample,

(1) $\mathit{M}=\left(m_B,m_{\mathit{CD}4},m_{\mathit{CD}8},m_{\mathit{granulocyte}},m_{\mathit{monocyte}},m_{\mathit{neutrophil}}\right)$(1)

(2) $M_{\left[\mathit{k},6\right]}=\left(\begin{array}{cccccc}{m}_{11},& m_{12},& m_{13},& m_{14},& m_{15},& m_{16}\\ m_{21},& m_{22},& m_{23},& m_{24},& m_{25},& m_{26}\\ \dots ,& \dots ,& \dots ,& \dots ,& \dots ,& \dots \\ m_{\mathit{k}1},& m_{\mathit{k}2},& m_{\mathit{k}3},& m_{\mathit{k}4},& m_{\mathit{k}5},& m_{\mathit{k}6}\end{array}\right)$(2)

(3) $m_{\mathit{ij}}=\frac{{\sum }_{\mathit{t}=1}^{\mathit{T}}{\beta }_{\mathit{it}}}{\mathit{T}}$(3)

Step 2. Estimation of tumor-infiltrative immune cells based on deconvolution

Nextly, we utilized matrix M to deconvolute immune cell infiltration for glioma patients. For the given number of cell types (C), each cell type has a DNA methylation profile ($m_c$). We used the robust partial correlation (RPC) algorithm to assess the intra-sample heterogeneity by the R package EpiDISH, where ${\omega }_c$ denotes the weight coefficient, and the model is based on the least-squares sense to estimate ${\omega }_c$ [37–39]. The deconvolution can be explained as (5),

(4) ${\mathit{y}}_{\left[\mathit{k},\mathit{\hspace{0.17em}}1\right]}=\left(\begin{array}{c}{\beta }_{11}\\ {\beta }_{21}\\ \dots \\ {\beta }_{\mathit{k}1}\end{array}\right)$(4)

(5) $\mathit{y}={\sum }_{\mathit{c}=1}^{\mathit{C}}{\omega }_c{m}_c+\mathit{\epsilon }$(5)

$\mathit{y}$ represents β-values of K CpG^ct in the given patient, expressed as (4). By performing the EpiDISH procedure, we derived the enrichment score matrix ($\mathit{W}$) of n patients explained as (6) or (7), where ${\omega }_{\mathit{ij}}$means infiltration levels in the j^th cell type of the i^th patient,

(6) $\mathit{W}=\left({\omega }_B,{\omega }_{\mathit{CD}4},{\omega }_{\mathit{CD}8},{\omega }_{\mathit{granulocyte}},{\omega }_{\mathit{monocyte}},{\omega }_{\mathit{neutrophil}}\right)$(6)

(7) $W_{\left[\mathit{n},6\right]}=\left(\begin{array}{cccccc}{\omega }_{11},& {\omega }_{12},& {\omega }_{13},& {\omega }_{14},& {\omega }_{15},& {\omega }_{16}\\ {\omega }_{21},& {\omega }_{22},& {\omega }_{23},& {\omega }_{24},& {\omega }_{25},& {\omega }_{26}\\ \dots ,& \dots ,& \dots ,& \dots ,& \dots ,& \dots \\ {\omega }_{\mathit{n}1},& {\omega }_{\mathit{n}2},& {\omega }_{\mathit{n}3},& {\omega }_{\mathit{n}4},& {\omega }_{\mathit{n}5},& {\omega }_{\mathit{n}6}\end{array}\right)$(7)

Step 3. Definition of CMHC score by performing multivariate linear regression

The methylation levels of CpG sites in tumour tissues are the average signal of CpG sites among cells [23]. We thus fitted a multivariate linear regression model to estimate the effect of different immune cell types on the β-values of CpG^ct. The equals were listed as (9),

(8) $Y_{\left[\mathit{n},1\right]}={\left(\begin{array}{cccc}{\mathit{\beta }}_{11},& {\beta }_{21},& \dots ,& {\beta }_{\mathit{n}1}\end{array}\right)}^{\mathit{T}}$(8)

${\beta }_{11}={\alpha }_0+{\omega }_{11}\cdot {\alpha }_1+{\omega }_{12}\cdot {\alpha }_2+{\omega }_{13}\cdot {\alpha }_3+{\omega }_{14}\cdot {\alpha }_4+{\omega }_{15}\cdot {\alpha }_5+{\omega }_{16}\cdot {\alpha }_6+\mathit{L}\cdot {\alpha }_l+\mathit{e}$

(9) ${\beta }_{21}={\alpha }_0+{\omega }_{21}\cdot {\alpha }_1+{\omega }_{22}\cdot {\alpha }_2+{\omega }_{23}\cdot {\alpha }_3+{\omega }_{24}\cdot {\alpha }_4+{\omega }_{25}\cdot {\alpha }_5+{\omega }_{26}\cdot {\alpha }_6+\mathit{L}\cdot {\alpha }_l+\mathit{e}$(9)

${\beta }_{\mathit{n}1}={\alpha }_0+{\omega }_{\mathit{n}1}\cdot {\alpha }_1+{\omega }_{\mathit{n}2}\cdot {\alpha }_2+{\omega }_{\mathit{n}3}\cdot {\alpha }_3+{\omega }_{\mathit{n}4}\cdot {\alpha }_4+{\omega }_{\mathit{n}5}\cdot {\alpha }_5+{\omega }_{\mathit{n}6}\cdot {\alpha }_6+\mathit{L}\cdot {\alpha }_l+\mathit{e}$

where $\mathit{Y}$ means the β-values of the given CpG^ct in n patients, explained as (8). L represents the latent factors estimated by the R package sva [40]. The coefficients ${\alpha }_1$to ${\alpha }_6$ (9) were considered as Cell-type-associated DNA Methylation Heterogeneity Contribution (CMHC) score. Hence, we derived the CMHC scores of K CpC^ct, expressed as (10) or (11). Where ${\alpha }_{\mathit{ij}}$ denotes the CMHC score of i^th CpG^ct corresponding to the j^th cell type, reflecting the influencing degree of cell type j on the β-values of i^th CpG^ct. The steps mentioned in this method were shown in Figure 4a.

(10) $\mathit{CMHC}=\left({\alpha }_B,{\alpha }_{\mathit{CD}4},{\alpha }_{\mathit{CD}8},{\alpha }_{\mathit{granulocyte}},{\alpha }_{\mathit{monocyte}},{\alpha }_{\mathit{neutrophil}}\right)$(10)

(11) $\mathbf{CMH}{\mathbf{C}}_{\left[\mathit{k},6\right]}=\left(\begin{array}{cccccc}{\alpha }_{11},& {\alpha }_{12},& {\alpha }_{13},& {\alpha }_{14},& {\alpha }_{15},& {\alpha }_{16}\\ {\alpha }_{21},& {\alpha }_{22},& {\alpha }_{23},& {\alpha }_{24},& {\alpha }_{25},& {\alpha }_{26}\\ \dots ,& \dots ,& \dots ,& \dots ,& \dots ,& \dots \\ {\alpha }_{\mathit{k}1},& {\alpha }_{\mathit{k}2},& {\alpha }_{\mathit{k}3},& {\alpha }_{\mathit{k}4},& {\alpha }_{\mathit{k}5},& {\alpha }_{\mathit{k}6}\end{array}\right)$(11)

Survival analysis

To examine the clinical application of CpG^ct in glioma, we used the matrix W and CMHC to calculate the CMHC_cell for each CpG^ct in individual glioma patients. The CMHC_cell score was displayed in (13), (14), and Figure 5a.

(12) $\mathit{CMH}{C}_{\mathrm{i}}=\left({\alpha }_{\mathit{i}1},{\alpha }_{\mathit{i}2},{\alpha }_{\mathit{i}3},{\alpha }_{\mathit{i}4},{\alpha }_{\mathit{i}5},{\alpha }_{\mathit{i}6}\right)$(12)

$W_{\left[\mathit{n},6\right]}=\left(\begin{array}{cccccc}{\omega }_{11},& {\omega }_{12},& {\omega }_{13},& {\omega }_{14},& {\omega }_{15},& {\omega }_{16}\\ {\omega }_{21},& {\omega }_{22},& {\omega }_{23},& {\omega }_{24},& {\omega }_{25},& {\omega }_{26}\\ \dots ,& \dots ,& \dots ,& \dots ,& \dots ,& \dots \\ {\omega }_{\mathit{n}1},& {\omega }_{\mathit{n}2},& {\omega }_{\mathit{n}3},& {\omega }_{\mathit{n}4},& {\omega }_{\mathit{n}5},& {\omega }_{\mathit{n}6}\end{array}\right)$

(13) $\mathit{CMH}{C}_{\mathit{cell}}=\left(\mathit{CMH}{C}_B,\mathit{CMH}{C}_{\mathit{CD}4},\mathit{CMH}{C}_{\mathit{CD}8},\mathit{CMH}{C}_{\mathit{gran}},\mathit{CMH}{C}_{\mathit{mono}},\mathit{CMH}{C}_{\mathit{neu}}\right)$(13)

(14) $\mathit{CMH}{C}_{\mathit{cell}}=\left(\mathit{CMH}{C}_B,\mathit{CMH}{C}_{\mathit{CD}4},\mathit{CMH}{C}_{\mathit{CD}8},\mathit{CMH}{C}_{\mathit{gran}},\mathit{CMH}{C}_{\mathit{mono}},\mathit{CMH}{C}_{\mathit{neu}}\right)$(14)

where $\mathit{CMH}{C}_{\mathrm{i}}$ represents the CMHC score of i^th CpG^ct (i from 1 to k), expressed as (12), $\mathit{W}$ denotes infiltration levels of six immune cell types in n glioma patients. We thus derived the CMHC_B, CMHC_CD4, CMHC_CD8, CMHC_gran, CMHC_mono, and CMHC_neu for each patient at the given CpG^ct, explained as (13) or (14). CMHC_cell score was used to measure the DNA methylation degree of a given CpG^ct in one specific cell type for different patients.

Nextly, we used GBM and LGG samples from TCGA as training datasets, performing the multivariate Cox regression model based on six CMHC_cell scores of each CpG^ct (the formula for each CpG^ct is: Surv(OS.time, OS)~ CMHC_B + CMHC_CD4 + CMHC_CD8 + CMHC_gran + CMHC_mono +CMHC_neu). The CpG^ct with P-values <0.05 were considered as candidate prognosis-relevant markers. We subsequently performed multivariate stepwise Cox regression on potential prognosis-related CpG^ct. We selected a model with the minimum Akaike Information Criterion value (AIC), producing eight CpG^ct as the prognosis signature. We constructed a risk score to examine the prognostic value of the eight CpG^ct, i.e., Cell-type-associated DNA Methylation Heterogeneity Risk (CMHR) score, the CMHR was calculated as (15),

(15) $\mathrm{CMHR}\mathrm{score}{\sum }_{\mathit{i}=1}^8\mathrm{CMH}{\mathrm{C}}_{\mathit{cell}}\mathit{\hspace{0.17em}}\left(\mathit{Cp}{G}_i^{\mathit{ct}}\right)\ast \mathrm{coef}\left(\mathit{Cp}{G}_i^{\mathit{ct}}\right)$(15)

where coef ($\mathit{Cp}{G}_i^{\mathit{ct}}$) is the regression coefficient of i^th prognosis-relevant CpG^ct. We further divided GBM and LGG patients from TCGA into high-risk and low-risk groups according to the median CMHR score, investigating the relationship between overall survival probability and the CMHR score by Kaplan-Meier (KM) analysis with the log-rank test (P-values <0.05), and tested in additional independent datasets GSE60274 [28] and GSE103659 [29]. Survival analysis was based on the R packages survival and survminer.

Definition of correlation degree

To assess the advantage of the CpG^ct for inferring the cellular distributions in the tumour immune microenvironment, we defined a correlation degree to quantify the association between marker gene expressions (Table S2) and cellular infiltration. The cellular infiltration was deconvolution by EpiDISH and methylCIBERSORT [41]. DNA methylation reference matrix M (containing CpG^ct) was used in EpiDISH. The methylation reference matrix of methyl-CIBERSORT (heterogeneity of CpG sites was not considered) was obtained from Chakravarthy et al. study [41]. The correlation degree was expressed as (16),

(16) $\mathit{correlation}\mathit{degree}=-log10{\left(\mathit{q}-\mathit{values}\right)}^{\ast }\left|\mathit{Rho}\right|$(16)

where q-values are P-values of the Spearman bivariate correlation analysis adjusted by BH, and |Rho| is the absolute value of Spearman’s rank correlation coefficients. The greater the correlation degree, the stronger the association between marker gene expressions and cellular infiltration.

Quantification of inter-patient differences and inter-subtype differences

To examine that patients within the same NMF subtype have similar epigenetic backgrounds and cellular compositions, we calculated the CV value [23], entropy [42], and Gini-index [43] of DNA methylation level, ssGSEA score, and the PIM score for TCGA glioma patients and NMF subtypes to measure differences between patients and between NMF subtypes, respectively. The larger the CV value, the greater the variation [23]. The formulas of entropy (17) and Gini-index (18) are following,

(17) $\mathit{Entropy}:\mathit{H}\left(\mathit{x}\right)=-\sum _{\mathit{i}=1}^n\mathit{p}\left(x_i\right)log\left(x_i\right)$(17)

(18) $\mathit{Gini}-\mathit{index}:\mathit{Gini}\left(\mathit{D}\right)=1-{\displaystyle \sum _{\left|\mathit{n}\right|}^{\mathit{k}=1}{\left(p_k\right)}^2}$(18)

where parameter $\mathit{p}\left(\mathit{x}\right)$ denotes the probability of occurrence of a given event $\mathit{x}$. We divided methylation values, ssGSEA scores, and PIM scores into different classes according to their respective quartiles so that we can derive $\mathit{p}\left(\mathit{x}\right)$ for each level.

Somatic mutation analysis

We obtained the MAF files of GBM and LGG samples from TCGA. The somatic alterations were detected for the patients in high-risk (n = 261) and low-risk (n = 268) by the R package maftools [44]. The top 25 driver genes with the highest mutation frequency were further analysed by Fisher’s exact test. Genes with P-values <0.001 were regarded as differentially mutated genes.

Identification of featured genes that related to cytotoxic T-lymphocyte infiltration and the survival of glioma patients

T-lymphocyte cytotoxicity was revealed to be associated with the survival of GBM patients [45]. We thus explored the expression differences of cytotoxicity-associated genes between high-risk and low-risk groups (TCGA samples). We collected gene sets related to T-lymphocyte infiltration and anti-tumour immunity from The Molecular Signatures Database (MsigDB), Kyoto Encyclopedia of Genes and Genomes (KEGG), and Gene Ontology (GO). The function of these genes includes (i) antigen presentation such as HLA-A [41]; (ii) as chemokines to participate in T-lymphocyte recruitment such as CXCL9 [46,47]; (iii) as effector molecules performing tumour killing such as granzyme and perforin [46,48].

We then filtered genes associated with the overall survival probability and tumour progression of TCGA patients (log-rank test, P-values <0.05; Cox-wald test, P-values <0.05). We performed the Spearman bivariate correlation analysis (q-values <0.05) to select genes related to the CHMR score and cytotoxic score (the geometric mean of GZMA and PRF1 expressions [41,46]). These genes were considered featured genes related to cytotoxic T-lymphocyte infiltration and the prognosis of glioma patients.

Results

Intertumoral DNA methylation heterogeneity characterizing glioma subtypes related to tumor immune microenvironment

Tumour masses are composed of multiple cell types, including malignant cells, immune cells, and stromal cells, known as the tumour immune microenvironment [49,50]. Several studies have shown that the intermediate methylation levels correlated with mixed cellular compositions in tumour tissues, reflecting DNA methylation heterogeneity in the tumour immune microenvironment [1,23,51]. We used the Proportion of sites with Intermediate Methylation (PIM) (see methods) to estimate DNA methylation heterogeneity. Higher PIM scores indicate more intermediate methylated CpG sites enriched in the sample. At different ranges of intermediate methylation levels, we calculated different PIM scores for glioma patients. The result showed when intermediate methylation levels fell within 0.2 to 0.6, the PIM score correlated more with immune cell infiltration, stromal cell infiltration, tumour purity, mutation count, and aneuploidy percentage (Figure 1a), and showed greater correlation degrees with cell marker genes (Figure 1b, Table S2), suggesting that DNA intermediate methylation (β-values ranges from 0.2 to 0.6) could be more indicative of the transcriptional diversity of cell marker genes. After comparing the methylation profiles of purified human immune cells with the tumour tissue (GBM sample from TCGA), we found the DNA methylation profiles of B cell, monocyte, macrophage, CD4⁺ T cell, natural killer cell, neutrophil, and two glioma cell lines, including U87 [52] and U251 [53] showed bimodal distributions. By contrast, more intermediate methylation sites were enriched in the tumour tissue (Figure 1c).

In previous studies, glioma subtypes were classified based on bulk gene expressions and DNA methylation levels to examine transcriptional and epigenetic heterogeneity [35,54,55]. Our analysis found that PIM score was positively associated with tumour purity and immune cell infiltration levels (Figure 2a, Fig.S1A), This result indicated that the DNA methylation heterogeneity increased with more complex cell compositions present in the tumour immune microenvironment. We then used Non-negative Matrix Factorization (NMF) algorithm (see methods) to perform clustering analysis based on the ssGSEA scores of TCGA samples (Table S2), the ssGSEA score reflects cellular infiltration levels. The clustering analysis stratified GBM and LGG samples from TCGA as three tumour immune microenvironmental subtypes (NMF subtypes) with different cellular compositions (see methods, Figure 2b). We found that the proportion of dendritic cells and M1 macrophages was higher in NMF-1, the infiltration levels of memory B cells and naïve B cells were higher in NMF-2, and T cell subsets such as Th1 cells and Th2 cells were enriched in NMF-3 (Figure 2b). The PIM score showed the lowest level in NMF-2 (Figure 2c). Additionally, we also discovered that DNA methylation level, tumour purity, and abundances of immune infiltration were lowest in NMF-2 (Fig.S1B, Fig.S1C). These results suggested DNA methylation heterogeneity is associated with NMF subtypes. Furthermore, we found the CV values of ssGSEA scores in NMF-2 were smallest in CD4 T-lymphocyte and CD8 T-lymphocyte (Figure 2d). This suggested that the weaker DNA methylation heterogeneity of NMF-2 might be attributed to the lower diversity of T cell types present in the tumour tissues. We then compared the differences in glioma patients before clustering and patients within the same NMF subtype, the CV values for PIM scores, mean DNA methylation levels, and ssGSEA score all showed larger values between glioma patients than between NMF subtypes (see methods, Figure 2e, Fig.S1D), revealing more similar epigenetic backgrounds and cellular compositions of glioma patients within the same NMF subtype. We discovered that NMF-2 aggregated the most significant amount of LGG, and there was a clear overlap in the IDH mutant subtype and the MGMT promoter methylation subtype (Figure 2f). Furthermore, NMF-2 showed better survival outcomes and slower disease progression than NMF-1 and NMF-3 (Fig.S1E, Fig.S1F). This finding is consistent with previous studies, suggesting that glioma patients with IDH mutations have better survival rates [56–59]. NMF subtype was a factor independent of the age and gender of glioma patients (Fig.S1G). These results indicated that DNA methylation heterogeneity may affect the survival of glioma patients.

Figure 2. Intertumoral DNA methylation heterogeneity characterizing NMF subtypes. (a). The correlation between DNA methylation heterogeneity and tumor purity (Spearman’s rank correlation coefficient). (b). Using non-negative matrix factorization, GBM and LGG samples from TCGA were classified as three tumor immune microenvironmental subtypes (NMF subtypes). (c). The DNA methylation heterogeneity differed between NMF subtypes (Wilcoxon signed-rank test, P-values*** < 0.001, P-values** < 0.01, P-values* < 0.05). (d). The bubble plot demonstrates the CV values of the cellular enrichment ssGSEA score across NMF subtypes. (e). The line chart shows the CV values in DNA methylation heterogeneity, DNA methylation level, and cellular enrichment score in all glioma patients before clustering and in patients within the same NMF subtypes. (f). 90.18% of NMF-2 samples were enriched in isocitrate dehydrogenase (IDH) mutant status and 89.47% of NMF-2 samples were enriched in O-6-methylguanine-DNA methyltransferase (MGMT) promoter methylation status (hypergeometric test-P < 0.001).

DNA methylation heterogeneity distinguishes glioma tumor progression and drug treatment response

We used GBM and LGG samples from TCGA to explore the association between clinical variables (Table S4) and DNA methylation heterogeneity. To examine the differences in DNA methylation heterogeneity intensity of glioma patients among NMF subtypes, we investigated the distribution of PIM scores, and we found that there were differences in the number of patients with different PIM scores across NMF subtypes (Figure 3a). The largest number of patients with low PIM scores were enriched in NMF-2 (Figure 3b). Based on Cox regression analysis, we found that the PIM score was a factor independent of the age, gender, tumour grade, MGMT promoter status, and IDH status in glioma patients (Figure 3c, Cox-wald test, P-value <0.05), which is consistent with previous studies that DNA methylation heterogeneity was related to the glioma prognosis [17]. The GBM and LGG samples of TCGA have different subtype classifications (Table S4), in contrast to GBM samples, the number of LGG samples overlaps more with other subtypes including IDH subtypes, MGMT subtypes, telomerase reverse transcriptase (TERT) subtypes, and histological subtypes (Figure 3d). We then compared the PIM differences across distinct glioma phenotypes, GBM with higher PIM scores than LGG. There were higher PIM scores in IDH wild-type (WT), MGMT promoter unmethylated, and TERT promoter mutated tumours. In addition, recurrent gliomas had higher PIM scores than primary tumours, and the tumours with higher grades (grade 4) had the highest PIM score than lower grades (Figure 3e). Taken together, these results suggested that DNA methylation heterogeneity is more pronounced in undifferentiated tumours with malignant behaviour and increases with higher tumour grade [60].

Figure 3. DNA methylation heterogeneity is related to clinical variables of GBM and LGG patients from TCGA. (a). The distribution of PIM scores differed among NMF subtypes (Kruskal-Wallis test, P-value*** < 0.001). (b). The line chart shows the number distribution of glioma patients with different PIM scores among NMF subtypes. (c). The forest plot shows the prognostic value of the PIM score based on cox regression analysis. (d). The numbers of GBM and LGG samples from TCGA correspond to IDH status, MGMT promoter status, TERT promoter status, histological subtype, original subtype, G-CIMP phenotype, tumor grade, and sample type. The box plots show the PIM scores in (e) different glioma phenotypes, tumor grades, (f). Patients with different disease progression after Temozolomide treatment, and (g). Different patients treated with different drug therapies (Wilcoxon signed-rank test, P-values*** < 0.001, P-values** < 0.01, P-values* < 0.05). (h), the bar plots display the distribution of PIM scores for patients who received chemotherapy. For the responder (complete response and partial response), the proportions of high-PIM, low-PIM, and medium-PIM are 23.81%, 33.33%, and 42.86%, for non-responders (progressive disease and stable disease), the proportions of high-PIM, low-PIM, and medium-PIM are 39.25%, 31.78%, and 28.97% (Wilcoxon signed-rank test, P-values*** < 0.001).

We further examined the association between drug therapies and the PIM score. Temozolomide, a major chemotherapy drug used for glioma treatment [61], in patients who received Temozolomide, we found that the PIM scores were higher in progressive disease (PD) patients than in stable disease (SD) patients (Figure 3f). Additionally, the PIM score was related to drug therapies. PIM scores were higher in patients treated with Temozolomide than those receiving Carmustine. In patients accepting targeted molecular drugs, PIM scores were higher in Avastin than in Bevacizumab (Figure 3g). Moreover, PIM scores were linked to patients’ treatment response to chemotherapy, there were differences in PIM scores among subgroups of high-PIM, low-PIM, and medium-PIM patients, respectively. (Wilcoxon signed-rank test, P-values*** < 0.001). The proportion of high PIM scores in responders is minimum and the percentage of high-PIM patients is maximum in non-responders (Figure 3h). These results suggested that DNA methylation heterogeneity might play a role in distinguishing glioma phenotypes, tumour progressions, and patients’ responses to drug therapy.

Cell-type-associated heterogeneous CpG sites correlate with genes regulating T-lymphocyte proliferation and differentiation

To dissect the cell-type-specific DNA methylation patterns in a complicated tumour immune microenvironment, we identified CpG sites related to B cell, CD4⁺ T cell, CD8⁺ T cell, granulocyte, monocyte, and neutrophil. These CpG sites were mainly distributed in the gene body regions (Table S5, Figure 4b, Fig.S2A), after filtering cell-type-associated heterogeneous CpG sites (CpG^ct), we constructed a Cell-type-associated DNA Methylation Heterogeneity Contribution (CMHC) score to assess the influence of six immune cell types on the DNA methylation levels of CpG^ct (see methods, Figure 4a). In addition, we defined a correlation degree (see methods) to evaluate the association between marker gene expressions and cellular infiltrations calculated by the CMHC model (see methods) and methylCIBERSORT [41]. Using the DNA methylation profiles and bulk gene expressions of GBM and LGG samples from TCGA to calculate the correlation degree for different cell types. The correlation degree for B cells, CD4⁺ T cells, CD8⁺ T cells, and neutrophils were higher in the CMHC model (Figure 4c). This suggested the heterogeneity of CpG sites could decipher the cellular distribution of tumour immune microenvironment more accurately.

Figure 4. Cell-type-associated heterogeneous CpG sites correlate with genes regulating the proliferation and differentiation of T-lymphocyte. (a). Identification of cell-type-associated heterogeneous CpG site (CpG^ct) and the construction of CMHC score. (b). The barplot shows the number and genome distribution of cell-type-associated CpG sites (hypermethylation refers to the mean DNA methylation level of the given CpG site in cell-type a is higher than in other cell types, while hypomethylation refers to the mean DNA methylation level of the given CpG site in cell-type a is lower than in other cell types). (c). The box plots show the correlation degrees in different cell types calculated based on CMHC and methylCIBERSORT. DNA methylation reference matrix (M) consists of CpG^ct, which was used in the deconvolution for the CMHC model. (Wilcoxon signed-rank test, P-values*** < 0.001, P-values** < 0.01, P-values* < 0.05). Enrichment of (d). GO terms and (e). KEGG pathways associated with DEGs. GO terms include biological process (BP), molecular function (MF), and cell component (CC).

Subsequently, we identified 178 differentially expressed genes (DEGs) related to CpG^ct between NMF subtypes (Table S6, Fig.S2B). These DEGs are involved in the activation and proliferation of T cells, the regulation of leukocyte-mediated immunity, interferon-gamma production, and cytokine receptor activity (Table S7, Figure 4d). Our analysis showed the KEGG pathways related to these DEGs including T cell receptor signalling pathway, T cell differentiation, cytokine-cytokine receptor interaction, and haematopoietic cell lineage (Table S7, Figure 4e), this result was consistent with previous gene studies, which also reported the immune response pathways mediated by T-cell in glioblastoma [62,63]. These findings suggested that CpG^ct-related genes play roles in the proliferation and differentiation of T-lymphocyte.

Cell-type-associated heterogeneous CpG sites in glioma tissues are related to neurogenesis and nervous system development

We evaluated the prognostic value of the CMHC score by calculating the CMHC_cell score (CMHC_B, CMHC_CD4, CMHC_CD8, CMHC_gran, CMHC_mono, and CMHC_neu) of each CpG^ct in individual patients (see methods, Figure 5a), The GBM and LGG samples from TCGA were used as training data, and GSE103659 and GSE60274 as validation data. Using the β-values of CpG sites, we performed univariate Cox regression and log-rank test and identified 94 Prognosis-Related CpG sites (PRS) (see methods, Figure 5b, Table S8). Parallelly, we did multivariate Cox regression and log-rank test based on the CMHC_cell score and selected 83 Prognosis-Related cell-type-associated heterogeneous CpG sites (PRH) (see methods, Figure 5b, Table S9). The DNA methylation levels of PRS and the CMHC_cell scores of PRH were associated with glioma patients’ overall survival (Figure 5e). We then checked the immune cell type distribution of PRH, there were 49 out of PRH (59.04%) were hypomethylated in immune cells, there were 44 out of PRH (53.01%) associated with the B cell, CD4⁺ T cell, and CD8⁺ T cell. B cell and CD4⁺ T cell-type-associated CpG^ct was distributed in more genomic locations including gene body, untranslated regions, and the first exon region (Figure 5c). We further analysed the genomic distribution of PRH, there were 30 out of PRH (36.14%) enriched in the gene body. However, there were only 24 out of PRS (25.53%) were enriched in cell-type-associated heterogeneous CpG sites (Figure 5d). As described in the method, each immune cell type has a CMHC_cell score, reflecting the effect of immune cell types on DNA methylation level for the given CpG^ct, we further investigated the levels of β-values and the CMHC_cell scores for PRH in different genomic regions. We found that mean β-values and CMHC_cell score of B cell both showed larger levels in the shore regions with lower CpG density, while the CMHC_cell score of granulocyte and monocyte showed lower levels in the CpG island (CGI) with higher CpG density (Figure 5f), suggesting the DNA methylation pattern in specific genomic regions differs across immune cell types.

Figure 5. Cell-type-associated heterogeneous CpG sites were enriched in neural-related biological processes. (a). The construction of CMHC_cell score. (b). Prognosis-related CpG sites and CpG^ct were identified by cox regression and Kaplan-Meier (KM) analysis with the log-rank test (P-values < 0.05). (c). The stack bar plots show different DNA methylation patterns of CpG^ct. (d). The stack bar plots show the distributions of different genome regions for CpG sites and CpG^ct. (e). Kaplan-Meier (KM) analysis shows the prognostic value of the CMHC_cell score in TCGA, GSE103659, and GSE60274. (f). The line charts display the distributions of β-values and the CMHC_cell score for 83 prognosis-related CpG^ct. (g). The negative correlation between gene expressions and DNA methylation levels of CpG^ct located in the gene promoter (Spearman’s rank correlation coefficients). GO enrichment of EWAS for (h). 94 CpG sites and (i). 83 CpG^ct.

We subsequently identified genes whose expression values were negatively correlated with the methylation levels of PRH located in the promoter region. Our findings showed that the expression of protein-coding gene SCAMP4 was negatively associated with the DNA methylation level of B cell-type-associated hypomethylation CpG^ct cg15018704, and the expression of protein-coding gene PALM was negatively linked to the DNA methylation level of granulocyte cell-type-associated hypomethylation CpG^ct cg16019612 (Figure 5g). SCAMP4 was revealed as a prognosis marker for GBM and LGG [64]. We also performed the correlation analysis between DNA methylation of PRS and gene expressions, however, we did not identify genes correlated with PRS methylation. These results suggested that taking into account the heterogeneity of CpG sites might be more useful in identifying prognostic markers for glioma patients.

To further understand the biological function of PRS and PRH, we did the GO enrichment analysis by Epigenome-Wide Association Study (EWAS) [65]. We found that PRS play a role in regulating cellular response, chromatin remodelling, and kinase activity (Figure 5h), while PRH were enriched in neuron development, neurogenesis, and nervous system development (Figure 5i). These results suggest that PRH considers the differences in DNA methylation patterns across cells, facilitating the identification of glioma-related prognostic markers and helping us understand the nervous system function of heterogeneous CpG sites.

CMHR score based on cell-type-associated heterogeneous CpG sites identifying high mutation frequency genes

We performed multivariate step Cox regression and selected a set of eight CpG^ct to construct the risk score, i.e., Cell-type-associated DNA Methylation Heterogeneity Risk (CMHR) score (see methods, Figure 6a). We then divided GBM and LGG patients from TCGA into high-risk and low-risk groups based on the median CMHR scores (Figure 6d). In TCGA, GSE103659, and GSE60274, low-risk patients had worse survival outcomes (Figure 6b). CMHR score was independent of the age, gender, tumour grade, MGMT promoter status, and IDH status of glioma patients (Figure 6c). These results indicated that the CMHR score could be a prognostic biomarker for glioma patients.

Figure 6. The prognostic value of CHMC score and CMHR score. (a). The construction of the CMHR score. (b). Kaplan-Meier (KM) analysis showed that low-risk patients had poor survival outcomes (training data: TCGA, and testing data: GSE103659 and GSE60274). (c). For TCGA samples, the forest plot shows the CMHR score is a prognostic factor independent of the age, gender, tumor grade, MGMT promoter status, and IDH status of glioma patients. For GSE103659 and GSE60274, the CHMR score is independent of the age and MGMT promoter status of glioma patients. (d). Oncoplots show the mutation genes differed in high-risk and low-risk patients, there were nine genes overlapped in both groups including IDH1, TP53, ATRX, TTN, RYR2, PIK3CA, MUC16, CIC, and FUBP1. High-risk patients include 221 LGG samples and 40 GBM samples, while low-risk patients contain 266 LGG samples and 2 GBM samples. (e). The expression differences of mutation genes between high- and low-risk patients. For genes mutated more frequently in the high-risk group, there were 4 out of 15 genes (26.67%) with greater expressions in high-risk patients, and there were 5 out of 15 genes (33.33%) expressed higher in low-risk patients. For genes mutated more frequently in the low-risk group, there were 11 out of 16 genes (68.75%) with higher expressions in low-risk patients, and there were 3 out of 16 genes (18.75%) expressed higher in high-risk patients. The protein expression levels of (f) EGFR and (g) PCLO were checked from the human protein atlas (https://www.Proteinatlas.org/humanproteome/pathology).

Several studies have shown that the poor survival of glioma patients was related to gene somatic mutations [66–68]. We identified the top 25 genes with high mutation frequencies in high- and low-risk patients. The mutation frequencies of IDH1, TP53, and ATRX all ranked in the top three positions in both high- and low-risk patients (Figure 6d). These genes were found to be associated with the survival of glioma patients [69]. Moreover, previous studies have shown that the mutation of proto-oncogene EGFR could be a clinical marker for GBM [70–72]. We then examined the expression differences of mutated genes between high- and low-risk patients. EGFR mutated more frequently (Fig.S3, Fisher’s test) in high-risk patients (better survival) and showed higher expressions in low-risk patients (poor survival). Additionally, the tumour-suppressor gene PTEN [73] had larger expression values in low-risk patients (Figure 6e). There has been revealed that the PTEN-related pathway involves in the regulation of glioma tumorigenesis [74]. We then explored the immunohistochemical staining characteristics of EGFR from The Human Protein Atlas (HPA, https://www.proteinatlas.org). We found that the protein expression of EGFR was higher in high-grade glioma (the high-risk group contains more high-grade glioma patients) than in low-grade glioma (Figure 6f). In addition, PCLO also showed greater protein expression in high-grade gliomas (Figure 6g). A past study revealed that the glioma phenotypes with PCLO mutation had poor prognosis [75]. These findings suggested a potential association between DNA methylation heterogeneity at the epigenetic level and gene mutations at the genomic level, indicating the implication of DNA methylation heterogeneity on the prognosis of glioma.

Enhanced T-lymphocytosis and infiltration contribute to improved prognosis in high-risk patients

CMHR score was created based on the CpG^ct of CD4⁺ T cell and CD8⁺ T cell, and the cytotoxic T-lymphocyte (CTL) infiltration was related to anti-tumour immunity and GBM prognosis [45,46,76]. GZMA and PRF1, as cytotoxic effector genes, and the cytotoxic score could be estimated by the geometric mean of their expression values [41,46]. We further calculated the cytotoxic score for high- and low-risk patients to delineate the correlation between CTL and the CMHR score. We found that high-risk patients had greater cytotoxic scores, indicating stronger cytotoxic activity of T-lymphocyte (Figure 7e). We then filtered 25 featured genes that affected the overall survival of GBM and LGG samples from TCGA (see methods, Table S10, Figure 7a). These featured genes are mainly involved in antitumor-immunity and cytotoxicity of T-lymphocyte (Figure 7b), such as effector molecules like GZMA and PRF1 that perform tumour killing [46,48,76], interferon-stimulated chemokines like CXCL9 and CXCL10 that promote T-lymphocyte recruitment [47], and transcription factors like EOMES that regulate CTL effector function [77,78]. We discovered that these featured genes all showed higher expressions in high-risk patients (Figure 7d), suggesting that the better survival of high-risk patients may be due to the stronger cytotoxic function and infiltration level of T-lymphocyte.

Figure 7. CMHR score was related to cytotoxic T-lymphocyte infiltration. (a). The featured genes associated with T-lymphocyte infiltration were collected from the molecular signatures database (MsigDB), Kyoto Encyclopedia of genes and genomes (KEGG), and gene ontology (GO). By performing the log-rank test and Cox-wald test, these featured genes were related to the overall survival of GBM and LGG samples from TCGA. These featured genes were positively correlated with the CMHR score as well as the cytotoxic score (Spearman’s rank correlation coefficients). (b). Featured genes were involved in anti-tumor immunity and tumor killing (MHC: major histocompatibility complex, TCR: T cell receptors). (c). The association between prognosis-related CpG^ct and featured genes in low-risk (upper) and high-risk (bottom) patients (Spearman’s rank correlation coefficients). (d). The box plots show the expressions of featured genes in high- and low-risk patients (Wilcoxon signed-rank test, P-values*** < 0.001). (e). The associations of the CMHR score with tumor purity, immune cell infiltration, stromal cell infiltration levels, and the percent of aneuploidy (Spearman’s rank correlation coefficients). The box plot shows the cytotoxic score in high- and low-risk patients (Wilcoxon signed-rank test, P-values*** < 0.001). (f). For GBM and LGG samples from TCGA, the receiver operating characteristic curves (ROC) show the CMHR score was more predictive for overall survival, IDH status, MGMT promoter status, and histological phenotypes. The predictive potential was measured by the area under curves (AUC). (g). For GBM and LGG samples from TCGA, the bar plots show the values of CMHR score, cytotoxic score, GZMA expression, and PRF1 expression in different histological phenotypes (Kruskal-Wallis test, P-value*** < 0.001), IDH status, and MGMT promoter status (Wilcoxon signed-rank test, P-values*** < 0.001). For GSE103659 and GSE60274, the CMHR score differed between MGMT promoter methylation and unmethylation status (Wilcoxon signed-rank test, P-values*** < 0.001).

It has been observed that in high-risk patients, there is a stronger positive correlation between cytotoxic effector gene (GZMA, GZMH, and PRF1), immune cell identity genes (CD3D and CD3G), and CTL-related pathway genes (ITGAL) with the DNA methylation level of T-lymphocyte-related CpG^ct (Figure 7c). This indicated that T-lymphocyte-associated CpG^ct might potentially play a role in CTL regulation. Furthermore, the CMHR score was found to be positively correlated with the cytotoxic score, tumour purity, and immune cell infiltration (Figure 7e), indicating that the CMHR score increased in a more complex tumour immune microenvironment.

Moreover, in GBM and LGG datasets from TCGA, the CMHR score showed better predictive performances in overall survival, histological phenotype, MGMT promoter status, and IDH status than the cytotoxic score and the GZMA expression (Figure 7f). We compared the predictive performance of the CMHC_cell score of eight CpG^ct with the β-values of CpG sites in the Kessler et al. study, which were related to the overall survival and MGMT promoter status of GBM patients [79]. Our markers displayed better AUC values for the survival and MGMT promoter status in TCGA, GSE103659, and GSE60274 (Fig.S4). In addition, the CMHR score, the cytotoxic score, GZMA expression, and PRF1 expression showed larger values in glioblastoma, IDH-WT subtypes, and MGMT promoter unmethylated subtypes (Figure 7g). These findings suggested that the CMHR score contributes to the identification of markers for distinguishing different glioma phenotypes.

DNA methylation alterations of prognosis-related CpG^ct correlate with the glioma drug therapy

Previous studies have investigated the correlation between CpG sites and the drug therapies of glioma patients to elucidate the clinical application and regulation functions of CpG biomarkers [80–82]. We aimed at exploring the association between prognosis-related CpG^ct and the treatment approaches for glioma. In TCGA patients treated without Temozolomide or radiotherapy, there were 3 out of 8 CpG^ct (37.5%) did not show DNA methylation changes between long- and short-survival patients. In GSE207426, there were 4 out of 8 CpG^ct (50%) had lower DNA methylation levels in long-survival patients receiving Temozolomide, and there were 3 out of CpG^ct (37.5%) showed increased DNA methylation levels in long-survival patients treated with radiation (Figure 8b), which is similar to previous studies that Temozolomide and radiation treatment could affect the DNA methylation levels of CpG sites [83,84]. DNA methylation levels of B cell-related cg03180426 increased more in long-survival patients who received Temozolomide than those who received radiation, while monocyte-related cg24770798 shows decreased DNA methylation in both long-survival patients treated with Temozolomide and radiation. Granulocyte-related cg16269755 showed decreased DNA methylation levels in long-survival patients accepting Temozolomide and increased DNA methylation levels in those receiving radiation (Figure 8a). These results suggested that prognosis-associated CpG^ct could be potential drug targets of Temozolomide and radiotherapy in glioma patients.

Figure 8. The correlation between DNA methylation level of prognosis-related CpG^ct and the drug therapies. (a). Heatmaps show the DNA methylation changes of prognosis-related CpG^ct in short-survival and long-survival glioma patients (TCGA and GSE207426). (b). The bar plots show the proportion of CpG^ct with different methylation changes between patients receiving Temozolomide and radiation as well as TCGA patients without treatment. (c). Line charts show the DNA methylation changes of prognosis-related CpG^ct for 19 GBM patients with primary and paired recurrent tumors (GSE85087), and 16 patients received different treatments. The pie plots show the proportion of 16 patients with different $\mathrm{\Delta \beta }$ values. $\mathrm{\Delta \beta }>0$, DNA methylation levels of CpG^ct were higher in recurrent tumors; $\mathrm{\Delta }\mathit{\beta }<0$, DNA methylation levels of CpG^ct were lower in recurrent tumors. TMZ represents the drug Temozolomide, and BVZ means the drug Bevacizumab. The more species of Δβ in the pie chart, the higher DNA methylation heterogeneity. (d). Kaplan-Meier (KM) analysis shows the prognosis value of B cell-related CpG^ct and monocyte-related CpG^ct in GBM patients accepting radiation (GSE195684). (e). The forest plot shows that the DNA methylation level of B cell-related CpG^ct was a prognostic factor independent of the age, gender, and DNA methylation levels of other CpG^ct (GSE195684).

We found that the CD4⁺ T cell, CD8⁺ T cell, monocyte, and neutrophil-related CpG^ct showed lower DNA methylation levels in recurrent tumours, while granulocyte and B cell-related CpG^ct showed higher DNA methylation levels in recurrent tumours (Figure 8c). The highest intertumoral DNA methylation heterogeneity was observed for CD4⁺ T cell-related CpG^ct cg04848682, while CD8⁺ T cell-related CpG^ct cg04359840 showed the lowest intertumoral DNA methylation heterogeneity (Figure 8c). This may be due to the Temozolomide treatment, which could change the number and the anti-tumour function of CD8⁺ T cells [85,86]. These results provided a new insight that the CpG^ct with lower intertumoral DNA methylation heterogeneity might serve as stable markers for Temozolomide therapy.

Out of 16 patients who received Temozolomide, 9 out of them (56.25%) had higher DNA methylation levels of B cell-related cg23757461 in their recurrent tumours compared with their primary tumours. According to the study, the protein-coding gene NEDD4L (cg23757461 was located in NEDD4L) affects the drug sensitivity of GBM patients to Temozolomide [82]. Moreover, the DNA methylation levels of B cell-related cg03180426 decreased more in recurrent tumours of patients treated with adjuvant Bevacizumab therapy compared with other CpG^ct (Figure 8c). These findings suggested that the DNA methylation alterations of prognosis-related CpG^ct in different immune cell types potentially distinguish the drug treatments for glioma.

Our research discovered a correlation between lower DNA methylation levels of B cell-related CpG^ct and improved survival rates of GBM patients who underwent radiation. Conversely, we found that lower DNA methylation level of monocyte-related CpG^ct was associated with poor prognosis in GBM patients receiving radiation (Figure 8d). Our multivariate Cox regression analysis revealed the DNA methylation level of B cell-related cg23757461 to be an independent prognostic factor for glioma patients, irrespective of their the age, gender, and other CpG^ct (Figure 8e). These findings indicated that the B cell-associated CpG^ct might be a useful clinical marker for evaluating the survival outcomes of glioma patients.

Discussion

In previous studies, researchers have investigated the relationship between DNA methylation heterogeneity and the phenotype and prognosis of glioma. Wenger et al. conducted a study where they stratified glioblastoma subtypes based on DNA methylation to better understand intertumoral epigenetic heterogeneity [1]. Klughammer et al. also looked at the correlation between the clinical presentation of glioblastoma patients and DNA methylation heterogeneity [17]. Additionally, Johnson et al. explored the influencing factors of DNA methylation heterogeneity in glioma tissues [51]. However, the regulatory mechanism of DNA methylation heterogeneity in the tumour immune microenvironment remains unclear. To address this, we delved into the different DNA methylation patterns across cells in the complex tumour immune microenvironment that gives rise to DNA methylation heterogeneity. We hope that our findings can contribute to the application of DNA methylation heterogeneity in subtype stratification and clinical care for glioma patients.

Our research revealed that there was a higher intermediate methylation level in tumour tissues compared to a single purified human immune cell. This finding suggests that there is DNA methylation heterogeneity among cells, which increases with higher abundances of CD4⁺ T cell and CD8⁺ T cell infiltrations in glioma patients. Additionally, we found that patients with fewer intermediate methylation sites had a better prognosis and slower disease progression. This was consistent with a previous study that DNA methylation heterogeneity affects the survival of glioma patients [17]. The PIM score was a prognostic factor independent of various patient factors such as age, gender, and tumour grade. We created the risk score CMHR to evaluate the prognostic value of CpG^ct. The low-risk patients had poorer survival outcomes and we identified differential mutation genes between high-risk and low-risk groups, some genes such as TP53 and ATRX were suggested to be prognostic biomarkers for glioma patients [69,71]. The CMHR score was positively correlated with T-lymphocyte infiltration, which has been indicated to be associated with GBM subtypes classified based on DNA methylation [87]. We found that the CMHR score was more predictive than the cytotoxic score for IDH status and histological phenotypes and showed higher levels in glioblastoma phenotype, IDH-WT subtypes, and MGMT promoter unmethylated subtypes. Our results suggested that the CMHR score may contribute to characterizing the glioma phenotypes and could be a potential biomarker for distinguishing the patients’ drug approaches. Previous studies have also investigated the association between the DNA methylation level of CpG sites and the drug treatment of glioma patients [80–82], adding to the potential clinical significance of our findings.

However, the study had some limitations. The prognostic model was created based on finite cohorts. So future studies need more datasets to test the prognostic value of DNA methylation heterogeneity score. Additionally, more datasets are needed to elucidate the association between prognosis-related CpG^ct and drug therapies for glioma patients. Finally, the potential function of prognosis-related CpG^ct needs to be further explored.

Conclusions

This study has shown that DNA methylation heterogeneity is linked to the tumour immune microenvironment, tumour progression, and treatment responses of glioma patients. The CMHR score is associated with the cytotoxicity and infiltration of T-lymphocyte, which could play a role in characterizing glioma phenotypes.

Authors’ contributions

Y.S.M. and X.P. designed the study, performed analysis, and wrote the manuscript. J.G., X.X.G., H.J.H., Y.H.H., and C.Y.W. revised the manuscript. H.Z. and W.S.N. supervised research and provided critical advice on the study. All authors approved the final manuscript.

Data availability statement

The datasets analysed during the current study are available in the TCGA database; https://www.cancer.gov/tcga/, GEO repository; https://www.ncbi.nlm.nih.gov/geo/, and Xena browser: https://gdc.xenahubs.net. Accession numbers in GEO have been described in the section ‘Materials and Methods’ and Table S1. And custom analysis scripts are available at https://github.com/ShuangyueMa/DMH_glioma.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Supplementary material

Supplemental data for this article can be accessed online at https://doi.org/10.1080/15592294.2024.2318506

Additional information

Funding

This work was supported by the National Natural Science Foundation of China [32170674].