Development and Validation of a Tumor Mutation Burden–Related Immune Prognostic Model for Lower-Grade Glioma

Tumor mutation burden (TMB) is a useful biomarker to predict prognosis and the efficacy of immune checkpoint inhibitors (ICIs). In this study, we aimed to explore the prognostic value of TMB and the potential association between TMB and immune infiltration in lower-grade gliomas (LGGs). Somatic mutation and RNA-sequencing (RNA-seq) data were downloaded from the Cancer Genome Atlas (TCGA) database. TMB was calculated and patients were divided into high- and low-TMB groups. After performing differential analysis between high- and low-risk groups, we identified six hub TMB and immune-related genes that were correlated with overall survival in LGGs. Then, Gene Set Enrichment Analysis was performed to screen significantly enriched GO terms between the two groups. Moreover, an immune-related risk score system was developed by LASSO Cox analysis based on the six hub genes and was validated with the Chinese Glioma Genome Atlas dataset. Using the TIMER database, we further systematically analyzed the relationships between mutants of the six hub genes and immune infiltration levels, as well as the relationships between the immune-related risk score system and the immune microenvironment in LGGs. The results showed that TMB was negatively correlated with OS and high TMB might inhibit immune infiltration in LGGs. Furthermore, the risk score system could effectively stratify patients into low- and high-risk groups in both the training and validation datasets. Multivariate Cox analysis demonstrated that TMB was not an independent prognostic factor, but the risk score was. Higher infiltration of immune cells (B cells, CD4⁺ T cells, CD8⁺ T cells, neutrophils, macrophages, and dendritic cells) and higher levels of immune checkpoints (PD-1, CTLA-4, LAG-3, and TIM-3) were found in patients in the high-risk group. Finally, a novel nomogram model was constructed and evaluated to estimate the overall survival of LGG patients. In summary, our study provided new insights into immune infiltration in the tumor microenvironment and immunotherapies for LGGs.

Introduction

Gliomas are the most common malignant tumors in the central nervous system (1). Traditionally, gliomas are divided into grades I to IV, including astrocytoma, oligoastrocytoma, oligodendroglioma, and glioblastoma (GBM) (2, 3). The phrase “low-grade glioma,” which refers to grade I and II gliomas, is being gradually replaced by term the “lower-grade glioma (LGGs).” Lower-grade gliomas comprise WHO grades II and III astrocytomas, oligodendrogliomas, and oligoastrocytomas (4, 5). Some LGGs tend to progress to WHO grade IV GBM within months, whereas others remain stable for years. The survival of patients with LGGs ranging from 1 to 15 years is closely related to therapeutic sensitivity (6). Because of high intraobserver and interobserver variability, the histopathological classification of LGGs has not been adequate for outcome prediction. As a result, genetic classification is also very important for guiding clinical decision-making (7). According to the updated classification system by the World Health Organization (WHO) in 2016, LGGs can be divided into three subtypes based on the mutation status of isocitrate dehydrogenase 1 (IDH1) and the codeletion status of 1p/19q, in which both tumor phenotypes and genotypes are considered (8). Although this classification system of LGGs has been adopted for molecular diagnosis, the known molecular markers are currently very limited for explaining the prognosis of LGGs. Thus, further exploration of the genetic mechanism and identification of new biomarkers to predict the prognosis of LGGs is important to develop precise treatments.

The tumor microenvironment plays an important role in tumor growth and development. Tumor-infiltrating immune cells are a critical part of the tumor microenvironment that regulate tumor growth and invasion. These immune cells include lymphocytes, natural killer (NK) cells, dendritic cells, macrophages, neutrophils, and so on (9). With the deep understanding of the tumor microenvironment, immunotherapy recently has been developed to be a more effective treatment for aggressive cancers (10–13). For example, multiple therapeutic antibodies that block immune checkpoints, such as cytotoxic T lymphocyte associated antigen 4 (CTLA4) and programmed cell death protein 1 (PD1), showed great effects in treating non-small-cell lung cancer, kidney cancer, and melanoma (14). In addition to solid tumors, immune checkpoint inhibitors (ICIs) have also shown remarkable efficacy in some refractory hematologic malignancies, such as leukemia and lymphoma (15). However, immunotherapy could only benefit a subset of current cancer patients, as some cancers are immunotherapy insensitive, some patients failed to respond at all, and some effective cases in the early stage achieved a limited response followed by tumor progression or recurrence (16). Thus, it is very important to find more immunotherapy targets and elucidate a more detailed molecular mechanism of immunotherapy responsiveness. Previous studies have shown that tumor mutation burden (TMB) has become a useful biomarker across many cancer types to predict the efficacy of immune checkpoint blockade (ICB) (16, 17). TMB is usually defined as the total number of somatic protein-coding base substitutions, but in some cases, it also includes insertion/deletion mutations. Theoretically, TMB should be determined by whole exome sequencing (WES) which is not routinely used as a clinical tool owing to its greater cost and complexity (17). With the development of next generation sequencing technology (NGS), large NGS-targeted panels are sufficient to substitute WES for TMB estimation in the clinic (17). Wang et al. found that TMB has a close relationship with immune infiltration and the prognosis of various cancers (18). However, few studies have focused on the relationship between TMB and immune infiltration in LGGs. Hence, we performed a comprehensive analysis to further explore the relationship between TMB and the immune response based on WES and RNA-sequencing (RNA-seq) data.

Currently, with the rapid development of sequencing technique, WES and RNA-seq data from many cancers are available in many public databases, such as The Cancer Genome Atlas (TCGA) database and Gene Expression Omnibus (GEO) database. In the present study, we downloaded somatic mutation and RNA-seq data of LGG patients from TCGA database. Then, we analyzed the influence of TMB on the immune microenvironment, and developed an immune-related risk score system based on six TMB-related immune genes to classify patients into high- and low-risk groups with distinct prognoses. Moreover, the risk score system was validated in the Chinese Glioma Genome Atlas (CGGA) dataset, and a reliable predictive nomogram model was constructed to evaluate overall survival (OS) for LGG patients. We believe that the immune-related risk score system has potential in patient management and that the selected hub genes can serve as potential therapeutic biomarkers for LGGs.

Materials and Methods

Somatic Mutation, RNA-seq Data, and Immune-Related Genes

The somatic mutation and RNA-seq expression data and corresponding clinical data sheets of LGGs were obtained from the TCGA database (https://cancergenome.nih.gov/) and used as the training dataset. The “maftools” R package was used to analyze and visualize the somatic mutation data (19). For RNA-seq data, only coding RNAs with an expression raw count value >10 in more than 75% of samples were retained for further analysis. We downloaded the RNA-seq data and corresponding clinical information from CGGA database (http://www.cgga.org.cn/) as the validation dataset. Furthermore, a comprehensive immune-related gene set was extracted from the Immunology Database and Analysis Portal (ImmPort) database (https://immport.niaid.nih.gov) (20).

TMB Scores and Prognostic Analysis

In our study, the TMB score of each sample was calculated as the number of mutations/length of exons (30 Mb). Then, LGG samples were divided into high and low TMB groups according the median data and Kaplan–Meier analysis was conducted between the high and low TMB groups. Moreover, TMB levels were also assessed according to the WHO grade, histopathological types, and molecular subtype of LGGs.

Differentially Expressed Genes (DEGs) and Gene Set Enrichment Analysis (GSEA)

Using the “limma,” “edgeR,” and “DEseq2” packages in the R software, the DEGs between the high and low TMB groups were obtained with the following thresholds: |Fold change| > 1 and False Discovery Rate (FDR) < 0.05. The expression profiles of genes were converted to log₂(x+1) for further analysis. GSEA was performed between the high and low TMB groups using the JAVA8 platform. The reference gene set (c5.bp.v6.2.symbols.gm) was obtained from the MSigDB database (http://software.broadinstitute.org/gsea/msigdb/). Only enrichment pathways with a p < 0.05 and FDR < 0.25 were considered significant. In addition, the intersection between the DEGs and a list of 1,811 immune-related genes from the Immport database was selected for further analysis.

Construction and Validation of Immune-Related Risk Score System

After excluding patients with missing mutation information and survival time <30 days, 474 samples were subjected to subsequent analysis. The clinical characteristics of these patients are shown in Table 1. Then, Kaplan–Meier analysis of the selected genes was performed to screen for prognostic genes in the TCGA dataset, which is then validated in the CGGA dataset. In total, six coding genes were significantly related to OS. To identify the best prognostic value of these genes, Cox analysis with least absolute shrinkage and selection operator (LASSO) L1-penalty was performed using the “glmnet” R package (21, 22). Finally, an immune-related risk score system was constructed utilizing Cox regression coefficients to multiply the expression values of immune genes in each patient. By applying the “survminer” R package, LGG patients were divided into low- and high-risk groups based on the optimal cutoff point of their risk score. Kaplan–Meier survival analysis and the log-rank test were employed to evaluate the prognostic value of this system. Using the “survival ROC” R package, we depicted time-dependent receiver operating characteristic (ROC) curves to evaluate the sensitivity and specificity of the system. The risk score system was also validated with CGGA database.

TABLE 1

Table 1. Clinical characteristics of 474 LGG patients from TCGA cohort included in this study.

TIMER Database Analysis

The TIMER database (https://cistrome.shinyapps.io/timer) is a web tool for the comprehensive analysis and visualization of immune cells infiltration among 10,897 tumors from 32 cancer types (23). Six tumor-infiltrating immune subsets (B cells, CD4 T cells, CD8 T cells, macrophages, neutrophils, and dendritic cells) are included in TIMER. The abundance of the six immune cell types in the tumor microenvironment is estimated by a novel statistical method. Tumor immunological, clinical, and genomic features can be comprehensively explored in the TIMER database (23). In the “SCNA,” “Survival,” and “Gene” modules, the association between immune infiltration and somatic CNVs, clinical outcome, and gene expression, respectively, can be analyzed. Based on the three modules in the TIMER database, hub immune-related gene mutation types and the Kaplan–Meier analysis of immune infiltration cells were evaluated. Moreover, data on the immune infiltration levels of LGG samples were extracted from the TIMER database to calculate the correlation with the risk score system.

Development and Evaluation of the Nomogram

To validate whether the risk score system has an independent prediction value, univariate and multivariate Cox regression analyses were performed together with traditional clinical features (gender, age, pathologic stage, IDH1 and 1p/19q status, and radiation status). Then, according to the results of multivariate Cox analysis using the “rms” R package, a nomogram was generated to predict the 1, 3, and 5-years survival probability. Calibrations and ROC analyses were used to predict the accuracy of the nomogram. A concordance index (C-index) was applied to evaluate the discrimination of the system.

Statistical Analysis

Statistical analyses were conducted using the R software (version 3.5.1) and GraphPad Prism (version 7.0.0), and a p < 0.05 was considered statistically significant. The log-rank test was used in the Kaplan–Meier survival analysis. Student's t-test and Kruskal–Wallis test were employed in the two-group comparisons.

Results

Landscape of the LGG Mutation Profiles

In total, we analyzed the somatic mutation profiles of 509 LGG patients in the VCF format using the “maftools” package. As shown in the waterfall plot, IDH1, TP53, and ATRX mutations are the top three mutated genes in LGG samples, and IDH1 mutations are found in more than 75% of LGG samples (Figure 1A). Moreover, missense mutations are the most common mutation classification, single nucleotide polymorphisms (SNPs) showed a higher fraction in the variant type than insertion or deletion, and C>T was the most common single nucleotide variant (SNV) in LGGs (Figure 1B). Furthermore, the number of variants in each sample was calculated, and the mutation types are also shown in box plot with different colors for LGGs (Figure 1B). The co-occurrence and exclusive associations between mutated genes are shown in Figure 1C.

FIGURE 1

Figure 1. Summary of the of LGG mutation information. (A) Landscape of mutation profiles in LGG samples. Waterfall plot showing the mutation information for each gene. Different colors represented the different mutation types. (B) Classification of mutation types according to different categories and tumor mutation burden in specific samples. (C) The coincident and exclusive associations across the top 25 mutated genes. SNP, single nucleotide polymorphism; SNV, single nucleotide variant.

TMB Correlated With OS, WHO Grades, and Histopathological Types of LGGs

After calculating the TMB value of each sample, all patients were divided into high and low TMB groups using the median TMB as the cutoff point. Interestingly, patients in high-TMB group have an obviously shorter OS than those in the low TMB group with p < 0.0001 (Figure 2A). We also found that the TMB levels are positively correlated with WHO grades (Figure 2B). Moreover, astrocytoma has higher TMB levels than oligoastrocytoma (p = 0.0121) and oligodendroglioma (p = 0.0301), but whereas TMB levels between oligoastrocytoma and oligodendroglioma showed no statistical significance (Figure 2C). Furthermore, LGG samples with an IDH1-mutant type have lower TMB levels than IDH1-wildtype samples; however, the TMB levels in IDH1-mutant and 1p/19q codeletion samples have shown no statistical significance when compared with IDH1-mutant and 1p/19q non-codeletion samples (Figure 2D).

FIGURE 2

Figure 2. Prognostic value of TMB and a comparison of the gene expression profiles in high and low TMB groups. (A) Kaplan–Meier analysis showing that higher TMB levels correlated with a poor prognosis (p < 0.0001). (B) A higher TMB level was found in advanced grade tumors (p < 0.0001). (C) Comparison of the TMB levels with the different pathological types. Astrocytoma has higher TMB levels than oligoastrocytoma and oligodendroglioma (p = 0.0121 and p = 0.0301). (D) Lower TMB levels were associated with IDH1 mutant LGGs (p < 0.01). (E) Identification of TMB-related immune genes. (F) GSEA showed immune-related biological processes between the high- and low-risk groups.

Identification of Immune-Related DEGs Between Low and High TMB Groups

In total, 14,848 coding genes were identified, according to the annotation information provided in the Ensembl database (http://asia.ensembl.org/index.html). Based on the limma, edgeR, and DEseq2 algorithms, 99 genes were identified as being differentially expressed between the low and high TMB groups, with thresholds of |Fold change| > 1 and FDR < 0.05. Among these genes, seven immune-related DEGs were identified by the ImmPort database for further analysis (Figure 2E). Moreover, GSEA analysis of the LGG samples in the high (n = 235) and low (n = 234) TMB groups was performed. The results showed that LGGs in the high TMB group were significantly enriched for 273 biological processes, and the following four immune-related biological processes were selected: GO_SOMATIC_DIVERSIFICATION_OF_IMMUNE_RECEPTORS (normalized enrichment score (NES) = 1.76, size = 56), GO_SOMATIC_DIVERSIFICATION_OF_IMMUNE_RECEPTORS_VIA_GERMLINE_RECOMBINATION_WITHIN_A_SINGLE_LOCUS (NES = 1.72, size = 50), GO_IMMUNOGLOBULIN_PRODUCTION_INVOLVED_IN_IMMUNOGLOBULIN_MEDIATED_IMMUNE_RESPONSE (NES = 1.60, size = 43), and GO_RESPONSE_TO_RADIATION (NES = 1.58, size = 385) (p < 0.05 and FDR < 0.25) (Figure 2F). In contrast, LGGs in the low TMB group did not enrich for any immune-related biological processes.

Associations of Hub TMB-Related Immune Genes With Immune Infiltration

The Kaplan–Meier analysis results showed that six hub immune genes (BIRC5, CRLF1, GDF15, LTF, PRLHR, and TNFRSF11B) were highly associated with OS in LGGs. Higher expression levels of BIRC5, GDF15, LTF, and TNFRSF11B were positively correlated with poor prognosis, whereas higher expression levels of CRLF1 and PRLHR were negatively correlated with poor prognosis (Figure 3). Furthermore, we analyzed the underlying relationships between mutants of these hub genes with the immune infiltration microenvironment in LGGs based on the TIMER database. The results showed that mutants of these hub immune genes were related to the immune infiltration microenvironment in LGGs. Among them, mutants of BIRC5, CRLF1, and GDF15 inhibit the infiltration of several immune cells; in contrast, mutants of PRLHR promote the infiltration of several immune cells (Figure 4). Moreover, we found that the infiltration levels of immune cells (B cells, CD4⁺ T cells, CD8⁺ T cells, neutrophils, macrophages, and dendritic cells) was negatively correlated with the OS of patients in TCGA database (Figure 5). In addition, the correlation between the six hub genes and the level of immune cell infiltration was also analyzed in TIMER, and the results showed that the expression of BIRC5, GDF15, LTF, and TNFRSF11B were positively correlated with the infiltrating levels of immune cells, whereas the expression of CRLF1 and PRLHR were negatively correlated (Figure 6).

FIGURE 3

Figure 3. Kaplan–Meier analysis of the six hub TMB-related genes (BIRC5, CRLF1, GDF15, LTF, PRLHR, and TNFRSF11B) in TCGA database (A) and CGGA database (B).

FIGURE 4

Figure 4. Immune cell infiltration levels of the six hub TMB-related gene mutants. (A) BIRC5, (B) CRLF1, (C) GDF15, (D) LTF, (E) PRLHR, (F) TNFRSF11B. *p < 0.05, **p < 0.01, and ***p < 0.001.

FIGURE 5

Figure 5. Kaplan–Meier analysis reveals that lower immune cell (B cells, CD4⁺ T cells, CD8⁺ T cells, neutrophils, macrophages, and dendritic cells) infiltration levels are correlated with poor survival outcomes in LGGs (p < 0.05).

FIGURE 6

Figure 6. Correlation of hub gene expressions with immune infiltration levels in LGGs. The expression levels of BIRC5, GDF15, LTF, and TNFRSF11B were positively correlated with infiltrating levels of immune cells (B cell, CD8⁺ T cell, CD4⁺ T cell, macrophage, neutrophil, and dendritic cell). In contrast, the expression levels of CRLF1 and PRLHR were negatively correlated with infiltrating levels of immune cells (B cell, CD8⁺ T cell, CD4⁺ T cell, macrophage, neutrophil, and dendritic cell).

Construction and Validation of the Immune-Related Risk Score System

To build an immune-related risk score system in the TCGA cohort, LASSO Cox analysis was performed to select genes (Figures 7A,B). Six hub immune genes were included in the immune-related risk score system. The risk scores formula is as follows: risk score = 0.222 × exp_BIRC5 – 0.088 × exp_CRLF1 + 0.264 × exp_GDF15 + 0.060 × exp_LTF – 0.063 × exp_PRLHR+ 0.357 × exp_TNFRSF11B (Figure 7C). Then, risk scores were calculated for each sample. Patients in the TCGA cohort were divided into high- and low-risk groups by the median value. The Kaplan–Meier analysis indicated that patients in the high-risk group have a poorer prognosis than patients in the low-risk group (p = 4e−13) (Figure 7D). Moreover, the time-dependent ROC curve analysis demonstrated a promising prognostic prediction (1-year AUC = 0.89, 3-years AUC = 0.87, 5-years AUC = 0.76) (Figure 7E).

FIGURE 7

Figure 7. Construction and validation of the immune-related risk score system. (A,B) The six hub genes were all selected by LASSO Cox analysis in TCGA dataset. (C) Coefficient values for each gene. (D) Kaplan-Meier curves of the OS, risk scores distribution and survival status of each patient in the training cohort (TCGA dataset). (E) Time-dependent ROC curve analysis of the immune-related risk score system in the training cohort (TCGA dataset). (F) Kaplan-Meier curves of the OS, risk scores distribution and survival status of each patient in the validation cohort (CGGA dataset). (G) Time-dependent ROC curve analysis of the immune-related risk score system in the validation cohort (CGGA dataset).

To validate that the system had robust prognostic prediction ability, the same risk score formula was applied to the CGGA dataset. Using the median value of the risk score as a cutoff point, 274 LGG patients were divided into the high- and low-risk groups. Consistently, patients in the high-risk group had significantly poorer prognoses than patients in the low-risk group (Figure 7F). The ROC analysis also indicated that the system showed promising prognostic prediction (1-year AUC = 0.72, 3-years AUC = 0.78, 5-years AUC = 0.76) (Figure 7G).

Immune Infiltration Landscape in the Low- and High-Risk Groups

To explore the potential relationship between our risk score system and the immune infiltration microenvironment, we analyzed the correlation between the risk score and infiltrating immune cells using the “TIMER” tool. The result showed that the risk score was positively correlated with infiltrating immune cells including B cells, CD4⁺ T cells, CD8⁺ T cells, neutrophils, macrophages, and dendritic cells (Figure 8A). Moreover, patients in the high-risk group had higher proportions of immune cell infiltration, with a p < 0.001 (Figure 8B).

FIGURE 8

Figure 8. Correlation of the risk score with infiltrating immune cell proportions and immune checkpoints. (A) Correlation of the risk score with infiltrating immune cell proportions. Pearson's correlation coefficient values are shown in the heatmap. (B) Violin plot showing the immune cell proportions between low- and high-risk patients. (C) Correlation of the risk score with the expression of crucial immune checkpoints. Pearson's correlation coefficient values and level of significance were also showed. ***p < 0.001 and **p < 0.01. (D) The expression of different immune checkpoints between the low- and high-risk patients is shown in the violin plot.

Immune checkpoints play important roles in immune regulation, and their inhibitors are promising strategies for cancer treatment. Next, we explored the relationship between the risk score and the expression of critical immune checkpoints (PD-1, CTLA-4, LAG-3, TIM-3, and TIGIT). We found that the risk score showed a significantly positive correlation with PD-1, CTLA-4, LAG-3, and TIM-3 expression (Figure 8C). Further, we analyzed the expression of these immune checkpoints in the high- and low-risk groups, finding that patients in the high-risk group had higher expressions of PD-1, CTLA-4, LAG-3, and TIM-3 (p < 0.001), but lower expression of TIGIT (p = 0.001) (Figure 8D).

Construction and Evaluation of the Nomogram Model

Next, univariate and multivariate Cox analyses were performed to comprehensively analyze whether the TMB and immune-related risk score system were independent prognosis factors for LGG patients. The univariate Cox analysis results showed that the TMB and risk score were significantly associated with OS. However, multivariate Cox analysis along with clinicopathological variables indicated that only the risk score can serve as an independent prognostic factor in the TCGA dataset (HR: 1.92, 95% CI: 1.92 (1.50–2.47), p = 2.33e−07; Table 2).

TABLE 2

Table 2. Univariate and multivariate Cox regression analysis in TCGA.

To systematically predict the prognosis of LGGs, we constructed a nomogram model based on the risk score and two independent prognostic factors (age and IDH1 and 1p/19q status) in the TCGA dataset (Figure 9A). The C-index for the nomogram was 0.862 indicating a high discrimination ability, and calibration plot showed excellent concordance for the 1-, 3-, and 5-years predicted and actual OS probabilities (Figure 9B). Moreover, the ROC curve analysis also demonstrated a satisfactory prediction for sensitivity and specificity with a 1-year AUC of 0.921, 3-years AUC of 0.89, and 5-years AUC of 0.80 (Figure 9C).

FIGURE 9

Figure 9. Construction of the nomogram based on TCGA dataset. (A) A nomogram for the quantitative prediction of 1-, 2-, and 3-years survival for LGG patients based on the RS, age and IDH1 and 1p/19q status. (B) Calibration curves for the nomogram. (C) Time-dependent ROC curves for the 1-, 2-, and 3-years survival used to assess the nomogram model.

Discussion

Genomic variations are considered a major cause of LGGs (24). Although great efforts have been made in neurosurgery, radiotherapy, and chemotherapy, the survival of LGG patients still ranges widely. Recently, immunotherapy has shown promising results in the treatment of advanced or aggressive cancers (25). Although many efforts have been made for glioma immunotherapy, there is still a lack of reliable molecular biomarkers to distinguish patients with potential sensitivity to immunotherapy (26). Hence, it is particularly important to identify more immune-related prognostic biomarkers that can be potential therapeutic targets or can be used to screen immunotherapy-sensitive patients.

TMB is a novel biomarker to predict cancer immunotherapeutic response, which has been shown to be effective for many tumors, such as lung cancer (11), melanoma (27), and so on. Wang et al. found that a high TMB may promote cancer-testis antigen expression and inflammatory response, and patients with a higher TMB could gain a more favorable prognosis if treated with immunotherapy in a variety of cancers (18). However, few studies have focused on the prognostic role of TMB and the association between TMB and immune cell infiltration in LGGs. Thus, we aimed to explore the prognostic role of TMB and its potential association with immune infiltration in LGGs in this study.

Interestingly, the results of the Kaplan–Meier analysis showed that patients in high TMB group had a poorer prognosis, higher tumor grades, and advanced pathological subtypes. A recent bioinformatics study showed that the high tumor proliferative activity in high TMB patients may lead to a shorter OS of glioma, but more experiments are still needed to validate their findings (28). As high TMB patients usually benefit from immunotherapy (17, 29), higher TMB glioma patients may be able to achieve a better prognosis once immunotherapy is widely utilized in the treatment of glioma. Moreover, GSEA analysis showed that more immune-related biological processes were enriched in the high TMB group, indicating that a high TMB enhanced the immune phenotype. Then, we identified six hub immune genes that were highly associated with OS in LGGs. Among them, the expression levels of four genes (BIRC5, GDF15, LTF, and TNFRSF11B) were negatively correlated with OS, whereas two genes (CRLF1 and PRLHR) were positively correlated with OS. Analysis with the TIMER database showed that high immune cells (B cells, CD4⁺ T cells, CD8⁺ T cells, neutrophils, macrophages, and dendritic cells) infiltration was proven to be associated with a poorer prognosis in LGGs. Indeed, the expression levels of BIRC5, GDF15, LTF, and TNFRSF11B were positively correlated with infiltration levels of immune cells, whereas the expression levels of CRLF1 and PRLHR were negatively correlated. Furthermore, TIMER database analysis also showed that the mutants of these hub immune genes were correlated with immune infiltration in LGGs.

BIRC5 (survivin), a member of the inhibitor of apoptosis proteins family, can suppress cell apoptosis and regulate cell proliferation. BIRC5 is overexpressed in various tumors and has been found as a prognostic marker in gastric cancer (30), renal cell carcinoma (31), and breast cancer (32). However, the function of BIRC5 in LGGs has never been reported. GDF15 is a member of the growth differentiation factors (GDFs) subfamily that belongs to transforming growth factor beta superfamily (33). Peng et al. found that GDF15 might be able to regulate the expression of PD-L1, and targeting the GDF15/PD-L1 pathway might be promising for the treatment of GBM patients (33). The lactoferrin (LTF) gene, an iron-binding protein that is involved in innate and adaptive immunity, acts as a tumor suppressor gene in several tumors, including nasopharyngeal carcinoma (34), prostate carcinogenesis (35), and so on. In contrast, it appears to be a cancer-promoting factor in LGGs, even though further studies are needed to clarify the underling mechanism. TNFRSF11B is a cytokine receptor and belongs to the tumor necrosis factor (TNF) receptor superfamily. Deng et al. demonstrated that TNFRSF11B was a prognostic biomarker and related to worse survival in LGGs for the first time (36). CRLF1, cytokine receptor-like factor 1, stimulates neuronal growth and differentiation and has been proven to be involved in neuroprotection (37). However, its function in LGGs remains unclear and our study found that elevated CRLF1 expression was related to worse survival in LGGs. PRLHR (prolactin releasing hormone receptor), namely G-protein-coupled receptor 10, is the receptor for prolactin releasing peptide (PrRP). Previous studies found that it was associated with the regulation of feeding and energy balance (38). Su et al. found that PRLHR gene variants are protective factors in colorectal cancer patients of Chinese Han population (39). More studies are needed to explore the function of the PRHLR gene in LGGs.

In this study, we comprehensively studied the role of TMB in the regulation of immune phenotype in LGGs. Then, an immune-related risk score system was developed based on the TCGA dataset and validated with the CGGA dataset. This risk score system has favorable prognostic prediction ability, which is independent of traditional prognostic factors, such as IDH1 and 1p/19q status, age, and WHO grade. More importantly, we constructed a novel nomogram model integrated risk score with age and IDH1 and 1p/19q status to predict the OS of LGG patients. According to the risk score system and nomogram, clinicians can calculate an individual score for a patient and then can predict the 1-, 3-, and 5-years OS. Using the TIMER database, we found that higher infiltrating levels of B cells, CD4⁺ T cells, CD8⁺ T cells, neutrophils, macrophages, and dendritic cells are all negatively correlated with the OS of LGGs. Furthermore, we explored the correlation between the risk score and immune cell infiltration in LGGs. The results showed that patients in the high-risk group had higher infiltrating levels of B cells, CD4⁺ T cells, CD8⁺ T cells, neutrophils, macrophages, and dendritic cells than patients in the low-risk group. This indicated that higher levels of immune cell infiltration in the high-risk group may contribute to the poorer prognosis. These results also suggested that immune cell infiltration plays a crucial role in the progression of LGGs. Thus, the risk score system could be used as a predictor for prognosis and immune cell infiltration, and has good prospects for clinical application.

Currently, many clinical trials are evaluating the effect of ICIs in glioma (40). Furthermore, we analyzed the association between the risk score and the expression of critical immune checkpoints. The result showed that patients in high-risk group had higher PD-1, CTLA-4, LAG-3, and TIM-3 expression. The immunosuppressive microenvironment may lead to a poor prognosis in these patients. Thus, the patients in the high-risk group were more likely to benefit from ICIs.

Although a previous study had developed an IDH1-associated immune prognostic signature for LGGs (36), no study has systematically explored the relationships between TMB and immune infiltration and constructed a TMB-related risk score system in LGGs. Thus, our study provides new insights into the immune cell infiltration of tumor microenvironment and immunotherapies for LGGs. However, there were some limitations to our study. First, this is a retrospective study, and the results should be further confirmed by prospective studies. Moreover, more experiments are needed to elucidate the underlying mechanisms of the selected genes on immune cell infiltration in the future.

In conclusion, TMB was negatively correlated with OS and a high TMB might inhibit immune infiltration in LGGs. The TMB-related immune-related risk score system can divide patients into low- and high-risk groups with different outcomes and immunophenotypes. Moreover, the patients in the high-risk group are more likely to benefit from ICI treatment in the future. These findings may aid clinicians in identifying patients who are most likely to benefit from ICIs and to develop valuable personalized immunotherapy regimens for LGG patients in the future.

Data Availability Statement

The datasets for this study can be found in TCGA (http://cancergenome.nih.gov/) and CGGA (http://www.cgga.org.cn/).

Author Contributions

XJ, CR, and GT conceived and designed the study. WY wrote the manuscript. JT, ZX, QZ, CZ, and XF analyzed the results. ZW, YG, and ZJ performed the image visualization. All the authors approved the final manuscript.

Funding

This study was supported by the National Natural Science Foundation of China (no. 81472355) and the Hunan Provincial Science and Technology Department (no. 2014FJ6006).

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.

References

1. Lapointe S, Perry A, Butowski NA. Primary brain tumours in adults. Lancet. (2018) 392:432–46. doi: 10.1016/S0140-6736(18)30990-5

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Chen R, Smith-Cohn M, Cohen AL, Colman H. Glioma subclassifications and their clinical significance. Neurotherapeutics. (2017) 14:284–97. doi: 10.1007/s13311-017-0519-x

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Ma Q, Long W, Xing C, Chu J, Luo M, Wang HY, et al. Cancer stem cells and immunosuppressive microenvironment in glioma. Front Immunol. (2018) 9:2924. doi: 10.3389/fimmu.2018.02924

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Zeng WJ, Yang YL, Liu ZZ, Wen ZP, Chen YH, Hu XL, et al. Integrative analysis of DNA methylation and gene expression identify a three-gene signature for predicting prognosis in lower-grade gliomas. Cell Physiol Biochem. (2018) 47:428–39. doi: 10.1159/000489954

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Eckel-Passow JE, Lachance DH, Molinaro AM, Walsh KM, Decker PA, Sicotte H, et al. Glioma groups based on 1p/19q, IDH, and TERT promoter mutations in tumors. N Engl J Med. (2015) 372:2499–508. doi: 10.1056/NEJMoa1407279

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Patil V, Mahalingam K. A four-protein expression prognostic signature predicts clinical outcome of lower-grade glioma. Gene. (2018) 679:57–64. doi: 10.1016/j.gene.2018.08.001

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Wesseling P, Capper D. WHO 2016 classification of gliomas. Neuropathol Appl Neurobiol. (2018) 44:139–50. doi: 10.1111/nan.12432

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Louis DN, Perry A, Reifenberger G, von Deimling A, Figarella-Branger D, Cavenee WK, et al. The 2016 World health organization classification of tumors of the central nervous system: a summary. Acta Neuropathol. (2016) 131:803–20. doi: 10.1007/s00401-016-1545-1

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Domingues P, González-Tablas M, Otero Á, Pascual D, Miranda D, Ruiz L, et al. Tumor infiltrating immune cells in gliomas and meningiomas. Brain Behav Immun. (2016) 53:1–15. doi: 10.1016/j.bbi.2015.07.019

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Kruger S, Ilmer M, Kobold S, Cadilha BL, Endres S, Ormanns S, et al. Advances in cancer immunotherapy 2019—latest trends. J Exp Clin Cancer Res. (2019) 38:268. doi: 10.1186/s13046-019-1266-0

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Calles A, Aguado G, Sandoval C, Álvarez R. The role of immunotherapy in small cell lung cancer. Clin Transl Oncol. (2019) 21:961–76. doi: 10.1007/s12094-018-02011-9

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Johnston MP, Khakoo SI. Immunotherapy for hepatocellular carcinoma: current and future. World J Gastroenterol. (2019) 25:2977–89. doi: 10.3748/wjg.v25.i24.2977

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Xu S, Tang L, Li X, Fan F, Liu Z. Immunotherapy for glioma: current management and future application. Cancer Lett. (2020) 476:1–12. doi: 10.1016/j.canlet.2020.02.002

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Topalian SL, Taube JM, Anders RA, Pardoll DM. Mechanism-driven biomarkers to guide immune checkpoint blockade in cancer therapy. Nat Rev Cancer. (2016) 16:275–87. doi: 10.1038/nrc.2016.36

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Thanarajasingam G, Thanarajasingam U, Ansell SM. Immune checkpoint blockade in lymphoid malignancies. FEBS J. (2016) 283:2233–44. doi: 10.1111/febs.13668

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Braun DA, Burke KP, van Allen EM. Genomic approaches to understanding response and resistance to immunotherapy. Clin Cancer Res. (2016) 22:5642–50. doi: 10.1158/1078-0432.CCR-16-0066

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Chan TA, Yarchoan M, Jaffee E, Swanton C, Quezada SA, Stenzinger A, et al. Development of tumor mutation burden as an immunotherapy biomarker: utility for the oncology clinic. Ann Oncol. (2019) 30:44–56. doi: 10.1093/annonc/mdy495

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Wang X, Li M. Correlate tumor mutation burden with immune signatures in human cancers, BMC Immunol. (2019) 20:4. doi: 10.1186/s12865-018-0285-5

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Mayakonda A, Lin D-C, Assenov Y, Plass C, Koeffler HP. Maftools: efficient and comprehensive analysis of somatic variants in cancer. Genome Res. (2018) 28:1747–56. doi: 10.1101/gr.239244.118

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Bhattacharya S, Andorf S, Gomes L, Dunn P, Schaefer H, Pontius J, et al. ImmPort: disseminating data to the public for the future of immunology. Immunol Res. (2014) 58:234–9. doi: 10.1007/s12026-014-8516-1

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Goeman JJ. L1 penalized estimation in the Cox proportional hazards model. Biom J. (2010) 52:70–84. doi: 10.1002/bimj.200900028

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Tang G, Yin W. Development of an immune infiltration-related prognostic scoring system based on the genomic landscape analysis of glioblastoma multiforme. Front Oncol. (2020) 10:154. doi: 10.3389/fonc.2020.00154

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Li T, Fan J, Wang B, Traugh N, Chen Q, Liu JS, et al. TIMER: a web server for comprehensive analysis of tumor-infiltrating immune cells. Cancer Res. (2017) 77:e108–10. doi: 10.1158/0008-5472.CAN-17-0307

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Zhang CM, Brat DJ. Genomic profiling of lower-grade gliomas uncovers cohesive disease groups: implications for diagnosis and treatment. Chin J Cancer. (2016) 35:12. doi: 10.1186/s40880-015-0071-1

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Yang Y. Cancer immunotherapy: harnessing the immune system to battle cancer. J Clin Invest. (2015) 125:3335–7. doi: 10.1172/JCI83871

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Kamran N, Alghamri MS, Nunez FJ, Shah D, Asad AS, Candolfi M, et al. Current state and future prospects of immunotherapy for glioma. Immunotherapy. (2018) 10:317–39. doi: 10.2217/imt-2017-0122

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Forschner A, Battke F, Hadaschik D, Schulze M, Weißgraeber S, Han CT, et al. Tumor mutation burden and circulating tumor DNA in combined CTLA-4 and PD-1 antibody therapy in metastatic melanoma—results of a prospective biomarker study. J Immunother Cancer. (2019) 7:180. doi: 10.1186/s40425-019-0659-0

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Wang L, Ge J, Lan Y, Shi Y, Luo Y, Tan Y, et al. Tumor mutational burden is associated with poor outcomes in diffuse glioma. BMC cancer. (2020) 20:213. doi: 10.1186/s12885-020-6658-1

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Stenzinger A, Allen JD, Maas J, Stewart MD, Merino DM, Wempe MM, et al. Tumor mutational burden standardization initiatives: recommendations for consistent tumor mutational burden assessment in clinical samples to guide immunotherapy treatment decisions. Genes Chromosomes Cancer. (2019) 58:578–88. doi: 10.1002/gcc.22733

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Zou J, Liao X, Zhang J, Wang L. Dysregulation of miR-195-5p/-218-5p/BIRC5 axis predicts a poor prognosis in patients with gastric cancer. J Biol Regul Homeostat Agents. (2019) 33:1377–85. doi: 10.23812/19-146-A

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Luo Y, Shen D, Chen L, Wang G, Liu X, Qian K, et al. Identification of 9 key genes and small molecule drugs in clear cell renal cell carcinoma. Aging. (2019) 11:6029–52. doi: 10.18632/aging.102161

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Dai JB, Zhu B, Lin WJ, Gao HY, Dai H, Zheng L, et al. Identification of prognostic significance of BIRC5 in breast cancer using integrative bioinformatics analysis. Biosci Rep. (2020) 40:BSR20193678. doi: 10.1042/BSR20193678

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Peng H, Li Z, Fu J, Zhou R. Growth and differentiation factor 15 regulates PD-L1 expression in glioblastoma. Cancer Manage Res. (2019) 11:2653–61. doi: 10.2147/CMAR.S192095

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Zhang H, Feng X, Liu W, Jiang X, Shan W, Huang C, et al. Underlying mechanisms for LTF inactivation and its functional analysis in nasopharyngeal carcinoma cell lines. J Cell Biochem. (2011) 112:1832–43. doi: 10.1002/jcb.23101

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Porter CM, Haffner MC, Kulac I, Maynard JP, Baena-Del Valle JA, Isaacs WB, et al. Lactoferrin CpG Island hypermethylation and decoupling of mRNA and protein expression in the early stages of prostate carcinogenesis. Am J Pathol. (2019) 189:2311–22. doi: 10.1016/j.ajpath.2019.07.016

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Deng X, Lin D, Chen B, Zhang X, Xu X, Yang Z, et al. Development and validation of an IDH1-associated immune prognostic signature for diffuse lower-grade glioma. Front Oncol. (2019) 9:1310. doi: 10.3389/fonc.2019.01310

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Looyenga BD, Resau J, MacKeigan JP. Cytokine receptor-like factor 1 (CRLF1) protects against 6-hydroxydopamine toxicity independent of the gp130/JAK signaling pathway. PLoS ONE. (2013) 8:e66548. doi: 10.1371/journal.pone.0066548

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Varghese BV, Koohestani F, McWilliams M, Colvin A, Gunewardena S, Kinsey WH, et al. Loss of the repressor REST in uterine fibroids promotes aberrant G protein-coupled receptor 10 expression and activates mammalian target of rapamycin pathway. Proc Natl Acad Sci USA. (2013) 110:2187–92. doi: 10.1073/pnas.1215759110

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Su Q, Wang Y, Zhao J, Ma C, Wu T, Jin T, et al. Polymorphisms of PRLHR and HSPA12A and risk of gastric and colorectal cancer in the Chinese Han population. BMC Gastroenterol. (2015) 15:107. doi: 10.1186/s12876-015-0336-9

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Eissa IR, Bustos-Villalobos I, Ichinose T, Matsumura S, Naoe Y, Miyajima N, et al. The current status and future prospects of oncolytic viruses in clinical trials against melanoma, glioma, pancreatic, breast cancers. Cancers. (2018) 10:356. doi: 10.3390/cancers10100356

PubMed Abstract | CrossRef Full Text | Google Scholar