Immune Infiltration-Related Signature Predicts Risk Stratification and Immunotherapy Efficacy in Grade II and III Gliomas

Background: Tumor microenvironment, especially infiltrating immune cell, is crucial for solid tumors including glioma. However, the hub genes as well as their effects on patient prognosis and immunotherapy efficacy remain obscure.

Methods: We employed a total of 952 lower grade glioma (LGG) patients from The Cancer Genome Atlas (TCGA) and Chinese Glioma Genome Atlas (CGGA) databases, and 24 samples in our hospital for subsequent analyses. Abundances of immune infiltrates were evaluated using CIBERSORT and ImmuCellAI. Their correlations with prognosis were assessed by log-rank test. Immune infiltration-related hub genes were obtained from overlapped differential expressed genes (DEGs) in various subsets of survival-related immune cell types. The risk signature was constructed by Least Absolute Shrinkage and Selection Operator (LASSO) Cox regression analysis. The functional analyses were estimated by GVSA and Gene Set Enrichment Analysis (GSEA) algorithms. And protein–protein interaction enrichment analysis was carried out with the Metascape database integrating STRING, BioGrid, OmniPath, and InWeb_IM.

Results: Among the 21 infiltrates, the abundances of five immune infiltrates were correlated with overall survival (OS) in LGG patients. Higher abundances of naïve CD4+ T cells (p = 0.002), activated mast cells (p = 0.015), and monocytes (p = 0.014) were correlated with better prognosis, while higher abundances of resting memory CD4+ T cells (p = 0.015) and M1 macrophages (p = 0.020) correlated with poorer OS. We finally obtained 44 hub genes and constructed an immune infiltration-related signature (IIRS). The IIRS correlates with clinicopathological characteristics and exhibited potential power in predicting the immunotherapy efficacy. The IRRS correlates with cancer related pathways, especially “epithelial-mesenchymal transition (EMT),” and cytotoxic T lymphocytes.

Conclusion: Our study constructed and validated a novel signature for risk stratification and prediction of immunotherapy response in grade II and III gliomas, which was closely associated with glioma immune microenvironment and could serve as a promising prognostic biomarker for glioma patients.

Introduction

Gliomas are the most common primary tumors of the central nervous system that arise from the intrinsic constituent cells of the brain (Sanai et al., 2005). They have historically been classified on the basis of their microscopic and immunohistochemical resemblance and have been graded according to histological features indicative of biological aggressiveness.

Genome-wide molecular-profiling studies have revealed comprehensive genomic landscapes for major types of gliomas (Suzuki et al., 2015). These developments have identified novel biomarkers for improved tumor classification and promising therapeutic targets. Predictive biomarkers recognized and used in clinics mainly included isocitrate dehydrogenase (IDH) mutation, the discovery of which constituted a key breakthrough in the understanding of WHO grade II/III gliomas (Yan et al., 2009). Besides, the presence of O6-methylguanine-DNA methyltransferase (MGMT) promoter methylation is predictive of response to temozolomide-based chemotherapy in patients with IDH-wild-type glioma (Wick et al., 2012). Furthermore, 1p/19q codeletion is suggested as a predictive marker of benefit from upfront combined radiotherapy and chemotherapy in two phase III trials (Bent et al., 2013; Cairncross et al., 2013). Novel pathogenesis-based treatments targeting oncogenic signaling pathways such as BRAF mutation (Robinson et al., 2014), epidermal growth factor receptor (EGFR) amplification (Phillips et al., 2016), and fibroblast growth factor receptor (FGFR)-TACC fusion demonstrate a potential for lower grade gliomas (LGG) elimination (Stefano et al., 2015). Despite a deconstruction at molecular level that furthered our understanding of tumorigenesis and personalized therapy, a certain LGG population acquired resistant to these targeted therapies.

For solid tumors, the tumor microenvironment, composed of extracellular matrix, stromal cells, and immune cells, plays a crucial role in the initiation and progression of cancer. The infiltrating immune cells there are very vital as they were associated with patient prognosis in various cancers (Yang et al., 2019; Huang et al., 2020; Zhang et al., 2020, 2021). Exploration of the glioma microenvironment will provide a better understanding of the occurrence and development of glioma. It is, however, worth noting that gliomas are not considered highly immunogenic since mutational loads are typically low, and gliomas are characterized by profound immunosuppression mediated by immune-inhibitory factors (Nduom et al., 2015). Therefore, we hope to decode the unique immune microenvironment and identify novel biomarkers to overcome immunosuppression, exploit antitumor immune responses, and guide individualized treatments.

We herein conducted a comprehensive analysis based on two independent cohorts, plus our own samples to explore the profile of infiltrating immune cells in gliomas in order to better understand its biological functions there. Furthermore, a risk signature based on immune infiltration (IIRS) was constructed to predict the prognosis of patients diagnosed with LGGs. Multifaceted performance of the IIRS was also examined to reveal its superior predictive ability for response to immunotherapy.

Materials and Methods

Data Extraction

All transcriptomic and clinical characteristics of enrolled samples were extracted from The Cancer Genome Atlas (TCGA) and Chinese Glioma Genome Atlas (CGGA) databases (Brat et al., 2015; Zhao et al., 2021). A total of 952 primary LGG samples with detailed clinical information were enrolled in our study, in which 508 samples extracted from the TCGA database were defined as the training set; 444 samples extracted from the CGGA database were defined as the validation set. Normal or glioblastoma (GBM, grade IV glioma) samples were excluded.

Immune Infiltration and Survival Analysis in Lower Grade Gliomas

CIBERSORT algorithm was employed for evaluating the percentage of 21 human hematopoietic cell phenotypes, including seven T cell types, naïve and memory B cells, plasma cells, natural killer cells, and myeloid subsets (Newman et al., 2015). The associations between infiltrating abundance and overall survival (OS) in TCGA LGG cohort were evaluated using univariate Cox analysis, and the survival curves were correspondingly established by Kaplan–Meier analysis.

Construction and Validation of the Immune Infiltration-Related Signature

Considering that the abundance of these infiltrating cells was mostly low, we used the mean abundance to divide the entire LGG population into high- and low-infiltrating groups for each infiltrate. For the immune infiltrates that were significantly correlated to the outcomes of LGG patients, we conducted a differential expression analysis using the R package “limma” and obtained fold change and p-value for each gene. Subsequently, for those infiltrating cells that were detrimental to patient prognosis, we selected genes ranked in the top 500 by fold change; whereas for those with a beneficial effect on patient prognosis, we selected genes whose fold change ranked in the bottom 500. Finally, we got five gene sets with 500 gene in each one.

Hub genes are derived from the intersection of the above gene sets, which were screened by using univariate Cox regression analysis. Thereafter, we used the R package “glmnet” to conduct Least Absolute Shrinkage and Selection Operator (LASSO) Cox regression analysis (with the penalty parameter estimated by 10-fold cross-validation), we developed an immune infiltration-related signature (IIRS) for the LGG patients. The risk score calculating formula is:

$$Riskscore=\sum _{i=1}^n\left({\mathrm{\beta }}_i\times x_i\right)$$

where “n” means the number of genes included in the model, “β_i” means the LASSO coefficients, “xi” is the expression level of each model gene.

Risk scores were subsequently computed for all patients included in our study. For both cohorts, the patients were divided into high- and low-risk groups according to the median risk score. Then risk plots, scatter diagrams, heatmaps, survival curves, and time-dependent receiver operating characteristic curves (ROC) were plotted using the R package ‘‘ggplot2.’’ The relationships between risk signature and survival as well as other clinicopathological characteristics were also assessed. We calculated the correlations among the signature genes, as well as the correlations between individual gene and OS and progression free survival (PFS) in LGG populations. A principal components analysis (PCA) analysis, based on Gene Expression Profiling Interactive Analysis (GEPIA2¹) webtool (Tang et al., 2019), was performed to examine the resolving power of the IRRS.

Sample Collection and RNA Sequencing

Twenty-four samples were collected and then immediately stored in liquid nitrogen. Total RNA was extracted from the tissues using TRIzol (Invitrogen, Carlsbad, CA, United States) following the instructions. The mRNA library was then constructed after quantification using NanoDrop and Agilent 2100 bioanalyzer (Thermo Fisher Scientific, MA, United States). Total RNA was purified and fragmented into small pieces for cDNA synthetization. The cDNA fragments were further amplified by polymerase chain reaction after incubating with A-tailing mix and RNA Adapter Index for end repair. The qualified double-stranded PCR products were then used to construct the final library. Eventually, the 24 qualified glioma samples were further sequenced on a BGISEQ-500 platform (BGI-Shenzhen, China). The gene expression levels were calculated using RSEM (v1.2.12). The sequencing and clinical data of these samples were summarized in Supplementary Table 1.

Prediction of Therapy Efficacy and Drug Response

The efficacies of four therapies (radiotherapy, chemotherapy, targeted therapy, and immunotherapy) in high-risk and low-risk groups were evaluated. We used Tumor Immune Dysfunction and Exclusion (TIDE²) algorithm to assess the ability of the IIRS in predicating response to the immunotherapy.

Two resources for therapeutic biomarker discovery in cancer cells, including Cancer Therapeutics Response Portal (CTRP) and Genomics of Drug Sensitivity in Cancer (GDSC) (Yang et al., 2012; Rees et al., 2016), were employed to evaluate the relationship between drug sensitivity (IC50) and mRNA expression.

Gene Set Variation Analysis, Cancer Related Pathway, and Infiltrating Immune Cells

The Gene Set Variation Analysis (GSVA) score represents the integrated level of the expression of model gene set. The GSVA score of each patient was calculated using R package “GSVA” (Hänzelmann et al., 2013).

Data of reverse phase protein array (RPPA), a high-throughput antibody-based technique with the procedures similar to that of western blots, were used to calculate pathway activity score of 10 cancer related pathways, including TSC/mTOR, RTK, RAS/MAPK, PI3K/AKT, Hormone ER, Hormone AR, EMT, DNA Damage Response, Cell Cycle, Apoptosis pathways.

The infiltrates of 24 immune cells were evaluated through Immune Cell Abundance Identifier (ImmuCellAI) (Miao et al., 2020).

Gene Set Enrichment Analysis, Mutational Profiles, and Protein–Protein Interaction

We calculated the degree to which the inputted gene set is overrepresented at the top or bottom of all genes ranked by gene expression fold change between high- and low-risk groups. The GSEA calculation was performed based on R package ‘‘fgsea.’’ The mutational profiles of model genes were assessed in cBioPortal website.³

For the differential expressed genes (DEGs) between high- and low-risk groups, protein–protein interaction enrichment analysis has been carried out with the Metascape databases integrating STRING (Szklarczyk et al., 2019), BioGrid (Stark et al., 2006), OmniPath (Türei et al., 2016), and InWeb_IM (Li et al., 2017). The resultant network contains the subset of proteins that form physical interactions with at least one other member in the list. The Molecular Complex Detection (MCODE) algorithm (Bader and Hogue, 2003) has been applied to identify densely connected network components.

Statistical Analysis

Kaplan–Meier curve and log-rank test were used to compare the survival between various subgroups. The Student’s t-test was used to compare the risk scores between pairs of subgroups based on the following clinicopathologic features: age at initial pathologic diagnosis (≤40 vs. > 40 years old), gender (male vs. female), WHO grade (II vs. III), histological type (astrocytoma, oligoastrocytoma, and oligodendroglioma), and IDH1 status (mutant vs. wild-type). Wilcoxon test and Kruskal–Wallis test were used for comparison between two groups, and for comparison among more than two groups, respectively. p < 0.05 was the significance threshold in most analyses. The statistical analyses were achieved by using R language (version 4.0.3).

Results

The Abundances of Five Immune Infiltrates Were Correlated With Overall Survival in Lower Grade Glioma Patients

We obtained data from two databases (TCGA, n = 508; CGGA, n = 444) to evaluate the abundances of immune infiltrates in LGGs. The characteristics of patients in the two cohorts were summarized in Table 1. The distribution and correlations of immune infiltrates in the TCGA cohort was displayed as Figures 1A,B. Among the 21 infiltrates, the abundances of five immune infiltrates were correlated with OS in grade II and III glioma patients (Figures 1C–G). Specifically, higher abundances of naïve CD4+ T cells (p = 0.002), activated mast cells (p = 0.015), and monocytes (p = 0.014) were related to better prognosis, while higher abundances of resting memory CD4+ T cells (p = 0.015) and M1 macrophages (p = 0.020) correlated with poorer OS.

TABLE 1

Table 1. Detailed characteristics of included patients.

FIGURE 1

Figure 1. Immune infiltration profile and its relationship to survival in gliomas. (A) Proportions of various immune infiltrates. (B) The correlations among various immune infiltrates (non-significant values were omitted). (C–G) Kaplan–Meier survival curves based on immune infiltrates that were significantly related to patient prognosis.

Identification of Candidate Genes and Construction of the Risk Signature

According to the mean value of abundances of the five immune infiltrates, DEGs were obtained, respectively (Figures 2A–E). We profiled 500 candidate genes in these five groups, respectively. And 44 hub genes were finally obtained when overlapping the five sets containing 2,500 genes (Figure 2F). All 44 genes were significantly associated with patient outcome, and 10 genes were profiled by LASSO regression analysis to construct the IIRS (Figures 2G,H). Detailed descriptions, LASSO coefficients, and hazard ratios for all model genes are summarized in Table 2. And we exhibited the detailed expression levels as well as prognostic curves of these genes in Supplementary Figure 1. Briefly, all ten genes except ACTN1 were differentially expressed between tumor and normal tissues. Among the remaining nine genes, FABP and PLAT were decreased in tumor tissues, while the others exhibited the opposite distribution (Supplementary Figure 1A). Regarding survival, these model genes were risk factors for OS, PFS, and disease specific survival (DSS) in patients with LGG (Supplementary Figure 1B). And the specific Kaplan–Meier survival curves (for OS) were displayed in Supplementary Figure 1C.

FIGURE 2

Figure 2. Identification of immune infiltration-related hub genes and LASSO analysis. (A–E) Volcano plot of DEGs based on the mean abundances of survival-related infiltrates in glioma. (F) Forty-four hub genes identified by overlapping the five gene sets. (G) LASSO coefficient profiles of the 44 hub genes. Each curve represents a coefficient, and the x-axis represents the regularization penalty parameter. As λ changes, a coefficient that becomes non-zero enters the LASSO regression model. (H) Cross-validation to select the optimal tuning parameter (λ). The left dotted vertical line crosses over the optimal log λ, which corresponds to the minimum value for multivariate Cox modeling.

TABLE 2

Table 2. Detailed information of model genes.

External and Subgroup Validation Demonstrates Stability of the Immune Infiltration-Related Signature

Risk plots, survival distributions, and model gene expressions were plotted in Figure 3A. Kaplan–Meier survival curve indicated that LGG patients with higher risk scores had significantly worse outcomes in the training set (p < 0.0001, Figure 3B). The time-dependent ROC curve demonstrated a promising ability of the model to predict OS in the training cohort (1-year AUC = 0.66, 3-year AUC = 0.69, 5-year AUC = 0.78; Figure 3C). The results were similar in the external CGGA cohort (Figure 3D). Higher risk scores also indicated poorer OS (p < 0.0001, Figure 3E). The risk model retained stable and high predication ability (1-year AUC = 0.69, 3-year AUC = 0.70, 5-year AUC = 0.74; Figure 3F). These results showed that the IIRS had a robust and stable OS-predictive ability for LGG patients. Furthermore, we performed a stratification analysis and found that the risk model maintained the ability to predict OS in most subgroups in both cohorts (Figures 3G,H).

FIGURE 3

Figure 3. Construction and validation of the IIRS. (A) Risk score, survival status, and expression pattern of model genes in each patient in the training cohort. (B) Kaplan–Meier analysis of patients in the high- and low-risk groups. (C) Time-dependent ROC analysis of the IIRS in predicting prognosis. (D–F) Validation of the IIRS in the CGGA cohort. (G,H) Kaplan–Meier analysis of patients in the high- and low-risk groups in stratified subgroups in the TCGA and CGGA cohorts.

The IRRS Correlates With Clinicopathological Characteristics and Predicts Immunotherapy Efficacy

Sankey diagrams were displayed showing the distribution of risk scores and clinicopathologic characteristics among LGG patients (Figures 4A,B). In the TCGA cohort, LGG patients with higher WHO grade had higher risk scores, while the risk score was not associated with gender. Besides, an individual patient would have a higher risk score if he had a pathologic type of astrocytoma (Figure 4C). In the CGGA cohort, risk scores were higher in patients with 1p/19q codeletion (Figure 4D). Importantly, LGG patients with wild-type IDH1 would have higher risk scores in both cohorts. Further Cox analyses showed that higher age, higher grade, pathological type (astrocytoma), and higher risk score were significantly associated with worse survival in both cohorts (Figures 4E,F).

FIGURE 4

Figure 4. Exploration the relationship between clinicopathologic characteristics and the IIRS. (A,B) Sankey diagrams exhibiting the distribution of clinicopathologic characteristics in high- and low-risk groups in the two cohorts. (C,D) The risk score in different subgroups stratified by age, gender, grade, histological type, IDH1 status, and 1p/19q codeletion status in the TCGA and CGGA cohorts. (E,F) Cox analyses examining different variables for LGG patient survival in the TCGA and CGGA cohorts.

Using samples from Xiangya hospital, we found the consistent results that glioma patients with histological type of astrocytoma, wild-type IDH1, and unmethylated MGMT had higher risk scores (Figure 5A). Although no differences in tumor purity were observed between low- and high-risk groups (Figure 5B), we found that, in our samples, high risk patients had higher abundances of M0 macrophages and activated mast cells, and lower abundances of naïve B cells, monocytes, M2 macrophages, resting mast cells, and neutrophils (Figure 5C).

FIGURE 5

Figure 5. Validation using collected samples and prediction of therapy efficacy. (A) The risk score in different subgroups stratified by age, gender, grade, histological type, IDH1 status, and MGMT status using own samples. (B) Stromal, immune, and ESTIMATE scores compared between high- and low-risk groups. (C) The abundances of various immune cells in the high- and low-risk groups. (D,E) The efficacies of chemotherapy, radiation therapy, targeted therapy, and immunotherapy in the low- and high-risk groups. (F) Comparison of the abilities in predicting the response to immunotherapy with recognized biomarkers or genes. *p < 0.05, **p < 0.01, ***p < 0.001, ns: not significant.

Furthermore, we investigated the efficacies of multiple treatments in low- and high- risk groups (Figures 5D,E). Chemotherapies demonstrated high efficacies in both groups, while radiation and targeted therapies did not improve patients’ outcome. Intriguingly, patients in the low-risk group uniquely responded well to immunotherapy. To better evaluate the potential of our IIRS in predicting patients’ responses to immunotherapies, we employed TIDE algorithm to compare our model with other published biomarkers in immunotherapy response prediction (Figure 5F). Compared with recognized signatures (or genes) including Merck18, TIDE score, microsatellite instability (MSI) score, tumor mutation burden (TMB), CD274, CD8, IFNG, T clonality, and B clonality, the IIRS showed robust ability in predicting the response to immunotherapies in patients with different cancers. Besides, the correlations between drug sensitivity and model gene mRNA expressions were exhibited in Supplementary Figures 2A,B.

The IRRS Correlates With Cancer Related Pathways and Cytotoxic T Lymphocytes

To better understand how the IRRS participates in oncologic processes, we analyzed our model both intrinsically and extrinsically. Strong and significant correlations existed among the model genes (Figure 6A), all of which were risk factors for the OS in the LGG cohort (Figure 6B). Besides, the expressions of all genes, except for LGALS3, were correlated with poorer PFS in patients with LGG (Figure 6B). The comprehensive investigation in the mutational profiles of model genes was displayed as Supplementary Figure 3, where we found BGN was most likely to mutate in LGGs. These mutations were significantly related to molecular status including IDH, 1p, and 19q. A subsequent PCA analysis suggested the IRRS could distinguish LGG from normal cerebral cortex tissues and glioblastomas (Figure 6C).

FIGURE 6

Figure 6. Comprehensive analysis of model genes in the IIRS. (A) Spearman correlations among the model genes. (B) Cox analyses exhibiting the relationships between the expressions of model genes and OS and PFS. (C) PCA analysis. (D) Involvement of individual model genes in cancer related pathways. (E) Correlation between GSVA score and cancer related pathways. (F) Correlation between GSVA score and immune cell types.

Next, we evaluated the involvement of individual genes in various cancer related pathways, and found that all genes participated in activated “epithelial-mesenchymal transition (EMT)” pathway while five of them were involved in inhibited “DNA damage” pathway (Figure 6D). Furthermore, GVSA score was significantly and positively correlated with EMT pathway (r = 0.27, p = 9.6e-9) while negatively correlated with “Hormone androgen/androgen receptor (AR)” pathway (r = −0.16, p = 5.1e-4) (Figure 6E).

We then focused on the association of the IRRS and multiple immune cell types (Figure 6F). Strong and positive associations were observed between GVSA scores and activated T lymphocytes including cytotoxic T lymphocytes (r = 0.61, p = 3.6e-60) and type 1 helper T cells (r = 0.51, p = 1.3e-39). In contrary, the GVSA score was strongly and negatively correlated with the abundance of neutrophil (r = −0.47, p = 2.4e-33).

Functional Analysis Reveals Deep Involvement of Immune Infiltration-Related Signature in the Glioma Immune Microenvironment

We obtained DEGs between low- and high-risk groups (Figure 7A), which were displayed as volcano plot (Figure 7B). Functional analysis indicated these DEGs were mainly located in extracellular matrix and were involved in binding functions such calcium ion, kinase, and cell adhesion molecule binding (Figure 7C). As for biological processes, they participated in responses to wounding and inflammatory responses (Figure 7D). Furthermore, the gene sets were mainly enriched in several Reactome pathways including “cytokine signaling in immune system,” “hemostasis,” and “neutrophil degranulation” (Figure 7E). To further capture the relationships between the terms, a subset of enriched terms has been selected and rendered as a network plot, where terms with a similarity above 0.3 are connected by edges. We select the terms with the best p-values from each of the 20 clusters and visualized the network as Supplementary Figures 2C,D (Shannon et al., 2003).

FIGURE 7

Figure 7. Functional analysis of the IRRS. (A) Heatmap of differential expressed genes between high- and low-risk groups. (B) Volcano plot of differential expressed genes (log₁₀FDR < 0.05, | log₂(fold change)| > 1). (C) GO and Reactome enrichment analysis. (D,E) GSEA analysis for KEGG and HALLMARK terms. (F,G) Protein–protein interaction displayed by enriched MCODE networks.

The MCODE networks identified for individual gene lists have been gathered and are shown in Figures 7F,G. Pathway and process enrichment analysis has been applied to each MCODE component independently, and the three best-scoring terms by p-value have been retained as the functional description of the corresponding components: (1) GO:0002399, MHC class II protein complex assembly (log₁₀p = −15.7); (2) GO:0002503, peptide antigen assembly with MHC class II protein complex (log₁₀p = −15.7); and (3) GO:0009611, response to wounding (log₁₀p = −15.4).

Discussion

Glioma is a common tumor in human central nervous system. Over the past decades, surgical section with radiotherapy and chemotherapy still represented the mainstream treatment against glioma. Due to the unique microenvironment, gliomas acquire immunosuppressive phenotypes and poorly response to established immunotherapies. Therefore, novel targets are in urgent need for the hope to predict response rate and improve patient outcome. We employed comprehensive bioinformatic analyses to build an IIRS, which helps clinicians to optimize the management of LGGs. We believe that the IRRS is a good predictor of outcomes in LGG, and targeting the model genes there will demonstrate encouraging efficacy in future preclinical and clinical practice.

Immune infiltration has gained widespread attention in the last decade, especially its unique involvement in malignant processes such as tumor progression and immunotherapy resistance, making it a promising target in tumor microenvironment (Hinshaw and Shevde, 2019). Previous studies have investigated the relationships between immune infiltration and prognosis in various cancer types. High abundance of M2 macrophages was reported to be related to poorer survival in patients with bladder cancers, while CD8+ T cells were related to improved prognosis in this cancer type (Zhang et al., 2020). In addition, CD4+ naive T cells, regulatory T cells, M2 macrophages, resting mast cells were identified as risk immune cells in digestive system cancers, while the abundances of naive B cells, CD8+ T cells, CD4 memory activated T cells, follicular helper T cells, and eosinophils were correlated with better relapse free survival (Yang et al., 2019). To the best our knowledge, this study is the first to comprehensively evaluate the effects of various immune infiltrates on glioma patients’ survival. We assessed the profile of immune infiltration, and found that naïve CD4+ T cells, activated mast cells, and monocytes were protective factors, while resting memory CD4+ T cells and M1 macrophages were risk factors for the prognosis of patients with grade II and III gliomas.

More recent studies have reported that myeloid cells and B cells in the meninges mainly originate from the calvaria bone marrow, rather than the peripheral circulation (Brioschi et al., 2021; Cugurra et al., 2021). Here, B cells mature and develop in the meninges rather than in the bone marrow as recognized (Brioschi et al., 2021). Immune cells may be directly transported through vessels present between the skull and dura. These findings collectively point to the fact that the brain, unlike peripheral organs, has a “self-sufficient” and relatively independent immune system. Further studies should be conducted to assess whether these cells elicit similar effects as we estimated. Meanwhile, studies should be conducted to examine whether these effects are mediated by cytotoxic T lymphocytes or helper T cells.

Forty-four hub genes were selected and used to establish the IIRS. The signature retained stability in different dataset and subgroups. For the model genes in the IIRS, previous studies had provided evidence for their effects in malignancies. Insulin like growth factor (IGF)-binding protein 2 (IGFBP2) is an IGF system regulator and a developmentally regulated gene. Accumulating evidence indicates that in solid tumors, IGFBP2 is upregulated and promotes several key oncogenic processes, such as epithelial mesenchymal transition, cell migration, invasion, angiogenesis, stemness, transcriptional activation, and epigenetic programming through signaling, thus being a hub of oncogenic networks and a potential therapeutic target for cancer treatment (Li T. et al., 2020). Furthermore, in a study including 2447 glioma samples with gene expression profiles, IGFBP2 was found to be involved in immunosuppressive activities and was an independent unfavorable prognostic biomarker (Cai et al., 2018). Actinin alpha 1 (ACTN1) has been identified as a glioma microenvironment-related gene with prognostic value in malignant gliomas (Li Y. et al., 2020). Moreover, lectin, galactoside-binding, soluble, 3 (LGALS3) is a poor prognostic factor in diffuse glioma (Hu et al., 2020), and it can promote therapeutic resistance there (Wang et al., 2019). Fibromodulin (FMOD) was upregulated in glioma and could promote glioma cell migration by inducing the formation of filamentous actin stress fibers. Both FMOD promoter methylation and transcript levels predict prognosis in gliomas (Mondal et al., 2017). Furthermore, fatty acid binding protein 5 (FABP5) was identified as one of the most enriched genes and its elevation revealed severe outcomes in malignant LGGs. And the malignant properties of LGGs were promoted by exogenous overexpression of FABP5 through tumor necrosis factor α-dependent NF-κB signaling (Wang et al., 2021). Importantly, a high-quality study linked metabolism and tumors, where it found that gliomas expressed high levels of branched chain amino-acid transaminase 1 (BCAT1). Inhibition of BCAT1 in glioma cell lines blocked glutamate excretion and resulted in decreased proliferation and invasiveness in vitro, as well as a significant decrease in tumor growth in a glioma xenograft model, indicating a central role of BCAT1 in glioma pathogenesis (Tönjes et al., 2013). Additionally, collagen, type I, alpha 2 (COL1A2) was identified as a hub gene in glioma in several studies (Yin et al., 2021). These findings comprehensively supported our IRRS in glioma, as most model genes were elaborated to be involved in the initiation, progression, or treatment resistance in gliomas.

A key finding in our study is that glioma patients in the low-risk group exhibited a unique response to immunotherapy, as those who received immunotherapy had significantly improved survival compared with those who did not. Taken together with the TIDE algorithm showing the excellent predictive ability of our model for immunotherapy response, we conclude that the IIRS can reflect the sensitivity of LGGs to immunotherapy and recommend this model to guide clinical decisions.

There were several limitations in our study. First, the size of samples for validation was too small, thus the accuracy of the validation is open to question. Second, we included 10 genes in the IIRS, proposing a great challenge for experimental validation. Although we used patients from the CGGA database (validation set) to validate the results obtained from the TCGA database (training set) and showed good concordance, multicenter cohorts with large sample size and complete clinical data are still needed to elaborate our conclusions. Third, there is no LGG cohort in the TIDE database, so the prediction of response to immunotherapy in LGG needs further validation.

Conclusion

In summary, our study establishes a model based on glioma immune infiltration profiles that accurately predicts patient prognosis and response to immunotherapy, with the expectation of aiding decision making in clinics.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: NCBI SRA BioProject, accession no: PRJNA767573.

Ethics Statement

The studies involving human participants were reviewed and approved by the Ethical Committee of Xiangya Hospital (No. 2017121019). The patients/participants provided their written informed consent to participate in this study.

Author Contributions

WY performed the data analysis and interpreted the data and prepared the draft. FL collected the samples in our cohort and they were responsible for the subsequent RNA sequencing of them. CL performed the visualization. ZL and FL revised the manuscript. CL and FL designed the research and supervised all the work. All authors have read and approved the final manuscript, and agreed to be accountable for the content of the work.

Funding

This research was funded by the National Natural Science Foundation of China under grant no. 82001223, the National Natural Science Foundation of China under grant no. 81901401, and the Natural Science Foundation for Young Scientist of Hunan Province, China (Grant No. 2019JJ50952).

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.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

We thank all reviewers for their comments, which have greatly improved the quality of our manuscript.

Supplementary Material

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

Supplementary Figure 1 | (A) Expression levels of model genes between tumor and normal tissues. (B) Survival map for OS, DSS, PFS, and disease free interval (DFI) with hazard ratios and cox p-values. (C) Kaplan–Meier survival curves of model genes.

Supplementary Figure 2 | (A,B) Correlations between drug sensitivity and model gene expression in the CTRP and GDSC databases. (C,D) Networks of enrichment terms based on similarity with best p-values.

Supplementary Figure 3 | Mutational profile of model genes. (A) Oncoprint mutation status, (B) lollipop charts, and (C) copy number variations of model genes. (D) Distinction of molecular features between mutant and wild-type groups (p < 0.01).

Footnotes

1. ^ http://gepia2.cancer-pku.cn
2. ^ http://tide.dfci.harvard.edu
3. ^ https://www.cbioportal.org

References

Bader, G. D., and Hogue, C. W. V. (2003). An automated method for finding molecular complexes in large protein interaction networks. BMC Bioinformatics 4:2. doi: 10.1186/1471-2105-4-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Bent, M. J. V. D., Brandes, A. A., Taphoorn, M. J. B., Kros, J. M., Kouwenhoven, M. C. M., Delattre, J.-Y., et al. (2013). Adjuvant procarbazine, lomustine, and vincristine chemotherapy in newly diagnosed anaplastic oligodendroglioma: long-term follow-up of EORTC brain tumor group study 26951. J. Clin. Oncol. 31, 344–350. doi: 10.1200/JCO.2012.43.2229

PubMed Abstract | CrossRef Full Text | Google Scholar

Brat, D. J., Verhaak, R. G. W., Aldape, K. D., Yung, W. K. A., Salama, S. R., Cooper, L. A. D., et al. (2015). Comprehensive, integrative genomic analysis of diffuse lower-grade gliomas. N. Engl. J. Med. 372, 2481–2498. doi: 10.1056/NEJMOA1402121

PubMed Abstract | CrossRef Full Text | Google Scholar

Brioschi, S., Wang, W.-L., Peng, V., Wang, M., Shchukina, I., Greenberg, Z. J., et al. (2021). Heterogeneity of meningeal B cells reveals a lymphopoietic niche at the CNS borders. Science 373:eabf9277. doi: 10.1126/SCIENCE.ABF9277

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Cairncross, G., Wang, M., Shaw, E., Jenkins, R., Brachman, D., Buckner, J., et al. (2013). Phase III trial of chemoradiotherapy for anaplastic oligodendroglioma: long-term results of RTOG 9402. J. Clin. Oncol. 31, 337–343. doi: 10.1200/JCO.2012.43.2674

PubMed Abstract | CrossRef Full Text | Google Scholar

Cugurra, A., Mamuladze, T., Rustenhoven, J., Dykstra, T., Beroshvili, G., Greenberg, Z. J., et al. (2021). Skull and vertebral bone marrow are myeloid cell reservoirs for the meninges and CNS parenchyma. Science 373:eabf7844. doi: 10.1126/SCIENCE.ABF7844

PubMed Abstract | CrossRef Full Text | Google Scholar

Hänzelmann, S., Castelo, R., and Guinney, J. (2013). GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics 14:7. doi: 10.1186/1471-2105-14-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Hinshaw, D. C., and Shevde, L. A. (2019). The tumor microenvironment innately modulates cancer progression. Cancer Res. 79, 4557–4566. doi: 10.1158/0008-5472.CAN-18-3962

PubMed Abstract | CrossRef Full Text | Google Scholar

Hu, W.-M., Yang, Y.-Z., Zhang, T.-Z., Qin, C.-F., and Li, X.-N. (2020). LGALS3 is a poor prognostic factor in diffusely infiltrating gliomas and is closely correlated with CD163+ tumor-associated macrophages. Front. Med. 7:182. doi: 10.3389/FMED.2020.00182

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, R., Mao, M., Lu, Y., Yu, Q., and Liao, L. (2020). A novel immune-related genes prognosis biomarker for melanoma: associated with tumor microenvironment. Aging 12, 6966–6980.

Google Scholar

Li, T., Forbes, M. E., Fuller, G. N., Li, J., Yang, X., and Zhang, W. (2020). IGFBP2: integrative hub of developmental and oncogenic signaling network. Oncogene 39, 2243–2257. doi: 10.1038/S41388-020-1154-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, T., Wernersson, R., Hansen, R. B., Horn, H., Mercer, J., Slodkowicz, G., et al. (2017). A scored human protein-protein interaction network to catalyze genomic interpretation. Nat. Methods 14, 61–64. doi: 10.1038/NMETH.4083

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, Y., Deng, G., Qi, Y., Zhang, H., Gao, L., Jiang, H., et al. (2020). Bioinformatic profiling of prognosis-related genes in malignant glioma microenvironment. Med. Sci. Moni. 26:e924054. doi: 10.12659/MSM.924054

PubMed Abstract | CrossRef Full Text | Google Scholar

Miao, Y.-R., Zhang, Q., Lei, Q., Luo, M., Xie, G.-Y., Wang, H., et al. (2020). ImmuCellAI: a unique method for comprehensive T-Cell subsets abundance prediction and its application in cancer immunotherapy. Adv. Sci. 7:1902880. doi: 10.1002/ADVS.201902880

PubMed Abstract | CrossRef Full Text | Google Scholar

Mondal, B., Patil, V., Shwetha, S. D., Sravani, K., Hegde, A. S., Arivazhagan, A., et al. (2017). Integrative functional genomic analysis identifies epigenetically regulated fibromodulin as an essential gene for glioma cell migration. Oncogene 36, 71–83. doi: 10.1038/ONC.2016.176

PubMed Abstract | CrossRef Full Text | Google Scholar

Nduom, E. K., Weller, M., and Heimberger, A. B. (2015). Immunosuppressive mechanisms in glioblastoma. Neuro Oncol. 17, vii9–vii14. doi: 10.1093/NEUONC/NOV151

PubMed Abstract | CrossRef Full Text | Google Scholar

Newman, A. M., Liu, C. L., Green, M. R., Gentles, A. J., Feng, W., Xu, Y., et al. (2015). Robust enumeration of cell subsets from tissue expression profiles. Nat. Methods 12, 453–457. doi: 10.1038/NMETH.3337

PubMed Abstract | CrossRef Full Text | Google Scholar

Phillips, A. C., Boghaert, E. R., Vaidya, K. S., Mitten, M. J., Norvell, S., Falls, H. D., et al. (2016). ABT-414, an antibody drug conjugate targeting a tumor-selective EGFR epitope. Mol. Cancer Ther. 15, 661–669. doi: 10.1158/1535-7163.MCT-15-0901

PubMed Abstract | CrossRef Full Text | Google Scholar

Rees, M. G., Seashore-Ludlow, B., Cheah, J. H., Adams, D. J., Price, E., Gill, S., et al. (2016). Correlating chemical sensitivity and basal gene expression reveals mechanism of action. Nat. Chem. Biol. 12, 109–116. doi: 10.1038/NCHEMBIO.1986

PubMed Abstract | CrossRef Full Text | Google Scholar

Robinson, G. W., Orr, B. A., and Gajjar, A. (2014). Complete clinical regression of a BRAF V600E-mutant pediatric glioblastoma multiforme after BRAF inhibitor therapy. BMC Cancer 14:258. doi: 10.1186/1471-2407-14-258

PubMed Abstract | CrossRef Full Text | Google Scholar

Sanai, N., Alvarez-Buylla, A., and Berger, M. S. (2005). Neural stem cells and the origin of gliomas. N. Engl. J. Med. 353, 811–822. doi: 10.1056/NEJMRA043666

PubMed Abstract | CrossRef Full Text | Google Scholar

Shannon, P., Markiel, A., Ozier, O., Baliga, N. S., Wang, J. T., Ramage, D., et al. (2003). Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 13, 2498–2504. doi: 10.1101/GR.1239303

PubMed Abstract | CrossRef Full Text | Google Scholar

Stark, C., Breitkreutz, B.-J., Reguly, T., Boucher, L., Breitkreutz, A., and Tyers, M. (2006). BioGRID: a general repository for interaction datasets. Nucleic Acids Res. 34, 535–539. doi: 10.1093/NAR/GKJ109

PubMed Abstract | CrossRef Full Text | Google Scholar

Stefano, A. L. D., Fucci, A., Frattini, V., Labussiere, M., Mokhtari, K., Zoppoli, P., et al. (2015). Detection, characterization and inhibition of FGFR-TACC fusions in IDH wild type glioma. Clin. Cancer Res. 21, 3307–3317. doi: 10.1158/1078-0432.CCR-14-2199

PubMed Abstract | CrossRef Full Text | Google Scholar

Suzuki, H., Aoki, K., Chiba, K., Sato, Y., Shiozawa, Y., Shiraishi, Y., et al. (2015). Mutational landscape and clonal architecture in grade II and III gliomas. Nat. Genet. 47, 458–468. doi: 10.1038/NG.3273

PubMed Abstract | CrossRef Full Text | Google Scholar

Szklarczyk, D., Gable, A. L., Lyon, D., Junge, A., Wyder, S., Huerta-Cepas, J., et al. (2019). STRING v11: protein-protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets. Nucleic Acids Res. 47, D607–D613. doi: 10.1093/NAR/GKY1131

PubMed Abstract | CrossRef Full Text | Google Scholar

Tang, Z., Kang, B., Li, C., Chen, T., and Zhang, Z. (2019). GEPIA2: an enhanced web server for large-scale expression profiling and interactive analysis. Nucleic Acids Res. 47, W556–W560. doi: 10.1093/NAR/GKZ430

PubMed Abstract | CrossRef Full Text | Google Scholar

Tönjes, M., Barbus, S., Park, Y. J., Wang, W., Schlotter, M., Lindroth, A. M., et al. (2013). BCAT1 promotes cell proliferation through amino acid catabolism in gliomas carrying wild-type IDH1. Nat. Med. 19, 901–908. doi: 10.1038/NM.3217

PubMed Abstract | CrossRef Full Text | Google Scholar

Türei, D., Korcsmáros, T., and Saez-Rodriguez, J. (2016). OmniPath: guidelines and gateway for literature-curated signaling pathway resources. Nat. Methods 13, 966–967. doi: 10.1038/NMETH.4077

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, H., Song, X., Huang, Q., Xu, T., Yun, D., Wang, Y., et al. (2019). LGALS3 promotes treatment resistance in glioblastoma and is associated with tumor risk and prognosis. Cancer Epidemiol. Biomarkers Prev. 28, 760–769. doi: 10.1158/1055-9965.EPI-18-0638

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, Y., Wahafu, A., Wu, W., Xiang, J., Huo, L., Ma, X., et al. (2021). FABP5 enhances malignancies of lower-grade gliomas via canonical activation of NF-κB signaling. J. Cell. Mol. Med. 25, 4487–4500. doi: 10.1111/JCMM.16536

PubMed Abstract | CrossRef Full Text | Google Scholar

Wick, W., Platten, M., Meisner, C., Felsberg, J., Tabatabai, G., Simon, M., et al. (2012). Temozolomide chemotherapy alone versus radiotherapy alone for malignant astrocytoma in the elderly: the NOA-08 randomised, phase 3 trial. Lancet Oncol. 13, 707–715. doi: 10.1016/S1470-2045(12)70164-X

CrossRef Full Text | Google Scholar

Yan, H., Parsons, D. W., Jin, G., McLendon, R., Rasheed, B. A., Yuan, W., et al. (2009). IDH1 and IDH2 mutations in gliomas. N. Engl. J. Med. 360, 765–773. doi: 10.1056/NEJMOA0808710

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, S., Liu, T., Cheng, Y., Bai, Y., and Liang, G. (2019). Immune cell infiltration as a biomarker for the diagnosis and prognosis of digestive system cancer. Cancer Sci. 110, 3639–3649. doi: 10.1111/CAS.14216

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, W., Soares, J., Greninger, P., Edelman, E. J., Lightfoot, H., Forbes, S., et al. (2012). Genomics of drug sensitivity in cancer (GDSC): a resource for therapeutic biomarker discovery in cancer cells. Nucleic Acids Res. 41, 955–961. doi: 10.1093/NAR/GKS1111

PubMed Abstract | CrossRef Full Text | Google Scholar

Yin, W., Zhu, H., Tan, J., Xin, Z., Zhou, Q., Cao, Y., et al. (2021). Identification of collagen genes related to immune infiltration and epithelial-mesenchymal transition in glioma. Cancer Cell Int. 21:276. doi: 10.1186/S12935-021-01982-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, C., Shen, L., Qi, F., Wang, J., and Luo, J. (2020). Multi-omics analysis of tumor mutation burden combined with immune infiltrates in bladder urothelial carcinoma. J. Cell. Physiol. 235, 3849–3863. doi: 10.1002/JCP.29279

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, N., Dai, Z., Wu, W., Wang, Z., Cao, H., Zhang, Y., et al. (2021). The predictive value of monocytes in immune microenvironment and prognosis of glioma patients based on machine learning. Front. Immunol. 12:656541. doi: 10.3389/fimmu.2021.656541

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, Z., Zhang, K.-N., Wang, Q., Li, G., Zeng, F., Zhang, Y., et al. (2021). Chinese glioma genome atlas (CGGA): a comprehensive resource with functional genomic data from chinese glioma patients. Genomics Proteomics Bioinformatics S1672-0229, 00045-0. doi: 10.1016/J.GPB.2020.10.005

PubMed Abstract | CrossRef Full Text | Google Scholar