Glioma-Associated Stromal Cells Stimulate Glioma Malignancy by Regulating the Tumor Immune Microenvironment

Background: The glioma-associated stromal cell (GASC) is a recently identified type of cell in the glioma microenvironment and may be a prognostic marker for glioma. However, the potential mechanisms of GASCs in the glioma microenvironment remain largely unknown. In this work, we aimed to explore the mechanisms of GASCs in gliomas, particularly in high-grade gliomas (HGG).

Methods: We used glioma datasets from The Cancer Genome Atlas (TCGA) and the Chinese Glioma Genome Atlas (CGGA). We utilized the Single-sample Gene Set Enrichment Analysis (ssGSEA) algorithm to discriminate between patients with high or low GASC composition. The xCELL and CIBERSORT algorithms were used to analyze the composition of stromal cells and immune cells. Risk score and a nomogram model were constructed for prognostic prediction of glioma.

Results: We observed for the first time that the levels of M2 macrophages and immune checkpoints (PD-1, PD-L1, PD-L2, TIM3, Galectin-9, CTLA-4, CD80, CD86, CD155, and CIITA) were significantly higher in the high GASC group and showed positive correlation with the GASC score in all glioma population and the HGG population. Copy number variations of DR3 and CIITA were higher in the high-GASC group. THY1, one of the GASC markers, exhibited lower methylation in the high GASC group. The constructed risk score was an independent predictor of glioma prognostics. Finally, a credible nomogram based on the risk score was established.

Conclusions: GASCs stimulate glioma malignancy through the M2 macrophage, and are associated with the level of immune checkpoints in the glioma microenvironment. The methylation of THY1 could be used as prognostic indicator and treatment target for glioma. However, further studies are required to verify these findings.

Introduction

Glioma is the most common primary malignant tumor of the central nervous system, and it generally has a poor prognosis. The World Health Organization (WHO) classified gliomas into grades I-IV, with grades III and IV indicating high-grade gliomas (HGG) (1). The current treatments for HGG involve tumor resection, radiotherapy (RT), and temozolomide (TMZ), but this strategy has not yielded optimal effects (2).

Immunotherapy has been extensively studied for human malignant tumors in the past few years (3). However, due to the “immune-cold” phenotype and inner complexity of glioma (4), only a minority of glioma patients benefit from immune checkpoint (ICP) inhibitors (5). Researchers are deepening our understanding of the complex interactions between glioma and the immune system and trying to maximize the effectiveness of immunotherapy for glioma (6).

The glioma-associated stromal cell (GASC) is a recently identified important stromal cell in the glioma microenvironment, with potential value for prognostic prediction and therapeutic perspectives (7). The available evidence indicates that GASCs facilitate angiogenesis, invasion, and tumor growth (7). However, the potential mechanisms of GASCs remain largely unknown.

We aimed to identify the underlying mechanisms of GASCs in a glioma microenvironment, particularly in HGG. Figure 1 illustrates the workflow of the study.

FIGURE 1

Figure 1 Flow diagram of this investigation. HGG, high-grade glioma; DEGs, differentially expressed genes; 4 scores, stemness, mesenchymal-EMT, tumorigenic cytokine, and angiogenic activity scores; SM, somatic mutations; CNV, copy number variations; ICPs, immune checkpoints.

Materials and Methods

Glioma Datasets

The Cancer Genome Atlas (TCGA, https://portal.gdc.cancer.gov/) and Chinese Glioma Genome Atlas (CGGA, www.cgga.org.cn/) are public databases. The mRNA sequencing data and clinical information data for 702 glioma samples from TCGA and 693 glioma samples from CGGA were downloaded. Among these samples, 393 samples from TCGA and 504 samples from CGGA were high-grade glioma (HGG). The somatic mutation data for 666 glioma samples from TCGA were downloaded. The copy number variation data for 692 samples from TCGA were downloaded from the UCSC Xena Project database (http://xena.ucsc.edu/). For methylation analysis, the methylation data for 34 samples from CGGA and mRNA sequencing data for 325 samples were downloaded from CGGA.

ssGSEA Analysis

A gene set of GASC markers (Table S1) was obtained from Clavreul et al. (7). Enrichment scores for GASCs were separately calculated for each sample with the Single-sample Gene Set Enrichment Analysis (ssGSEA) algorithm. We also used the ssGSEA algorithm to calculate the stemness score (8), mesenchymal-epithelial-to-mesenchymal transition (EMT) score (9), tumorigenic cytokine score (10) and angiogenic activity score (11) based on the corresponding gene sets (Table S1). The “GSVA” R package (version 1.34.0) was applied to conduct an ssGSEA analysis.

Principle Component Analysis (PCA)

PCA was used to show the differentiation of high- and low-GASC groups and was visualized with the “ggfortify” R package (version 0.4.11).

Differential Analysis of Expressed Gene

We used Morpheus (https://software.broadinstitute.org/morpheus) to identify significantly differentially expressed genes (DEGs) between the high- and low-GASC groups. P < 0.05 and |log2 FC (fold-change)| ≥ 1 were selected as the cutoff values for statistically significant DEGs. A heatmap of DEG expression was produced by the “pheatmap” R package (version 1.0.12).

Functional Annotation

To reveal the probable biofunctions and signaling pathways that were correlated with the DEGs, we performed Gene Ontology (GO) annotations enrichment analysis, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis and enrichment analysis, and Gene Set Enrichment Analysis (GSEA) using the “clusterProfiler” (12) package (version 3.14.3) in R. Adjusted p < 0.05 was selected as the cutoff criterion.

xCELL Analysis and Cell Type Identification by Estimating Relative Subsets of RNA Transcripts (CIBERSORT) Analysis

xCell is an R package (version 1.1.0) that estimates the comprehensive levels of 64 cell types, which include 14 stromal cells. CIBERSORT can accurately quantify the abundance scores of 22 types of immune cells for each sample. We applied xCELL and CIBERSORT to separately calculate the abundance scores for stromal cells and immune cells in glioma samples.

Analysis of Somatic Mutations and Copy Number Variations

The somatic mutations of glioma samples from TCGA were calculated and visualized by the “Maftools” R package (version 2.2.10) (13). The copy number variations were visualized by the “ComplexHeatmap” R package (version 2.2.0).

Prediction of the Immunotherapy Response

The Tumor Immune Dysfunction and Exclusion (TIDE) algorithm (14) was employed to predict the clinical response of immune checkpoint inhibitors for each glioma sample.

Construction of Prognostic Model

The glioma datasets from TCGA and CGGA were used separately as a training dataset and validation dataset during the construction of the prognostic model. In the filtering process, least absolute shrinkage and selection operator (LASSO) regression analysis was applied to filter input parameters with p < 0.05. The input parameters included GASC score, GASC markers, immune checkpoints, stemness score, mesenchymal-EMT score, tumorigenic cytokine score, angiogenic activity score, stromal cell scores, and immune cell scores. Then, multivariate Cox regression analysis was conducted, and the risk score for glioma was computed via this formula: risk = score ${\sum }_{i=1}^n{\beta }_i\times X_i.X_i$ indicates the input parameter of multivariate regression analysis, and β_i represents the coefficient of X_i. Risk score and clinicopathological features were used to construct a prognostic model with uni- and multivariate Cox regression analysis. A nomogram was built to show the prognostic model. Receiver operating characteristic (ROC) curve analysis was conducted to evaluate the effect of the prognostic model in the training and validation datasets. “glmnet” (version 4.1), “rms” (version 6.1.0), and “timeROC” (version 0.4) R packages were used for the construction of prognostic model.

Statistical Analysis

All statistical analyses were performed utilizing R software (version 3.5.1), and statistical significance was set at p < 0.05. Comparisons between 2 continuous variables were evaluated by Student’s t-test and one-way ANOVA with ≥ 3 variables. Boxplots and bar charts were utilized to display these comparisons using the “ggplot” R package (version 3.3.3). The chi-square test or Fisher’s exact test were used for comparisons of categorized variables. The Kaplan-Meier approach was conducted for survival analysis, and the log-rank test was used to compare the overall survival (OS). Spearman correlation analysis was applied to evaluate two continuous variables, and the data were visualized with “ggplot” and “corrgram” (version 1.13) R packages. Univariate Cox regression was applied to identify potential predictors of survival, and the data were displayed with “forestplot” (version 1.10.1) R package.

Results

Identification of High- and Low-GASC Groups With ssGSEA

To analyze the potential mechanisms of GASCs in the glioma microenvironment, we obtained mRNA sequencing data for 702 samples from TCGA and 693 samples from CGGA, and then calculated the GASC score for each sample using the ssGSEA algorithm (Figure 1). Samples from TCGA and CGGA were classified separately into high- or low-GASC groups according to the median of the GASC score. Information for the high- and low-GASC groups is shown in Figures 2A, B and Table 1. Separate classification was also performed for 393 HGG samples from TCGA and 504 HGG from CGGA into high- and low-GASC groups using the same method. K-M curves were drawn, and the results revealed that a higher GASC score was associated with worse OS in all glioma population and the HGG population (p < 0.0001; Figures 2C–F). PCA showed robust differences in the expression portraits of the GASC markers between the high- and low-GASC groups (Figures 2G–J).

FIGURE 2

Figure 2 Identification of high- and low-GASC groups. (A, B) Heatmap of GASC markers in all glioma population (A for TCGA and B for CGGA). (C–F) Kaplan-Meier overall survival (OS) curves for samples in high- and low-GASC groups from all glioma population (C for TCGA and D for CGGA) and high-grade glioma population (E for TCGA and F for CGGA). (G–J) Principal component analysis of high- and low-GASC groups from all glioma population (G for TCGA and H for CGGA) and high-grade glioma population (I for TCGA and J for CGGA).

TABLE 1

Table 1 Correlations between GASC groups and clinical characteristics in glioma patients.

Enrichment Analysis of DEGs Between the High- and Low-GASC Groups

DEGs between the high- and low-GASC groups were identified with the Morpheus webtool (Figure S1). Functional enrichment revealed a significant association between DEGs and immune-related terms. Biological process (BP) terms enriched in the GO analysis included “lymphocyte chemotaxis” and “neutrophil activation” in all glioma population (Figures 3A, D and Table S2) and the HGG population (Figures 4A, D and Table S2). The “JAK-STAT signaling pathway” and “IL-17 signaling pathway” were enriched in the KEGG analysis in all glioma population (Figures 3B, E and Table S3) and the HGG population (Figures 4B, E and Table S3). GSEA analysis revealed immune-related terms such as “Antigen processing and presentation” and “PD-L1 expression and PD-1 checkpoint pathway” in all glioma population (Figures 3C, F and Table S6 ) and the HGG population (Figures 4C, F and Table S6).

FIGURE 3

Figure 3 Functional annotation of upregulated DEGs between the high- and low-GASC groups from all glioma population. (A, B) GO analysis (A) and KEGG pathway analysis (B) of up-regulated DEGs from TCGA data. (C) GSEA analysis of genes from TCGA data. (D, E) GO analysis (D) and KEGG pathway analysis (E) of up-regulated DEGs from CGGA data. (F) GSEA analysis of genes from CGGA data.

FIGURE 4

Figure 4 Functional annotation of upregulated DEGs between the high- and low-GASC groups from the high-grade glioma population. (A, B) GO analysis (A) and KEGG pathway analysis (B) of up-regulated DEGs from TCGA data. (C) GSEA analysis of genes from TCGA data. (D, E) GO analysis (D) and KEGG pathway analysis (E) of up-regulated DEGs from CGGA data. (F) GSEA analysis of genes from CGGA data.

Correlation of Stemness, Mesenchymal-EMT, Tumorigenic Cytokine, and Angiogenic Activity Scores With GASCs

To explore the potential mechanisms of GASCs in glioma, we also calculated stemness, mesenchymal-EMT, tumorigenic cytokine, and angiogenic activity scores for each glioma sample using the ssGSEA algorithm (Table S1). The results showed that mesenchymal-EMT, tumorigenic cytokine, and angiogenic activity scores were significantly higher in the high-GASC group (Figure 5A), and were positively correlated with the GASC score in all glioma population (Figures 5B, C) and the HGG population (Figures 5D, E).

FIGURE 5

Figure 5 Correlations between GASCs and Stemness, Mesenchymal-EMT, Tumorigenic cytokine, and Angiogenic activity scores. (A) Correlations between GASC groups and stemness, mesenchymal-EMT, tumorigenic cytokine and angiogenic activity scores. (B–E) Scatterplot of GASC score and 4 scores in all glioma population (B for TCGA and C for CGGA) and high-grade glioma population (D for TCGA and E for CGGA). *Statistical significance.

Associations Between GASCs and Stromal Cells

To discover the relationship between GASCs and other stromal cells, we computed the levels of 14 stromal cells using the xCELL algorithm. Bar charts showed that endothelial cells, lymphatic endothelial cells, and microvascular endothelial cells were higher in the high-GASC group in all glioma population from the TCGA (Figure 6A) and CGGA (Figure 6D) databases. Univariate Cox regression revealed that the level of mesenchymal stem cells is a protective factor for glioma (Figures 6B, E). A coefficient matrix showed that the GASC score was positively correlated with the levels of endothelial cells, lymphatic endothelial cells, and microvascular endothelial cells (Figures 6C, F). Similar results were found in the HGG population (Figure S3).

FIGURE 6

Figure 6 Associations between GASCs and stromal cells in all glioma population. (A, D) Bar charts illustrating the differences of xCELL scores between high- and low-GASC groups (A for TCGA and D for CGGA). ns: p > 0.05, *p ≤ 0.05, **p ≤ 0.01, ***p ≤ 0.001, ****p ≤ 0.0001. (B, E) Forest plots of univariate Cox regression analysis of stromal cells (B for TCGA and E for CGGA). (C, F) Correlograms of GASC score and stromal cells intercorrelation (C for TCGA and F for CGGA).

Immune Landscape of the High- and Low-GASC Groups

Because some immune-related terms were enriched in the functional annotation analysis, we explored the relationship between GASCs and the immune microenvironment. The CIBERSORT algorithm computed the relative abundance of 22 types of immune cells, which are shown in Figure 7. Overall, the adaptive immunity was at a relatively lower level in the high-GASC group compared with that in the low-GASC group. Notably, the level of M2 macrophages was significantly higher in the high-GASC group (Figures 8A, C) and was positively correlated with GASC score (Figure 9A) in all glioma population. Univariate Cox regression also revealed that the level of M2 macrophages is a risk factor for glioma (Figures 8B, D). Similar results were found in the HGG population (Figures 8E–H, 9B).

FIGURE 7

Figure 7 The relative abundances of the 22 types of immune cells. (A) Results from TCAG data. (B) Results from CGGA data.

FIGURE 8

Figure 8 Associations between GASCs and immune cells. (A, C, E, G) Bar charts illustrating the differences in CIBERSORT scores between high- and low-GASC groups in all glioma population (A for TCGA and C for CGGA) and high-grade glioma population (E for TCGA and G for CGGA). (B, D, F, H) Forest plots of univariate Cox regression analysis of immune cells in all glioma population (B for TCGA and D for CGGA) and high-grade glioma population (F for TCGA and H for CGGA). ns: p ≥ 0.05, *p ≤ 0.05, **p ≤ 0.01, ***p ≤ 0.001, ****p ≤ 0.0001.

We also analyzed the correlation of GASCs and 14 important ICPs. As shown in Figure 10, the expression levels of most ICPs were statistically higher in the high-GASC group. The univariate Cox regression showed that the expression levels of PD-L2, TIM3, CD80, CD86, CD155, and CIITA were risk factors for glioma in all glioma population and the HGG population. Correlation analysis indicated strong positive correlations within ICPs. The GASC score was positively correlated with PD-1, PD-L1, PD-L2, TL1A, TIM3, Galactin-9, CTLA-4, CD80, CD86, CD155, LAG3, and CIITA in all glioma population and the HGG population (Figure 9).

FIGURE 9

Figure 9 Correlogram of GASC score, immune cells, and expression of ICP intercorrelation. (A) Correlogram of data in all glioma population. (B) Correlogram of data in high-grade glioma population.

FIGURE 10

Figure 10 Associations between GASCs and expression of ICPs. (A, C, E, G) Bar charts illustrating the differences of ICP’ expressions between high- and low-GASC groups in all glioma population (A for TCGA and C for CGGA) and high-grade glioma population (E for TCGA and G for CGGA). (B, D, F, H) Forest plots of univariate Cox regression analysis of ICP’ expressions in all glioma population (B for TCGA and D for CGGA) and high-grade glioma population (F for TCGA and H for CGGA). ns: p > 0.05, *p ≤ 0.05, ***p ≤ 0.001, ****p ≤ 0.0001.

Copy Number Variations (CNVs) in DR3 and CIITA Indicated Worse OS

We also downloaded the somatic mutation and CNV data for glioma to analyze the difference in genomic alterations between the high- and low-GASC groups in all glioma population. The 20 genes with the greatest amounts of somatic mutations and CNVs are shown in Figure S4. We also compared the somatic mutations and CNVs of GASC markers between the high- and low-GASC groups (Figure S5), but found no significant difference in genomic alterations. However, in the comparison of somatic mutations and CNVs of ICPs between the high- and low-GASC groups (Figure S6), the results showed that the CNVs of DR3 and CIITA were significantly higher in the high-GASC group (Figure S6E). Survival analysis indicated that the CNVs of DR3 and CIITA significantly decreased the OS of glioma patients (Figures S6F, G).

Higher THY1 and CD80 Methylation Indicated Better OS

In the search for a possible treatment target for glioma, we conducted methylation analysis of GASC markers and ICPs. Because there was not a satisfactory match between samples with methylation data in the CGGA database and samples in the CGGA_693 mRNA dataset, we also downloaded the CGGA_325 mRNA dataset and separated these samples into high- and low-GASC groups with the previously mentioned method. Overall, 26 glioma samples (6 in the high-GASC group and 20 in the low-GASC group) and 8 normal samples with methylation data were included.

The GASC markers indicated that the methylation levels of THY1, CD9, CD14, CD44, ITGAM, and ACTA1 were significantly different among the high-GASC, low-GASC, and normal groups (Figure 11A and Figure S7). Then, we divided the glioma samples into high- and low-methylation groups according to the median of the gene methylation level. Survival analysis indicated that statistical difference was only observed between high- and low-THY1 methylation groups (p = 0.018; Figure 11B). High THY1 methylation suggested greater patient OS. In the methylation analysis of ICPs, significant differences were detected in Galactin-9, CD80, CD155, and LAG3 (Figure 11C and Figure S8). However, only high- and low-CD80 methylation groups showed statistical difference in the survival analysis (p = 0.031; Figure 11D), and high CD80 methylation indicated better OS.

FIGURE 11

Figure 11 Methylation analysis of GASC markers and ICPs. (A) Box plots illustrating the differences in THY1, CD9, CD14, CD44, ITGAM and ACTA1 methylation levels across high-GASC, low-GASC and normal groups. (B) Kaplan-Meier overall survival (OS) curves for samples in high- and low-methylation groups of THY1, CD9, CD14, CD44, ITGAM and ACTA1. (C) Box plots illustrating the differences in Galectin-9, CD80, CD155 and LAG3 methylation levels across high-GASC, low-GASC and normal groups. (D) Kaplan-Meier OS curves for samples in high- and low-methylation groups of Galectin-9, CD80, CD155 and LAG3.

Predicted Potential Immunotherapy Responses Between the High- and Low-GASC Groups

The TIDE webtool was applied to predict the likelihood of immune response for each sample. The results showed that in all glioma population, the low-GASC group (56%, 197/351 in TCGA; 40%, 138/347 in CGGA) was more likely to respond to immunotherapy than the high-GASC group (40%, 141/351 in TCGA; 29%, 102/346 in CGGA). However, in the HGG population, difference was found only in the TCGA dataset (high-GASC vs. low-GASC: 31% vs. 44%; Figures S9A–E).

In order to further analyze the immune infiltration between the responder and no responder groups, we compared the levels of immune cells between these two groups. The results showed that in all glioma population, the responder group had lower “T cells CD8” and “Macrophages M0” and higher “Mast cells activated” (Figures S9F, G). For the HGG population, “T cells CD8” was lower in the responder group, and “Mast cells activated” was higher in the responder group (Figures S9H, I).

Construction of a Risk Score System and Establishment and Validation of a Nomogram Survival Model

The mRNA sequencing data from TCGA (702 samples) was used as training dataset, and the data from CGGA (693 samples) was set as an independent validation dataset. For the training dataset, GASC score, GASC markers, immune checkpoints, stemness score, mesenchymal-EMT score, tumorigenic cytokine score, angiogenic activity score, stromal cell scores, and immune cell scores were filtered using LASSO regression with the “glmnet” R package. The change in trajectory of each variable was plotted in Figure 12A. We utilized 10-fold cross-validation to construct the model, and Figure 12B shows the confidence interval under each lambda. When lambda equaled 0.03431609, the model reached the optimal value, and 19 variables were selected for the next analysis. In the multivariable Cox regression analysis, the number of variables was reduced to 9, and the final 9-variable signature formula was: Risk score = 0.20 × CSPG4 – 0.30 × ALCAM + 35.16 × Adipocytes – 11.11 × Osteoblast – 7.05 × Pericytes – 2.50 × Plasma cells + 0.22 × CD274 + 0.16 × CD80 + 13.45 × angiogenesis. The risk score was calculated for each sample in the training and validation datasets. Thus, we divided samples into high- and low-risk groups according to the median risk score. Survival analysis revealed that in the training dataset, glioma patients in high-risk group have worse OS (p < 0.0001; Figure 12C), which was also confirmed in the validation dataset (p < 0.0001; Figure 12D).

FIGURE 12

Figure 12 Construction of the risk score system and establishment and validation of the nomogram survival model. (A, B) LASSO Cox regression analysis of training dataset. (C, D) Kaplan-Meier overall survival (OS) curves for samples in high- and low-risk score groups (C for training dataset and D for validation dataset). (E) Uni- and multi-variate Cox regression analysis for prognostic model. *Statistical significance. (F) The nomogram for predicting 1-, 3‐, or 5‐year OS. (G, H) Time-dependent ROC curves of the nomogram prediction on the 1-, 3-, and 5-year survival rates (G for training dataset and H for validation dataset). (I, J) Time-dependent ROC curves of the IDH mutation status prediction on the 1-, 3-, and 5-year survival rates (I for training dataset and J for validation dataset).

Then, we constructed a nomogram model to predict the prognosis of glioma, which included the risk score and clinicopathologic features. The uni- and multivariate Cox regression analysis (Figure 12E) indicated that risk score was an independent predictor for glioma prognostics. We finally included four features (age, WHO grade, isocitrate dehydrogenase (IDH) mutation status, and risk score) in the nomogram model (Figure 12F). Time-dependent ROC analysis further indicated that the area under the curve (AUC) for 1-, 3-, and 5-year OS were 0.902, 0.948, and 0.911, respectively, in the training dataset (Figure 12G). These AUCs were better compared with IDH mutation status, which is a traditional indicator, and were 0.842, 0.862, and 0.813 at the 1-, 3-, and 5-year marks (Figure 12I). Similar results were obtained in the validation datasets (Figures 12H, J).

Discussion

The GASC is a recently identified particular type of cell in the glioma microenvironment, with various names, e.g., glioma-associated human MSCs (GA-hMSCs) (7). The phenotypic and functional properties of GASCs are similar to those of cancer-associated fibroblasts and mesenchymal stem cells. The mechanism of GASCs in the glioma microenvironment is still largely unknown. In this work, we explored the potential mechanisms of GASCs in the glioma microenvironment, and discovered that GASCs may upregulate the level of M2 macrophages and ICPs. We also found that the CNVs of DR3 and CIITA were higher in the high-GASC group, and the methylation level of THY1 was lower in the high-GASC group, which could be a potential treatment target for glioma, particularly in HGG.

The tumor microenvironment determines the invasiveness of glioma. The EMT regulates this invasive state of glioma, particularly in HGG (9). Studies reported that GASCs drive cell invasion through HA synthase-2 (HAS2) induction (15), the UCA1/miR-182/PFKFB2 axis (16), the C5a/p38/ZEB1 axis (17) and CCL2/JAK1/MLC2 signaling (18). In the current work, we also discovered that the GASC score was positively correlated with the mesenchymal-EMT score in all glioma population and the HGG population (Figure 5). Terms from functional annotation include adhesion-related terms (Tables S2, S6), e.g., “Cell adhesion molecules”.

In the current study, we found a strong correlation between the GASC score and tumorigenic cytokine score, indicating the tumor-supporting function of GASCs. Studies reported that GASCs have tumor-promoting effects in vitro and in vivo (19–21). Additionally, Figueroa et al. also suggested that the tumor-supporting role of GASCs is mediated by the exosomal delivery of specific oncogenic miRNAs (21).

Although GASCs infiltrate into the glioma stroma, they are predominantly located around blood vessels (22), particularly abnormal vessels (23). Previous studies indicated that GASCs increase the angiogenesis of glioma (24, 25). Zhang et al. suggested that CD90^low (THY1) GASCs stimulate angiogenesis via vascular endothelial cells (25). In the current work, we detected a high correlation between the GASC score and the angiogenesis score (Figure 5) in all glioma population and the HGG population. We also found that in addition to endothelial cells, the GASC score also was positively correlated with microvascular endothelial cells (Figures 6 and Figure S3). The levels of endothelial cells and microvascular endothelial cells were higher in the high-GASC group (Figures 6 and Figure S3). These results indicated that GASCs may promote angiogenesis of glioma by stimulating the growth of both blood vessels and microvessels, which requires further verification.

Tumor-associated macrophages (TAMs) play an emerging role in glioma progression and are found in high proportions in the immune landscape of malignant glioma (26, 27). There are continuous phenotypes in the activation state of TAMs, in which M1 and M2 represent two extreme phenotypes (28). M2 has an anti-inflammatory phenotype, which leads to downregulation of immune responses, and thus prevents tissue damage and supports healing processes (27). In this work, our results suggest for the first time that GASCs are highly correlated with M2 macrophages in the glioma microenvironment. Based on our results, the level of M2 macrophages in the high-GASC group is statistically higher than that in the low-GASC group in all glioma population and the HGG population (p ≤ 0.0001; Figure 8). We also found high correlation coefficients between GASC scores and M2 macrophages in all glioma population (R = 0.46 (TCGA); R = 0.30 (CGGA); Figure 9A) and the HGG population (R = 0.44 (TCGA); R = 0.26 (CGGA); Figure 9B). These results indicated that TAMs may be phenotypically polarized to M2 macrophages by GASCs, which may further depress the immunity of the microenvironment and stimulate malignant progression of glioma. Conversely, the M2 macrophages may also upregulate the level of GASCs and further increase the malignant properties of glioma, e.g., invasion and angiogenesis.

Immune checkpoint blockade is the most developed immunotherapy in clinical use (4), but its efficiency still remains doubtful. We analyzed the expression levels of 14 important ICPs and found that the expression levels of most ICPs were higher in the high-GASC group. Although this result suggested that the high-GASC group may have a more optimal immunotherapy response, the results from the TIDE prediction were puzzling because they showed a contrary tendency (Figure S9). These contradictory results reflect the inner complexity of glioma, and in response to this, further high-quality studies of immunotherapy in glioma are required.

McDonald et al. reported that deletion of DR3 (Tumor Necrosis Factor Receptor Superfamily Members 25, TNFRSF25) was found in oligodendroglioma (29). The results from Qian et al.’s work suggested that suppression of CIITA (class II transactivator) downregulates the expression of MHC class II molecules in glioma (30). In the current study, we discovered that the CNVs of DR3 (P < 0.001) and CIITA (p = 0.015) were significantly higher in the high-GASC group (Figure S6E). The glioma patients with amplified/deleted DR3 or amplified CIITA had worse OS compared with wild-type glioma patients (Figures S6F, G). These results indicate that the CNVs of DR3 and CIITA may be potential prognostic indicators for glioma, and further studies are expected to verify their efficiency.

THY1 (CD90) is a surrogate marker for a variety of stem cells, including glioblastoma stem cells (GSC) (31) and GASC (7). Svensson et al. detected CD90^- and CD90⁺ GASC subpopulations by cell sorting and discovered that the CD90^- subpopulation exhibited greater tumor vascularization and immunosuppression activity than the CD90⁺ subpopulation (32). Zhang et al. further investigated these two subpopulations. They found that CD90^high GASCs drove glioma progression via increasing proliferation, migration, and adhesion. However, CD90^low GASCs contributed to glioma progression through the stimulation of vascular formation via vascular endothelial cells (25). In the current work, we discovered that the methylation levels were different among high-GASC, low-GASC, and normal groups (Figure 11A), and the high THY1 methylation group had better OS compared with the low THY1 methylation group (p = 0.018; Figure 11B). These results suggested for the first time that the methylation of THY1 could be a potential prognostic indicator of glioma as well as a treatment target.

To create a comprehensive risk score, we included the following features, produced from mRNA sequencing data, in the filtering process: GASC score, stemness score, mesenchymal-EMT score, tumorigenic cytokine score, angiogenic activity score, stromal cell scores, and immune cell scores. The results showed that the risk score could be used to differentiate patients with high or low risk (Figures 12C, D), and the risk score was an independent prognostic indicator for glioma (Figure 12E). The validation results from the CGGA validation dataset verified the robustness of our nomogram model (Figure 12H). We also compared the efficiency of the nomogram model with a traditional prognostic indicator, IDH mutation status. The AUCs for the nomogram model were better than those for IDH mutation status in the training dataset (Figures 12G, I) and the validation dataset (Figures 12H, J). To verify the credibility of this nomogram model, further high-quality clinical studies are required.

There are limitations in our study. First, because of a lack of public mRNA resources with proportions of GASCs, stromal cells, and immune cells, we selected the ssGEAS algorithm to compute these data, because it has been widely used with proven reliability. Second, because this is a retrospective study, the efficiency of our risk score and nomogram model needs to be verified in further high-quality prospective cohorts. In addition, our predicted results for immunotherapy response were contradictory with the expression level of ICPs. Glioma, particularly HGG, is characterized by remarkably high tumor heterogeneity and an “immune-cold” phenotype, denoting an immunosuppressive microenvironment. The GASC is an important cell type in the microenvironment and might influence immunotherapy responses. Based on previous research, the overexpression of some ICPs on malignant cells may increase the anti-tumor immune responses (33). Nevertheless, although our data suggested that higher ICP expression occurs in the high-GASC group, TIME prediction revealed worse immunotherapy responses in the high-GASC group. This contradiction may result from the inner complexity of the glioma microenvironment. Further prospective clinical trials to test immunotherapy for glioma with the GASC proportion data are required to produce a reliable conclusion regarding the relationship between GASCs and immunotherapy responses.

Conclusion

We found potential mechanisms of GASCs in the glioma microenvironment, particularly HGG, and we also discovered that GASCs are positively correlated with the level of M2 macrophages and ICPs. The methylation of THY1 decreased in the high-GASC group, which could be a prognostic indicator and treatment target for glioma. We also developed a prognostic nomogram for glioma. Further studies to verify these findings and the performance of our model with large prospective cohorts are warranted.

Data Availability Statement

Publicly available datasets were analyzed in this study. This data can be found here: TCGA database (https://portal.gdc.cancer.gov/) and CGGA database (www.cgga.org.cn/).

Author Contributions

CM conceived and designed the investigation. XC and FY analyzed the data and drafted the manuscript. JZ, JY, CT, and ZC conducted statistical analyses. All authors contributed to the article and approved the submitted version.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

Special thanks to Sheng Zhao for his support.

Supplementary Material

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

Supplementary Figure 1 | Differential analysis of expressed genes. (A, B) Heatmap of DEGs in all glioma population (A for TCGA and B for CGGA). (C, D) Heatmap of DEGs in high-grade glioma population (C for TCGA and D for CGGA).

Supplementary Figure 2 | Functional annotation of downregulated DEGs between the high- and low-GASC groups. (A, B) GO analysis (A) and KEGG pathway analysis (B) of down-regulated DEGs in all glioma population from TCGA data. (C, D) GO analysis (C) and KEGG pathway analysis (D) of down-regulated DEGs in all glioma population from CGGA data. (E, F) GO analysis (E) and KEGG pathway analysis (F) of down-regulated DEGs in high-grade glioma population from TCGA data. (G, H) GO analysis (G) and KEGG pathway analysis (H) of down-regulated DEGs in high-grade glioma population from CGGA data.

Supplementary Figure 3 | Associations between GASCs and stromal cells in the high-grade glioma population. (A, D) Bar chart illustrating the differences in xCELL scores between high- and low-GASC groups (A for TCGA and D for CGGA). ns: p > 0.05, *p ≤ 0.05, **p ≤ 0.01, ***p ≤ 0.001, ****p ≤ 0.0001. (B,E) Forest plot of univariate Cox regression analysis of stromal cells (B for TCGA and E for CGGA). (C,F) Correlogram of GASC score and stromal cells intercorrelation (C for TCGA and F for CGGA).

Supplementary Figure 4 | Comparison of genomic alterations between the high- and low-GASC groups in the TCGA dataset. (A, B) Differential copy number variation analysis between high- (A) and low- (B) GASC groups. (C, D) Differential somatic mutation analysis between high- (C) and low- (D) GASC groups.

Supplementary Figure 5 | Comparison of genomic alterations of GASC markers between the high- and low-GASC groups in the TCGA dataset. (A, B) Differential copy number variation analysis between high- (A) and low- (B) GASC groups. (C, D) Differential somatic mutation analysis between high- (C) and low- (D) GASC groups. (E) Correlations between GASC groups and somatic mutation of GASC markers.

Supplementary Figure 6 | Comparison of genomic alterations of ICPs between the high- and low-GASC groups in the TCGA dataset. (A, B) Differential copy number variation analysis between high- (A) and low- (B) GASC groups. (C, D) Differential somatic mutation analysis between high- (C) and low- (D) GASC groups. (E) Correlations between GASC groups and somatic mutation of ICPs. (F,G) Kaplan-Meier overall survival (OS) curves for samples of amplified, deleted and wildtype DR3 groups (F) and CIITA groups (G).

Supplementary Figure 7 | Methylation analysis of GASC markers. (A) Box plots illustrating the differences in CD34, ALCAM, CSPG4, ENG, GFAP, S100A4, NT5E, PDGFRB, PECAM1 and PTPRC methylation levels across high-GASC, low-GASC and normal groups.

Supplementary Figure 8 | Methylation analysis of immune checkpoints. (A) Box plots illustrating the differences in PD-1, PD-L1, PD-L2, DR3, TL1A, CTLA-4, CD86, TIM3 and CIITA methylation levels across high-GASC, low-GASC and normal groups.

Supplementary Figure 9 | Predicted potential immunotherapy responses between the high- and low-GASC groups. (A–D) Predicted potential immunotherapy responses of samples from all glioma population (A for TCGA and B for CGGA) and high-grade glioma population (C for TCGA and D for CGGA). (E) Correlation of GASC and Predicted immunotherapy responses. (F–I) Bar chart illustrating the differences in immune cell scores between responder and no responder groups in all glioma population (F for TCGA and G for CGGA) and high-grade glioma population (H for TCGA and I for CGGA).

References

1. Louis DN, Ohgaki H, Wiestler OD, Cavenee WK, Burger PC, Jouvet A, et al. The 2007 WHO Classification of Tumours of the Central Nervous System. Acta Neuropathol (2007) 114:97–109. doi: 10.1007/s00401-007-0243-4

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Liu S, Zhao Q, Shi W, Zheng Z, Liu Z, Meng L, et al. Advances in Radiotherapy and Comprehensive Treatment of High-Grade Glioma: Immunotherapy and Tumor-Treating Fields. J Cancer (2021) 12:1094–104. doi: 10.7150/jca.51107

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Hegde PS, Chen DS. Top 10 Challenges in Cancer Immunotherapy. Immunity (2020) 52:17–35. doi: 10.1016/j.immuni.2019.12.011

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Qi Y, Liu B, Sun Q, Xiong X, Chen Q. Immune Checkpoint Targeted Therapy in Glioma: Status and Hopes. Front Immunol (2020) 11:578877. doi: 10.3389/fimmu.2020.578877

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Luo H, Tao C, Wang P, Li J, Huang K, Zhu X. Development of a Prognostic Index Based on Immunogenomic Landscape Analysis in Glioma. Immun Inflammation Dis (2021) 1–13. doi: 10.1002/iid3.407

CrossRef Full Text | Google Scholar

6. Sonabend AM, Stupp R, Lee-Chang C, Okada H. Glioma Immunoediting, a Driver of Tumor Evolution, and the Next Battle for Immunotherapy. Oncotarget (2021) 12:8–9. doi: 10.18632/oncotarget.27865

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Clavreul A, Menei P. Mesenchymal Stromal-Like Cells in the Glioma Microenvironment: What are These Cells? Cancers (Basel) (2020) 12(9):2628. doi: 10.3390/cancers12092628

CrossRef Full Text | Google Scholar

8. 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

9. Iser IC, Pereira MB, Lenz G, Wink MR. The Epithelial-to-Mesenchymal Transition-Like Process in Glioblastoma: An Updated Systematic Review and In Silico Investigation. Med Res Rev (2017) 37:271–313. doi: 10.1002/med.21408

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Sheu BC, Chang WC, Cheng CY, Lin HH, Chang DY, Huang SC. Cytokine Regulation Networks in the Cancer Microenvironment. Front Biosci (2008) 13:6255–68. doi: 10.2741/3152

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Chang HY, Sneddon JB, Alizadeh AA, Sood R, West RB, Montgomery K, et al. Gene Expression Signature of Fibroblast Serum Response Predicts Human Cancer Progression: Similarities Between Tumors and Wounds. PloS Biol (2004) 2:E7. doi: 10.1371/journal.pbio.0020007

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Yu G, Wang LG, Han Y, He QY. clusterProfiler: An R Package for Comparing Biological Themes Among Gene Clusters. OMICS (2012) 16:284–7. doi: 10.1089/omi.2011.0118

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Mayakonda A, Lin DC, 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

14. Fu J, Li K, Zhang W, Wan C, Zhang J, Jiang P, et al. Large-Scale Public Data Reuse to Model Immunotherapy Response and Resistance. Genome Med (2020) 12:21. doi: 10.1186/s13073-020-0721-z

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Lim EJ, Suh Y, Yoo KC, Lee JH, Kim IG, Kim MJ, et al. Tumor-Associated Mesenchymal Stem-Like Cells Provide Extracellular Signaling Cue for Invasiveness of Glioblastoma Cells. Oncotarget (2017) 8:1438–48. doi: 10.18632/oncotarget.13638

PubMed Abstract | CrossRef Full Text | Google Scholar

16. He Z, You C, Zhao D. Long non-Coding RNA Uca1/miR-182/PFKFB2 Axis Modulates Glioblastoma-Associated Stromal Cells-Mediated Glycolysis and Invasion of Glioma Cells. Biochem Biophys Res Commun (2018) 500:569–76. doi: 10.1016/j.bbrc.2018.04.091

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Lim E-J, Kim S, Oh Y, Suh Y, Kaushik N, Lee J-H, et al. Corrigendum to: Crosstalk Between GBM Cells and Mesenchymal Stemlike Cells Promotes the Invasiveness of GBM Through the C5a/p38/ZEB1 Axis. Neuro-Oncology (2020) 22:1229. doi: 10.1093/neuonc/noaa125

CrossRef Full Text | Google Scholar

18. Lim EJ, Suh Y, Kim S, Kang SG, Lee SJ. Force-Mediated Proinvasive Matrix Remodeling Driven by Tumor-Associated Mesenchymal Stem-Like Cells in Glioblastoma. BMB Rep (2018) 51:182–7. doi: 10.5483/bmbrep.2018.51.4.185

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Behnan J, Isakson P, Joel M, Cilio C, Langmoen IA, Vik-Mo EO, et al. Recruited Brain Tumor-Derived Mesenchymal Stem Cells Contribute to Brain Tumor Progression. Stem Cells (2014) 32:1110–23. doi: 10.1002/stem.1614

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Hossain A, Gumin J, Gao F, Figueroa J, Shinojima N, Takezaki T, et al. Mesenchymal Stem Cells Isolated From Human Gliomas Increase Proliferation and Maintain Stemness of Glioma Stem Cells Through the IL-6/gp130/STAT3 Pathway. Stem Cells (2015) 33:2400–15. doi: 10.1002/stem.2053

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Figueroa J, Phillips LM, Shahar T, Hossain A, Gumin J, Kim H, et al. Exosomes From Glioma-Associated Mesenchymal Stem Cells Increase the Tumorigenicity of Glioma Stem-Like Cells Via Transfer of Mir-1587. Cancer Res (2017) 77:5808–19. doi: 10.1158/0008-5472.CAN-16-2524

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Kim SM, Kang SG, Park NR, Mok HS, Huh YM, Lee SJ, et al. Presence of Glioma Stroma Mesenchymal Stem Cells in a Murine Orthotopic Glioma Model. Childs Nerv Syst (2011) 27:911–22. doi: 10.1007/s00381-011-1396-y

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Trylcova J, Busek P, Smetana K Jr., Balaziova E, Dvorankova B, Mifkova A, et al. Effect of Cancer-Associated Fibroblasts on the Migration of Glioma Cells In Vitro. Tumour Biol (2015) 36:5873–9. doi: 10.1007/s13277-015-3259-8

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Clavreul A, Guette C, Faguer R, Tétaud C, Boissard A, Lemaire L, et al. Glioblastoma-Associated Stromal Cells (Gascs) From Histologically Normal Surgical Margins Have a Myofibroblast Phenotype and Angiogenic Properties. J Pathol (2014) 233:74–88. doi: 10.1002/path.4332

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Zhang Q, Yi D-Y, Xue B-Z, Wen W-W, Lu Y-P, Abdelmaksou A, et al. CD90 Determined Two Subpopulations of Glioma-Associated Mesenchymal Stem Cells With Different Roles in Tumour Progression. Cell Death Disease (2018) 9:1101. doi: 10.1038/s41419-018-1140-6

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Akins EA, Aghi MK, Kumar S. Incorporating Tumor-Associated Macrophages Into Engineered Models of Glioma. iScience (2020) 23:101770–. doi: 10.1016/j.isci.2020.101770

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Roesch S, Rapp C, Dettling S, Herold-Mende C. When Immune Cells Turn Bad-Tumor-Associated Microglia/Macrophages in Glioma. Int J Mol Sci (2018) 19:436. doi: 10.3390/ijms19020436

CrossRef Full Text | Google Scholar

28. Murray PJ. Macrophage Polarization. Annu Rev Physiol (2017) 79:541–66. doi: 10.1146/annurev-physiol-022516-034339

PubMed Abstract | CrossRef Full Text | Google Scholar

29. McDonald JM, Dunmire V, Taylor E, Sawaya R, Bruner J, Fuller GN, et al. Attenuated Expression of DFFB is a Hallmark of Oligodendrogliomas With 1p-Allelic Loss. Mol Cancer (2005) 4:35. doi: 10.1186/1476-4598-4-35

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Qian J, Luo F, Yang J, Liu J, Liu R, Wang L, et al. Tlr2 Promotes Glioma Immune Evasion by Downregulating Mhc Class Ii Molecules in Microglia. Cancer Immunol Res (2018) 6:1220–33. doi: 10.1158/2326-6066.Cir-18-0020

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Avril T, Etcheverry A, Pineau R, Obacz J, Jegou G, Jouan F, et al. Cd90 Expression Controls Migration and Predicts Dasatinib Response in Glioblastoma. Clin Cancer Res (2017) 23:7360–74. doi: 10.1158/1078-0432.Ccr-17-1549

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Svensson A, Ramos-Moreno T, Eberstål S, Scheding S, Bengzon J. Identification of Two Distinct Mesenchymal Stromal Cell Populations in Human Malignant Glioma. J Neurooncol (2017) 131:245–54. doi: 10.1007/s11060-016-2302-y

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Yearley JH, Gibson C, Yu N, Moon C, Murphy E, Juco J, et al. Pd-L2 Expression in Human Tumors: Relevance to Anti-PD-1 Therapy in Cancer. Clin Cancer Res (2017) 23:3158–67. doi: 10.1158/1078-0432.Ccr-16-1761

PubMed Abstract | CrossRef Full Text | Google Scholar