A novel defined risk signature of interferon response genes predicts the prognosis and correlates with immune infiltration in glioblastoma

Abbreviations: AUC: area under the ROC curve; CGGA: the Chinese Glioma Genome Atlas; Coef: coefficient; EMT: epithelial-mesenchymal transition; GBM: glioblastoma; GSVA: gene set variation analysis; HR: hazard ratio; IFN: interferon; IHC: immunohistochemical; LASSO: the least absolute shrinkage and selection operator; OS: overall survival; PCA: principal components analysis; PPI: protein-protein interaction; SOCS3: suppressor of cytokine signaling 3; STAT2: signaling transducer and activator of transcription 2; STAT3: signaling transducer and activator of transcription 3; TCGA: the Cancer Genome Atlas; TMZ: temozolomide; TSC: tumor stem-like cell

1.   Introduction

Glioblastoma (GBM) is the most prevalent and lethal primary tumor in the central nervous system, whose 5-year survival rate is only 5.5% ^[1]. Despite GBM patients receive standard therapy including surgical resection and subsequent radiotherapy and chemotherapy, their median survival period is only 15 months ^[2]. Therefore, it is urgent to investigate and introduce novel treatment modalities to cure this malignancy. Chemotherapy is a vital process in treating GBM ^[3]. Temozolomide (TMZ), an alkylating agent, is the standard GBM chemotherapy, but response remains unsatisfactory owing to resistance ^[4]. Recent reports reveal interferons (IFNs) can re-sensitize glioma cells which are resistant to TMZ ^[5,6,7], thus combination of TMZ and IFNs could be an optimal method to cure GBM.

The IFNs are cytokines with diverse functions such as regulation of cell proliferation and immunity; therefore, IFNs have been studied as antitumor agents in many cancers, including glioma. IFNs family contains type I (IFN-α and IFN-β) and type II (IFN-γ) members. Each IFNs type has distinct antitumor effects: type I IFNs inhibit malignant cells or tumor stem-like cells (TSCs) proliferation by inducing apoptosis or interrupting cell cycle ^[8,9,10], and can also modulate tumor immune system through inactivating regulatory T cells or increasing survival of memory cytotoxic T cells ^[11,12]; type II IFNs disrupt cell cycle of tumor cells and inhibit angiogenesis ^[13]. Although IFNs have been proved to decrease the survival of glioma cell lines, only a small number of GBM patients can benefit from it as reported by clinical trials ^[13,14,15]. This situation suggests the effect of IFNs in GBM maybe controversial.

To data, increasing researches declare IFNs related genes can predict patient prognosis and response to IFN therapy ^[16,17,18]. The pleiotropic effects of IFNs are mediated by cell-specific expression of nearly 500 IFNs-stimulated genes. Nonetheless, no existing study has comprehensively screened and studied IFNs response genes as the risk signature for GBM prognosis. And the availability of high-throughput transcriptome GBM datasets bring the opportunity to develop and validate robust risk model. In present study, 3 IFNs-stimulated gene [signaling transducer and activator of transcription 3 (STAT3), signaling transducer and activator of transcription 2 (STAT2) and suppressor of cytokine signaling 3 (SOCS3)] signature was an independent biomarker for GBM prognosis. Pathway enrichment analysis and immune infiltration analysis revealed several immune modules and immune cells were enriched in higher risk GBM patients. These findings of this study may provide novel references for further investigating the usage of IFNs in GBM patients.

2.   Materials and methods

2.1.   Datasets and interferon response genes acquisition

Transcriptional profiles with clinical survival information of GBM patients were retrieved from the Chinese Glioma Genome Atlas (CGGA) ^[19] and the Cancer Genome Atlas (TCGA) ^[20]. CGGA GBM cohort 1, DataSet ID: mRNAseq_693, containing 237 GBM patients were used as the training cohort, and 152 GBM samples from TCGA and 137 GBM patients from CGGA cohort 2, DataSet ID: mRNAseq_325, were utilized as the validation cohorts. We downloaded IFN response gene sets from the molecular signature database of GSEA (http://www.gsea-msigdb.org/gsea/msigdb/). Ten IFN response gene sets (summarized in Table S1) and 491 IFN response genes in total were used in this study.

2.2.   Role of IFN response gene expression in GBM patient survival

Uni-variate Cox analysis was done and 1951 genes had significant hazard ratio (HR) related to GBM patient overall survival (OS). Then, we intersected these genes with 491 IFN response genes. Based on this analysis, 54 IFN response genes were related to patient OS and were used in subsequent study. Survival curves were performed by Kaplan–Meier analysis with Mantel-Cox log-rank test for statistical examination. Meta-analysis was done through PrognoScan (http://dna00.bio.kyutech.ac.jp/PrognoScan/) to further confirm the prognostic value of chosen candidate IFN response genes. PrognoScan is a web tool evaluating potential tumor markers and therapeutic targets ^[21]. A p-value less than 0.05 was regarded as statistically significant.

2.3.   Protein-protein interaction (PPI) network

The 54 survival-related IFN response genes were applied to construct PPI network by using the String website (https://string-db.org/). All parameters were default set. This raw protein-protein interaction network was further analyzed by Cytoscape (v.3.8.2). We adopted the cytoHubba plugin in Cytoscape to explore important nodes in the network. cytoHubba predicts and explores important nodes and subnetworks in a given PPI network by several topological algorithms. All nodes were ranked by degree, and the top 10 nodes were picked out and used in further analysis.

2.4.   Identification and validation of interferon response subtypes

The former 10 candidate IFN response genes were selected for consensus clustering by using R package ConsensusClusterPlus (v.1.56.0) to reveal robust clusters. The consensus heatmap and cumulative distribution function were applied to access the optimal K. Furthermore, principal components analysis (PCA) conducted by R package FactoMineR (v.2.4) and expression heatmap of these candidate genes were used to validate the difference between subtypes. Finally, the OS was also compared to confirm the clinical difference between these subtypes.

2.5.   Construction and validation of interferon response gene signature

These 10 candidate IFN response genes were further screened by the least absolute shrinkage and selection operator (LASSO) regression analysis using R package glmnet (v.4.1-2). Three genes with regression coefficient (Coef) were selected, and the risk score for each sample was calculated by the formula: Risk score$= {\sum }_{i}^{n}Coef \ i\times Expression \ i$. According to the median risk score, GBM patients were divided into higher and lower risk group. The uni-variate and multi-variate Cox Analysis were used to evaluating the prognosis predictor effect of this IFN response gene risk model. Additionally, we estimated the specificity and sensitivity of the model by calculating area under the ROC curve (AUC). CGGA GBM cohort 1 was the training cohort, and TCGA GBM dataset and CGGA GBM cohort 2 were the validation cohort to test our risk model.

2.6.   Correlation between DNA methylation and GBM patient survival

The relationship between OS and the methylation level of STAT2, STAT3 and SOCS3 in GBM was retrieved from MethSurv (https://biit.cs.ut.ee/methsurv/). It is a web tool to perform multivariable survival analysis using DNA methylation data from TCGA database ^[22].

2.7.   Functional enrichment analysis

To explore the difference on biological pathways between higher and lower risk group, gene set variation analysis (GSVA) was carried out by R package GSVA (v.1.40.1), and method = 'ssgsea' was set. The pathway related gene sets are in the "h.all.v7.2.symbols" document from the MSigDB database ^[23].

2.8.   Immune infiltration analysis

We used the TIMER2.0 (http://timer.cistrome.org/) platform and CIBERSORTx portal for immune infiltration analysis ^[24,25]. In this study, the "Immune Association" module in TIMER2.0 was taken to evaluate relationship between immune infiltrates and gene expression of STAT2, STAT3 and SOCS3 respectively. Then, the relationship between immune infiltration and risk scores was also investigated. CIBERSORTx web tool was utilized to evaluate the enrichment levels of 22 immune cells infiltration in GBM samples, and the enrichment scores predicted relative abundance of these immune cells. Then, the difference in immune cell infiltrates between higher and lower risk subgroup was checked by t test. In addition, expression level of immune checkpoints was also compared between these two risk groups using t test. The p-value < 0.05 was regarded as statistically significant.

2.9.   Immunohistochemical (IHC) staining data

The IHC data of normal brain tissues and GBM samples used in this study were downloaded from the Human Protein Atlas website (https://www.proteinatlas.org/) ^[26], which allows exploration of proteins, as well as global expression patterns, in major tissues, organs and tumors in the human body.

3.   Results

3.1.   Selection of candidate interferon response genes

The workflow diagraph of this study was portrayed in Figure 1. The first thing we had to do was selecting out candidate genes after obtaining the IFN response genes. Based on initial uni-variate Cox survival analysis, there were 1951 survival-related genes. Then, 54 survival-related IFN response genes were acquired by intersecting all 491 IFN response genes with the former survival genes (Figure 2A and Table S2). To get the most significant signatures among 54 genes, we constructed the PPI network and explored the core nodes ranked by score from high to low (Table S3). The top 10 genes, namely STAT3, CD44, CCL2, ICAM1, IRF1, NFKB1, MYD88, TNFAIP3, SOCS3 and STAT2, were deemed as the hub signatures and were used in downstream analysis (Figure 2B, C). The spearman correlation in gene expression was calculated to better understand the interactions among these 10 hub genes (Figure 2D). Majority of them were positively correlated with each other. STAT3 and STAT2 belong to signaling transducer and activator of transcription (STAT) family, which plays an important role in the progression of multiple cancers ^[27]. CD44, a complex trans-membrane adhesion glycoprotein, is a well-known marker of glioma TSCs and its mRNA expression level can predict low grade glioma (LGG) patient prognosis ^[28,29]. CCL2, a cytokine, is involved in immunoregulatory and inflammatory processes ^[30]. ICAM1, a transmembrane glycoprotein, can serve as a diagnostic tool for tumor prognosis ^[31]. IRF1, interferon regulatory factor 1, suppresses tumor cell growth and stimulates immune response against cancer ^[16]. NFKB1 is a member of the NFKB family and can suppress cancer progress ^[32]. MYD88 behaves as a vital signaling transducer in the interleukin-1 and Toll-like receptor signaling pathways, and has predictive prognostic value in glioma patients ^[33]. TNFAIP3 inhibits NFKB activation as well as TNF-mediated apoptosis, and is an independent factor for poor survival in carcinoma ^[34]. SOCS3 is one cytokine-inducible negative regulator of cytokine signaling, and is associated with glioma progression ^[35]. In general, these 10 selected IFN response genes are involved in immunity or related to tumor progression.

Figure 1.  Workflow diagraph of this study.

DownLoad: Full-Size Img PowerPoint

Figure 2.  Selection of candidate interferon response genes. (A) Venn diagraph of IFN response genes and GBM survival-related genes. Fifty-four genes were shared among them. (B) Node degree scores of top10 genes among 54 survival-related IFN response genes were ranked from high to low, which were calculated based on PPI network. (C) PPI network of top10 survival-related IFN response genes. Genes were colored by node degree scores as shown in (B). (D) Spearman correlation analysis of expression levels among these 10 survival-related IFN response genes from CGGA GBM cohort 1. The proportion of circle represented the spearman correlation coefficient.

DownLoad: Full-Size Img PowerPoint

3.2.   Consensus clustering identifies two subtypes in GBM

Before constructing the survival risk model, we had to evaluate the survival significance of these 10 selected core genes. Firstly, we did meta-analysis of prognostic value of these genes through PrognoScan, and the major genes were independent risk factors for glioma patients with HR > 1 and COX p-value < 0.05 (Table S4). Next, we clustered the GBM patients into different subgroups based on expression levels of these 10 genes, and compared patient survival difference in different subgroups. Consensus clustering analysis revealed that GBM patients from CGGA cohort 1 could be divided into two different clusters with the optimal k = 2 (Figure 3A). We further used the PCA to check the stability of consensus clustering result. The 1st principal component (PC) could distinguish cluster1 from cluster 2 well, which confirmed the solidity of clustering (Figure 3B). The clinical features between two clusters were shown in heatmap (Figure 3C). CGGA GBM cohort 1 had two fully distinct expression pattern of 10 IFN response genes, and cluster 2 was characterized by the increased expression of STAT3, CD44, CCL2, ICAM1, IRF1, NFKB1, MYD88, TNFAIP3, SOCS3 and STAT2 (Figure 3C). In addition, cluster 2 GBM patients had shorter OS period than cluster 1 (Figure 3D). We got the same results in the validation cohort TCGA GBM dataset (Figure 3E–H). The above results confirmed that these 10 selected core signatures harbored the survival predicting effect in GBM, and could be utilized in following risk model construction.

Figure 3.  Two subtypes in GBM patients with distinct clinical outcomes classified by 10 IFN response genes. Consensus clustering matrix at optimal k = 2 (A. CGGA GBM cohort 1, E. TCGA GBM cohort). PCA (PC1 vs. PC2) of 10 survival-related IFN response gene expression in GBM patients, and they can be divided into two clusters by PC1 (B. CGGA GBM cohort 1, F. TCGA GBM cohort). Heatmap of 10 survival-related IFN response gene expression in two clusters of GBM patients together with clinicopathological features, and cluster 2 GBM patients had high expression of these 10 genes (C. CGGA GBM cohort 1, G. TCGA GBM cohort). Survival analysis of two clusters in GBM patients (D. CGGA GBM cohort 1, H. TCGA GBM cohort). GBM patients in cluster 2 had shorter survival time than those in cluster 1.

DownLoad: Full-Size Img PowerPoint

3.3.   Construction of a risk model including 3 interferon response genes in GBM

For a more precise prediction of GBM prognosis by IFN response genes, the Cox regression analysis penalized by LASSO was applied (Figure 4A, B). Three genes including STAT3, STAT2 and SOCS3 were selected for constructing a risk model after cross validation.

Figure 4.  Analysis of three genes in the risk signature. (A) LASSO regression analysis was applied to calculate the coefficient of IFN response genes. (B) Three genes were chosen as active covariates to construct the risk model. (C) RNA expression difference between GBM and normal brain tissue. T test; ****p < 0.0001; ***p < 0.001. (D) The protein expression of 3 risk signatures in GBM and normal brain tissues using the Human Protein Atlas database. Survival analysis of 3 risk signatures in CGGA GBM cohort 1 (E) and TCGA GBM cohort (F).

DownLoad: Full-Size Img PowerPoint

Prior to constructing a risk model, we checked these chosen risk signatures in mRNA, protein and DNA methylation levels respectively. It found all of them were up-regulated in GBM compared to normal brain tissues in both mRNA and protein level (Figure 4C, D). In addition, GBM patients with high RNA expression of these 3 risk IFN response genes had bad clinical outcomes (Figure 4E, F). DNA methylation also plays a vital role in prognostic assessment and can serve as potential biomarker in glioma treatment ^[36]. We used MethSurv to determine the prognostic patterns of single CpG methylation of STAT2, STAT3 and SOCS3 in GBM respectively. It revealed that CpG site methylations of STAT2, STAT3 and SOCS3 indicated better prognosis with HR < 1, excluding cg25857307 and cg04439252 in STAT3 (Figure 5). As DNA methylation regulates gene expression by inhibiting the binding of transcription factors to DNA or by recruiting proteins involved in gene repression ^[37], survival effects of three risk genes from mRNA and DNA methylation level coincided with each other. Thus, all these 3 risk genes were positively related to malignancy and were associated with poor prognosis in GBM.

Figure 5.  The Kaplan-Meier survival of the promoter methylation of STAT2, STAT3 and SOCS3 in GBM patients.

DownLoad: Full-Size Img PowerPoint

Risk score of each sample was calculated with LASSO regression coefficient and corresponding mRNA expression level following the formulation: Risk score = 0.106993101 * STAT3 + 0.007796376 * STAT2 + 0.058969299 * SOCS3. GBM patients in CGGA cohort 1 were divided into higher and lower risk subgroups based on the median risk score of 0.82176. The difference of clinicopathological information and risk signatures mRNA expression between higher and lower group was detailed in Figure 6A, and higher risk group had high expression of STAT3, STAT2 and SOCS3. The risk score distribution and associated survival status revealed that higher risk GBM patients had more chance to death (Figure 6B). Furthermore, Kaplan–Meier survival curve confirmed that GBM patients in higher risk group had much lower overall survival probability than that in the lower risk ones (p = 0.00167, Figure 6C). Meanwhile, HR of risk score in GBM was 3.87 (95% confidence interval (CI): 1.88–7.95, p = 0) based on uni-variate Cox proportional hazards regression analysis (Figure 6D). And only when the factor passed uni-variate Cox analysis (p < 0.05), it would be brought into further multi-variate Cox proportional hazards regression analysis. The HR of risk score from multi-variate Cox analysis was 2.6 (95% CI: 1.18–5.9, p = 0.019; Figure 6E). Finally, the prognosis prediction of risk score alone had an AUC value of 0.646 (three years) based on ROC curve analysis, which was higher than Age and IDH state (Figure 6F, G). These results revealed that the risk model could accurately predict the prognosis of GBM patients in CGGA GBM cohort 1.

Figure 6.  Construction of risk model for GBM patients in CGGA cohort 1. (A) Heatmap of 3 risk signatures expression among the higher and lower risk GBM subgroups together with clinicopathological features. Higher risk group had relative higher expression of STAT2, STAT3 and SOCS3. (B) Risk score distribution together with survival status for GBM patients. (C) Survival analysis for the OS of GBM patients in the higher and lower risk subgroups based on risk score. Higher group confronted poor prognosis than the lower one. Uni-variate (D) and multi-variate (E) Cox regression analysis of the association between risk score, clinicopathological features and patient OS revealed the risk score was an independent factor for GBM patient survival. (F) ROC curves presented the predictive efficiency. (G) Multi-indicator ROC curves for risk score, age and IDH state. Risk score was a better predictor than others for GBM OS.

DownLoad: Full-Size Img PowerPoint

3.4.   Validation of the risk model using independent GBM cohorts

To verify the robustness of the risk model independently, we examined the three-gene prognostic model in two independent GBM data sets (TCGA GBM cohort and CGGA GBM cohort2). The TCGA GBM patients were classified into higher and lower risk subgroup based on median risk score = 1.99871. The higher risk scores were associated with more death in GBM patients (Figure 7A), and patients in the higher risk subgroup had an obvious reduced survival probability than the ones in the lower risk group (p = 0.0165, Figure 7B), implying the higher risk score carried more chance of death in GBM patients. The clinical features and expression levels of 3 genes were also compared between higher and lower risk group, detailed in heatmap (Figure 7C). STAT3, STAT2 and SOCS3 mRNA expression increased in higher risk GBM patients in the TCGA GBM cohort. Uni-variate and multi-variate Cox regression analysis were done to estimate whether the risk score could serve as an independent prognostic factor when compared with other clinicopathological features, including Gender, Age, IDH type and MGMT methylation state. In uni-variate Cox analysis, HR of the risk score was 15.82 in TCGA GBM cohort (95% CI: 2.79–89.72, p = 0.002; Figure 7D); HR of the risk score was 11.0 (95% CI: 1.35–90.2, p = 0.025; Figure 7E) from multi-variate Cox analysis. The ROC curve analysis was also carried out to evaluate the prognostic accuracy of the risk model. The prognostic prediction had the AUC value of 0.721 based on the risk score, which was higher than other clinicopathological features, including Age (0.643) and IDH wildtype (0.61) (Figure 7F, G). These analyses confirmed the risk model could serve as an independent prognostic factor in TCGA GBM cohort. In addition, the other independent cohort, namely CGGA cohort 2, was applied to validate the reliability of our risk model, and the detailed results were shown in Figure S1, which also confirmed our risk model was solid and effective.

Figure 7.  Validation of risk model for GBM patients in TCGA cohort. (A) Risk score distribution together with survival status for GBM patients. (B) Survival analysis for the OS of GBM patients in the higher and lower risk subgroups based on risk score. The lower risk group had better clinical outcome than higher risk group. (C) Heatmap of 3 risk signatures expression among the higher and lower risk GBM subgroups together with clinicopathological characteristics. Higher risk group had relative higher expression of STAT2, STAT3 and SOCS3. Uni-variate (D) and multi-variate (E) Cox regression analysis of the association between risk score, clinicopathological features and patient OS confirmed the risk score as an independent factor for GBM patient survival. (F) ROC curves presented the predictive efficiency. (G) Multi-indicator ROC curves for risk score, age and IDH state. Risk score was a better predictor than others for GBM survival.

DownLoad: Full-Size Img PowerPoint

3.5.   Functional analysis of higher and lower risk group

To further explore the biological difference between higher and lower risk subgroups, GSVA was done with Hallmark gene sets. Functional analysis results revealed that the higher risk group GBM samples were enriched in P53 pathway, hypoxia, angiogenesis and epithelial-mesenchymal transition (EMT), suggesting the malignancy phenotype of higher risk subgroup GBM patients (Figure 8). Tumor cells are marked with aerobic glycolysis, known as the Warburg effect; and glycolysis can promote tumor progression ^[38]. While lower risk GBM patients had low HALLMARK_Glycolysis pathway enrichment levels, the higher risk ones were obviously enriched in glycolysis pathway. Meanwhile, the most abundant pathways in higher risk GBM subgroup were correlated with immunity (e.g., IL2/STAT5 signaling, TNFA signaling via NFKB, IL6/JAK/STAT3 signaling, Complement, Inflammatory response), implying tumor immune microenvironment played vital role in the poor prognosis of higher risk GBM patients.

Figure 8.  Functional analysis between the two-risk subgroup of GBM. Gene set variation analysis of "h.all.v7.2.symbols" gene sets was carried out in higher and lower risk GBM groups. EMT, epithelial-mesenchymal transition.

DownLoad: Full-Size Img PowerPoint

3.6.   The risk model defines specific immune cell infiltration and immune checkpoint gene expression in GBM

Immune cells play a critical role in immune process, therefore the connection between risk subgroups and infiltrated immune cells in GBM were studied. Because TCGA-Mesenchymal subtype is featured as abundant immune cell infiltration, we firstly explored the TCGA-subtype presence of risk core genes and risk scores. The mRNA expressions of STAT2, STAT3 and SOCS3 were quite different in the three different GBM characteristic subtypes. While the Proneural subtype had the least expression level of STAT2 and STAT3, the Mesenchymal subtype expressed SOCS3 higher than the other two subtypes (Figure 9A). Furthermore, GBM patients with high risk score were mainly concentrated in the Mesenchymal subtype (Figure 9B). These results revealed our risk model were associated with specific GBM immune subtype.

Figure 9.  GBM TCGA-subtype presence of risk core genes and risk scores. (A) STAT2 and STAT3 expressed lower in the Proneural subtype, and SOCS3 was highly upregulated in the Mesenchymal subtype in TCGA and CGGA GBM samples. (B) Higher risk GBM samples were mainly concentrated in the Mesenchymal subtype. T test; ****p < 0.0001; ***p < 0.001; **p < 0.01; *p < 0.05; ^nsp ≥ 0.05.

DownLoad: Full-Size Img PowerPoint

Next, we explored how STAT2, STAT3 and SOCS3 correlate with immune cell infiltration by using the TIMER2.0 database. The immune infiltration trend of STAT2 and STAT3 was similar. They were all positively correlated with CD4+ T cell, B cell, neutrophil, macrophage and dendritic cell (Cor > 0.1 and p < 0.05) (Figure 10A, B). While STAT3 was negatively related to CD8+ T cell (Cor = -0.204 and p < 0.05), there was no obvious correlation between the expression of STAT2 and CD8+ T cell infiltration in GBM (Figure 10A, B). SOCS3 was positively correlated with neutrophil, macrophage and dendritic cell (Cor > 0.2 and p < 0.05), but negatively correlated with B cells (Cor = -0.174 and p < 0.05) (Figure 10C). Furthermore, CIBERSORTx tool was utilized to evaluate the relationship between the risk model and 22 immune cells infiltration. Among these infiltrating immune cells, we found the higher risk GBM samples were characterized with more Macrophage M0, Macrophage M2 and Neutrophils infiltration than lower risk ones (Figure 11A, B); there were also positive correlations between GBM risk scores and Macrophage M0, Macrophage M2 and Neutrophils infiltration respectively (Figure 11C, D), which was consistent with the former results (Figure 10). As we all known, Macrophage M2 in cancers accelerate tumor cell invasion and proliferation, and suppress antitumor immunity ^[39], thus higher risk GBM samples would confront short OS owing to more Macrophage M2 infiltration.

Figure 10.  Correlation between immune cell infiltrations and STAT2, STAT3 and SOCS3 expression in GBM patients. The correlation between the abundance of immune cells and the RNA expression of (A) STAT2, (B) STAT3 and (C) SOCS3 in GBM.

DownLoad: Full-Size Img PowerPoint

Figure 11.  Comparison of immune cell infiltration between the higher and lower risk groups in GBM. Histogram portrayed the distribution of infiltration ratio of 22 immune cells in higher and lower risk subgroups in GBM (A. CGGA GBM cohort 1, B. TCGA GBM cohort). T test; ****p < 0.0001; ***p < 0.001; **p < 0.01; *p < 0.05; ^nsp ≥ 0.05. Spearman correlation analysis between risk score and infiltration ratio of Macrophage M0, Macrophage M2 and Neutrophils in GBM (C. CGGA GBM cohort 1, D. TCGA GBM cohort).

DownLoad: Full-Size Img PowerPoint

It is promising to cure cancers by targeting immune checkpoints in immune cells ^[40]. We further compared the differential expression of 18 immune checkpoint molecules, including CD40, CD70, CD86, CD274, PDCD1LG2, HAVCR2, ICOSLG, ITGAM, ITGB2, LGALS9, LTBR, SIGLEC10, TNFRSF9, TNFRSF14, TNFRSF18, TNFSF4, TNFSF14 and TRAF3, between higher and lower groups (Figure 12). The results revealed that higher risk group had high expression of these 18 immune checkpoints in both CGGA and TCGA GBM datasets (p < 0.05), suggesting higher risk GBM patients had distinctive immune checkpoint molecules expression pattern.

Figure 12.  Differentially expressed immune checkpoint genes between the higher and lower risk GBM patients. The expression level of immune checkpoints was compared between the two-risk subgroup of GBM. T test; ****p < 0.0001; ***p < 0.001; **p < 0.01; *p < 0.05.

DownLoad: Full-Size Img PowerPoint

4.   Discussion

The development of sequencing technology enables us understand cancer clearer than before. With the expanding scope of studied tumors and genes, bioinformatics is becoming increasing essential. In this study, we explored the biological functions of a prognostic IFNs response gene risk model using bioinformatics analysis. Firstly, we identified 54 GBM survival-related IFNs response genes. Top 10 hub genes were selected out by calculating the connectivity degree in protein-protein interaction network; and based on these top 10 genes, GBM samples could be classified into 2 subgroups. GBM clinical data and mRNA expression profiles were combined and were further analyzed using LASSO regression model to build up a survival risk model, which was then validated by two independent GBM cohorts. Furthermore, enrichment analysis, infiltrating immune cells profiling and immune checkpoints expression analysis revealed the immunosuppression in higher-risk GBM patients.

In the past few decades, clinical features such as age, IDH state and MGMT state are widely used to predict the clinical outcome of GBM patients. However, these factors are usually with low predictive abilities by themselves solely owing to tumor heterogeneity ^[41]. As the explosion of high-throughput sequencing technology, mRNA expression level of numerous genes is deemed as predictive biomarkers in malignancies ^[29,42,43]. Because gene expression level can be regulated by multiple signaling pathways, a single gene is not accurate enough. Adopting the key genes that act in the same biological activity or signaling pathway to set up multi-gene-based models is vital to improve the predictive accuracy and may help to uncover new therapy targets ^[44,45]. Recently, IFNs have been reported to show a potent immunomodulatory capacity in GBM clinical treatment and their use is relatively safe ^[13,14], but the connections between their stimulated genes and GBM prognosis need to be further clarified. Our work constructed a 3-IFNs-response signature (STAT3, STAT2 and SOCS3) and revealed this risk model could predict GBM prognosis well.

Among the 3 genes in our risk model, STAT2 and STAT3 are both STAT family members. STAT2 is a vital mediator firstly identified in IFNs signaling in mammalian cells and is regarded as a pivotal component of adaptive immunity and antiviral response ^[46]. And it is required in signals transmission triggered by IFNs and in inducing T cells activation in host defense progress ^[47]. In present study, we found that expression of STAT2 was higher in GBM tumor samples than normal brain tissues (Figure 4C, D), and higher STAT2 expression was related to poor prognosis in GBM (Figure 4E, F). However, the function of STAT2 in glioma has rarely been studied ^[48]. STAT3 has been verified as an oncogene in glioma, and its activation is associated with clinically more aggressive behavior and bad outcomes ^[49]. Recent evidences prove the crucial role of STAT3 in tumor-induced immunosuppression ^[50]. The active STAT3 is the key transcription factor mediating tumor-associated macrophages/microglia (TAMs) M2 type polarization which is the immunosuppressive phenotype and accelerates tumorigenicity ^[51]. Previous studies also show aberrant STAT3 signaling in GBM promotes cell proliferation, angiogenesis and resistance to apoptosis ^[49]. The last signature gene was SOCS3. It is the endogenous inhibitor of JAK/STAT3 pathway and has anti-tumor activities in many malignancies ^[52]. But SOCS3 was revealed up-regulated in GBM compared to normal brain tissues and high SOCS3 RNA expression predicted poor prognosis in GBM (Figure 4C–F). Inhibition of SOCS3 can decrease the proliferation of GBM cells ^[53]. SOCS3 also plays a role in immunosuppression. Loss of SOCS3 in macrophages or microglia induces decreased infiltration of M2-TAMs, and results in increased CD8+ T cells and decreased regulatory T cells infiltration in glioma orthotopic transplantation mouse model ^[54]. In addition, SOCS3 takes part in radiation resistance in GBM and blocking its expression can sensitize tumor cells to radiotherapy ^[55].

Using the constructed risk model, GBM patients from three different cohorts could be divided into the higher and lower risk subgroups, and the lower risk patients had significantly better clinical outcome than the higher risk ones. GSVA results indicated the higher risk group was positively related to hypoxia, angiogenesis, EMT, aerobic glycolysis and several immune profiling signatures (including IL2/STAT5 signaling, TNFA signaling via NFKB, IL6/JAK/STAT3 signaling, Complement, Inflammatory response), which suggested higher risk group was strongly related to immunity. And our work further explored the relationship between risk score and immune cell infiltration because infiltrating immune cells are important regulator of tumor progression ^[56]. The higher risk GBM samples were enriched with Macrophage M0, Macrophage M2 and Neutrophils infiltration (Figure 11). Furthermore, the higher risk group had the higher expression of CD274 (namely PD-L1) and PDCD1LG2 (namely PD-L2). PD-L1 inhibit antitumor immunity by connecting with PDCD1 receptor on T cells ^[57]. Therefore, the higher risk score may represent an immunosuppressive microenvironment in GBM.

5.   Conclusions

In summary, this work uncovered that STAT3, STAT2 and SOCS3 expression levels increased in GBM samples compared to normal brain tissues and derived a risk model which could accurately predict the prognosis of GBM patients. Our work would further our understanding of IFNs stimulated genes and their role in signaling pathways, immune cells infiltration, and affecting immune checkpoints expression in GBM. The present study may help to guide more effective immunotherapy treatments against GBM.

Acknowledgements

This work was supported by the National Natural Science Foundation of China, Grant/Award Number: 81972350; Jiangsu Science and Education Strengthening Engineering Innovation Team Project, Grant/Award Number: CXTDA2017050; Medical Research Foundation of Jiangsu Health Commission, Grant/Award Number: H2019059; Medical Science, Technology Development Foundation of Nanjing, Grant/Award Number: ZDX16011; Postgraduate Research & Practice Innovation Program of Jiangsu Province, Grant/Award Number: SJCX19_0331 and KYCX20_1419. We thank Dr. Yanxiang Deng from the Department of Biomedical Engineering at Yale University and Dr. Liangtao Zheng at Peking University for preparing the manuscript.

Conflict of interest

The authors declare that they have no conflict of interest.