A Novel DNA Methylation-Based Signature Can Predict the Responses of MGMT Promoter Unmethylated Glioblastomas to Temozolomide

Glioblastoma (GBM) is the most malignant glioma, with a median overall survival (OS) of 14–16 months. Temozolomide (TMZ) is the first-line chemotherapy drug for glioma, but whether TMZ should be withheld from patients with GBMs that lack O6-methylguanine-DNA methyltransferase (MGMT) promoter methylation is still under debate. DNA methylation profiling holds great promise for further stratifying the responses of MGMT promoter unmethylated GBMs to TMZ. In this study, we studied 147 TMZ-treated MGMT promoter unmethylated GBM, whose methylation information was obtained from the HumanMethylation27 (HM-27K) BeadChips (n = 107) and the HumanMethylation450 (HM-450K) BeadChips (n = 40) for training and validation, respectively. In the training set, we performed univariate Cox regression and identified that 3,565 CpGs were significantly associated with the OS of the TMZ-treated MGMT promoter unmethylated GBMs. Functional analysis indicated that the genes corresponding to these CpGs were enriched in the biological processes or pathways of mitochondrial translation, cell cycle, and DNA repair. Based on these CpGs, we developed a 31-CpGs methylation signature utilizing the least absolute shrinkage and selection operator (LASSO) Cox regression algorithm. In both training and validation datasets, the signature identified the TMZ-sensitive GBMs in the MGMT promoter unmethylated GBMs, and only the patients in the low-risk group appear to benefit from the TMZ treatment. Furthermore, these identified TMZ-sensitive MGMT promoter unmethylated GBMs have a similar OS when compared with the MGMT promoter methylated GBMs after TMZ treatment in both two datasets. Multivariate Cox regression demonstrated the independent prognostic value of the signature in TMZ-treated MGMT promoter unmethylated GBMs. Moreover, we also noticed that the hallmark of epithelial–mesenchymal transition, ECM related biological processes and pathways were highly enriched in the MGMT unmethylated GBMs with the high-risk score, indicating that enhanced ECM activities could be involved in the TMZ-resistance of GBM. In conclusion, our findings promote our understanding of the roles of DNA methylation in MGMT umethylated GBMs and offer a very promising TMZ-sensitivity predictive signature for these GBMs that could be tested prospectively.

Introduction

Glioma is the most common type of malignant brain tumor in adults (Jiang et al., 2016; Chai et al., 2019b). Glioblastoma (GBM, WHO IV) is the most malignant glioma, accounting for 50–60% of total glioma (Louis et al., 2016). Currently, the prognosis for patients with GBM is still dismal, with a median overall survival (OS) of 14–16 months (Jiang et al., 2016; Louis et al., 2016; Chai et al., 2019a). The alkylating agent temozolomide (TMZ) is the first-line chemotherapy drug for glioma. TMZ is used concurrently with radiation and then provided as monotherapy during adjuvant treatment. The promoter methylation level of the O6-methylguanine-DNA methyltransferase (MGMT), a ubiquitous DNA repair enzyme which can rapidly reverse alkylation at the O6 position, has been acknowledged as a predictive marker for TMZ sensitivity (Hegi et al., 2005; Chai et al., 2019a; Chai et al., 2019e). MGMT promoter methylated GBM displays higher sensitivity to TMZ treatment than MGMT promoter unmethylated GBM (Hegi et al., 2005; Wick et al., 2014; Chai et al., 2019a). However, we noticed that the prognosis for TMZ treated MGMT promoter unmethylated GBM is still largely heterogeneous, indicating that some other factors may also affect the sensitivity of MGMT promoter unmethylated GBM to TMZ treatment. Thus, further stratification of these GBM is urgently needed.

In the central nervous system, DNA methylation profiling has been used as a robust and reproducible method to further stratify the tumors into different subgroups (Sturm et al., 2012; Pajtler et al., 2015; Sturm et al., 2016). Moreover, general DNA methylation or a group of CpGs methylation profiling could also serve as biomarkers to evaluate drug- or radio-therapeutic sensitivity in various diseases, including tumors (Kumar et al., 2018; Zhao et al., 2018b; Chen et al., 2019b). In a recent study, a five-CpG DNA methylation score has shown its value in predicting metastatic-lethal outcomes of males suffering localized prostate cancer, treated with radical prostatectomy (Zhao et al., 2018b). The rapid accumulation of DNA methylation datasets makes it also possible to further stratify the glioma and may uncover novel biomarkers for management of gliomas. Recently, DNA methylation profiling of 23 DNA damage response (DDR) genes was shown to be associated with benefit from RT or TMZ therapy in IDH mutant low-grade glioma (Bady et al., 2018). Nevertheless, whether a group of CpGs DNA methylation profiling can predict the TMZ sensitivity of MGMT promoter unmethylated GBM remains unclear.

Here, we aimed to identify TMZ-sensitive GBMs in the entity of MGMT promoter unmethylated GBMs, using DNA methylation profiling. We adopted 107 and 40 TMZ treated MGMT promoter unmethylated GBMs as the training set and the validation set, respectively. We identified a list of CpGs whose methylation levels are significantly associated with the OS of TMZ-treated MGMT promoter unmethylated GBMs by univariate Cox regression analyses. Based on this, we developed a 31-CpGs TMZ therapeutic prognosis risk signature in the MGMT promoter unmethylated GBMs. This risk signature could successfully identify a subgroup of TMZ treated MGMT promoter unmethylated GBMs which have a similar prognosis when compared with the TMZ treated MGMT promoter methylated GBMs.

Materials and Methods

Samples Information

A total of 376 cases were enrolled in this study according to the following criteria: (a) diagnosed with GBM; (b) the DNA methylation data could be obtained; (c) the TMZ treatment option is available. The DNA methylation data and corresponding clinicopathological features for these cases were obtained from The Cancer Genome Atlas (TCGA) (http://cancergenome.nih.gov/). Within the 376 cases, the DNA methylation information of 279 cases (the 27K cohort) was collected from the HumanMethylation27 (HM-27K) BeadChips dataset, and the other 97 cases (the 450K cohort) were obtained from the HumanMethylation450 (HM-450K) BeadChips dataset. Clinicopathological information for all cases is summarized in Supplementary Table 1.

Of all 279 cases in the 27K cohort, 107 cases who received TMZ treatment and also with unmethylated MGMT were used to investigate the TMZ therapeutic prognosis value of CpGs methylation levels, and we also developed a risk signature using these cases. Of the 97 cases in the 450K cohort, 40 TMZ treated cases with unmethylated MGMT were used as the validation cohort. Clinicopathological information for these 147 cases is summarized in Table 1. There is no statistically significant difference for the clinicopathological features between the training and validation cohorts.

TABLE 1

Table 1 Clinicopathological characteristics for MGMT unmethylated GBM patients who received TMZ.

Analytical Approach

The approach and workflow for the selection of TMZ therapeutic prognosis associated CpGs, functional annotation for the genes corresponding to these CpGs, development and validation of a TMZ therapeutic prognostic risk signature, analysis of the correlation between the risk signature and other clinicopathological features, and the functional analysis of genes associated with the risk signature are summarized in Figure 1.

FIGURE 1

Figure 1 The workflow for this study. The workflow for the selection of TMZ therapeutic prognosis related CpGs, development and validation of a TMZ therapeutic prognostic risk signature, and the functional analysis of genes that are correlated with the signature risk score.

Identification of the Risk Signature

We performed univariate Cox regression analysis of the CpGs methylation to identify CpGs significantly correlated with the prognosis of TMZ treated MGMT unmethylated GBM in the 27K cohort. Then, we used the least absolute shrinkage and selection operator (LASSO) Cox regression algorithm to develop an optimal risk signature with the minimum number of CpGs (Dai et al., 2018; Zhou et al., 2018; Chai et al., 2019d). Finally, a set of 31 CpGs and their coefficients were determined by the minimum criteria, which involves selecting the best penalty parameter λ associated with the smallest 10-fold cross validation within the training set. The risk score for the risk signature was calculated using the formula:

$$\text{Risk\hspace{0.17em}score\hspace{0.17em}}={\Sigma }_{i=1}^{n}Coe{f}_i\ast x_i$$

where Coef is the coefficient and x_i is the beta-value of each selected CpGs. In both groups (cohorts), we used the beta-value [beta-value = the methylated signal/(methylated signal + unmethylated signal)] to represent the methylation level of each CpGs. Since the Risk score was calculated as a weighted sum of the methylation level of all selected CpG sites (Chai et al., 2019d; Chen et al., 2019a), we just used the original beta value of each CpG sites to calculate the risk scores.

We did not directly compare the samples in two different groups (cohorts). In order to avoid the bias caused by the different arrays, we only compared the methylation levels among samples in the same cohort. We first developed the risk signature in 107 samples used HumanMethylation27 (HM-27K) BeadChips. Then, we used another 40 samples to validate the prognostic value of the proposed signature. Patients were divided into “high-risk” and “low-risk” groups using the respective median risk score as the cutoff value in both the training and validation datasets.

Bioinformatic Analysis

Significance analysis of microarray (SAM) was performed to identify differentially expressed genes within the low- and high- risk scores. We performed Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses with the Database for Annotation, Visualization, and Integrated Discovery (http://david.abcc.ncifcrf.gov/home.jsp) to functionally annotate genes corresponding to the CpGs with prognosis of TMZ treated MGMT unmethylated GBM and genes that were differentially expressed between the low- and high-risk groups in the 27K cohort. Gene Set Enrichment Analysis (GSEA) was performed to investigate the functions of genes that were differentially expressed between the low- and high-risk groups in the 27K cohort.

Statistical Analysis

We used the nonparametric test to compare the distribution of age between the low- and high-risk groups, and Chi-square tests were used to compare the distribution of other clinicopathological features. A one-way analysis of variance was performed to compare the risk scores in patients grouped by the TCGA defined subtypes. Student’s t test was performed to compare the risk scores in patients grouped by other clinical or molecular-pathological characteristics.

Univariate and multivariate Cox regression analysis was performed to determine the prognostic value of the risk score and various clinical and molecular–pathological characteristics.

The Kaplan–Meier method with a two-sided log-rank test was used to compare the OS of patients stratified by the risk scores or other clinicopathological features. All statistical analyses were conducted using R v3.4.1 (https://www.r-project.org/), SPSS 16.0 (SPSS, Inc., Chicago, IL, USA) and Prism 7 (GraphPad Software, Inc., La Jolla, CA, USA).

Results

A Set of CpGs’ Methylation Profile Could Predict the TMZ Therapeutic Response of MGMT Unmethylated GBMs

To assess the TMZ therapeutic response value of the methylation of CpGs, we performed univariate Cox regression analysis of all CpGs methylation levels in the 107 TMZ treated MGMT unmethylated GBMs of the 27K cohort. We found that the methylation levels of 3,565 CpGs were significantly correlated with the OS of these GBMs (Supplementary Table 2). Based on the methylation profile of these genes, we could divide the 107 TMZ treated MGMT unmethylated GBMs into 3 clusters (Cluster A–C) in the heatmap (Figure 2A). We observed that patients in the Cluster A had significantly shorter survival than patients in the Cluster B and C, and the patients in the Cluster B and C had a similar OS with the TMZ treated MGMT methylated GBM patients (Figure 2B).

FIGURE 2

Figure 2 TMZ therapeutic prognosis associated CpGs’ methylation profile in MGMT unmethylated GBMs. (A) Heatmap showing the methylation levels of the 3,565 GpGs associated with the overall survival of TMZ treated patients with MGMT unmethylated GBMs. The MGMT unmethylated GBMs could be clustered into 3 clusters (Cluster A–C) according to the CpGs methylation levels. (B) Kaplan–Meier overall survival (OS) curves of TMZ treated MGMT unmethylated GBM patients (stratified by Cluster A–C) and TMZ treated MGMT methylated GBM patients. (C, D) GO biological process terms (C) and KEGG pathways (D) enriched among the genes positively and negatively corresponding to the 3,565 GpGs.

We also investigated the functions of the respective genes for the 3,565 CpGs. Three thousand one hundred eighty-two of these CpGs methylation levels were found to have a HR < 1 and were considered protective-associated, and the remaining 383 CpGs methylation levels with a HR >1 were considered risk-associated. GO terms of biological progress (BP) and KEGG pathway analysis indicated that the genes corresponding to the protective-associated CpGs were enriched in the processes of mitochondrial translation, protein modification, cell cycle, DNA repair, others, and pathways in cancer (Figures 2C, D). In contrast, the genes corresponding to the risk-associated CpGs were mainly enriched in the cellular membrane-associated biological processes and pathways (Figures 2C, D).

Identification of a 31-CpGs Panel as a TMZ Therapeutic Prognosis Risk Signature in MGMT Unmethylated GBMs

We next sought to develop a representative “risk signature” with a small number of CpGs to predict the TMZ therapeutic responses of the MGMT unmethylated GBMs. We applied the LASSO Cox regression algorithm to the 3,565 CpGs in 107 GBMs of the 27K cohort (Figure 3A). Finally, a total of 31 CpGs were contained in the risk signature, and the respective genes and the coefficients of these CpGs were also shown (Figure 3B and Supplementary Table 3). Twenty-four of the 31 CpGs are located in the CpG islands of prospective genes, and 5 of the other 7 CpGs are located within 200 bp of the transcription start site of the prospective genes (Supplementary Table 3). Most of the genes corresponding to the 31 CpGs have been reported to be involved in the tumorigenesis or prognosis of cancer, including ATOH1, ATPIG1, ELL3, RBM15B, GATA4, TXN, DLX5, THSD4, Polr2d, LGALS3BP, HIST1H3D, FLRT1, IFI35 and OSBPL5. Among these genes, hypermethylation of THSD4 has been reported to be associated with the prediction of prognosis in GBM (Ma et al., 2015), and Polr2d expression is associated with the therapy response of GBM (Serao et al., 2011).

FIGURE 3

Figure 3 Identification of the risk signature could stratify the TMZ therapeutic prognosis of the MGMT unmethylated GBM. (A) Ten-fold cross validation for tuning parameter selection in the LASSO model. The minimum criterion was indicated by the dashed vertical line (left). (B) Heatmap shows the association of risk scores and clinicopathological features based on the methylation profile of the 31 CpGs in the signature. The coefficients were calculated by multivariate Cox regression analysis using LASSO. (C–D) Kaplan–Meier overall survival (OS) curves for TMZ treated patients with MGMT methylated GBMs, TMZ treated patients with MGMT unmethylated GBMs with low- or high-risk significance scores in the training set (C) and validation set (D), respectively.

We divided patients into high-risk and low-risk groups using their median risk-score as the cutoffs. We observed significant differences between the low- and high- risk groups with respect to IDH status (P = 0.0431), age (P = 0.0069) and TCGA defined subtype (0.0047), but no differences in gender or chromosome 7 gain combined with chromosome 10 loss (chr 7 gain and chr 10 loss) (Figure 3B and Supplementary Table 4).

Then we investigated the relationship between the risk signature and OS of TMZ treated MGMT unmethylated GBM patients. The data showed that patients with low-risk-scores had significantly longer OS than patients with high-risk-scores in both the training (P < 0.0001) and validation (P = 0.0331) datasets (Figures 3C, D). In addition, although the OS of MGMT methylated GBM patients was significantly longer than that of MGMT unmethylated GBM patients (Supplementary Figure 1), we noticed that the OS of MGMT unmethylated GBM patients in the low-risk group was similar to that of MGMT methylated GBM patients in both the training and validation datasets (Figures 3C, D).

Association of the Risk Signature and Other Clinicopathological Features

Considering that the TMZ therapeutic prognosis value of the risk signature may be associated with other known clinicopathological features, we examined this in the MGMT unmethylated GBMs. We observed that the risk scores were only significantly different between patients stratified by age (P < 0.05), rather than gender, chr 7 gain and chr 10 loss, and the TCGA defined subtypes (Supplementary Figure 2). We did not compare the risk scores in patients with different IDH status, as there were only four IDH-mutant patients.

We also performed univariate and multivariate Cox regression analyses in the TMZ treated MGMT unmethylated GBMs of both the training and validation datasets. By univariate analysis, the risk score [hazard ratio (HR) = 12.674 (7.661–20.968) in the training set; HR = 1.685 (1.058–2.682) in the validation set] and age [HR = 1.029 (1.009–1.048) in the training set; HR = 1.075 (1.023–1.13) in the validation set] were significantly correlated with the OS in both two datasets (Table 2). When including these factors into the multivariate Cox regression analysis, the risk score remained significantly associated with the OS in the training [HR = 12.748 (7.767–21.173)] and validation [HR = 2.157 (1.139–4.086)] datasets (Table 2). These results indicated that the risk score can independently predict the TMZ therapeutic prognosis of patients with MGMT unmethylated GBMs.

TABLE 2

Table 2 Univariate and multivariate Cox regression analyses for the risk score in the training and validation set, respectively.

We also investigated the association of risk scores and clinicopathological features in all GBM. We found that the risk scores were only significantly different between patients with different IDH status (P < 0.0001) or between Proneural subtype and Mesenchymal subtype (P < 0.01), but not between patients stratified by age, gender, MGMT promoter methylation status, chr 7 gain and chr 10 loss, or treated with or without TMZ (Supplementary Figure 3).

Prognosis Value of the Risk Signature in Stratified GBMs

To further understand the TMZ therapeutic prognostic value of the risk signature in MGMT unmethylated GBMs, we compared the OS of MGMT unmethylated GBMs patients stratified by TMZ treatment status in the low-risk and high-risk groups respectively. The results indicated that patients with TMZ treatment had longer OS than that of patients without TMZ treatment in the low-risk group of both the training set (P < 0.0001, Figure 4A) and the validation set (P = 0.0456, Figure 4F). In contrast, there was no significant difference between patients with or without TMZ treatment in the high-risk group (Figures 4B, G).

FIGURE 4

Figure 4 Clinical outcomes prediction of the signature in patients with stratified GBMs. (A–B) Kaplan–Meier overall survival (OS) curves for MGMT unmethylated GBM patients with or without TMZ treatment in the low-risk group (A) and high-risk groups (B) of the training set. (C–E) Kaplan–Meier overall survival (OS) curves for stratified GBM patients (C) MGMT unmethylated GBM without TMZ; (D) MGMT methylated GBM with TMZ; (E) MGMT methylated GBM without TMZ) with low- or high-risk scores in the training set. (F–J) Kaplan–Meier overall survival (OS) curves for stratified GBM patients in the validation set.

We also investigated the prognostic value of the risk signature in other stratified GBMs. We respectively stratify the GBM patients into four subgroups according to MGMT status and TMZ treatment option. In the training set, the risk signature could not stratify the prognosis of three subgroups (TMZ non-treated MGMT unmethylated GBM, TMZ treated MGMT methylated GBM, and TMZ non-treated MGMT methylated GBM) (Figures 4C–E). Similar results could also be observed in the validation set except TMZ non-treated MGMT unmethylated GBM (Figures 4G, H).

The Potential Functions Underlying the TMZ Therapeutic Prognostic Value of the Risk Signature

To determine the functional differences between the high-risk and low-risk cases of the TMZ treated MGMT unmethylated GBM in the 27K cohort, we identified the differentially (P < 0.05) expressed genes by SAM (Figure 5A). GO and KEGG analyses revealed that extracellular matrix related biological processes and signaling pathways were significantly enriched in the high-risk group (Figures 5B, C). In contrast, the biological processes of T cell differentiation, nervous system development, and transcription were significantly enriched in the low-risk group (Figure 5B). Meanwhile, GSEA also indicated that the high-risk cases showed enrichment of “regulation of endothelial cell apoptotic process,” “extracellular structure organization,” “aminoglycan metabolic process,” and “extracellular matrix disassembly biological progresses” (Figure 5D). Moreover, the hallmarks of “epithelial–mesenchymal transition,” “PI3K-AKT-mTOR signaling,” “glycolysis”, and “angiogenesis” also enriched in the high-risk cases (Figure 5E). The results indicated that the extracellular matrix related functions and mesenchymal phenotype could contribute to the TMZ-resistant of glioma.

FIGURE 5

Figure 5 Functional annotation for genes differentially expressed between low- and high-risk groups. (A) The differential genes between low- and high-risk groups are shown by green (enriched in the low-risk group) and red (enriched in the high-risk group) dots. (B–C) Go analysis (B) and KEGG analysis (C) are used to evaluate differential genes between low-and high-risk groups. (D and E) GSEA analysis reveals the biological processes (D) and cancer hallmarks (E) enriched in the high-risk groups.

Discussion

Undoubtedly, MGMT promoter methylation status is critical for the chemotherapeutic management of glioma, especially for GBM (Hegi et al., 2005; Chai et al., 2019a; Chai et al., 2019e). However, whether TMZ should be withheld from patients with GBMs that lack MGMT promoter methylation is still under debate, and some of these patients indeed benefit from the treatment (Wick et al., 2014). Thus, it is critical to uncover novel biomarkers to identify TMZ-sensitive individuals with MGMT promoter unmethylated GBMs. In this study, we successfully developed a 31-CpG methylation signature which could identify the TMZ-sensitive GBMs in the MGMT promoter unmethylated GBMs from both the training and validation datasets, and OS of these TMZ-sensitive GBMs is similar to that of the MGMT promoter methylated GBMs after TMZ treatment in both two datasets. Considering the robust and reproducible nature of DNA methylation in the classification of brain tumors, this signature has great value in predicting the TMZ sensitivity of the GBMs that lack MGMT promoter methylation.

In this study, we systematically investigated 107 MGMT promoter unmethylated GBMs to obtain the TMZ therapeutic prognostic value of each of the CpGs that were included in the HM-27K BeadChip, and we identified that 3,565 CpGs are significantly associated with the OS of these GBMs. Previous studies have indicated that abnormal metabolism could alter the response of tumor cells to chemotherapy through inhibiting the activities of DNA repair enzymes (Gusyatiner and Hegi, 2018). DNA instability and DNA injury repair have been linked to the chemo-resistance of cancer cells (Kanai et al., 2012; Roos et al., 2018; Zhao et al., 2018a; Ha Thi et al., 2019; Zhang et al., 2019). Here we also investigated the functions of genes corresponding to the 3,565 CpGs, and the results indicated that biological processes or pathways of mitochondrial translation, cell cycle and DNA repair could be involved in the TMZ-sensitivity of MGMT promoter unmethylated GBMs. Given that DNA proliferation rate is positively correlated to the sensitivity to chemotherapy (Li et al., 2017; Krell et al., 2019; Qiang et al., 2018), our finding supports that transcriptional activities of genes enriched in mitochondria, DNA injury and repair, and cell cycle processes could be important in the sensitivity of GBM cells to TMZ chemotherapy.

The extracellular matrix (ECM) components and their partners, including the glycosaminoglycans, glycoproteins, and proteoglycans, play a crucial role in the glioma invasion through promoting tumor cell migration and angiogenesis (Ferrer et al., 2018). The up-regulation of ECM partners, such as CD44, has been acknowledged as a marker for the “proneural–mesenchymal transition” of GBM cells (Yang et al., 2017). Here, we noticed that not only the hallmark of epithelial–mesenchymal transition but also ECM related biological processes and pathways were highly enriched in the MGMT unmethylated GBMs with the high-risk score, indicating that enhanced ECM activities could be involved in the TMZ-resistance of GBMs. This may be associated with the roles of ECM in regulating the extracellular microenvironments and intracellular signaling pathways (Wang et al., 2018). Chemokine (C-X-C motif) ligand 12 (CXCL12) and its receptors CXCR4 and CXCR7, which are stored in or attached to the ECM, are extremely important in forming a more invasive and resistant phenotype of glioma (Gatti et al., 2013; Zhao et al., 2018a). Recently, we also identified that the glycoprotein ADAMTS4, which is important for the upregulation of integrins, is also a novel immune-related biomarker for the primary GBM (Zhao et al., 2019). Transforming growth factor-beta (TGF-β), an ECM-bound bioactive factor, is involved in both the activation of NF-κB signaling and mesenchymal transition of GBM (Song et al., 2018; Batlle and Massague, 2019). Both of these two processes have been involved in the TMZ-resistance of GBM (Ming et al., 2017; Yang et al., 2017; Chai et al., 2019c; Chai et al., 2019d). All of these emphasize the value of the ECM in glioma TMZ sensitivity. Thus, the ECM and microenvironment should not be neglected in drug development, especially in developing an ideal in vitro drug screening model for glioma.

Chr 7 gain and Chr 10 loss is quite common in GBM (Bady et al., 2016; Chai et al., 2019a). Patients with high-grade gliomas harboring deletions of chromosomes 9p and 10q may benefit more from TMZ treatment (Wemmert et al., 2005), and the MGMT resides on chromosome 10q. Here, we also investigated the association between the risk signature and deletion of one copy of chromosome 10, and the results indicated that the predictive value of the risk signature was not affected by the status of Chr 7 gain and Chr 10 loss. This finding excludes the possibility that the predictive value of the risk signature may be caused by the unbalanced MGMT expression between GBM with or without Chr 7 gain and Chr 10 loss. Moreover, we have reported that chromosome 10/10q deletion does not significantly affect MGMT expression of GBM in the TCGA dataset (Chai et al., 2019a).

In conclusion, our findings reveal the predictive value of DNA methylation profiling in GBMs with an unmethylated MGMT promoter. The developed 31-CpG methylation signature could accurately predict the TMZ-sensitivity of MGMT promoter unmethylated GBMs. Though the risk signature still needs to be confirmed in future prospective studies with specific test kits, our current findings can promote our understanding of the roles of DNA methylation in GBMs with an unmethylated MGMT promoter and also offer a very promising TMZ-sensitivity predictive signature for these GBMs.

Data Availability Statement

All methylation and clinical data used in this study were available from the TCGA database (http://cancergenome.nih.gov). Other information is available through contacting the corresponding authors.

Author Contributions

R-CC conceived and designed the study, R-CC, Y-ZC, and Q-WW crafted the literature search, figures, and tables and were responsible for the writing and critical reading of the manuscript. K-NZ, J-JL, and HH contributed to the data analysis and the critical reading of the manuscript. FW, Y-QL, and Y-ZW supervised the analysis and contributed to the critical reading of the manuscript.

Funding

This work was supported by the National Key Research and Development Program of China (2018YFC0115604), the National Natural Science Foundation of China (81773208, 81802994), the National Natural Science Foundation of China (NSFC)/Research Grants Council (RGC) Joint Research Scheme (81761168038), Beijing Municipal Administration of Hospitals’ Mission Plan (SML20180501, 2018.03-2022.02).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

The authors gratefully acknowledge contributions from the TCGA Network.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2019.00910/full#supplementary-material

References

Bady, P., Delorenzi, M., Hegi, M. E. (2016). Sensitivity analysis of the MGMT-STP27 Model and impact of genetic and epigenetic context to predict the MGMT methylation status in gliomas and other tumors. J. Mol. Diagn. 18 (3), 350–361. doi: 10.1016/j.jmoldx.2015.11.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Bady, P., Kurscheid, S., Delorenzi, M., Gorlia, T., van den Bent, M. J., Hoang-Xuan, K., et al. (2018). The DNA methylome of DDR genes and benefit from RT or TMZ in IDH mutant low-grade glioma treated in EORTC 22033. Acta Neuropathol. 135 (4), 601–615. doi: 10.1007/s00401-018-1810-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Batlle, E., Massague, J. (2019). Transforming growth factor-beta signaling in immunity and cancer. Immunity 50 (4), 924–940. doi: 10.1016/j.immuni.2019.03.024

PubMed Abstract | CrossRef Full Text | Google Scholar

Chai, R. C., Liu, Y. Q., Zhang, K. N., Wu, F., Zhao, Z., Wang, K. Y., et al. (2019a). A novel analytical model of MGMT methylation pyrosequencing offers improved predictive performance in patients with gliomas. Mod. Pathol. 32 (1), 4–15. doi: 10.1038/s41379-018-0143-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Chai, R. C., Wang, N., Chang, Y. Z., Zhang, K. N., Li, J. J., Niu, J. J., et al. (2019b). Systematically profiling the expression of eIF3 subunits in glioma reveals the expression of eIF3i has prognostic value in IDH-mutant lower grade glioma. Cancer Cell Int. 19, 155. doi: 10.1186/s12935-019-0867-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Chai, R. C., Wu, F., Wang, Q. X., Zhang, S., Zhang, K. N., Liu, Y. Q., et al. (2019c). m(6)A RNA methylation regulators contribute to malignant progression and have clinical prognostic impact in gliomas. Aging (Albany N.Y.) 11 (4), 1204–1225. doi: 10.18632/aging.101829

PubMed Abstract | CrossRef Full Text | Google Scholar

Chai, R. C., Zhang, K. N., Chang, Y. Z., Wu, F., Liu, Y. Q., Zhao, Z., et al. (2019d). Systematically characterize the clinical and biological significances of 1p19q genes in 1p/19q non-codeletion glioma. Carcinogenesis. doi: 10.1093/carcin/bgz102

CrossRef Full Text | Google Scholar

Chai, R. C., Zhang, K. N., Liu, Y. Q., Wu, F., Zhao, Z., Wang, K. Y., et al. (2019e). Combinations of four or more CpGs methylation present equivalent predictive value for MGMT expression and temozolomide therapeutic prognosis in gliomas. CNS Neurosci. Ther. 25 (3), 314–322. doi: 10.1111/cns.13040

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, D., Song, Y., Zhang, F., Wang, X., Zhu, E., Zhang, X., et al. (2019a). Genome-wide analysis of lung adenocarcinoma identifies novel prognostic factors and a prognostic score. Front. Genet. 10, 493. doi: 10.3389/fgene.2019.00493

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, S., Pu, W., Guo, S., Jin, L., He, D., Wang, J. (2019b). Genome-wide DNA methylation profiles reveal common epigenetic patterns of interferon-related genes in multiple autoimmune diseases. Front.Genet. 10 (223). doi: 10.3389/fgene.2019.00223

CrossRef Full Text | Google Scholar

Dai, W., Feng, Y., Mo, S., Xiang, W., Li, Q., Wang, R., et al. (2018). Transcriptome profiling reveals an integrated mRNA-lncRNA signature with predictive value of early relapse in colon cancer. Carcinogenesis 39 (10), 1235–1244. doi: 10.1093/carcin/bgy087

PubMed Abstract | CrossRef Full Text | Google Scholar

Ferrer, V. P., Moura Neto, V., Mentlein, R. (2018). Glioma infiltration and extracellular matrix: key players and modulators. Glia 66 (8), 1542–1565. doi: 10.1002/glia.23309

PubMed Abstract | CrossRef Full Text | Google Scholar

Gatti, M., Pattarozzi, A., Bajetto, A., Wurth, R., Daga, A., Fiaschi, P., et al. (2013). Inhibition of CXCL12/CXCR4 autocrine/paracrine loop reduces viability of human glioblastoma stem-like cells affecting self-renewal activity. Toxicology 314 (2-3), 209–220. doi: 10.1016/j.tox.2013.10.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Gusyatiner, O., Hegi, M. E. (2018). Glioma epigenetics: from subclassification to novel treatment options. Semin. Cancer Biol. 51, 50–58. doi: 10.1016/j.semcancer.2017.11.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Ha Thi, H. T., Kim, H. Y., Lee, Y. J., Kim, S. J., Hong, S. (2019). SMAD7 in keratinocytes promotes skin carcinogenesis by activating ATM-dependent DNA repair and an EGFR-mediated cell proliferation pathway. Carcinogenesis 40 (1), 112–120. doi: 10.1093/carcin/bgy121

PubMed Abstract | CrossRef Full Text | Google Scholar

Hegi, M. E., Diserens, A. C., Gorlia, T., Hamou, M. F., de Tribolet, N., Weller, M., et al. (2005). MGMT gene silencing and benefit from temozolomide in glioblastoma. N. Engl. J. Med. 352 (10), 997–1003. doi: 10.1056/NEJMoa043331

PubMed Abstract | CrossRef Full Text | Google Scholar

Jiang, T., Mao, Y., Ma, W., Mao, Q., You, Y., Yang, X., et al. (2016). CGCG clinical practice guidelines for the management of adult diffuse gliomas. Cancer Lett. 375 (2), 263–273. doi: 10.1016/j.canlet.2016.01.024

PubMed Abstract | CrossRef Full Text | Google Scholar

Kanai, R., Rabkin, S. D., Yip, S., Sgubin, D., Zaupa, C. M., Hirose, Y., et al. (2012). Oncolytic virus-mediated manipulation of DNA damage responses: synergy with chemotherapy in killing glioblastoma stem cells. J. Natl. Cancer Inst. 104 (1), 42–55. doi: 10.1093/jnci/djr509

PubMed Abstract | CrossRef Full Text | Google Scholar

Krell, A., Wolter, M., Stojcheva, N., Hertler, C., Liesenberg, F., Zapatka, M., et al. (2019). MiR-16-5p is frequently down-regulated in astrocytic gliomas and modulates glioma cell proliferation, apoptosis and response to cytotoxic therapy. Neuropathol. Appl. Neurobiol. 45 (5), 441–458 doi: 10.1111/nan.12532

PubMed Abstract | CrossRef Full Text | Google Scholar

Kumar, R., Liu, A. P. Y., Orr, B. A., Northcott, P. A., Robinson, G. W. (2018). Advances in the classification of pediatric brain tumors through DNA methylation profiling: from research tool to frontline diagnostic. Cancer 124 (21), 4168–4180. doi: 10.1002/cncr.31583

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, T., Li, S., Chen, D., Chen, B., Yu, T., Zhao, F., et al. (2017). Transcriptomic analyses of RNA-binding proteins reveal eIF3c promotes cell proliferation in hepatocellular carcinoma. Cancer Sci. 108 (5), 877–885. doi: 10.1111/cas.13209

PubMed Abstract | CrossRef Full Text | Google Scholar

Louis, D. N., Perry, A., Reifenberger, G., von Deimling, A., Figarella-Branger, D., Cavenee, W. K., et al. (2016). The 2016 World Health Organization classification of tumors of the central nervous system: a summary. Acta Neuropathol. 131 (6), 803–820. doi: 10.1007/s00401-016-1545-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Ma, J., Hou, X., Li, M., Ren, H., Fang, S., Wang, X., et al. (2015). Genome-wide methylation profiling reveals new biomarkers for prognosis prediction of glioblastoma. J. Cancer Res. Ther. 11 Suppl 2, C212–C215. doi: 10.4103/0973-1482.168188

PubMed Abstract | CrossRef Full Text | Google Scholar

Ming, J., Sun, B., Li, Z., Lin, L., Meng, X., Han, B., et al. (2017). Aspirin inhibits the SHH/GLI1 signaling pathway and sensitizes malignant glioma cells to temozolomide therapy. Aging (Albany N.Y.) 9 (4), 1233–1247. doi: 10.18632/aging.101224

CrossRef Full Text | Google Scholar

Pajtler, K. W., Witt, H., Sill, M., Jones, D. T., Hovestadt, V., Kratochwil, F., et al. (2015). Molecular classification of ependymal tumors across All CNS compartments, histopathological grades, and age groups. Cancer Cell 27 (5), 728–743. doi: 10.1016/j.ccell.2015.04.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Qiang, Z., Jun-Jie, L., Hai, W., Hong, L., Bing-Xi, L., Lei, C., et al. (2018). TPD52L2 impacts proliferation, invasiveness and apoptosis of glioblastoma cells via modulation of wnt/beta-catenin/snail signaling. Carcinogenesis 39 (2), 214–224. doi: 10.1093/carcin/bgx125

PubMed Abstract | CrossRef Full Text | Google Scholar

Roos, W. P., Frohnapfel, L., Quiros, S., Ringel, F., Kaina, B. (2018). XRCC3 contributes to temozolomide resistance of glioblastoma cells by promoting DNA double-strand break repair. Cancer Lett. 424, 119–126. doi: 10.1016/j.canlet.2018.03.025

PubMed Abstract | CrossRef Full Text | Google Scholar

Serao, N. V., Delfino, K. R., Southey, B. R., Beever, J. E., Rodriguez-Zas, S. L. (2011). Cell cycle and aging, morphogenesis, and response to stimuli genes are individualized biomarkers of glioblastoma progression and survival. BMC Med. Genomics 4, 49. doi: 10.1186/1755-8794-4-49

PubMed Abstract | CrossRef Full Text | Google Scholar

Song, Y., Chen, Y., Li, Y., Lyu, X., Cui, J., Cheng, Y., et al. (2018). Metformin inhibits TGF-beta1-induced epithelial-to-mesenchymal transition-like process and stem-like properties in GBM via AKT/mTOR/ZEB1 pathway. Oncotarget 9 (6), 7023–7035. doi: 10.18632/oncotarget.23317

PubMed Abstract | CrossRef Full Text | Google Scholar

Sturm, D., Orr, B. A., Toprak, U. H., Hovestadt, V., Jones, D. T. W., Capper, D., et al. (2016). New brain tumor entities emerge from molecular classification of CNS-PNETs. Cell 164 (5), 1060–1072. doi: 10.1016/j.cell.2016.01.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Sturm, D., Witt, H., Hovestadt, V., Khuong-Quang, D. A., Jones, D. T., Konermann, C., et al. (2012). Hotspot mutations in H3F3A and IDH1 define distinct epigenetic and biological subgroups of glioblastoma. Cancer Cell 22 (4), 425–437. doi: 10.1016/j.ccr.2012.08.024

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, X., Li, X., Dai, X., Zhang, X., Zhang, J., Xu, T., et al. (2018). Coaxial extrusion bioprinted shell-core hydrogel microfibers mimic glioma microenvironment and enhance the drug resistance of cancer cells. Colloids Surf. B Biointerfaces 171, 291–299. doi: 10.1016/j.colsurfb.2018.07.042

PubMed Abstract | CrossRef Full Text | Google Scholar

Wemmert, S., Ketter, R., Rahnenfuhrer, J., Beerenwinkel, N., Strowitzki, M., Feiden, W., et al. (2005). Patients with high-grade gliomas harboring deletions of chromosomes 9p and 10q benefit from temozolomide treatment. Neoplasia 7 (10), 883–893. doi: 10.1593/neo.05307

PubMed Abstract | CrossRef Full Text | Google Scholar

Wick, W., Weller, M., van den Bent, M., Sanson, M., Weiler, M., von Deimling, A., et al. (2014). MGMT testing—the challenges for biomarker-based glioma treatment. Nat. Rev. Neurol. 10 (7), 372–385. doi: 10.1038/nrneurol.2014.100

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, F., Liu, X., Liu, Y., Liu, Y., Zhang, C., Wang, Z., et al. (2017). miR-181d/MALT1 regulatory axis attenuates mesenchymal phenotype through NF-kappaB pathways in glioblastoma. Cancer Lett. 396, 1–9. doi: 10.1016/j.canlet.2017.03.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, Y., Li, J., Yi, K., Feng, J., Cong, Z., Wang, Z., et al. (2019). Elevated signature of a gene module coexpressed with CDC20 marks genomic instability in glioma. Proc. Natl. Acad. Sci. U S A 116 (14), 6975–6984. doi: 10.1073/pnas.1814060116

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, K., Yao, Y., Luo, X., Lin, B., Huang, Y., Zhou, Y., et al. (2018a). LYG-202 inhibits activation of endothelial cells and angiogenesis through CXCL12/CXCR7 pathway in breast cancer. Carcinogenesis 39 (4), 588–600. doi: 10.1093/carcin/bgy007

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, S., Leonardson, A., Geybels, M. S., McDaniel, A. S., Yu, M., Kolb, S., et al. (2018b). A five-CpG DNA methylation score to predict metastatic-lethal outcomes in men treated with radical prostatectomy for localized prostate cancer. Prostate. doi: 10.1002/pros.23667

CrossRef Full Text | Google Scholar

Zhao, Z., Zhang, K. N., Chai, R. C., Wang, K. Y., Huang, R. Y., Li, G. Z., et al. (2019). ADAMTSL4, a secreted glycoprotein, is a novel immune-related biomarker for primary glioblastoma multiforme. Dis. Markers 2019, 1802620. doi: 10.1155/2019/1802620

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, Z., Huang, R., Chai, R., Zhou, X., Hu, Z., Wang, W., et al. (2018). Identification of an energy metabolism-related signature associated with clinical prognosis in diffuse glioma. Aging (Albany N.Y.) 10 (11), 3185–3209. doi: 10.18632/aging.101625

CrossRef Full Text | Google Scholar