DNA methylation regulator-mediated modification patterns and tumor microenvironment characterization in glioma

Landscape of genetic variations in the DNA methylation regulatory genes in glioma

We selected 20 DNA methylation regulators, containing 3 writers, 3 erasers, and 14 readers after literature survey [22, 23]. The DNA methylation regulatory proteins that regulate the dynamic process of DNA methylation and demethylation and the underlying molecular mechanisms are shown in Figure 1A. The transcriptome and clinicopathological data of the 2228 patients from the TCGA, CGGA1, CGGA2, GSE16011, and GSE108474 glioma datasets, is summarized in Supplementary Table 1.

We then analyzed tumor somatic mutations and copy number variations (CNV) in these DNA methylation regulatory genes in the glioma samples. Somatic mutations in the DNA methylation regulatory genes were observed in only 32 patients; tumor mutation burden (TMB) rate of the 20 regulators was < 1% in the glioma patients (Figure 1B). We observed co-occurrence of mutations in some DNA methylation regulatory genes (MBD4 and TET2; MBD4 and DNMT1; TET3 and DNMT3A; ZBTB33 and DNMT1), mutation-exclusive phenomenon was not observed in any of the 20 DNA methylation regulators (Supplementary Figure 1B). The somatic mutation frequency of glioma patients among DNA methylation regulators was low (<1%), which revealed the different functions of these regulators were not caused by genetic alterations. CNV analysis in the glioma samples demonstrated high frequency of amplifications in MBD3, TDG, UHRF1, DNMT1, and MECP2, and copy number deletions in UHRF2, SMUG1, ZBTB38, MBD4, MBD2, TET2, MBD1, and ZBTB4 (Figure 1C). Figure 1E shows the chromosomal locations of the CNV alternations in these DNA methylation regulatory genes. PPI network analysis showed co-expression of DNA methylation regulatory genes belonging to the same functional group as well as association between writers, erasers, and readers (Figure 1D, Supplementary Table 2). Univariate Cox regression analysis showed that DNMT1, DNMT3A, DNMT3B, MBD2, MBD4, UHRF2, MECP2, UNG, and SMUG1 were potential risk genes, whereas TET1, TET2, TET3, ZBTB4, and MECP2 were potential protective genes in glioma (Supplementary Figure 1C). The expression of 11 DNA methylation regulators (DNMT1, DNMT3A, DNMT3B, MBD2, MBD3, MBD4, UHRF1, UHRF2, UNG, TDG, and SMUG1) was significantly higher in the glioblastoma (GBM) compared to the low-grade glioma (LGG), whereas the expression of 7 DNA methylation regulators (TET1, TET2, TET3, ZBTB33, ZBTB38, ZBTB4, and MECPE) was significantly lower in the GBM (Figure 1F). We performed the immunohistochemistry assay to validate the bioinformatics results, and the assay revealed that three selected DNA methylation regulators (DNMT3A, TET3, and UNG) were expressed to different extents in normal brain tissue (NBT) and LGG (Supplementary Figure 2A–2C), in accord with the bioinformatics results.

Our results indicated most DNA methylation regulators (DNMT1, MBD3, UHRF1, and TDG) with amplificated CNV demonstrated markedly higher expression in GBM, which revealed that the alterations of CNV could be an important element resulting in perturbations on the DNA methylation regulators expression in glioma. This suggested crucial roles for these 20 DNA methylation regulatory genes in glioma progression.

The relationship between DNA methylation regulator expression, prognosis, and tumor immune characteristics in gliomas

Next, we investigated the role of these DNA methylation regulators in the TME. Spearman’s correlation analysis showed that the expression of UNG, ZBTB33, MECP2, and DNMT3A genes was positively associated with the proportion of several immune cell types in the glioma (Supplementary Figure 3A). Because UNG expression was positively associated with naïve B cells, we systematically investigated the biological functions of UNG in the TME. Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis showed that the high UNG expression group was enriched in several tumor-related pathways, such as p53 signaling pathway, TNF signaling pathway, DNA replication, focal adhesion, NF-kB signaling pathway, and Human T-cell leukemia virus 1 infection (Supplementary Figure 3B). Moreover, gliomas with high UNG expression showed significantly higher infiltration of immune cell types such as naïve B cells, delta gamma T cells, M0 macrophages, activated DCs, and activated mast cells compared to gliomas with low UNG expression (Supplementary Figure 3C). Furthermore, correlation between expression of UNG and clinical response to immunotherapy was evaluated in two immunotherapy cohorts, IMvigor210 and GSE78220. Glioma patients with high UNG high expression in the two immunotherapy cohorts showed better prognosis and clinical response (Supplementary Figure 3D–3G).

These results demonstrated association between the expression of DNA methylation regulatory genes and the infiltration of immune cell types into the TME. Our analysis also suggested that UNG regulated the tumor immune microenvironment and is a potential marker for evaluating the efficacy of immunotherapy in glioma patients.

Characterization of three DNA methylation regulator patterns in the gathered glioma datasets

DNA methylation regulators might play a crucial role in the heterogeneity of gliomas and the tumor immune microenvironment status [24]. Therefore, to further determine the functions of the DNA methylation regulator, we used the package of ConsensusClusterPlus to identify distinct glioma patient subgroups based on the specific expression patterns of the 20 DNA methylation regulator in the gathered glioma cohort. We performed unsupervised clustering with a K value of 3 and classified glioma samples into three subgroups or modification patterns, namely, pattern A (n = 1007), pattern B (n = 694), and pattern C (n = 527, Figure 2A, Supplementary Figure 4A–4E). The glioma samples with patterns A and C showed high expression of ZBTB4 and NTHL, respectively, whereas, glioma samples with pattern B showed high expression of most DNA methylation regulators (Figure 2B, Supplementary Figure 4F). PCA showed three distinct expression patterns of all 20 DNA methylation regulators, thereby confirming the unsupervised clustering results (Figure 2C). Kaplan-Meier survival curve analysis showed that OS of glioma patients with pattern B were significantly lower compared to those with patterns A and C, probably because majority of the DNA methylation regulators were highly expressed in the pattern B (P < 0.001; Figure 2D). Then, we investigated the biological behaviors of the gliomas with these three distinct patterns using GSVA. Patterns A and C demonstrated an immune activation phenotype with enrichment of immune pathways such as the chemokine signaling pathway, cytokine-cytokine receptor interaction, and arachidonic acid metabolism; glioma in pattern B were enriched in stromal activation-related pathways, such as p53 signaling pathway, RNA degradation, cell cycle, and DNA replication (Figure 2E, 2F; Supplementary Table 3).

Stromal activation is associated with inhibition of immune cell infiltration and poor prognosis of cancer patients. Therefore, we analyzed the correlation between the three DNA infiltration and several immune cell infiltration and stromal activation signatures. Gliomas with pattern C were significantly rich in immune cell types such as activated DCs, activated NK cells, monocytes, DCs, activated mast cells, co-stimulation T cell, CD8+ T cells, and Tfh (Figure 3A, 3C). These results were in accordance with favorable prognosis of glioma patients with pattern C. Gliomas with patterns B showed significant stromal activation including upregulation of epithelial-mesenchymal transition (EMT)-related genes and pan-fibroblast-THF-β response signature (pan-F-TBRS; Figure 3B). Moreover, pattern B also showed positive correlation with TMB (Figure 3D).

These results showed that all three DNA methylation regulator patterns were associated with significantly distinct immune cell infiltration phenotypes. Pattern A demonstrated immune-excluded tumor phenotype characterized by infiltration of innate immune cell types and stromal activation. Pattern B demonstrated immune-desert tumor phenotype characterized by low immune cell infiltration and immune repression status. Pattern C demonstrated inflamed immune phenotype characterized by high immune cell infiltration into the tumor microenvironment and activation of the adaptive immune response.

Validation analysis of the three DNA methylation regulator patterns in the TCGA glioma dataset

Next, we classified the patients in the TCGA glioma dataset (n = 649) into three clusters based on the expression of the 20 DNA methylation regulators using the package of ConsensusClusterPlus (Supplementary Figure 5A, 5B, 5D). The results for the TCGA dataset were similar to the entire glioma cohort in this study (Supplementary Figure 5C). PCA demonstrated that the 20 DNA methylation regulators formed three distinct patterns (A, B, and C) in the TCGA dataset (Supplementary Figure 5E). Kaplan-Meier survival analysis showed that OS was significantly lower for glioma patients in pattern C compared to those with patterns A and B (P < 0.001; Supplementary Figure 5G), probably due to high expression of most DNA methylation regulator in pattern C (Supplementary Figure 5C). Furthermore, glioma in pattern B showed significantly higher proportions of various immune cells types such as activated DCs, CD8+ T cells, activated mast cells, monocytes, activated NK cells, and plasma cells, which was in accordance with significantly increased OS of glioma patients with pattern B (Supplementary Figure 5F). Moreover, higher stromal activity was observed in pattern C including activation of pan-F-TBRS and EMT-related markers (Supplementary Figure 5H).

This integrated analysis demonstrated that glioma patients could be sorted into immune-inflamed, immune-excluded, and immune-desert phenotypes based on the expression levels of 20 DNA methylation regulators both in the gathered glioma cohort and the TCGA glioma dataset.

Characterization of differences in the clinicopathological parameters and gene expression patterns in the three DNA methylation related gene clusters and estimation of DNA methylation score (DMS)

We then analyzed the DEGs between the glioma belonging to the three distinct patterns of DNA methylation regulator patterns using the “limma” R package and identified 8291 pattern-related DEGs (Supplementary Figure 6D). GO analysis showed that these DEGs were enriched in some immune response and DNA methylation-related pathways such as asparagine N-linked glycosylation, cell cycle, DNA-dependent DNA replication, and Class I MHC-mediated antigen processing and presentation (Figure 3E). Moreover, univariate Cox regression analysis demonstrated that 5679 DEGs were associated with prognosis of the glioma patients. We then classified the glioma patients into different genomic phenotypes based on the expression levels of the 5679 prognosis-related DEGs. Furthermore, we performed unsupervised cluster analysis of the glioma patients based on the expression patterns of 5679 prognosis-related DEGs and identified three distinct gene clusters, A, B, and C (Supplementary Figure 6A–6C). These results confirmed three distinct DNA methylation-related gene clusters in the gliomas. Patients in the DNA methylation-related gene cluster B were associated with WHO IV grade, old age, non-codeletion status of 1p/19q, and wild-type IDH gene status (Figure 4A). Kaplan-Meier survival curves analysis showed that glioma patients in the DNA methylation-related gene cluster B were associated with poorer OS compared to patients in gene clusters A and C (P < 0.05; Figure 4B). Moreover, consistent with the clinicopathological features, glioma samples in gene cluster B showed higher proportions of resting CD4 memory T cells, resting NK cells, and resting mast cells in the TME and were significant positively correlated with the activation of EMT and pan-F-TBRS, thereby demonstrating the immune-desert phenotype (Figure 4C, 4D). Glioma in DNA methylation-related gene cluster C showed higher proportions of activated NK cells, activated DCs, and activated mast cells, and were significantly associated with activation of immune-response and immune-checkpoint related genes, thereby demonstrating immune-inflamed phenotype and activated tumor immune microenvironment (Figure 4C; Supplementary Figure 6E).

These results demonstrated that DNA methylation changes contributed to three distinct immune phenotypes in the glioma. However, we could not correctly forecast the DNA methylation status of individual glioma samples using unsupervised cluster analysis. Therefore, we calculated the DNA methylation score (DMS) of each glioma, based on the expression of 5679 prognosis-related DEGs. Our results showed significant difference in the DMS between the three DNA methylation modification patterns and DNA methylation-related gene clusters, with pattern B and gene cluster B showing the lowest median DMS (Figure 4E, 4F). WHO grade II gliomas showed the highest median DMS and WHO grade IV glioma showed the lowest median DMS (Figure 4G). The association between the DNA methylation modification patterns, tumor grade, DNA methylation-related gene clusters, and DMS categories are presented in the Sankey diagram (Figure 5D, Supplementary Table 4).

The glioma samples were classified into low- and high-DMS categories based on a cutoff value of 7.705903. Kaplan-Meier survival analyses showed that OS was significantly shorter in glioma patients with low-DMS compared to those with high-DMS across all grades in the gathered glioma cohort (Figure 5A). Moreover, LGG and GBM patients with low DMS showed significantly shorted OS compared to those with high DMS (P < 0.001; Figure 5B, 5C). Further, we analyzed the correlation between DMS and clinicopathological characteristics of the glioma patients. We observed significant differences in the DMS between glioma patient groups stratified by WHO grade (P < 0.001), survival status (P < 0.001), age (P < 0.05), 1p/19q status (P < 0.001), and isocitrate dehydrogenase (IDH) status (P < 0.01; Supplementary Figure 7A).

Univariate Cox regression analysis showed that WHO grade, age, gender, IDH status, 1p/19q status, and DMS were strongly associated with the OS of glioma patients in our study. Multivariate Cox regression analysis revealed that WHO grade, age, IDH status, 1p/19q status, and DMS were positively correlated with the OS of patients with gliomas (Supplementary Figure 7C). These results confirmed that DMS was an independent prognostic marker in glioma.

Characterization of DMS in the TCGA dataset and its association with tumor somatic mutations in glioma samples

Next, we further validated the prognostic value of DMS in the glioma samples using the TCGA dataset. Unsupervised cluster analysis showed 5312 DEGs between the three DNA methylation modification patterns and different DNA methylation-related gene clusters (Supplementary Figure 8A–8D). Furthermore, Kaplan-Meier survival curves showed that OS was significantly lower in the glioma patients belonging to DNA methylation gene cluster C compared to those from DNA methylation-related gene clusters A and B (Supplementary Figure 8F); cluster C glioma samples showed high expression levels of most DNA methylation regulatory genes (Supplementary Figure 8E) and positive correlation with activated EMT 1-3 and pan-F-TBRS (Supplementary Figure 8G), thereby representing the immune-desert phenotype.

We classified glioma samples in the TCGA dataset into low- and high-DMS categories using the cutoff value of −28.66371. Kaplan-Meier survival curve showed that OS was significantly shorter for glioma patients with low DMS compared to those with high DMS across different grades based on the TCGA dataset (P < 0.001; Supplementary Figure 9A–9C). This data was in agreement with the results for the entire glioma cohort. Our results also validated that DMS was an independent prognostic biomarker for predicting OS of glioma patients in the TCGA cohort (Supplementary Figure 7D). Furthermore, the high DMS group showed significantly OS in glioma patients belonging to different grades (ALL, LGG, and GBM) in the CGGA1 (P < 0.001; Supplementary Figure 9D–9F), CGGA2 (P < 0.001; Supplementary Figure 9G–9I), GSE180474 (P < 0.001; Supplementary Figure 9J–9L), and GSE16011 (P < 0.001; Supplementary Figure 9M–9O) datasets. We also observed positive correlation between DMS and other clinicopathological characteristics of the glioma patients based on the TCGA dataset (P < 0.01; Supplementary Figure 7B).

We then investigated differences in the TMB between the low DMS and high DMS subgroups of the TCGA glioma dataset. Low-DMS subgroup demonstrated more extensive TMB than the high DMS subgroup (Figure 5G, 5H). TMB quantification analysis showed that gliomas with low DMS were positively associated with higher TMB (Figure 5E, 5F). Our results also showed that gliomas with high DMS positively correlated with high TMB. Previous studies have shown that TMB was negatively associated with OS of glioma patients and high TMB was associated with reduced infiltration of immune cells into the glioma [25, 26]. Therefore, our results would provide novel perspective for exploring the mechanisms of DNA methylation modification in TMB, and shaping of TME landing.

Predictive value of DMS in immunotherapy

CTLA-4 and PD-1 blockades therapies significantly improve survival rates in many cancer types [27, 28]. Tumor mutation load (TML), PD-L1 expression levels, microsatellite instability (MSI) status, Tumor Immune Dysfunction and Exclusion (TIDE) scores, and IPS, are used to determine the status of tumor immune response [29, 30]. Because our study demonstrated strong association between DMS and tumor immune response, we further investigated the correlation between TIDE scores and DMS. We observed significantly reduced TIDE scores in the high-DMS subgroup of the gather glioma cohort and five glioma datasets (Figure 6A–6F). These results demonstrated that the expression levels of DNA methylation regulator significantly influenced the tumor immune response in glioma.

Therefore, we investigated whether DMS affected the clinical response to anti-PD-L1/PD-1 immunotherapy in the IMvigor210 and GSE78220 cohorts. The results showed that tumors with high DMS were associated with significantly increased OS in both immunotherapy cohorts (P < 0.05; Figure 6G, 6K). We further demonstrated significant differences in the clinical response of the high DMS subgroup patients to anti-PD-1/PD-L1 immunotherapy compared to the low DMS subgroup (Figure 6H, 6I, 6L, 6M). Furthermore, patients in the high DMS subgroup were positively associated with significantly high expression of PD-L1 in the two immunotherapy cohorts, thereby demonstrating the potential molecular mechanism underling the clinical response to immunotherapy (P < 0.05; Figure 6J, 6N).

In summary, our study showed that DNA methylation modification patterns were significantly associated with the immune response of the glioma. Moreover, the DNA methylation-related gene signature accurately predicted clinical responses to anti-PD-1/PD-L1 immunotherapy.

Pre-processing of patient data from multiple glioma datasets and two immune-checkpoint blockade cohorts

The flow chart of our research strategy is shown in Supplementary Figure 1A. We evaluated mRNA expression data and clinicopathological information of glioma patients in The Cancer Genome Atlas (TCGA), GEO, and Chinese Glioma Genome Atlas (CGGA) databases. We excluded samples with missing prognostic information. Finally, five glioma datasets including 2228 glioma samples ((TCGA (n = 648), CGGA1 (n = 309), CGGA2 (n = 593), and GEO: GSE16011 (n = 263), GSE108474 (n = 415)) were included for further analysis. The mRNA expression and clinicopathological information for the TCGA glioma datasets were downloaded from the University of California Santa Cruz (UCSC) Xena browser. The mRNA expression and clinicopathological information for the glioma samples in the GSE16011 and GSE108474 datasets were downloaded from the GEO database (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?) and those from the CGGA1 and CGGA2 datasets were downloaded from CGGA database (http://www.cgga.org.cn/). We also downloaded mRNA expression and clinicopathological data of patients in the two independent anti-PD-L1 immunotherapy cohorts, namely, IMvigor210, which included 298 urothelial cancer patients that underwent atezolizumab treatment (n = 298; http://research-pub.gene.com/IMvigor210CoreBiologies/); and GSE78220, included 26 metastatic melanoma patients that underwent treatment with pembrolizumab (n = 26; https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?) [11, 43]. Then the gene expression data of glioma samples was transformed into transcripts per kilobase million (TPM) values, and we merged and removed batch effect from non-biological technical biases by the “ComBat” algorithm of sva package. The information of the five gathered glioma databases was shown in Supplementary Table 1. The genomic mutation data was obtained from the UCSC Xena browse. The RNA-seq data of all samples was normalized with the following formula: log₂ (N+1).

Use of immunohistochemistry assay to validate bioinformatics results

We collected 10 NBTs and 10 LGG tissues from the Second Affiliated Hospital of Nanchang University from June 2020 to April 2022. Our study was approved by the Ethics Committee of this hospital. We performed the immunohistochemistry assay on human tissues by methods described preciously.

Protein-protein interaction network between 20 DNA methylation regulators

We retrieved the literatures databases and identified 20 DNA methylation regulatory genes for analysis in this study, containing 3 writers, 3 erasers, and 14 readers. The STRING and Cytoscape databases were used to construct the protein-protein interaction (PPI) network between these 20 DNA methylation regulatory proteins [44, 45].

Unsupervised clustering of DNA methylation regulatory genes

The samples in the five integrated glioma datasets were classified according to distinct DNA methylation modification patterns based on the expression levels of various DNA methylation regulatory genes using the R package “ConsensusClusterPlus” [46]. The Euclidean distance was used to calculate similarity distance between the glioma patients. The k-means algorithm was used for validation with 50 iterations and 1000 times repetitions. Principal component analysis (PCA) was performed to verify the classification into multiple patterns based on DNA methylation regulatory genes expression.

Gene set variation analysis and annotation of gene ontology terms

The gene set variation analysis (GSVA) R package was used to identify molecular functions related to different patterns of DNA methylation modifications. The gene signatures were obtained from the Molecular Signature Database (MSigDB) using the gene set “c2.cp.kegg.v6.2 symbols” [47]. The differentially expressed genes (DEGs) between the different patterns of DNA methylation modifications were obtained using the R package of “limma” [48]. Geno Ontology (GO) analysis was performed to determine the biological functions related to the DNA methylation regulatory related genes using P < 0.05 as the cut-off value [49].

Estimation of the relative abundance of immune cell types in the TME of gliomas

The single-sample gene-set enrichment analysis (ssGSEA) formula was performed using the “GSVA” R package to determine the relative abundance of the immune cell types in the glioma TME. A list of representative marker genes that represent different immune cells types were acquired from Charoentong’s study (Supplementary Table 5) that summarized 366 microarrays of several immune cell types containing activated mast cells, plasma cells, activated dendritic cells (DCs), activated NK cells, natural killer T cells, activated CD8+ T cells, eosinophils, activated CD4+ T cells, macrophages, and others [50, 51]. The enrichment scores were computed using the ssGSEA algorithm to determine the relative abundance of tumor-infiltrating immune cells in the glioma patients.

Generation of the DNA methylation score (DMS)

We developed a risk score system based on the DNA methylation modification in the glioma samples. The DNA methylation score (DMS) was calculated by: first identifying survival-associated DEGs using the “survival” R package, PCA was used to evaluate DMS by selecting PC1 and PC2, DMS was based on the largest block of highly correlated survival-associated DEGs and was calculated as follows: DMS = ∑ (PC1_i + PC2_i), where i, represents the expression of survival-associated DEGs expression in the cohort of glioma patients [52, 53].

Statistical analysis

The Shapiro-Wilk normality test was used to analyze normality of the variables [54]. The expression levels of DNA methylation regulators between two normally distributed groups were analyzed using unpaired t tests, whereas, the Wilcoxon rank-sum test was used to analyze non-normally distributed variables. Kruskal-Wallis or one-way ANOVA tests were used to compare differential gene expression between the three subgroups [55]. Spearman and distance correlation analyses were performed to determine the relationships between different subgroups. We computed the best cutoff point for the overall survival (OS) rates using R package of “survminer”. Then prognostic analysis was performed using the Kaplan-Meier survival curves [56]. Log-rank tests were used to estimate the differences between pairs of glioma subtypes. Univariate and multivariate Cox regression analysis was performed to determine the hazard ratios (HRs) for the 20 DNA methylation regulators and to determine if DMS was an independent prognostic biomarker. The mutation landscape of the DNA methylation regulatory genes in the TCGA dataset was determined using “maftools” R package, the landscape of copy number variation (CNVs) in the DNA methylation regulatory genes was determined using the “RCircos” R package [57]. Statistical analysis was performed using the R programming language v3.6.3, P < 0.05 was considered statistically significant.