Construction of Molecular Subtypes and Related Prognostic and Immune Response Models Based on M2 Macrophages in Glioblastoma

Introduction

Glioblastoma multiforme (GBM) is the most aggressive and deadliest primary brain tumor of adults.¹ Although many treatments, including surgical resection with chemotherapy and radiotherapy, may improve the outcome, the median survival time is still only 14–16 months² and the 5-year survival rate is just 9.8%.³ GBM is a biologically heterogeneous intracranial neoplasm and it has most of the classical hallmarks of a tumor. The poor prognosis in GBM is due to significant differences between patients and the wide dissemination of cancer cells before diagnosis.⁴ Even if resections are repeatedly performed, recurrence can be predicted with certainty.⁵ Recent studies have focused on genetic and epigenetic research to discover the subtypes of GBM with different prognoses, and to guide personalized therapies targeting specific molecules or pathogenies. However, the challenge remains that GBM is still incurable, while any progress and achievements have been limited.

Over the past decade, immunotherapy has revolutionized the field of oncology and developed into one of the most effective cancer therapies for many fatal cancers, such as melanoma and colorectal cancer, which largely explains its success. The main breakthrough came with immune checkpoint blockade (ICB) therapy. Among the immune check inhibitors, monoclonal antibodies, including PD-1/PD-L1 and CTLA-4 inhibitors, have provided promising therapeutic outcomes. However, GBM has benefited little from immune checkpoint inhibitors so far. The poor results may be caused by several mechanisms, including the heterogeneous nature of GBM and the immunosuppressive status of the tumor microenvironment.

Tumor-associated macrophages (TAMs), accounting for approximately 30–50% of the bulk in GBM, play an important role in cancer maintenance and progression. TAMs have different phenotypes owing to their activation state. M0 is a precursor “pre-activation” state, which can polarize into “classical” (M1) with seemingly anti-tumor functions and “alternative” (M2) with pro-tumor effects due to chemokines and cytokines.⁶ M2 macrophages play a pivotal role in immune suppression and pro-tumor functions during GBM progression through the secretion of immunosuppressive factors.⁷ A high ratio of M2/M1 in TAMs results in poor clinical prognosis, and inhibiting the polarization of M0 to M2 could dramatically reduce the growth of glioma stem cells. Gordon et al found that PD-1 expression by TAMs had anti-phagocytic and tumor immunity inhibition effects, implying that PD-1 may also function through a direct effect on macrophages, other than modulating T cells.⁸ The underlying mechanism of M2 macrophages in GBM has great significance for estimating the prognosis of GBM patients and improving the therapeutic outcome.

According to the expression patterns, GBM tumors were originally grouped into four subtypes: neural, proneural, classical and mesenchymal.⁹ The revised classification system then eliminated the neural subtype owing to normal brain tissue contamination.¹⁰ Although many other classification systems have been developed,^11–14 no M2 macrophage-based molecular subtypes have been reported. We hypothesized that analyzing the M2 macrophage-related prognostic gene expression pattern in GBM could enable us to differentiate the subtypes efficiently. This new classification system would provide us with a new direction for GBM prognosis and response to therapies.

Methods

Dataset Acquisition from GlioVis Database

GlioVis is a web-based database for data visualization and analysis to explore brain tumor expression datasets.¹⁵ This user-friendly database offers the research community an unprecedented fast and intuitive portal to molecular profiles.¹⁵ The Cancer Genome Atlas–Glioblastoma Multiforme (TCGA-GBM) and a GBM dataset reported by Gravendeel et al,¹⁶ including corresponding clinical data, were obtained for downstream study. The relevant CEL files were downloaded and the robust multiarray average (RMA) normalized analysis was applied. Then, the probe was converted into gene symbols by the R program and the batch effects were removed before downstream analysis.¹⁷ GBM patients with complete outcome data and expression data were included in the subsequent analysis. Somatic mutation datasets of TCGA-GBM and The Cancer Genome Atlas Low-Grade Glioma (TCGA-LGG) were downloaded from TCGA, and copy number variation (CNV) data were obtained from the University of California, Santa Cruz (UCSC) Xena website (http://xena.ucsc.edu/). No additional ethical approval was applicable, owing to the data being publicly available.

Tumor Immunity Analyses and WGCNA

The immune and stromal cell expression profiles of all GBM samples were qualified using xCell tools (https://xcell.ucsf.edu). xCell is a simple web tool that calculates independent enrichment scores of 64 kinds of immune and stromal cells based on enrichment.¹⁸

Weighted gene co-expression network analysis (WGCNA) was applied to explore the correlations between immune-related genes and immune cell types by constructing significant modules. WGCNA was conducted using the “WGCNA” package for R (version 3.6.2).¹⁹ The immune-related gene list was downloaded from the immunology database and analysis portal (ImmPort) system. First, both the average linkage method and Pearson’s correlation matrices were applied to establish the co-expression similarity matrix. To enhance the matrix similarity and construct a co-expression network, appropriate β-values were selected. Next, the adjacency matrix was transformed into a topological overlap matrix (TOM) to determine immune cell connectivity. Finally, an average linkage hierarchical clustering tree was established according to TOM-based dissimilarity, and immune cells with similar expression profiles were clustered to form modules.²⁰ The associations between immune cell type and immune-related genes were evaluated by Pearson’s correlation analysis.

Discovering Subtypes of GBM

To distinguish cancer subtypes based on M2 macrophages, we repeatedly perturbed the data and partitioned the patients using a classical clustering algorithm.²¹ A range of potential cluster numbers K was explored to evaluate the stability and consistency of the clustering, determine the clustering numbers of samples and establish the consistent clustering algorithm. Two principles were applied to select the number K: first, the growth rate of the cumulative distribution function (CDF) value should be low; second, the correlation in the group should be high, while there should be no small clusters in the clustering group.²² Classification of TCGA-GBM patients into various clinically significant subtypes was performed using the “ConsensusClusterPlus” package in this process. For the purpose of elucidating the molecular mechanism of patients in the different clusters, gene set enrichment analysis (GSEA) was conducted by comparing the patients in the different clusters.

Identification of Differentially Expressed Genes in Different Clusters

The mRNA expression profiles extracted from patients in different clusters were subjected to differential expression analysis by the R software package. The differentially expressed genes (DEGs) were selected according to the criteria: P<0.05, |log2FC)| >1. The “limma” package of R was conducted to find DEGs.²³ The Gene Ontology (GO) pathway enrichment of DEGs was performed, and a P-value less than 0.05 was set to indicate a significant difference.

Establishment of the M2 Macrophage-Related Gene Signature and Prognostic Model Construction

Based on the M2 macrophage-related genes, principal component analysis (PCA) was used to construct the M2 macrophage-related prognostic model (M2score). PCA1, PCA2, PCA3, PCA4 and PCA5 were extracted to serve as the M2score, since PCA1–PCA5 made independent contributions to more than 70% of the overall explained variation. This method focused the score on the set with the largest block of well-correlated genes in the set, while the contributions from other genes, which did not track with other set members, were decreased.²⁴ After obtaining the prognostic value of each gene signature score, we applied an algorithm similar to GGI²⁵ to define the M2score of each patient:

where j is the expression of genes.

Protein–Protein Interaction (PPI) Network Construction and Identification of Hub Genes

The PPI network was constructed using the Search Tool for Retrieval of Interacting Genes (STRING) (https://string-db.org) and Cytoscape software.²⁶ MCODE, a Cytoscape plugin,²⁷ was applied to analyze significant modules in the PPI network. We used the default settings, with degree cutoff at 2, node score cutoff at 0.2, K core at 2 and a maximum depth of 100 to infer modules. The PPI network topological structure and relationship characteristics, including MCC, DMNC, MNC, Degree, EPC, BottleNeck, EcCentricity, Closeness, Radiality, Betweenness, ClusteringCoefficient and Stress, were calculated using the CytoHubba plugin, and the top ranked genes can be considered as the hub genes.

Prediction of Immunotherapy Response

The IMvigor210 cohort is an urothelial carcinoma cohort treated with atezolizumab, an anti-PD-L1 antibody approved by the FDA. This cohort, which is based on the Creative Commons 3.0 License, was used to predict the patients’ response to immunotherapy.²⁸ The clinical information and complete expression data were obtained from http://research-pub.Gene.com/imvigor210-corebiologies. The raw data were normalized and the count value was transformed into the TPM value.

Tumor Immune Dysfunction and Exclusion (TIDE) is a web platform (http://tide.dfci.harvard.edu) that infers gene functions in modulating tumor immunity and evaluates biomarkers to predict the ICB clinical response.²⁹ The TIDE algorithm is a calculation method and can evaluate two different tumor-immune escape mechanisms, namely, the dysfunction of tumor-infiltrating cytotoxic T lymphocytes (CTLs) and the rejection of CTLs by immunosuppressive factors.³⁰ Patients with a higher TIDE score have a higher chance of anti-tumor immune escape, while patients with lower TIDE score show a lower ICB treatment response rate.³¹

Statistical Analysis

The expression of immune-related genes in tumor and normal tissues was compared by one-way ANOVA. The clinical characteristics between the two groups were analyzed by the chi-squared test. The Kaplan–Meier method was applied for the logarithmic rank test in survival analysis. The waterfall function in the “maftools” package was used to visualize the mutational landscape in TCGA-GBM and TCGA-LGG patients. The area under the curve (AUC) of the receiver operating characteristics (ROC) curve was applied to evaluate the performance of the prognostic signature using the R package. All statistical analyses were performed using the R program (version 3.6.1, https://www.r-project.org/). A P-value of 5% or lower was considered to be statistically significant.

Results

Low Fractions of M2 Macrophage Are Associated with Better Overall Survival in GBM

After downloading and processing the data from TCGA and Gravendeel, the gene expression profiles were obtained and then uploaded to the xCell website to analyze the different immune cell fractions. The fractions of M2 macrophages were significantly different in the GBM subtypes, with the mesenchymal subtype having the highest fraction (Figure 1A and B). According to the M2 macrophage fractions, patients were divided into high- and low-fraction groups. The overall survival was analyzed by the Kaplan–Meier survival curve and the results demonstrated that the groups with low fractions of M2 macrophages in both TCGA (Figure 1C) and Gravendeel datasets (Figure 1D) were significantly associated with better outcomes in GBM.

Figure 1 Distribution of M2 macrophages in different subtypes and the prognostic value of M2 macrophages. (A) Distribution of M2 macrophages in GBM subtype – mesenchymal, proneural and classical – in TCGA-GBM dataset. (B) Distribution of M2 macrophages in the Gravendeel dataset. (C) Kaplan–Meier survival curve for proportion of M2 macrophages in TCGA-GBM dataset. (D) Kaplan–Meier survival curve for proportion of M2 macrophages in the Gravendeel dataset. Statistical significance at ***P<0.001 and ****P<0.0001. Abbreviation: ns, not significant.

Acquisition of the M2 Macrophage-Related Gene Module

To identify the M2 macrophage-related genes, we downloaded 4723 immune genes from the ImmPort system. A WGCNA network was constructed by the R programming WGCNA package for this purpose. After removing outlier samples, the soft threshold (β=4, scale-free R²=0.9 in TCGA; β=3 scale-free R²=0.9) was used to guarantee a scale-free network, ultimately identifying 16 distinct modules in TCGA-GBM and five distinct modules in the Gravendeel dataset (Figure 2A–D). The immune genes in the same color module have similar gene expression patterns. Correlations between modules and the fractions of M2 macrophages were calculated, and the blue module, including 424 genes in TCGA-GBM and 461 genes in the Gravendeel dataset, was found to be correlated positively with M2 macrophages (Figure 2E and F). Thus, the genes in the blue module were selected as potential prognostic genes in the following univariate Cox analyses. Finally, 68 genes were confirmed to be significantly associated with overall survival by univariate Cox analyses (Figure 2G). So, these 68 genes were defined as M2 macrophage-related prognostic genes and were included in the subtype analysis.

Figure 2 WGCNA network and module detection. (A) Selection of the soft-thresholding powers in TCGA-GBM. Power 4 was chosen because the fit index curve flattened out upon reaching a high value (>0.8). (B) Selection of the soft-thresholding powers in the Gravendeel dataset. Power 3 was chosen because the fit index curve flattened out upon reaching a high value (>0.9). (C) Correlation matrix for eigengene and M2 macrophages in TCGA-GBM dataset. Each cell includes the corresponding correlations and the P-value. (D) Correlation matrix for eigengene and M2 macrophages in the Gravendeel dataset. (E) Scatter plot of the correlation between the blue module eigengene and M2 macrophages in TCGA-GBM dataset. The correlation coefficient and P-value are indicated in the plot. (F) Scatter plot of the correlation between the blue module eigengene and M2 macrophages in the Gravendeel dataset. (G) Venn diagram showing prognostic immune genes and M2 macrophage-related genes in TCGA-GBM dataset and Gravendeel dataset.

Clustering of GBM by M2 Macrophage-Related Prognostic Genes

The consistent clustering analysis of GBM was performed based on the expression of M2 macrophage-related prognostic genes in TCGA-GBM dataset. According to the cumulative distribution function (CDF), the area under CDF curve analyses and the correlation between groups, the patients were clustered into two subgroups, which we defined as the M2clusterA group and M2clusterB group (Figure 3A). The M2clusterA group included 287 GBM patients while the M2clusterB group included 238 GBM patients. Kaplan–Meier survival analysis confirmed that the M2clusterA group had a better overall survival rate (Figure 3B). Furthermore, the M2 macrophage-related prognostic gene expression, cluster information and clinical features are displayed in the heatmap, which illustrates that the mesenchymal subtype was mainly grouped in M2clusterB while MGMT promoter methylation status and CpG island methylator phenotype (G-CIMP) were significantly associated with the M2clusterA group (Figure 3C). The GO enrichment analysis indicated that the M2 macrophage-related prognostic genes were significantly enriched in immune response, immune regulation and immune receptor activity (Figure 3D). Based on these results, these two cluster groups had different clinical features and molecular characteristics.

Figure 3 Identification of M2 macrophage-related subtypes of GBM in TCGA-GBM dataset. (A) Consensus score matrix of TCGA-GBM samples when k=2. (B) Kaplan–Meier survival curves derived from TCGA-GBM survival data for M2clusterA and M2clusterB groups. (C) Heatmap of the M2 macrophage-related prognostic genes ordered by subtype, with annotations associated with MGMT status, CIMP status, M2cluster and M2score. (D) Signaling pathways associated with M2 macrophage-related prognostic genes. The enrichment of Gene Ontology (GO) pathways related to immune response, immune regulation and immune receptor activity. (E and F) GSEA results on the pathway of M2clusterA. E2F targets, G2M checkpoint and myc targets signaling pathways were significantly upregulated (E), while angiogenesis, epithelial–mesenchymal transition and IL-6/JAK/STAT1 signaling were downregulated (F). (G) Histogram showing differences in the composition of immune cells between M2clusterA and M2clusterB. *P<0.05, **P<0.01, ***P<0.001 and ****P<0.0001. Abbreviation: ns, not significant.

First, for further study of the potential molecular mechanisms between M2clusterA and M2clusterB groups, GSEA was conducted by comparing these two cluster groups, and the result showed that E2F targets, G2M checkpoint and myc targets signaling pathways were significantly upregulated in the M2clusterA group (Figure 3E), while angiogenesis, epithelial–mesenchymal transition and IL-6/JAK/STAT1 signaling were significantly downregulated in the M2clusterA group (Figure 3F).

Second, the immune infiltration in the M2clusterA and M2clusterB groups was analyzed by CIBERSORTx, which is another platform used to verify the different tumor immune microenvironment in the two groups. TCGA-GBM gene expression profiles were uploaded to the website and the results suggested that the fractions of immune cells, such as M0 macrophages, M2 macrophages, activated mast cells and resting NK cells, were different in the M2clusterA and M2clusterB groups, indicating a varied tumor immune microenvironment (Figure 3G).

Somatic Mutation Frequency Landscape and Protein–Protein Interaction Network Construction

To explore the mutations of M2 macrophage-related prognostic genes in different glioma types, we constructed waterfall plots for M2 macrophage-related prognostic genes in GBM and in low-grade glioma (LGG). The graphs showed that 91 out of 393 GBM patients (23.16%) had alterations (Figure 4A), while 43 out of 508 LGG patients (8.46%) had alterations (Figure 4B). Thus, more somatic mutation was observed in M2 macrophage-related prognostic genes in GBM than in LGG.

Figure 4 Somatic mutation frequency landscape and protein–protein interaction network construction. (A) Somatic mutation of M2 macrophage-related prognostic genes in TCGA-GBM dataset. (B) Somatic mutation of M2 macrophage-related prognostic genes in TCGA-LGG dataset. (C) PPI networks of M2 macrophage-related prognostic genes using the STRING tool. Genes are represented as nodes in the plot, and their interactions are denoted by lines. (D) Top 10 hub genes using the MCC algorithm. The color of a node in the PPI network reflects the P-value.

The STRING tool was applied and then a PPI network was constructed to better understand the interplay among the 68 M2 macrophage-related prognostic genes. There were 68 nodes and 473 edges in this network (Figure 4C). CytoHubba is a Cytoscape plugin that allows the use of several topological analysis algorithms to identify the hub genes. In the present study, the MCC method was applied, and 10 hub genes (FCGR2A, CCR1, CD68, ITGAM, C3AR1, TLR2, CYBB, CXCL8, CCR5 and ITGB2) of the whole network are visualized in Figure 4D. The gradual color change from red to yellow represents the connection score, from high to low.

Clinical Application of the PCA-Derived M2score and Relationship to Immunotherapy

Based on PCA and the calculation method, the M2score was obtained (Figure 5A). First, we grouped the IMvigor210 cohort into high- and low-M2score groups; although the Kaplan–Meier survival curve between the two groups did not differ significantly, there was a trend toward statistical significance (P=0.068) (Figure 5B). Then, patients in both TCGA-GBM and the Gravendeel dataset were divided into high- and low-M2score groups. The Kaplan–Meier survival results indicated that the high-M2score group was strongly associated with better overall survival (Figure 5C and D). A dendrogram, including the M2score, M2cluster, CIMP status, survival status, MGMT status, subtype, isocitrate dehydrogenase-1 (IDH1) status, treatment plan and clinical features, is illustrated in Supplementary Figure 1A (Figure S1A). The figure demonstrates that the IDH mutations were significantly associated with higher M2scores. Meanwhile, the IDH mutations were more likely to be grouped in clusterA and had better overall survival. A prognostic nomogram was then developed by combining prognostic factors, including IDH mutation, gender, patient age and M2score (Figure S1B). The predicted probabilities corresponded well with the actual 1-, 3- and 5-year overall survival rates of GBM patients (Figure S1C).

Figure 5 Clinical application of the PCA-derived M2score and its relationship to immunotherapy. (A) Principal component analysis (PCA) for GBM. (B) Kaplan–Meier survival curve for high- and low-M2score groups in the IMvigor210 cohort. (C) Kaplan–Meier survival curve for high- and low-M2score groups in the Gravendeel dataset. (D) Kaplan–Meier survival curve for high- and low-M2score groups in TCGA-GBM dataset. (E) Stacked barplot for immunotherapy response of high- and low-M2score groups based on the IMvigor 210 cohort. (F) ROC curve AUC statistics assessing the predictive power of the M2score for immunotherapy response in TCGA-GBM dataset and Gravendeel dataset based on the TIDE tools.

Anti-PD-1/PDL1 and anti-CTLA4 have provided breakthrough immunotherapy in cancer treatment. To determine whether our immune classification and M2score could predict GBM patients’ response to ICB therapy, we verified the high- and low-M2score based on the IMvigor 210 cohort. The stacked column chart showed that were higher percentages of complete and partial responses (CR and PR) in the high-M2score group (Figure 5E).

TIDE is a computational method which integrates T-cell dysfunction and removal characteristics and simulates tumor immune escape at the level of tumor-infiltrating CTLs to predict the ICB response, which has an advantage in predicting the efficacy of anti-PD1/PDL1 and anti-CTLA4 therapy. After uploading the expression data to the website, the results indicated that 35 GBM patients had immunotherapy efficacy, while 490 patients failed on immunotherapy in the TCGA-GBM cohort. Using the ROC algorithm, we constructed the ROC curve and found that M2score had the capacity to predict the responsiveness of immunotherapy (AUC 0.657, 95% CI: 0.556–0.758) (Figure 5F). In addition, the AUC was 0.643 (95% CI: 0.502–0.643) to evaluate the responsiveness to immunotherapy of M2score in the Gravendeel dataset (Figure 5F), which also demonstrated an appropriate predictive accuracy. Therefore, we concluded that M2score could effectively predict the survival outcome of GBM. Moreover, M2score was a potential biomarker for evaluating the clinical response to immunotherapy, to a certain extent.

Discussion

Based on morphological similarities, GBM has traditionally been classified into gliosarcoma, giant cell GBM and epithelioid GBM by the WHO Classification of Tumors of the Central Nervous System (CNS).³² With features containing both neoplastic cells and stroma, GBM has histological heterogeneity, with disappointing prognosis. Even for patients who have similar histopathological features, the molecular biological characteristics are variable, which makes it challenging to effectively target patients using different treatment strategies. Genome-wide transcriptome analysis led to the classification of GBM into four distinct molecular subtypes – mesenchymal, classical, proneural and neural – based on the bulk tumor gene expression profile.⁹ However, the clinical relevance of subtyping is still controversial.¹³ Gill et al reported that the neural subtype contained a variable mixture of infiltrating glioma cells and non-neoplastic brain cells, providing a particularly robust and clinically relevant dataset, as the samples were taken from the margins of the tumor.¹⁰ Therefore, only three common subtypes that related to the center of the tumor were investigated in the current study.

Although many achievements have been made in targeted therapies and immunotherapies, the results are still unsatisfactory owing to complex pathophysiology. Therefore, there is an urgent need to establish an individualized risk stratification approach, and to classify GBM patients into appropriate risk groups to guide individualized treatment.³³ In the present study, high-throughput gene and transcript-level data from TCGA-GBM and Gravendeel’s microarray data were analyzed to screen immune-related DEGs. In order to efficiently estimate the level of immune cell infiltration, xCell and CIBERSORTx tools, which can analyze traditional bulk expression profiles and then calculate the immune infiltration scores, were applied. The fractions of M2 macrophages and M2 macrophage-related immune genes were identified, and then clustering of the GBM by M2 macrophages was conducted and the M2score was established successfully. Finally, we found that the M2score could effectively predict the survival outcome of GBM and this is a potential biomarker for evaluating the clinical response to immunotherapy.

In the past few decades, the methods of uncovering tumor and immune cell interactions have improved the available treatment options for GBM.³⁴ Mounting evidence suggests that immune suppression induced by cancer cells is the major factor responsible for poor outcomes.³⁵ Macrophages are one of the most abundant innate immune cell types in the microenvironment of various tumors.³⁶ M2 macrophages play a pivotal role in immune suppression and pro-tumor functions during GBM progression because of the secretion of immunosuppressive factors. Our results verified that GBM patients with a low fraction of M2 macrophages had significantly better outcomes in both TCGA and Gravendeel datasets. Sørensen et al also found that M2 macrophages had an unfavorable prognostic value in high-grade gliomas.³⁷ Furthermore, a study performed by Zhu et al³⁸ showed that PD-L1-mediated immunosuppression attributed to the infiltration and M2-polarization and correlated with poor survival. Taking these findings into account, the presence of M2 macrophages is important for the progression of GBM and has been well studied. Assessing the mechanism of M2 macrophages in GBM may provide more information on treatment options.

The WGCNA network was constructed and then 68 immune genes were identified as M2 macrophage-related prognostic genes. Based on these 68 gene expression profiles, the patients in TCGA-GBM were divided into M2clusterA and M2clusterB groups. The M2clusterA group had higher M2score, a smaller fraction of M2 macrophages, more MGMT promoter methylation and CpG island methylator phenotype, and consequently a better overall survival rate. Many trials have reported that MGMT promoter methylation is associated with longer survival in GBM patients treated with TMZ, both in the overall population and in older patients.^39–42 As for CpG island methylation, researchers have found that it is a relevant prognostic factor, independently of the known predictors, such as grade and age.⁴³ The favorable prognostic value of CpG island methylator phenotype in GBM patients has been reported in many other studies.^43–45 Compared with the non-carrying CpG island methylator phenotype, the authors found that the carrying CpG island methylator phenotype was more enriched among the proneural subtypes, which had a better prognosis.^9,45 The GSEA results showed that E2F targets, G2M checkpoint and myc targets signaling pathways were significantly upregulated in the M2clusterA group, while angiogenesis, epithelial–mesenchymal transition and IL-6/JAK/STAT1 signaling were significantly downregulated in the M2clusterA group. Conversely, Yu et al found that expression of E2Fs was significantly upregulated and linked with grade progression.⁴⁶ Meanwhile, chromatin remodeler HELLS strongly correlated with targets of E2F3 and MYC transcriptional activity in GBM patients.⁴⁷ These results imply that E2F targets and myc targets signaling pathway could also affect the prognosis, independently of the M2 macrophage or immune response. Little is known about the role of the G2M checkpoint in the pathophysiology of GBM.

Ten hub genes, namely FCGR2A, CCR1, CD68, ITGAM, C3AR1, TLR2, CYBB, CXCL8, CCR5 and ITGB2, were identified by CytoHubba, all of which, except for CXCL8, were associated with the plasma membrane. CXCL8, also known as IL-8, was predominantly secreted by macrophages and contributed to the immunosuppressive microenvironment.⁴⁸ Seven genes were enriched in the inflammatory response (CCR1, CCR5, CXCL8, C3AR1, CYBB, ITGB2 and TLR2), while four genes each were enriched in the immune response (CCR1, CCR5, CXCL8 and TLR2), positive regulation of angiogenesis (CXCL8, C3AR1, CYBB and ITGB2) or chemotaxis (CCR1, CCR5, CXCL8 and C3AR1). TLR2, C3AR1, ITGAM, CD68 and CCR5 have been widely reported as a marker of macrophages in GBM, and studies have shown that TLR2, CD68 and CCR5 act as immune suppressors and contribute toward GBM progression in the tumor microenvironment.^49–53 FCGR2A, CYBB, CXCL8 and ITGB2 were significantly upregulated in GBM and associated with shorter overall survival of GBM patients.^54–58 As for the role of CCR1 in GBM, the results are still controversial. Zhang et al found that CCR1 and CCR5 were two key receptors of CCL8 and they triggered CCL8-induced invasion of GBM cells,⁵⁹ while Pham et al reported that individual deletion of CCR1 or CCR5 had little or no effect on the survival of tumor-bearing mice, or the number of GBM-infiltrated microglia/macrophages and lymphocytes.⁶⁰ Consequently, these 10 hub genes are associated with a higher degree of malignancy and poorer prognosis in gliomas.

Furthermore, we investigated the relationship between the M2score and immunotherapy to guide personalized treatment decisions. Nowadays, immune checkpoint inhibitors have been approved for gastric cancer, breast cancer, colon cancer and other cancers.⁶¹ However, they have demonstrated limited efficacy in GBM.⁶² Hence, for the future development of immunotherapy in GBM, it is necessary to construct an efficient immune classification or prognostic model to predict the response of GBM patients to ICB therapy.^63,64 Based on the results presented in this report, we successfully verified M2score as an appropriate predictive score to predict the response to ICB therapy. Both the IMvigor210 cohort and the TIDE score showed us that the M2score is a potential biomarker for evaluating the clinical response to immunotherapy.

we still have some limitations in the current study. First, this is a retrospective study that may display statistical bias, and the bulk sequence transcriptome data were used, in which there is a lack of comprehensive exploration for intratumoral heterogeneity. Second, the tumor immunity analysis was calculated by xCell, which cannot represent real immune cell infiltration. Finally, the urothelial cancer cohort, instead of a GBM cohort, was applied to explore the connection between gene expression and the immunotherapy response. Therefore, prospective studies and in vitro and in vivo experiments regarding these genes are required to validate our in silico results, and the response of immunotherapy should be further verified in randomized clinical trials.

Conclusion

In summary, by detecting distinct fractions of the M2 macrophage, this study successfully identified immune-related prognostic genes and constructed the M2score, which could be employed to quantitatively predict the prognosis of GBM patients. The activation pathway, immune infiltration, mutation features and drug response underlying the M2score were also identified. Our findings provide a basis for understanding the role of M2 macrophages and immune-related prognostic genes, and indicate the potential clinical implications of M2score in GBM.