Liang Zhao1†
Peng Zhao- 1Department of Neurosurgery, The First Affiliated Hospital of Nanjing Medical University, Nanjing, China
- 2Department of Respiratory and Critical Care Medicine, The First Affiliated Hospital of Nanjing Medical University, Nanjing, China
Background: O6-methylguanine-DNA methyltransferase (MGMT) methylation status affects tumor chemo-resistance and the prognosis of glioblastoma (GBM) patients. We aimed to investigate the role of MGMT methylation in the regulation of GBM immunophenotype and discover an effective biomarker to improve prognosis prediction of GBM patients.
Methods: A total of 769 GBM patients with clinical information from five independent cohorts were enrolled in the present study. Samples from the Cancer Genome Atlas (TCGA) dataset were used as the training set, whereas transcriptome data from the Chinese Glioma Genome Atlas (CGGA) RNA-seq, CGGA microarray, GSE16011, and the Repository for Molecular Brain Neoplasia (REMBRANDT) cohort were used for validation. A series of bioinformatics approaches were carried out to construct a prognostic signature based on immune-related genes, which were tightly related to the MGMT methylation status. In silico analyses were performed to investigate the influence of the signature on immunosuppression and remodeling of the tumor microenvironment. Then, the utility of this immune gene signature was analyzed by the development and evaluation of a nomogram. In vitro experiments were further used to verify the immunologic function of the genes in the signature.
Results: We found that MGMT unmethylation was closely associated with immune-related biological processes in GBM. Sixty-five immune genes were more highly expressed in the MGMT unmethylated than the MGMT-methylated group. An immune gene-based risk model was further established to divide patients into high and low-risk groups, and the prognostic value of this signature was validated in several GBM cohorts. Functional analyses manifested a universal up-regulation of immune-related pathways in the high-risk group. Furthermore, the risk score was highly correlated to the immune cell infiltration, immunosuppression, inflammatory activities, as well as the expression levels of immune checkpoints. A nomogram was developed for clinical application. Knockdown of the five genes in the signature remodeled the immunosuppressive microenvironment by restraining M2 macrophage polarization and suppressing immunosuppressive cytokines production.
Conclusions: MGMT methylation is strongly related to the immune responses in GBM. The immune gene-based signature we identified may have potential implications in predicting the prognosis of GBM patients and mechanisms underlying the role of MGMT methylation.
Introduction
Glioma is a type of central nervous system (CNS) neoplasms and accounts for the majority of intracranial malignant tumors in adults. Glioblastoma (GBM), which is defined as the World Health Organization (WHO) grade IV glioma, shows a highly aggressive, heterogeneous, and lethal phenotype. Currently, the etiology of GBM is still largely unknown. Although surgical resection combined with chemotherapy or radiotherapy has been widely used as a routine clinical treatment, the prognosis of patients has even not improved significantly, with a median overall survival fewer than 15 months (Stupp et al., 2005). Due to the persistent proliferation of tumor cells and penetration into the surrounding tissues, total tumor resection is seldom possible (Zhang et al., 2017). Alkylating agents, such as temozolomide (TMZ) and procarbazine, are applied as first-line chemotherapeutic drugs (Hegi et al., 2008). However, many GBM patients do not respond to alkylating agents or develop chemo-resistance after a period of chemotherapy (Wick et al., 2014). Therefore, there is an urgent need to explore the underlying mechanisms of gliomagenesis and develop reliable therapeutic approaches.
The intracellular DNA repair enzyme O6-methylguanine-DNA methyltransferase (MGMT) catalyzes the DNA repair process by removing the alkylation of the O6 position of guanine (Pegg, 1990). By generating methylguanine adducts, DNA-alkylating drugs can trigger base mismatching during DNA replication, affecting the cell cycle and promoting the death of tumor cells (Karran and Bignami, 1992; Dolan and Pegg, 1997; Liu and Gerson, 2006). As a most striking modification form, methylation of the MGMT promoter can induce a decrease of MGMT expression and further improve the curative effect of chemotherapy with alkylating agents (Stupp et al., 2005). In the case of GBM, clinical trials revealed that patients with MGMT methylation benefit more than those with unmethylated tumors when received combined radiochemotherapy (Hegi et al., 2005). Furthermore, patients without MGMT methylation have a shorter overall survival time (Molenaar et al., 2014). So far, these underlying molecular mechanisms have not been exhaustively described, especially about tumor immune microenvironment.
Recently, breakthroughs in tumor immunotherapy have revolutionized the treatment of several cancer types (Fukumura et al., 2018). The specific intratumoral immune microenvironment can facilitate immune escape, drive tumorigenesis, and promote the malignant progression of the lesion (Heimberger et al., 2008; Nduom et al., 2015). Immune checkpoints play an intermediary role in the tumor-immune system dynamics, and inhibitory immune checkpoints have been studied intensively as therapeutic targets for multiple malignant tumors. In non-small cell lung cancer (NSCLC) and melanoma, the application of immunosuppressive agents has produced remarkable antitumor effects (Hodi et al., 2010; Rizvi et al., 2015). The discovery of lymphatic system vessels in the CNS offered new hope for developing novel immunotherapeutic methods for GBM (Louveau et al., 2015b). The first large phase III clinical trial of nivolumab combined with ipilimumab (PD-1 inhibitors) in recurrent GBM (NCT02017717) was launched in 2014. A recent clinical trial indicated that the neoadjuvant PD-1 blockade could upregulate the amount of T cells, promote the expression of interferon-γ-related genes, and further offer a promising survival benefit in recurrent GBM (Cloughesy et al., 2019). The previous study has demonstrated that isocitrate dehydrogenase (IDH) mutation status was closely associated with immune response in GBM, and a prognostic model based on IDH mutation was further developed to predict the clinical outcomes of patients (Qian et al., 2018). However, the role of MGMT methylation status in regulating tumor immunity in GBM remains to be discussed, and an MGMT methylation-related biomarker should be exploited.
Here, we utilized data from a large number of GBM samples from multiple publicly available datasets to comprehensively investigate the potential relationship between MGMT status and the immunological tumor microenvironment. We hypothesized that the shorter survival of patients with MGMT-unmethylated GBM is partly due to the pro-tumor immune response. More importantly, we developed a robust prognostic model based on MGMT methylation status to predict the clinical outcomes of GBM patients.
Materials and Methods
GBM Patient Datasets
We retrospectively collected a large-scale profile composed of five independent datasets involving 769 GBM patients (Supplementary Table 1). The RNA sequencing (RNA-seq) gene expression data for 165 GBM samples were obtained from The Cancer Genome Atlas (TCGA) database. The microarray and RNA-seq transcriptome data for 112 and 113 GBM patients were downloaded from the Chinese Glioma Genome Atlas (CGGA) database (www.cgga.org.cn). The microarray gene expression profile matrix files of the GSE16011 project (including 159 GBM samples) and the Repository for Molecular Brain Neoplasia (REMBRANDT) database (including 220 GBM samples) were downloaded from the Gene Expression Omnibus (GEO) database (www.ncbi.nlm.nih.gov/geo/). The corresponding profiles comprising detailed clinicopathological characteristics were also generated from each data source. The TCGA RNA-seq dataset was used as the training set and the CGGA RNA-seq dataset as the validation set. Furthermore, the CGGA microarray profile, as well as the GSE16011 and REMBRANDT datasets, were used to test the prognostic value of the risk model developed in the TCGA cohort.
Immune-Related Genes
The 1,039 known immunologically relevant genes were downloaded from the ImmPort database (www.immport.org), and 1,100 immune-related genes were extracted from the “IMMUNE_RESPONSE” gene set (Molecular Signatures Database V7.0, http://software.broadinstitute.org/gsea/msigdb/index.jsp). Finally, after combining these two independent gene sets, 1,763 immune-related genes were selected for the next analyses (Supplementary Table 2).
Identification of Prognostic MGMT Methylation-Related Immune Genes
Hundred and thirty-seven GBM patients with data on gene expression and MGMT methylation status was selected in the TCGA cohort, and the differentially expressed genes (DEGs) between MGMT unmethylated and methylated samples were identified using the “DESeq2” R package (Love et al., 2014). Before performing the differential gene expression analysis, low-abundance genes with raw counts <10 in more than 75% of samples were removed (Bullard et al., 2010). Then, the differentially expressed immune genes were filtered using the criterion of log2 fold change > 1 and FDR-adjusted P < 0.05. Similarly, DEGs (log2 fold change > 0.8 and adjusted P < 0.05) between different MGMT methylation status in the CGGA RNA-seq cohort were identified using the “limma” R package (Ritchie et al., 2015). The prognostic value of the identified genes in the TCGA cohort was further investigated using a multivariate Cox regression model integrating prognostic clinicopathological information, including age, Karnofsky performance score (KPS), IDH mutation, and MGMT methylation status. Then, genes with significant prognostic value (P < 0.05) were screened for further analyses.
Construction of a Risk Model Based on Immune-Related Genes
The least absolute shrinkage and selection operator (LASSO) regression model with 10-fold cross-validation was adopted to select the most useful prognostic factors for GBM patients. The optimal lambda value was estimated based on the minimum criteria. Then, the risk score for each sample was calculated using the expression levels of genes and corresponding regression coefficients obtained from the multivariable Cox regression analysis. The formula of the risk score model was as follows:
where N is the number of genes in the signature, βgenei is the regression coefficient, and expgenei represents the expression value of a specific gene in a single sample. The cutoff for classifying GBM patients into high and low-risk groups was determined using X-title 3.6.1 software (Camp et al., 2004). Kaplan–Meier survival analysis was used to evaluate the predictive value of the risk model, and time-dependent receiver operating characteristic (ROC) curves were drawn to assess the efficiency of the risk model in predicting the clinical outcomes at different times. This procedure was completed using the “timeROC” R package.
Functional Enrichment Analyses
Metascape, an easy-to-use web portal that provides a comprehensive analysis for the functional annotation of lists of genes (Zhou et al., 2019), was used to perform gene ontology (GO) term enrichment and the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis for genes of interest. Gene set enrichment analysis (GSEA) analysis was conducted using the R package “fgsea” (Sergushichev, 2016). GSEA gene set (C5 collection) was downloaded from the Molecular Signatures Database (v7.0). To obtain robust results, we performed gene-set permutations for 10,000 times during GSEA. Enriched terms with FDR-adjusted P < 0.05 were deemed as statistically significant.
Quantification of Immune-Cell Infiltration by Single-Sample Gene Set Enrichment Analysis (ssGSEA)
The estimation of overall immune cell type fractions in GBM samples was performed using the ssGSEA algorithm in the “GSVA” R package (Subramanian et al., 2005), which uses 24 different types of immune cell marker genes derived from previous research to analyze each GBM sample (Bindea et al., 2013). Markers associated with cells of the innate immune system, including natural killer (NK) cells, NK CD56dim cells, NK CD56bright cells, dendritic cells (DCs), immature DCs (iDCs), activated DCs (aDCs), neutrophils, mast cells, eosinophils, and macrophages, as well as those associated with cells of the adaptive immune system, including B, T central memory (Tcm), CD8+ T, T effector memory (Tem), T follicular helper (Tfh), Tγδ, Th1, Th2, Th17, and Treg cells, were included in the gene list. A numeric matrix consists of the enrichment score of each immune cell type in a single sample was obtained via ssGSEA. For the specific macrophage subtype (M2 macrophage), the xCell algorithm, a novel gene signature-based method, was performed to infer the contents (Aran et al., 2017).
Prediction of Clinical Responses to Immune Checkpoint Blockade
We used the computational method of Tumor Immune Dysfunction and Exclusion (TIDE) (Jiang et al., 2018) to predict the likelihood of response to immune checkpoint blockade therapy for GBM patients. Due to a lack of sufficient data from normal brain tissues in the TCGA dataset, the TCGA TARGET GTEx dataset from UCSC Toil RNA-seq Recompute project (Vivian et al., 2017), which includes a conjoint analysis of three independent datasets in order to make the gene expression of samples from different datasets or batches comparable, was adopted for subsequent studies. The gene expression values of GBM patients were normalized toward the corresponding average values of normal brain tissues [cortex, frontal cortex (Ba9), and anterior cingulate cortex (Ba24)] from the GTEx project before subjecting the data to the TIDE algorithm. Next, the subclass mapping method (SubMap) (Hoshida et al., 2007) was used to compare the gene expression matrices of different risk groups with the expression profile of melanoma patients who received immunotherapy targeting CTLA-4 and PD-1 (Roh et al., 2017). This step was performed using the default parameters in the SubMap module on the GenePattern website (http://genepattern.broadinstitute.org/).
Construction and Validation of a Nomogram
Multivariable Cox proportional hazards regression analysis was applied with the following clinical-relevant covariates: gender, age, MGMT methylation status, Karnofsky Performance Status (KPS) score, IDH mutation status, and risk score. A combined nomogram was generated as a quantitative tool for predicting the likelihood to die of each patient using the “regplot” R package. The concordance index (C-index) was calculated to assess the consistency between model prediction and actual clinical outcomes of patients. The calibration plot was generated to evaluate the accuracy of the prediction for 1-, 2-, and 3-year overall survival using this nomogram by the “rms” R package.
Clinical Specimens
Archival paraffin-embedded GBM tissues were collected from 15 patients who underwent surgery at the First Affiliated Hospital of Nanjing Medical University. The diagnosis of GBM was confirmed by two experienced pathologists. This study was approved by the ethics committee of the First Affiliated Hospital of Nanjing Medical University. All samples were collected under protocols approved by the institutional review boards of Nanjing Medical University, and all donors provided informed consent.
The methylation status of the MGMT promoter of these samples was assessed using the pyrosequencing (PSQ) assay. Briefly, genomic DNA was extracted from the tumor tissues using the QIAamp DNA Mini Kit according to the manufacturer's instructions (Qiagen, Germany), and then DNA was treated with the EpiTect bisulfite kit to convert cytosine to uracil (Qiagen, Germany). Bisulfite-treated DNA was amplified using the specific PSQ primer 5′-YGTTTTGYGTTTYGAYGTTYGTAGGTTTTYGYGGTGYGTA-3′ and then treated with the PyroMark Q24 MGMT kit (Qiagen, Germany). The methylation status was determined by the mean value of the methylated alleles percentage, and the value over 10% was regarded as MGMT methylation (Xie et al., 2015).
Cell Culture
Human monocyte THP-1 cells and glioblastoma cell lines U87 and U251 were purchased from the National Collection of Authenticated Cell Cultures (Shanghai, China). U87 and U251 cells were grown in Dulbecco's modified Eagle's medium (DMEM; Gibco, United States) supplemented with 10% fetal bovine serum (FBS; Gibco, United States), 2 mM L-glutamine, and 1% penicillin–streptomycin. THP-1 cells were maintained in RPMI-1640 (Gibco, United States) supplemented with 10% FBS. All cells were cultured with 5% CO2 at 37°C in a humidified atmosphere.
Macrophage Induction and M2 Phenotype Polarization
THP-1 cells were treated with 100 ng/mL phorbol 12-myristate 13-acetate (PMA; Sigma-Aldrich, United States) for 24 h to simulate a macrophage-like (M0) phenotype. Following treatment, the M0 macrophages were polarized into M2 macrophages by incubating with 20 ng/mL IL-4 (PeroTech, United States) for 24 h.
Quantitative Real-Time Polymerase Chain Reaction (qRT-PCR)
Total RNA from cells was extracted using TRIzol reagent (Invitrogen, United States), and complementary DNA (cDNA) was synthesized using a reverse transcription-PCR Kit (Roche, Switzerland). RT-PCR was performed using SYBR qPCR Master Mix (Vazyme, China) and operated with LightCycler 480 System (Roche, Switzerland). The primer sequences used in this study were listed as follows: GAPDH primers, forward: 5′-AGAACATCATCCCTGCCTCTACTG-3′, reverse: 5′-ACGCCTGCTTCACCACCTTC-3′; CD68 primers, forward: 5′-GGAAATGCCACGGTTCATCCA-3′, reverse: 5′-TGGGGTTCAGTACAGAGATGC-3′; CD14 primers, forward: 5′-ACGCCAGAACCTTGTGAGC-3′, reverse: 5′-GCATGGATCTCCACCTCTACTG-3′; CD206 primers, forward: 5′-TCCGGGTGCTGTTCTCCTA-3′, reverse: 5′-CCAGTCTGTTTTTGATGGCACT-3′; CD163 primers, forward: 5′-TTTGTCAACTTGAGTCCCTTCAC-3′, reverse: 5′-TCCCGCTACACTTGTTTTCAC-3′; IL-10 primers, forward: 5′-GACTTTAAGGGTTACCTGGGTTG-3′, reverse: 5′-TCACATGCGCCTTGATGTCTG-3′; VDR primers, forward: 5′-GTGGACATCGGCATGATGAAG-3′, reverse: 5′-GGTCGTAGGTCTTATGGTGGG-3′; GATA3 primers, forward: 5′-GCGGGCTCTATCACAAAATGA-3′, reverse: 5′-GCCTTCGCTTGGGCTTAAT-3′; TNFSF9 primers, forward: 5′-GGCTGGAGTCTACTATGTCTTCT-3′, reverse: 5′-ACCTCGGTGAAGGGAGTCC-3′; TNFRSF9 primers, forward: 5′-AGCTGTTACAACATAGTAGCCAC-3′, reverse: 5′-GGACAGGGACTGCAAATCTGAT-3′; LILRA5 primers, forward: 5′-TCACGGCTGAGATTCGACAG-3′, reverse: 5′-CCTGCGAGAGCCATAGCATC-3′.
RNA Interference
The VDR, GATA3, TNFSF9, TNFRSF9 siRNA reagents, and a pool of three LILRA5-targeting siRNAs (GenePharma, China), were transfected into U87, U251 cells, and polarized M2 macrophage at the concentration of 50 nM with Lipofectamine 3000 reagent (Invitrogen, United States). The target sequences were as follows: GATA3, 5′-GGGCUCUACUACAAGCUUCTT-3′; VDR, 5′-CCCACCUGGCUGAUCUUGUCAGUUA-3′; TNFSF9, 5′-UAUUCCGACCUCGGUGAAGGG-3′; TNFRSF9, 5′-AAGCAGTTACTACAAGGATCC-3′; LILRA5-siRNA1, 5′-GUCCUUGGGAUUCUGAUAUTT-3′; LILRA5-siRNA2, 5′-GGAAUACCGUCUGGUUAAATT-3′; LILRA5-siRNA3, 5′-GUGACAGGAUUCUACAACATT-3′. Cells were harvested for the next assay after incubated for 48 h.
Immunofluorescence (IF) and Immunohistochemistry (IHC)
The polarized M2 macrophages were treated with VDR, GATA3, TNFSF9, TNFRSF9, and LILRA5 siRNA for 48 h, blocked with 5% BSA and 0.01% Triton-X 100 before incubation with primary antibodies. The cells were treated with primary antibodies at 4°C overnight and then incubated with fluorescence-conjugated secondary antibodies. DAPI (Invitrogen) was applied to stain the nuclei before capturing by a fluorescence confocal microscope (LSM-710; CarlZeiss Microscopy, Germany). The primary antibodies used in the IF assay were listed as follows: CD68 (Abcam, United Kingdom, 1:100 dilution), CD163 (Proteintech, China, 1:200), CD206 (Proteintech, China, 1:200). For IHC staining, the assay was conducted by following the manufacturer's protocol as previously published (Zhao et al., 2020a). Antibodies used in the IHC assay were list as follows: GATA3 (Proteintech, China, 1:100), VDR (Santa Cruz, United States, 1:100), TNFSF9 (R&D Systems, United States, 10 μg/mL), TNFRSF9 (R&D Systems, United States, 15 μg/mL), LILRA5 (Biorbyt, United Kingdom, 1:100). IHC staining of GBM tissues for individual target protein was assessed by measuring the ratio of the integrated optical density (IOD) using Image-pro plus software (version 6.0). Images were captured from six random fields, and the final score was set as an average per field.
Enzyme-Linked Immunosorbent Assay (ELISA)
U87 and U251 cells were seeded in six-well plates and treated with siRNAs-targeting GATA3, VDR, TNFSF9, TNFRSF9, and LILRA5 for 48 h. The culture supernatants were collected, and then the protein expression levels of VEGFA, TGF-β, and CCL2 were detected using the ELISA kits. All the ELISA kits were purchased from Proteintech. The procedures were performed according to the supplier's protocol.
Statistical Analysis
The LASSO Cox regression analysis was carried out using the “glmnet” R package. Restricted mean survival (RMS) represents the loss in average life expectancy for patients, and the RMS, as well as time ratio for the risk model in each dataset, were calculated using the “survRM” R package. Tumor purity, which represents the heterogeneity of each tumor sample, was estimated by the “ESTIMATE” R package (Yoshihara et al., 2013). Univariate and multivariate Cox regression analyses were used to identify the risk-score-based signature as an independent prognostic factor for GBM patients. R software (version 3.6.1, www.r-project.org) was used for all statistical analyses.
Fisher's exact test was used to compare the differences of clinicopathologic features as well as immunotherapy response between high- and low-risk groups, and was implemented using the “fisher.test” function in R software. Mann–Whitney U test was used to assess the statistical significance of differences between the means of continuous data. All hypothetical tests were two-sided, and P values <0.05 were considered to indicate statistical significance.
Results
Association Between MGMT Methylation and the Immunological Phenotype of GBM
There is a large body of work implicating that MGMT methylation can affect the prognosis as well as the effectiveness of chemotherapy in GBM. However, the influence of MGMT promoter status on the immunological phenotype of GBM had, to our best knowledge, not been investigated prior to this study. Here, we utilized the RNA-seq gene expression profiles and clinical information of GBM patients from both the TCGA and CGGA RNA-seq databases to explore the molecular differences in underlying immunological mechanisms according to the MGMT methylation status. Patients were divided into MGMTU (MGMT-unmethylated) and MGMTM (MGMT methylated) groups in these two cohorts, respectively. GSEA analysis showed that immune-related biological processes, such as humoral immune response [normalized enrichment score (NES) = 1.95, FDR < 0.001], T cell activation involved in immune response (NES = 1.75, FDR < 0.001), adaptive immune response (NES = 1.72, FDR < 0.001), cytokine production involved in immune response (NES = 1.65, FDR <0.001), and regulation of innate immune response (NES = 1.30, FDR < 0.01), were significantly enriched in the MGMTU group in the TCGA cohort (Figure 1A). Notably, all immune-related GO terms with NES > 1 were only enriched in the MGMTU group in the TCGA cohort (Supplementary Table 3). In another independent GBM dataset (CGGA RNA-seq cohort), immune-related biological processes were also highly enriched in the MGMTU samples (Figure 1B, Supplementary Table 3). Top-enriched gene terms were all implicated in regulating immune activities in these two GBM datasets. GO and KEGG analyses showed that MGMT-unmethylation-specific upregulation of genes associated with the immune-related biological processes and signaling pathways (Figures 1C,D). Genes with roles in the activation and migration of different immune cells, cytokine and chemokine signaling axes, and IL-17 pathway were characterized by enrichment of genes associated with immunological functions. These findings suggested that the MGMT promoter methylation status may affect the prognosis of GBM patients by mediating the immune response in the malignant tumor environment.
Figure 1. Functional annotation of MGMT methylation status in the TCGA and CGGA RNA-seq datasets. Gene set enrichment analyses showing genes upregulated in MGMT-unmethylated samples evaluated in the gene sets representative for immune-related biological processes and pathways in both TCGA (A) and CGGA RNA-seq (B) GBM cohorts. Normalized enrichment score (NES) and FDR were generated for each term. GO and KEGG analyses showing the enriched biological processes and signaling pathways in MGMT-unmethylated samples of the TCGA (C) and CGGA RNA-seq (D) datasets.
Construction of an Immune Risk Model Based on MGMT Methylation Status
To further investigate the potential immune pathways and key immune-related genes driving this phenomenon, differentially expressed genes between the MGMTU and MGMTM groups in the TCGA database were identified (Figure 2A, Supplementary Table 4). Sixty-five immune-related genes were found to be more highly expressed in the MGMTU than the MGMTM group (log2 fold change > 1, FDR-adjusted P < 0.05), such as CXC motif chemokine ligand (CXCL1, 2, 6, 12, 13), immunoglobulin kappa variable family (IGKV1-27, IGKV2-28, IGKV3-11, and IGKV3-15), and immunoglobulin heavy constant gamma cluster (IGHG1, 2, 3).
Figure 2. Construction and prognostic analysis of the immune gene-based signature. (A) Volcano plot of fold changes of genes in groups with different MGMT methylation status in the TCGA dataset. Overexpressed genes in MGMT-unmethylated group were plotted as red dots while blue dots represent downregulated genes. The x- and y-axis represent fold-change and the P-values estimated by DESeq2. (B) Results of multivariate Cox regression analysis using differentially expressed immune genes in different MGMT methylation status. Genes with p-values < 0.05 were obtained. (C) The partial likelihood deviance was plotted using vertical lines with the red dot, and the dotted vertical lines represent values based on minimum criteria and 1-SE criteria, respectively. (D) LASSO coefficient profiles of the candidate prognostic immune genes by 10-fold cross-validation. (E) The relationship between the expression patterns of five candidate immune genes and the increasing risk score, as well as other clinicopathological features. Kaplan–Meier analysis of overall survival based on high vs. low risk in different cohorts (F–J). The prognostic value of this model was further validated in subgroups with MGMT methylation (K) and unmethylation (L) from the TCGA dataset. P values were obtained from the log-rank test. Time-dependent ROC curves analysis of this signature in the TCGA (M), CGGA RNA-seq (N), CGGA microarray (O), GSE16011 (P), and REMBRANDT (Q) datasets.
The prognostic effect of the differentially expressed immune genes was investigated by applying multivariate Cox regression adjusted for other survival-related clinicopathological covariates, including gender, age, MGMT methylation status, KPS score, and IDH mutation status. TNF superfamily member 9 (TNFSF9), TNF receptor superfamily member 9 (TNFRSF9), interleukin 1 receptor type 1 (IL1R1), vitamin D receptor (VDR), CD70, leukocyte immunoglobulin-like receptor A5 (LILRA5), GATA-binding protein 3 (GATA3), and CXC motif chemokine ligand 13 (CXCL13) were selected as independent factors associated with the overall survival of GBM patients (Figure 2B, Supplementary Table 5). All these immune genes were regarded as risk factors with a hazard ratio (HR) > 1. Among these candidates, LASSO Cox regression identified five highly relevant immune-related genes: LILRA5, TNFSF9, GATA3, VDR, and TNFRSF9 (Figures 2C,D). The robustness of the result was confirmed using 10-fold cross-validation. Circos plots (Gu et al., 2014) were generated to delineate the correlations between these genes in both the TCGA and CGGA RNA-seq datasets (Supplementary Figures 1A,B). It is worth noting that the relationships of these genes were quite similar in different datasets indicating that the identified relevant immune genes had a stable positive correlation in GBM.
A schematic view of MGMT methylation related immune gene selection and prognostic gene signature development is delineated in Supplementary Figure 2. The risk score of each GBM patient was calculated by combining the expression level of genes and corresponding coefficients derived from the multivariate Cox regression analysis. Patients were split into high- and low-risk populations according to the optimal cutoff value estimated using X-title software. The relationship between the expression patterns of these five prognostic immune genes and the increasing risk score in the TCGA cohort were shown (Figure 2E). More older and IDH wild-type patients were found in the high-risk group compared with those with low-risk score, while no difference in the proportion of MGMT methylation between these groups was found. Kaplan–Meier survival analysis further showed that patients from the high-risk group had shorter overall survival than those with low risk in the TCGA cohort [HR = 2.6, 95% confidence interval (CI): 1.75–3.87, P < 0.001] (Figure 2F). We also applied this computational formula to the CGGA RNA-seq, CGGA microarray, GSE16011, and REMBRANDT datasets to explore the prognostic value of the identified immune-genetic signature. Consistent with the results of the TCGA cohort, high-risk patients had significantly worse outcomes than those with low risk (Figures 2G–J). Significant RMS time ratios ranging from 1.465 to 1.854 were observed in the five datasets (P < 0.05, Supplementary Table 6). In the TCGA dataset, we also checked the efficiency of this model in patients with or without MGMT methylation (Figures 2K,L). The results showed that low-risk patients had a significant survival advantage over the high-risk group in these two counterparts (P < 0.05). Time-dependent ROC curves were generated to evaluate the efficiency of this risk signature in predicting 1-, 2-, and 3-year survival of GBM patients, and the risk model showed considerable predictive potential in different cohorts (Figures 2M,Q).
To confirm whether the risk model can be used as an independent prognostic tool for GBM, univariate and multivariate Cox regression analyses were applied to the TCGA dataset. The risk score signature was adjusted using prognostic information of clinical factors that had been deemed statistically significant in univariate analysis (P < 0.05). The HRs for the signature in the univariate and multivariate analyses were 3.107 (P < 0.001, 95% CI: 2.005–4.816) and 2.476 (P < 0.001, 95% CI: 1.581–3.878), respectively. The risk model was further validated in the CGGA RNA-seq cohort, with HRs of 1.896 (P = 0.003, 95% CI: 1.246–2.886) and 1.644 (P = 0.03, 95% CI: 1.05–2.573), respectively (Table 1). Thus, the gene signature based on immune-related molecules can be used independently to predict the overall survival of GBM patients.
Table 1. Univariate and multivariate analyses of clinicopathological characteristics and immune gene-based signature with overall survival in GBM cohorts.
Functional Analysis of the MGMT Methylation-Based Signature
To further identify the potential biological processes and signaling pathways associated with the immune gene signature, GO and KEGG functional analyses were performed, and the results were graphed using Metascape online tool. Firstly, we selected a list of genes that were strongly correlated with the risk score [|Pearson correlation coefficient| (|R|) > 0.4, P < 0.05] in the TCGA and CGGA RNA-seq cohorts (Supplementary Table 7). Heatmaps for these genes and corresponding detailed clinicopathological features were plotted, as depicted in Supplementary Figures 3A,B. Functional analysis demonstrated that these genes were significantly enriched in immune processes, such as lymphocyte activation, cytokine production, cytokine-mediated signaling pathway, leukocyte migration, activation of the immune response, regulation of the inflammatory response, and interferon-gamma production (P < 0.05) (Figures 3A,B). Similarly, immune-related terms such as cytokine-mediated signaling pathway, myeloid leukocyte activation, cytokine production, leukocyte migration, TNF signaling pathway, and NF-kappa B signaling pathway were also enriched in genes derived from the CGGA RNA-seq cohort (P < 0.05) (Supplementary Figures 4A,B).
Figure 3. Biological functions of the immune gene signature in GBM in the TCGA dataset. (A) Bar plot showing the top 20 terms derived from the gene set enrichment of risk score. The x-axis represents statistical significance. (B) The enrichment network plot visualizing the relationship between a set of representative terms. Each term is assigned with a unique color.
Association of the Risk Score With Tumor Purity and Immune Cell Infiltration in the GBM Microenvironment
The intratumoral niche is a complex microenvironment that consists of the tumor and non-tumor cells, such as stromal and immune cells (Golebiewska et al., 2013). It is now clear that a mass of non-tumor cells plays an important role in the initiation and progression of various human cancers (Mao et al., 2013; Bussard et al., 2016; Wang M. et al., 2017). Tumor purity is defined as the proportion of cancer cells in the mixed tumor tissues (Zhang et al., 2015) and has been identified as an effective prognostic index in glioma (Zhang et al., 2017). To clarify the relationship between the immune gene-based risk score and the tumor purity, we used the ESTIMATE algorithm to calculate the purity of GBM samples based on the transcriptomic gene expression data. In both the TCGA and CGGA RNA-seq cohorts, the risk score was negatively correlated with the tumor purity (R = −0.6, P < 0.0001 in TCGA, and R = −0.4, P < 0.0001 in CGGA RNA-seq) (Figures 4A,B).
Figure 4. The risk signature is associated with tumor purity and immune infiltration in the microenvironment of GBM. Correlation of risk score and tumor purity in GBM samples in the TCGA (A) and CGGA RNA-seq (B) cohorts. The coefficient and P-value were obtained from Pearson correlation analysis. Heatmaps delineate the relationship between risk score and the contents of 24 types of immune cells infiltrated in the tumor tissues in the TCGA (C) and CGGA RNA-seq (D) cohorts. Box plots showing significantly different immune cells among high- and low-risk patients in the TCGA (E) and CGGA RNA-seq (F) datasets. *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001, with Mann–Whitney U test.
We then used the ssGSEA algorithm to evaluate the immune cell populations in the tumors from the two independent datasets. The contents of different types of immune cells varied as the risk score increased. The results showed that samples from the high-risk group in the TCGA cohort exhibited a higher abundance of pro-tumor immune cells, such as plasmacytoid dendritic cells (pDCs), immature dendritic cells (iDCs), CD56dim natural killer cells, macrophages, neutrophils, regulatory T cells (Tregs), and mast cells (P < 0.05) (Figures 4C,E). At the same time, B cells, which indicate a favorable prognosis of GBM patients (Li et al., 2016), were highly abundant in the tumor samples from the low-risk group (P < 0.05). Similarly, high-risk patients in the CGGA RNA-seq cohort had significantly higher quantities of pro-tumor immune cells, except for pDCs and mast cells (P < 0.05) (Figures 4D,F). Therefore, we concluded that the immune-related biological processes and pathways associated with the risk score might result from the observed significant differences in differentiation or recruitment of various immune cell types.
The Risk Score Is Correlated With Immunosuppressive Processes and the Inflammatory Responses
To investigate the potential mechanisms underlying the correlation between the genes related to the immune response and the risk score, we assessed the relationship between the risk score and glioma-associated immunosuppressors. The gene list of glioma-related immunosuppressors, including immunosuppressive cytokines and checkpoints, tumor-supportive macrophage chemotactic and skewing molecules, immunosuppressive signaling pathways, and immunosuppressors were obtained from the published literature (Doucette et al., 2013; Cai et al., 2018). The expression of immunosuppressive genes in tumors from the two GBM cohorts was plotted against the increasing risk score (Supplementary Figures 5A,B). Following ssGSEA analysis, the score of each immunosuppression-related term was calculated for each tumor sample. Interestingly, the risk score was positively correlated with the identified immunosuppression-related terms in both the TCGA and CGGA RNA-seq cohorts (Figures 5A,B). This finding further corroborated that the immune gene signature is related to immunological suppression in GBM via changes in the expression of immunosuppressors.
Figure 5. The risk signature is associated with immunosuppression and inflammatory activities in GBM. Scatterplots showing the correlation between risk score and each immunosuppressive activity in the TCGA (A) and CGGA RNA-seq (B) cohorts. The coefficients and P-values were calculated from Pearson correlation analyses. Corrplots showing the correlation between risk score and seven types of inflammatory activities in the TCGA (C) and CGGA RNA-seq (D) cohorts. The black cross indicates no statistical significance. Corrplots showing the correlation between risk score and the expression levels of several immune checkpoint molecules in the TCGA (E) and CGGA RNA-seq (F) cohorts. The black cross indicates no statistical significance.
As demonstrated above, genes correlated with the risk score were found to be significantly enriched in the category inflammatory response. Metagenes, i.e., expression patterns of various inflammation-related genes, were used to depict the relationship between the risk score and inflammatory activities in GBM. Except for interferon, most metagenes (HCK, IgG, LCK, MHC-I, MHC-II, and STAT1) were found to be positively associated with the risk score in both the TCGA and CGGA RNA-seq cohorts (P < 0.05) (Figures 5C,D, Supplementary Figures 6A,B). These results revealed that the risk score was intimately tied to the B- and T-cell mediated immune response, as well as macrophage activation, but not to genes regulated by interferon during the remodeling of the immune microenvironment.
Pioneering investigations revealed that immunotherapy targeting immune checkpoints offered great hope for the clinical treatment of human cancers (Sharma and Allison, 2015). We investigated the correlation between the risk score and the expression levels of several well-known immune checkpoints, both in the TCGA and CGGA RNA-seq datasets. Following Pearson correlation analysis, the risk score was found to be tightly associated with the expression of immune checkpoint molecules, including TIM-3, TIGIT (except in the CGGA RNA-seq cohort), LAG-3, PD-1, CTLA-4, CD80, PD-L1, PD-L2, B7-H3, and IDO1 (P < 0.05; Figures 5E,F). In conclusion, these findings implied a potential correlation of the risk score with the regulation of immunosuppression in GBM.
High-Risk Patients Are Likely to Be More Sensitive to Immunotherapy, Chemotherapy, and Radiotherapy
Given that the immune checkpoint blockade has not been used as a routine treatment for GBM, we modeled the potential response of patients in the different risk groups to immunotherapy, including antibodies against PD-1 and CTLA-4. Following the TIDE instructions, normal specimens from the GTEx database were used to adjust the gene expression matrix of GBM. The analysis indicated that the high-risk group (53.164%, 42/79) was significantly more likely to benefit from treatment with immune checkpoint inhibitors than the low-risk group (20.690%, 12/58) (P < 0.001) (Figure 6A, Supplementary Table 8). SubMap is an unsupervised clustering method that can match subclasses in two independent gene expression data sets (Hoshida et al., 2007). Here, SubMap analysis further indicated that patients in the high-risk group are more likely to respond to anti-PD-1 immunotherapy (Bonferroni corrected P = 0.024) (Figure 6B).
Figure 6. The risk signature may be tightly associated with the responses of immunotherapy, chemotherapy, and radiotherapy. (A) Comparison of the predicted immunotherapeutic responses in different risk groups using the TIDE algorithm. P-value was obtained from Fisher's exact test. (B) The SubMap analysis concluded that high-risk patients might be more sensitive to anti-PD-1 immunotherapy. Patients from the high-risk group can benefit more from TMZ-based chemotherapy (C,D) and radiotherapy (G,H) in the TCGA GBM cohort, compared with those of the low-risk group. Similar responses to chemotherapy and radiotherapy were found in the CGGA RNA-seq GBM cohort (E,F,I,J).
There is a large body of work implicating the influence of MGMT methylation status on the chemosensitivity to TMZ in GBM patients. Since the immune-gene-based signature was built based on the MGMT methylation status, we investigated whether this risk model can be used as an effective indicator for stratifying patients into different subgroups which may show distinct sensitivities to TMZ treatment. In both the TCGA and CGGA RNA-seq cohorts, whole GBM samples were divided into two groups according to their TMZ-based chemotherapy history, and then survival analyses for patients with high and low risk were carried out. In the TCGA cohort, patients who received TMZ treatment showed significantly better prognosis than those without chemotherapy in the high-risk group (HR = 0.47, 95% CI: 0.26–0.85, P = 0.003) (Figure 6C), while no significant difference was found in the low-risk group (P > 0.05) (Figure 6D). As we mentioned above, no difference in the proportion of MGMT methylation status between high- and low-risk groups suggested that the immune gene signature may serve as an effective biomarker for predicting the response of chemotherapy independent of MGMT methylation. The similar finding was also observed in the CGGA RNA-seq cohort that patients showed satisfactory responses to TMZ treatment in the high-risk group (HR = 0.51, 95% CI: 0.25–1.04, P = 0.027) (Figure 6E), while no difference was found between TMZ treated and non-TMZ treated subgroups in the low-risk group (P > 0.05) (Figure 6F). As for the role of this immune gene signature in regulating the sensitivity to radiotherapy in GBM patients, the same methods were performed as we described above. Patients with high risk can benefit more from radiotherapy both in the TCGA (HR = 0.4, 95% CI: 0.2–0.82, P < 0.001) and CGGA RNA-seq (HR = 0.21, 95% CI: 0.05–0.98, P < 0.001) cohorts (Figures 6G,I), while no difference existed in the low-risk group (Figures 6H,J). These findings indicated that the five-gene signature might serve as an effective biomarker for directing the applications of immunotherapy as well as routine treatment, including chemotherapy and radiotherapy.
Development of a Nomogram Based on MGMT Methylation-Related Signature
To develop a quantitative tool for predicting the prognosis of GBM patients in the TCGA cohort, we established a nomogram by integrating clinicopathological risk factors and immune gene signature based on the multivariable Cox proportional hazards model (Figure 7A). The point scale in the nomogram was utilized to generate point to these variables, and the risk of death of each GBM patient was qualified by accumulating total points of all variables. The risk score was found to have the most excellent weight among all these variables, which was consistent with the result of the previous multivariable Cox regression analysis. The C-index of this nomogram reached 0.704 (95% CI: 0.672–0.736). The development of the calibration plots further confirmed the significant consistency between predicted and observed actual clinical outcomes of GBM patients (Figure 7B). We additionally constructed a nomogram for the CGGA GBM cohort based on the risk model and other available clinicopathological features, and the feasibility of this nomogram was confirmed using the calibration plots (Figures 7C,D).
Figure 7. Developed nomogram to predict the probability of survival in GBM patients. Nomogram built with clinicopathological factors incorporated estimating 1-, 2- and 3-year overall survival for GBM patients in the TCGA (A) and CGGA RNA-seq (C) cohorts. The calibration curves describing the consistency between predicted and observed overall survival at different time points in the TCGA (B) and CGGA RNA-seq (D) cohorts. The estimated survival was plotted on the x-axes, and the actual outcomes were plotted on the y-axes. The gray 45-degree dotted line represents an ideal calibration mode.
Validation of the Five Key Genes in Clinical Specimens
In silico analysis revealed that all these five key genes in the risk signature were more highly expressed in the MGMTU GBM samples. Here, we further validated the protein expression levels of these genes in our clinical specimens. Fifteen GBM samples were divided into MGMTM (n = 8) and MGMTU (n = 7) groups according to the MGMT promoter methylation status. The detailed PSQ results of the samples were illustrated in Supplementary Figure 7. All these proteins were significantly upregulated in MGMTU groups compared with the MGMTM samples (Figures 8A,B).
Figure 8. The protein expression levels of the genes in the risk signature were higher in MGMT-unmethylated GBM tissues. (A) Representative photographs of IHC staining for GATA3, VDR, TNFSF9, TNFRSF9, and LILRA5 in clinical GBM samples with MGMT methylation and unmethylation. Scale bar = 50 μm. (B) Boxplots showing the quantification of IHC staining intensity for each protein in the tumor samples. U represents MGMT-unmethylation, and M represents MGMT-methylation. *p < 0.05, **p < 0.01.
Genes in the Signature Promote Immunosuppressive Microenvironment by Facilitating M2 Macrophage Polarization and Immunosuppressive Cytokines Production
The M2-type macrophage has been proved to promote the progression of glioma by reducing antitumor functions and inducing the immunosuppressive mediators (Kennedy et al., 2013). We further explored the content of M2 macrophages in the GBM tumors using xCell analysis. M2 macrophage score of the high-risk group was significantly higher compared with the low-risk group in TCGA-GBM cohort (Figure 9A). Meanwhile, we found a positive correlation between the risk scores and M2 macrophage scores (R = 0.4, P < 0.001) (Figure 9B). These results suggested that the genes in the risk signature may mediate the immunosuppression of the tumor microenvironment by directly regulating the M2 macrophage polarization. The expression levels of VDR, TNFSF9, TNFRSF9, and LILRA5 were also positively correlated with the M2 macrophage scores (Supplementary Figure 8).
Figure 9. The immune-gene based signature promotes the polarization of M2 macrophage in the tumor microenvironment of GBM. (A) A higher proportion of M2 macrophages was found in the high-risk group compared with the low-risk group. (B) The risk scores of TCGA-GBM samples were positively correlated with the contents of the M2 macrophage. The coefficient and P-value were calculated from the Pearson correlation analysis. (C) The morphological change implies the THP-1 cells differentiated into M0 macrophages after treated with PMA for 24 h. Magnification, ×100, and ×200. (D) The mRNA expression of macrophage markers (CD68, CD14) increased after PMA treatment. (E) The mRNA expression of M2 macrophage markers (CD163, CD206, IL-10, and TGF-β) increased after the induced macrophage treated with IL-4 for 24 h. (F) ELISA showing the level changes of the secreted immunosuppressive proteins (CCL2, VEGFA, and TGF-β) in the supernatant from GBM cells (U87, U251) after knocking down GATA3, VDR, TNFSF9, TNFRSF9, and LILRA5. (G) Immunofluorescence confocal microscopy of CD68, CD206, and CD163 in M2 macrophages. Red fluorescence indicates CD206 and CD163, green fluorescence indicates CD68, and blue fluorescence indicates DAPI staining of nuclear DNA. Scale bar = 50 μm. The data were presented as the mean ± SD; *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001.
Human THP-1 cells were differentiated into M0 macrophages by incubating with PMA reagent. The induced macrophage was confirmed by the cell morphology and the expression of macrophage markers, such as CD68 and CD14 (Figures 9C,D). Then, M0 macrophages were polarized into M2 macrophages using IL-4, and this process was evaluated by detecting the expression changes of M2 macrophage markers, including CD206, CD163, IL-10, and TGF-β. The M2 markers expression increased when macrophages were treated with IL-4 compared with the control group (Figure 9E). The silencing efficiency of the siRNAs for the five genes was determined by qRT-PCR (Supplementary Figure 9). The ELISA revealed that the immunosuppressive cytokines secreted by siRNA-treated GBM cells decreased when compared with untreated cells (Figure 9F). The IF assays also showed that macrophages cultured with siRNAs expressed less M2-type markers (CD163 and CD206) than those of the control group (Figure 9G). These findings demonstrated that the genes in the risk signature promote immunosuppressive microenvironment by facilitating M2 macrophage polarization and immunosuppressive cytokines production.
Discussion
In GBM, MGMT promoter methylation status is an essential determinant of the aggressiveness of the tumor. Malignant gliomas often show resistance to alkylating drugs due to the increased expression of MGMT, which remains a significant obstacle to clinical treatment (Sarkaria et al., 2008). Previously, considerable efforts have been made to exploit the mechanisms of chemo-resistance induced by MGMT promoter demethylation. Over the past few years, the promising outcomes of immunotherapy in malignant tumors have attracted enormous attention. However, whether tumor immunity can be affected by the methylation status of the MGMT promoter has, to our best knowledge, not been explored before. Here, we sought to investigate the potential immunological mechanisms associated with the methylation status of MGMT in GBM by a comprehensive bioinformatic analysis. We addressed the critical role of MGMT methylation in regulating the remodeling of the tumor microenvironment of GBM.
The solid tumor mass of GBM is composed of both cancer cells and non-cancerous cells, such as immune cells, endothelial cells, fibroblasts, as well as some CNS-specific cell types such as astrocytes and microglia (Quail and Joyce, 2017; Wang Q. et al., 2017). The discovery that multiple immune cell types can infiltrate the blood-brain barrier and the presence of reactivity against antigens derived from the CNS has made the view of the brain as an immunologically privileged site obsolete (Weiss et al., 2009; Louveau et al., 2015a). The most abundant immune cells, macrophages, comprise up to 30% of the total cells inside a glioma tumor (Graeber et al., 2002). Tumor-associated macrophages suppress T cell responses due to insufficient expression of molecules associated with T cell co-stimulation (Hussain et al., 2006). M2 polarization of microglia in the tumor niche can remodel the GBM immunosuppressive microenvironment (Meng et al., 2019). Neutrophil infiltration and neutrophil extracellular traps formation mediate the crosstalk between glioma progression and the tumor microenvironment by promoting IL-8 secretion (Zha et al., 2020). Dendritic cells are responsible for presenting antigens to T cells to induce immune responses in the CNS. The application of dendritic cell-based vaccination in GBM patients has produced survival benefits due to enhanced stimulated T-cell responses (Prins et al., 2011). Additionally, several studies demonstrated that reprogramming immunosuppressive T-cell subsets can stimulate antitumor immunity in glioma (Fecci et al., 2006; Hussain et al., 2006). In this study, we discovered that genetic signatures of immune-related biological processes and signaling pathways were highly enriched in MGMT-unmethylated specimens. Then, we filtered differentially expressed immune-related genes according to the MGMT promoter status and generated a signature based on six genes that can predict the clinical outcomes of GBM patients. At the same time, we investigated the contents of immune cells in the bulk tumor to predict immune infiltration using a marker gene-based approach based on transcriptomic data. Notably, we found that the infiltrating immune cell types differed significantly between the high- and low-risk groups. In agreement with our hypothesis, the tumors in the high-risk group were severely infiltrated by immune cells characterized as belonging to pro-tumor subsets. Thus, this risk model provides an immunological interpretation of the effect of MGMT methylation status on the prognosis of GBM patients.
There is a large body of work that confirms the immunity-related roles of the genes identified in this study. For example, the transcription factor GATA3 has been characterized as a critical regulator in the development of T, B, NK, and innate-like lymphoid cells (Wan, 2014). GATA3 is responsible for maintaining the function of mature T and Treg cells via the IL-7 and IL-2 signaling pathways, respectively (Ma et al., 2006). Additionally, ectopic expression of GATA3 was found to promote the function of human type-2 innate lymphoid cells (Mjösberg et al., 2012). TNSRSF9, also known as CD137, is a member of the TNF-receptor family that can activate cytokine production in cytotoxic T lymphocytes (CTLs) and Th1 cells by binding to its native ligand TNFSF9 (Ascierto et al., 2010). Recently, clinical trials with antibodies targeting CD137 have been launched with the expectation to improve cancer immunotherapy (Yonezawa et al., 2015). VDR negatively regulates Forkhead box M1 (FOXM1) signaling and further suppresses the development of pancreatic ductal adenocarcinoma (Li et al., 2015). In the immune system, VDR expression regulates the proliferation, differentiation, and function of CTLs (Sarkar et al., 2016). Stromal VDR activation reduces tumor-associated fibrosis and inhibits tumor-supportive signaling events in pancreatic cancer (Sherman et al., 2014).
Immune suppression has been recognized as a hindrance to successful antitumor immunotherapy. Tumors create a microenvironment that sustains tumor growth as well as reduces adaptive immunity to tumor-derived antigens (Hurwitz and Watkins, 2012). As for glioma, several genes were proved as critical molecules to contribute to the immune heterogeneity and immunosuppression in the tumor microenvironment (Cai et al., 2018; Chen et al., 2019). We found a strong correlation between the risk score and immunosuppressive cytokines, tumor-supportive macrophages, as well as immunosuppressive signaling pathways. This finding implied that the risk score reflects the immunosuppressive microenvironment of GBM. Both in silico analysis and in vitro experiments confirmed that the five genes in the signature could promote the immunosuppression by mediating the M2 macrophage polarization and infiltration. Furthermore, monoclonal antibodies that block immune checkpoint signaling have shown clinical benefits in several malignant tumors, including melanoma, NSCLC, Hodgkin lymphoma, and urothelial carcinoma. Such antibodies include pembrolizumab, nivolumab, and pidilizumab, which target PD-1, ipilimumab and tremelimumab for CTLA-4, or atezolizumab and durvalumab for PD-L1. Additionally, many clinical trials have been launched to evaluate the application of checkpoint inhibitors in GBM. A randomized phase III study (NCT02017717) was carried out to test the efficacy and safety of nivolumab alone vs. bevacizumab (Sampson et al., 2014). Clinical treatment of newly diagnosed GBM patients with ipilimumab and nivolumab is still under investigation (NCT02311920) (Binder et al., 2016). Therefore, we explored the relationship between the risk score and the expression levels of several well-known immune checkpoint molecules. The strong positive correlation implied that the identified signature mediates immune escape mechanisms in GBM, which provides a foundation for identifying novel therapeutic targets.
It is noteworthy that not all probands will respond to immune checkpoint blockade effectively in most cancer types (Sharma et al., 2017). Various factors were reported as predictors for evaluating immunotherapy effectiveness, such as the mutation burden (Snyder et al., 2014), tumor infiltration by cytotoxic T cells (Van Allen et al., 2015), PD-L1 expression (Nishino et al., 2017), and tumor aneuploidy (Davoli et al., 2017). The known mechanisms of tumor immune evasion include infiltration by CTLs and the prevention of T cell infiltration (Gajewski et al., 2013; Joyce and Fearon, 2015). The TIDE method underscores these two distinct aspects and can predict the outcome of melanoma patients receiving immunotherapy more accurately than other factors (Jiang et al., 2018). We demonstrated that more patients could benefit from immune checkpoint blockade in the high than in the low-risk group. Under a strict criterion of adjusted P-values, SubMap analysis further confirmed that anti-PD-1 therapy might produce favorable outcomes in high-risk patients.
In summary, our study revealed that MGMT promoter methylation status is correlated with immunological processes and the remodeling of the tumor microenvironment in GBM. The developed immune-related gene signature could serve as a useful prognostic tool for GBM patients as well as for predicting which patients would benefit from immunotherapy. These findings may provide new insights into fundamental aspects of the critical role of MGMT methylation status in GBM.
Data Availability Statement
Publicly available datasets were analyzed in this study. These data can be found here: TCGA (https://portal.gdc.cancer.gov/), CGGA (http://www.cgga.org.cn/), and GEO database (https://www.ncbi.nlm.nih.gov/geo/).
Ethics Statement
The studies involving human participants were reviewed and approved by the ethics committee of the First Affiliated Hospital of Nanjing Medical University. The patients/participants provided their written informed consent to participate in this study.
Author Contributions
PZ, LZ, and JZ: conception and design. PZ: administrative support. LZ, SX, and YW: provision of study materials or patients. LZ, SX, JZ, and ZL: collection and assembly of data. LZ and JZ: data analysis and interpretation. All authors: contributed manuscript writing and final approval of manuscript.
Funding
This work was supported by the National Natural Science Foundation of China under (Grant Number 81673210).
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
All authors would like to thank all contributors to the TCGA and CGGA project, especially. This manuscript has been released as a pre-print at biorxiv (Zhao et al., 2020b).
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fcell.2021.600506/full#supplementary-material
References
Aran, D., Hu, Z., and Butte, A. J. (2017). xCell: digitally portraying the tissue cellular heterogeneity landscape. Genome Biol. 18, 1–14. doi: 10.1186/s13059-017-1349-1
Ascierto, P. A., Simeone, E., Sznol, M., Fu, Y.-X., and Melero, I. (2010). Clinical experiences with anti-CD137 and anti-PD1 therapeutic antibodies. Semin. Oncol. 37, 508–516. doi: 10.1053/j.seminoncol.2010.09.008
Bindea, G., Mlecnik, B., Tosolini, M., Kirilovsky, A., Waldner, M., Obenauf, A. C., et al. (2013). Spatiotemporal dynamics of intratumoral immune cells reveal the immune landscape in human cancer. Immunity 39, 782–795. doi: 10.1016/j.immuni.2013.10.003
Binder, D. C., Davis, A. A., and Wainwright, D. A. (2016). Immunotherapy for cancer in the central nervous system: current and future directions. Oncoimmunology 5:e1082027. doi: 10.1080/2162402X.2015.1082027
Bullard, J. H., Purdom, E., Hansen, K. D., and Dudoit, S. (2010). Evaluation of statistical methods for normalization and differential expression in mRNA-Seq experiments. BMC Bioinformatics 11, 1–13. doi: 10.1186/1471-2105-11-94
Bussard, K. M., Mutkus, L., Stumpf, K., Gomez-Manzano, C., and Marini, F. C. (2016). Tumor-associated stromal cells as key contributors to the tumor microenvironment. Breast Cancer Res. 18:84. doi: 10.1186/s13058-016-0740-2
Cai, J., Chen, Q., Cui, Y., Dong, J., Chen, M., Wu, P., et al. (2018). Immune heterogeneity and clinicopathologic characterization of IGFBP2 in 2447 glioma samples. Oncoimmunology 7:e1426516. doi: 10.1080/2162402X.2018.1426516
Camp, R. L., Dolled-Filhart, M., and Rimm, D. L. (2004). X-tile: a new bio-informatics tool for biomarker assessment and outcome-based cut-point optimization. Clin. Cancer Res. 10, 7252–7259. doi: 10.1158/1078-0432.CCR-04-0713
Chen, Q., Han, B., Meng, X., Duan, C., Yang, C., Wu, Z., et al. (2019). Immunogenomic analysis reveals LGALS1 contributes to the immune heterogeneity and immunosuppression in glioma. Int. J. Cancer 145, 517–530. doi: 10.1002/ijc.32102
Cloughesy, T. F., Mochizuki, A. Y., Orpilla, J. R., Hugo, W., Lee, A. H., Davidson, T. B., et al. (2019). Neoadjuvant anti-PD-1 immunotherapy promotes a survival benefit with intratumoral and systemic immune responses in recurrent glioblastoma. Nat. Med. 25, 477–486. doi: 10.1038/s41591-018-0337-7
Davoli, T., Uno, H., Wooten, E. C., and Elledge, S. J. (2017). Tumor aneuploidy correlates with markers of immune evasion and with reduced response to immunotherapy. Science 355:eaaf8399. doi: 10.1126/science.aaf8399
Dolan, M. E., and Pegg, A. E. (1997). O6-benzylguanine and its role in chemotherapy. Clin. Cancer Res. 3, 837–847.
Doucette, T., Rao, G., Rao, A., Shen, L., Aldape, K., Wei, J., et al. (2013). Immune heterogeneity of glioblastoma subtypes: extrapolation from the cancer genome atlas. Cancer Immunol. Res. 1, 112–122. doi: 10.1158/2326-6066.CIR-13-0028
Fecci, P. E., Sweeney, A. E., Grossi, P. M., Nair, S. K., Learn, C. A., Mitchell, D. A., et al. (2006). Systemic anti-CD25 monoclonal antibody administration safely enhances immunity in murine glioma without eliminating regulatory T cells. Clin. Cancer Res.12, 4294–4305. doi: 10.1158/1078-0432.CCR-06-0053
Fukumura, D., Kloepper, J., Amoozgar, Z., Duda, D. G., and Jain, R. K. (2018). Enhancing cancer immunotherapy using antiangiogenics: opportunities and challenges. Nat. Rev. Clin. Oncol. 15:325. doi: 10.1038/nrclinonc.2018.29
Gajewski, T. F., Schreiber, H., and Fu, Y.-X. (2013). Innate and adaptive immune cells in the tumor microenvironment. Nat. Immunol. 14:1014. doi: 10.1038/ni.2703
Golebiewska, A., Bougnaud, S., Stieber, D., Brons, N. H., Vallar, L., Hertel, F., et al. (2013). Side population in human glioblastoma is non-tumorigenic and characterizes brain endothelial cells. Brain 136, 1462–1475. doi: 10.1093/brain/awt025
Graeber, M. B., Scheithauer, B. W., and Kreutzberg, G. W. (2002). Microglia in brain tumors. Glia 40, 52–259. doi: 10.1002/glia.10147
Gu, Z., Gu, L., Eils, R., Schlesner, M., and Brors, B. (2014). circlize implements and enhances circular visualization in R. Bioinformatics 30, 2811–2812. doi: 10.1093/bioinformatics/btu393
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, 997–1003. doi: 10.1056/NEJMoa043331
Hegi, M. E., Liu, L., Herman, J. G., Stupp, R., Wick, W., Weller, M., et al. (2008). Correlation of O6-methylguanine methyltransferase (MGMT) promoter methylation with clinical outcomes in glioblastoma and clinical strategies to modulate MGMT activity. J. Clin. Oncol. 26, 4189–4199. doi: 10.1200/JCO.2007.11.5964
Heimberger, A. B., Abou-Ghazal, M., Reina-Ortiz, C., Yang, D. S., Sun, W., Qiao, W., et al. (2008). Incidence and prognostic impact of FoxP3+ regulatory T cells in human gliomas. Clin. Cancer Res. 14, 5166–5172. doi: 10.1158/1078-0432.CCR-08-0320
Hodi, F. S., O'Day, S. J., McDermott, D. F., Weber, R. W., Sosman, J. A., Haanen, J. B., et al. (2010). Improved survival with ipilimumab in patients with metastatic melanoma. N. Engl. J. Med. 363, 711–723. doi: 10.1056/NEJMoa1003466
Hoshida, Y., Brunet, J.-P., Tamayo, P., Golub, T. R., and Mesirov, J. P. (2007). Subclass mapping: identifying common subtypes in independent disease data sets. PLoS ONE 2:e1195. doi: 10.1371/journal.pone.0001195
Hurwitz, A. A., and Watkins, S. K. (2012). Immune suppression in the tumor microenvironment: a role for dendritic cell-mediated tolerization of T cells. Cancer Immunol. Immunother. 61, 289–293. doi: 10.1007/s00262-011-1181-5
Hussain, S. F., Yang, D., Suki, D., Aldape, K., Grimm, E., and Heimberger, A. B. (2006). The role of human glioma-infiltrating microglia/macrophages in mediating antitumor immune responses. Neuro-oncology 8, 261–279. doi: 10.1215/15228517-2006-008
Jiang, P., Gu, S., Pan, D., Fu, J., Sahu, A., Hu, X., et al. (2018). Signatures of T cell dysfunction and exclusion predict cancer immunotherapy response. Nat. Med. 24, 1550–1558. doi: 10.1158/2326-6074.CRICIMTEATIAACR18-B077
Joyce, J. A., and Fearon, D. T. (2015). T cell exclusion, immune privilege, and the tumor microenvironment. Science 348:74–80. doi: 10.1126/science.aaa6204
Karran, P., and Bignami, M. (1992). Self-destruction and tolerance in resistance of mammalian cells to alkylation damage. Nucleic Acids Res. 20:2933. doi: 10.1093/nar/20.12.2933
Kennedy, B. C., Showers, C. R., Anderson, D. E., Anderson, L., Canoll, P., Bruce, J. N., et al. (2013). Tumor-associated macrophages in glioma: friend or foe? J. Oncol. 2013:486912. doi: 10.1155/2013/486912
Li, B., Severson, E., Pignon, J. C., Zhao, H., Li, T., Novak, J., et al. (2016). Comprehensive analyses of tumor immunity: implications for cancer immunotherapy. Genome Biol. 17:74. doi: 10.1186/s13059-016-1028-7
Li, Z., Jia, Z., Gao, Y., Xie, D., Wei, D., Cui, J., et al. (2015). Activation of vitamin D receptor signaling downregulates the expression of nuclear FOXM1 protein and suppresses pancreatic cancer cell stemness. Clin. Cancer Res. 21, 844–853. doi: 10.1158/1078-0432.CCR-14-2437
Liu, L., and Gerson, S. L. (2006). Targeted modulation of MGMT: clinical implications. Clin. Cancer Res. 12, 328–331. doi: 10.1158/1078-0432.CCR-05-2543
Louveau, A., Harris, T. H., and Kipnis, J. (2015a). Revisiting the mechanisms of CNS immune privilege. Trends Immunol. 36, 569–577. doi: 10.1016/j.it.2015.08.006
Louveau, A., Smirnov, I., Keyes, T. J., Eccles, J. D., Rouhani, S. J., Peske, J. D., et al. (2015b). Structural and functional features of central nervous system lymphatic vessels. Nature 523, 337–341. doi: 10.1038/nature14432
Love, M. I., Huber, W., and Anders, S. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 15:550. doi: 10.1186/s13059-014-0550-8
Ma, A., Koka, R., and Burkett, P. (2006). Diverse functions of IL-2, IL-15, and IL-7 in lymphoid homeostasis. Annu. Rev. Immunol. 24, 657–679. doi: 10.1146/annurev.immunol.24.021605.090727
Mao, Y., Keller, E. T., Garfield, D. H., Shen, K., and Wang, J. (2013). Stromal cells in tumor microenvironment and breast cancer. Cancer Metastasis Rev. 32, 303–315. doi: 10.1007/s10555-012-9415-3
Meng, X., Duan, C., Pang, H., Chen, Q., Han, B., Zha, C., et al. (2019). DNA damage repair alterations modulate M2 polarization of microglia to remodel the tumor microenvironment via the p53-mediated MDK expression in glioma. EBio Med. 41, 185–199. doi: 10.1016/j.ebiom.2019.01.067
Mjösberg, J., Bernink, J., Golebski, K., Karrich, J. J., Peters, C. P., Blom, B., et al. (2012). The transcription factor GATA3 is essential for the function of human type 2 innate lymphoid cells. Immunity 37, 649–659. doi: 10.1016/j.immuni.2012.08.015
Molenaar, R. J., Verbaan, D., Lamba, S., Zanon, C., Jeuken, J. W., Boots-Sprenger, S. H., et al. (2014). The combination of IDH1 mutations and MGMT methylation status predicts survival in glioblastoma better than either IDH1 or MGMT alone. Neuro-oncology 16, 1263–1273. doi: 10.1093/neuonc/nou005
Nduom, E. K., Weller, M., and Heimberger, A. B. (2015). Immunosuppressive mechanisms in glioblastoma. Neuro-oncology 17(Suppl. 7), vii9–vii14. doi: 10.1093/neuonc/nov151
Nishino, M., Ramaiya, N. H., Hatabu, H., and Hodi, F. S. (2017). Monitoring immune-checkpoint blockade: response evaluation and biomarker development. Nat. Rev. Clin. Oncol. 14:655. doi: 10.1038/nrclinonc.2017.88
Pegg, A. E. (1990). Mammalian O6-alkylguanine-DNA alkyltransferase: regulation and importance in response to alkylating carcinogenic and therapeutic agents. Cancer Res. 50, 6119–6129.
Prins, R. M., Soto, H., Konkankit, V., Odesa, S. K., Eskin, A., Yong, W. H., et al. (2011). Gene expression profile correlates with T-cell infiltration and relative survival in glioblastoma patients vaccinated with dendritic cell immunotherapy. Clin. Cancer Res. 17, 1603–1615. doi: 10.1158/1078-0432.CCR-10-2563
Qian, Z., Li, Y., Fan, X., Zhang, C., Wang, Y., Jiang, T., et al. (2018). Molecular and clinical characterization of IDH associated immune signature in lower-grade gliomas. Oncoimmunology 7:e1434466. doi: 10.1080/2162402X.2018.1434466
Quail, D. F., and Joyce, J. A. (2017). The Microenvironmental Landscape of Brain Tumors. Cancer Cell 31, 326–341. doi: 10.1016/j.ccell.2017.02.009
Ritchie, M. E., Phipson, B., Wu, D., Hu, Y., Law, C. W., Shi, W., et al. (2015). limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 43:e47. doi: 10.1093/nar/gkv007
Rizvi, N. A., Mazières, J., Planchard, D., Stinchcombe, T. E., Dy, G. K., Antonia, S. J., et al. (2015). Activity and safety of nivolumab, an anti-PD-1 immune checkpoint inhibitor, for patients with advanced, refractory squamous non-small-cell lung cancer (CheckMate 063): a phase 2, single-arm trial. Lancet Oncol. 16, 257–265. doi: 10.1016/S1470-2045(15)70054-9
Roh, W., Chen, P.-L., Reuben, A., Spencer, C. N., Prieto, P. A., Miller, J. P., et al. (2017). Integrated molecular analysis of tumor biopsies on sequential CTLA-4 and PD-1 blockade reveals markers of response and resistance. Sci. Transl. Med. 9:eaah3560. doi: 10.1126/scitranslmed.aah3560
Sampson, J. H., Vlahovic, G., Desjardins, A., Friedman, H. S., Baehring, J. M., Hafler, D., et al. (2014). Randomized phase IIb study of nivolumab (anti-PD-1; BMS-936558, ONO-4538) alone or in combination with ipilimumab versus bevacizumab in patients (pts) with recurrent glioblastoma (GBM). Am. Soc. Clin. Oncol. 37. doi: 10.1200/jco.2014.32.15_suppl.tps2101
Sarkar, S., Hewison, M., Studzinski, G. P., Li, Y. C., and Kalia, V. (2016). Role of vitamin D in cytotoxic T lymphocyte immunity to pathogens and cancer. Crit. Rev. Clin. Lab. Sci. 53, 132–145. doi: 10.3109/10408363.2015.1094443
Sarkaria, J. N., Kitange, G. J., James, C. D., Plummer, R., Calvert, H., Weller, M., et al. (2008). Mechanisms of chemoresistance to alkylating agents in malignant glioma. Clin. Cancer Res. 14, 2900–2908. doi: 10.1158/1078-0432.CCR-07-1719
Sergushichev, A. (2016). An algorithm for fast preranked gene set enrichment analysis using cumulative statistic calculation. bioRxiv 060012. doi: 10.1101/060012
Sharma, P., and Allison, J. P. (2015). The future of immune checkpoint therapy. Science 348, 56–61. doi: 10.1126/science.aaa8172
Sharma, P., Hu-Lieskovan, S., Wargo, J. A., and Ribas, A. (2017). Primary, adaptive, and acquired resistance to cancer immunotherapy. Cell 168, 707–723. doi: 10.1016/j.cell.2017.01.017
Sherman, M. H., Ruth, T. Y., Engle, D. D., Ding, N., Atkins, A. R., Tiriac, H., et al. (2014). Vitamin D receptor-mediated stromal reprogramming suppresses pancreatitis and enhances pancreatic cancer therapy. Cell 159, 80–93. doi: 10.1016/j.cell.2014.08.007
Snyder, A., Makarov, V., Merghoub, T., Yuan, J., Zaretsky, J. M., Desrichard, A., et al. (2014). Genetic basis for clinical response to CTLA-4 blockade in melanoma. N. Engl. J. Med. 371, 2189–2199. doi: 10.1056/NEJMoa1406498
Stupp, R., Mason, W. P., van den Bent, M. J., Weller, M., Fisher, B., Taphoorn, M. J., et al. (2005). Radiotherapy plus concomitant and adjuvant temozolomide for glioblastoma. N. Engl. J. Med. 352, 987–996. doi: 10.1056/NEJMoa043330
Subramanian, A., Tamayo, P., Mootha, V. K., Mukherjee, S., Ebert, B. L., Gillette, M. A., et al. (2005). Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc. Natl. Acad. Sci. 102, 15545–15550. doi: 10.1073/pnas.0506580102
Van Allen, E. M., Miao, D., Schilling, B., Shukla, S. A., Blank, C., Zimmer, L., et al. (2015). Genomic correlates of response to CTLA-4 blockade in metastatic melanoma. Science 350, 207–211. doi: 10.1126/science.aad0095
Vivian, J., Rao, A. A., Nothaft, F. A., Ketchum, C., Armstrong, J., Novak, A., et al. (2017). Toil enables reproducible, open source, big biomedical data analyses. Nat. Biotechnol. 35, 314–316. doi: 10.1038/nbt.3772
Wan, Y. Y. (2014). GATA3: a master of many trades in immune regulation. Trends Immunol. 35, 233–242. doi: 10.1016/j.it.2014.04.002
Wang, M., Zhao, J., Zhang, L., Wei, F., Lian, Y., Wu, Y., et al. (2017). Role of tumor microenvironment in tumorigenesis. J. Cancer 8:761. doi: 10.7150/jca.17648
Wang, Q., Hu, B., Hu, X., Kim, H., Squatrito, M., Scarpace, L., et al. (2017). Tumor evolution of glioma-intrinsic gene expression subtypes associates with immunological changes in the microenvironment. Cancer Cell 32, 42–56.e6. doi: 10.1016/j.ccell.2017.06.003
Weiss, N., Miller, F., Cazaubon, S., and Couraud, P.-O. (2009). The blood-brain barrier in brain homeostasis and neurological diseases. Biochim. Biophys. Acta Biomembr. 1788, 842–857. doi: 10.1016/j.bbamem.2008.10.022
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:372. doi: 10.1038/nrneurol.2014.100
Xie, H., Tubbs, R., and Yang, B. (2015). Detection of MGMT promoter methylation in glioblastoma using pyrosequencing. Int. J. Clin. Exp. Pathol. 8, 636–642.
Yonezawa, A., Dutt, S., Chester, C., Kim, J., and Kohrt, H. E. (2015). Boosting cancer immunotherapy with anti-CD137 antibody therapy. Clin. Cancer Res. 21, 3113–3120. doi: 10.1158/1078-0432.CCR-15-0263
Yoshihara, K., Shahmoradgoli, M., Martínez, E., Vegesna, R., Kim, H., Torres-Garcia, W., et al. (2013). Inferring tumour purity and stromal and immune cell admixture from expression data. Nat. Commun. 4:2612. doi: 10.1038/ncomms3612
Zha, C., Meng, X., Li, L., Mi, S., Qian, D., Li, Z., et al. (2020). Neutrophil extracellular traps mediate the crosstalk between glioma progression and the tumor microenvironment via the HMGB1/RAGE/IL-8 axis. Cancer Biol. Med. 17:154–168. doi: 10.20892/j.issn.2095-3941.2019.0353
Zhang, C., Cheng, W., Ren, X., Wang, Z., Liu, X., Li, G., et al. (2017). Tumor purity as an underlying key factor in glioma. Clin. Cancer Res. 23, 6279–6291. doi: 10.1158/1078-0432.CCR-16-2598
Zhang, N., Wu, H.-J., Zhang, W., Wang, J., Wu, H., and Zheng, X. (2015). Predicting tumor purity from methylation microarray data. Bioinformatics 31, 3401–3405. doi: 10.1093/bioinformatics/btv370
Zhao, L., Zhang, J., Liu, Z., and Zhao, P. (2020a). Identification of biomarkers for the transition from low-grade glioma to secondary glioblastoma by an integrated bioinformatic analysis. Am. J. Transl. Res. 12:1222.
Zhao, L., Zhang, J., Xuan, S., Liu, Z., Wang, Y., and Zhao, P. (2020b). Molecular and clinicopathological characterization of a prognostic immune gene signature associated with MGMT methylation in glioblastoma. bioRxiv:2020.07.16.206318. doi: 10.1101/2020.07.16.206318
Keywords: glioblastoma, MGMT methylation, gene signature, prognosis, immune microenvironment, biomarker
Citation: Zhao L, Zhang J, Xuan S, Liu Z, Wang Y and Zhao P (2021) Molecular and Clinicopathological Characterization of a Prognostic Immune Gene Signature Associated With MGMT Methylation in Glioblastoma. Front. Cell Dev. Biol. 9:600506. doi: 10.3389/fcell.2021.600506
Received: 30 August 2020; Accepted: 08 January 2021;
Published: 05 February 2021.
Edited by:
César López-Camarillo, Universidad Autónoma de la Ciudad de México, MexicoReviewed by:
Tao Xu, Shanghai Changzheng Hospital, ChinaChuanlu Jiang, Second Affiliated Hospital of Harbin Medical University, China
Copyright © 2021 Zhao, Zhang, Xuan, Liu, Wang and Zhao. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Peng Zhao, emhhb3BlbmcmI3gwMDA0MDtuam11LmVkdS5jbg==
†These authors have contributed equally to this work