High expression of stromal signatures correlated with macrophage infiltration, angiogenesis and poor prognosis in glioma microenvironment

Glioma is one of the most fatal tumors in central nervous system. Previous studies gradually revealed the association between tumor microenvironment and the prognosis of gliomas patients. However, the correlation between tumor-infiltrating immune cell and stromal signatures are unknown. In our study, we obtained gliomas samples from the Chinese Glioma Genome Atlas (CGGA) and The Cancer Genome Atlas (TCGA). The landscape of tumor infiltrating immune cell subtypes in gliomas was calculated by CIBERSORT. As a result, we found high infiltration of macrophages was correlated with poor outcome (P < 0.05). Then functional enrichment analysis of high/low macrophage-infiltrating groups was performed by GSEA. The results showed three gene sets includes 102 core genes about angiogenesis were detected in high macrophage-infiltrating group. Next, we constructed PPI network and analyzed prognostic value of 102 core genes. We found that five stromal signatures indicated poor prognosis which including HSPG2, FOXF1, KDR, COL3A1, SRPX2 (P < 0.05). Five stromal signatures were adopted to construct a classifier. The classifier showed powerful predictive ability (AUC = 0.748). Patients with a high risk score showed poor survival. Finally, we validated this classifier in TCGA and the result was consistent with CGGA. Our investigation of tumor microenvironment in gliomas may stimulate the new strategy in immunotherapy. Five stromal signature correlated with poor prognosis also provide a strong predator of gliomas patient outcome.


INTRODUCTION
Glioma is the most deadly central nervous system (CNS) tumor. Patients with gliomas always show poor outcome, and their median overall survival remains from 14.6 to 16.8 months. The prognosis of gliomas patients is correlated with tumor subtype, age and sex (Mortazavi, Mortazavi & Paknahad, 2018). Based on tumor morphology and METHODS Collecting RNA-seq and matched clinical data from CGGA and TCGA We built two independent cohorts for training set and testing set from CGGA and TCGA. The gene expression profiles were downloaded from The Chinese Glioma Genome Atlas (CGGA, http://www.cgga.org.cn/). A total of 695 RNA-seq data of Chinese gliomas patients were obtained from the CGGA database. A series of measures were taken: (1) Genes with a variance of 0 will be filtered out.
(2) Complete follow-up information and samples with a follow-up time >30 days will be served. Finally, 693 samples meeting the inclusion criteria were included. And each patient matched specific clinical data, which includes: histology, WHO grades, age, gender, chemotherapy, radiotherapy, survival status, survival duration in days. To build the testing cohort, The RNA-seq FPKM of gliomas including corresponding outcome data were downloaded from The Cancer Genome Atlas (TCGA, https://cancergenome.nih.gov/). The following measures were taken: (1) The IDs were annotated on the basis of hg38 reference genome. (2) Genes with a variance of 0 will be filtered out. (3) For the same gene corresponding to multiple IDs or a patient with multiple tumor samples, the average will be taken. (4) Samples with a follow-up time >30 days were remained. Eventually, a total of 668 TCGA samples fulfilled our criteria. The details of patients information was showed in Table 1.

Analysis of tumor-infiltrating immune cells
CIBERSORT is a computational method for quantifying immune cell subtypes fractions from normalized gene expression profiles based on deconvolution algorithm. The leukocyte gene signature matrix (LM22) is generated for the calculation of 22 human hematopoietic subsets. LM22 has been verified on Affymetrix HGU133 and Illumina Beadchip platforms. In our study, we used LM22 gene file to calculate 22 immune cells subtypes of 693 gliomas patients. These immune cell subtypes included CD8+ T cells, naïve CD4+ T cells, resting CD4+ memory T cells, activated CD4+ memory T cells, follicular helper T cells (Tfh), Tregs, γδ T cells, naïve B cells, memory B cells, plasma cells, resting NK cells, activated NK cells, resting NK cells, activated NK cells, resting

Statistical analysis
The differences of immune cell subtypes between normal brain tissues and gliomas tissues were assessed using the unpaired t-test, and the differences of tumor-infiltrating immune cells among different phenotypes were performed by t-test and one One-way analysis of variance (ANOVA). Kaplan-Meier curves were used to perform the correlations of 22 immune cell subtypes or gene signature and corresponding clinical follow-ups (log-rank test), and high vs low groups all based on median cut off. Statistical analyses were calculated by SPSS statistical package software or GraphPad Prism. P values < 0.05 were considered significant.

Construction of classifier
To generate and optimize the prognostic classifier, the multivariate Cox regression analysis was performed in CGGA and TCGA cohort. We calculated the RS of each sample based on the multivariate COX coefficient, and the low/high risk groups were defined according to the median cutoff RS. The receiver operator characteristics (ROC) curve analysis was applied to assess the classifier's ability to distinguish samples with a high or low RS, and it was draw by survivalROC package. The area under the curve (AUC) of the ROC curve was calculated and compared to examine the performance of the classifier in both training and testing cohorts. The median RS was determined to separate the genes into the high-risk or low-risk groups. KM curves and independent testify were performed by survival package to assess the effective of classifier in gliomas patients (log-rank test). All analysis were carried out by R version 3.6.1 and corresponding packages.

Proteomics and histology
To verify the infiltration of M2 macrophage in glioma patients, we downloaded proteomics data of 110 glioma patients from The National Cancer Institute's Clinical Proteomic Tumor Analysis Consortium (CPTAC; https://cptac-data-portal.georgetown.edu/). And the histological level research of glioma patients was performed in the human protein atlas (https://www.proteinatlas.org/).

RESULT The landscape of tumor infiltrating immune cell subtypes in gliomas
Based on the 693 RNA-seq from CGGA database, the different infiltration of 22 immune cell subtypes between normal brain tissue and gliomas were analyzed by t-test. The fraction of M0 macrophages, M2 macrophages, resting NK cells, plasma cells, CD8+ T cells was significantly higher in gliomas tissue than normal brain tissue (Fig. 1A). On the contrary, the fraction of naïve B cells, resting dendritic cells, Eosinophils, Neutrophils, M1 macrophages, activated mast cells, follicular helper T cells (Tfh), γδ T cells and resting CD4+ memory T cells was significantly lower in gliomas tissue than normal brain tissue (Fig. 1B). A total of nine clinical parameters were analyzed, which includes: age, gender, histology, WHO grades, primary or recurrence, chemotherapy, radiotherapy, mutations in IDH and 1p/19q co-deletion. As a result, gender, primary or recurrence, chemotherapy showed no significant difference of immune cells. Monocytes were decreased in elderly patients ( Fig. 2A). M0 macrophages, Tregs and activated dendritic cells were increased in high grade gliomas, whereas monocytes were decreased (Fig. 2B). The fraction of M0 macrophages, Tregs and γδ T cells was higher in glioblastoma (GBM) than astrocytoma (AOA), whereas monocytes and activated mast cells was lower (Fig. 2C). The fraction of activated mast cells, monocytes and resting CD4+ memory T cells was higher in IDH mutant than wildtype, while M0 macrophages, Tregs, γδ T cells and Tfh was lower

Prognostic value of tumor-infiltrating immune cells in gliomas
In this study, we divided all 22 immune cell subtypes into high/low infiltration group based on median cutoff. And we analyzed the correlation between 22 immune cell subtypes and patients prognosis by generating Kaplan-Meier survival curve. Result showed that high infiltration of M0 macrophages (P = 0.0104) and resting CD4+ memory T cells (P = 0.027) had an unfavorable prognosis. (Figs. 3A and 3B) Based on the previous studies, M2 macrophages showed pro-tumor function as well as M1 macrophages showed anti-tumor function in gliomas. We constructed a macrophage score: (M0 + M2)/M1. And results showed gliomas patients with high macrophage score correlated with poor outcome. (Fig. 3C).

Functional enrichment analysis of gene signature of macrophages
As the infiltration of macrophages play a significant role in glioma patients' outcomes, we further investigated the potential mechanism and gene signature. First, we divided gliomas patients into high/low macrophage-infiltrating groups. Functional enrichment analysis was performed by Gene Set Enrichment Analysis Gene Set Enrichment Analysis (GSEA). GO categories include BP, MF, CC and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway were analyzed. As a result, a total of 5 GO terms of biological process, 2 GO terms of cellular component, 8 GO terms of molecular function and 4 KEGG pathway terms were identified to be significant (FDR < 0.1 and P < 0.001) (Fig. 4). Interestingly, the 3 enrichment GO terms with high macrophage infiltration were correlated with angiogenesis, so we obtained the 102 core gene signatures for next investigation.

Investigation of stromal signatures correlated with macrophage infiltration and angiogenesis
We built PPI networks to further investigate the 102 core gene signatures which correlated with high macrophage infiltration and angiogenesis by STRING tool. Cytoscape and Cytohubba were used to analyze hub genes (Fig. S3). And the gene set which contain

Construction and verification of prognostic classifier based on stromal signatures
According to univariate COX regressions, all of five stromal signatures were significantly correlated with poor outcome in gliomas patients (HR > 1) (Fig. 6A). We calculated the RS Next, we generated univariate and multivariate analyses which enrolled clinical features and stromal signatures in the overall set. The result showed that stromal signature classifier was an independent factor for gliomas patients (Figs. 6C and 6D). Furthermore, we found all of five stromal signatures were correlated with WHO grades (P < 0.05) (Fig. 6F). As for histology, the expression of stromal signatures in GBM was significantly higher than AOA and oilgodendroglioma. The expression of COL3A1, HSPG2 and SRPX2 in patients with chemotherapy was higher than patients without chemotherapy. And the expression of COL3A1, KDR and SRPX2 was higher than patients with radiotherapy than patient without radiotherapy (Fig. S1). Principal component analysis showed significant different between groups of low and high RS (Fig. 6G). A receiver operating characteristic (ROC) curve was plotted to assess the effect of the classifier. The AUC of the classifier was 0.748. Although stromal classifier was less efficient than the WHO grade, it was much superior to histological type (Fig. 6H). To verify the classifier, we generated an independent cohort based on TCGA dataset. And the classifier also showed strong predictive ability in TCGA cohort (AUC = 0.726) (Fig. 6I).

Verifying the infiltration of M2 macrophage in glioma patients
CD163 is the marker protein of M2 macrophage. By measuring the expression of CD163, the infiltration degree of M2 macrophage in glioma can be determined. We demonstrated high expression of CD163 in glioma patients through the CPTAC proteomics database. Meanwhile, we compared the immunohistochemistry of CD163 in normal brain tissue and glioma through the human protein atlas, and found that the expression of CD163 in glioma was significantly higher than that in normal brain tissue (Fig. S2).

DISCUSSION
First of all, based on the RNA-seq data from CGGA and the CIBERSORT, we calculated 22 subtypes of tumor-infiltrating immune cells in gliomas. We comprehensively analyzed the tumor-infiltrating immune cells present in gliomas and revealed the prognosis value.
The landscape of tumor-infiltrating immune cells in gliomas as follows: 1. B cell lines: Naïve B cells can differentiate into antibody-secreting plasma cells after being stimulated by antigen (Liao et al., 2016). Our study found that naïve B cell was a significant decrease in gliomas, meanwhile, the plasma cells were significant increased. It was illustrated that naïve B cell may be differentiated into plasma cells by gliomas antigens.
2. T cell lines: Tregs are correlated with unfavorable prognosis in the several kinds of tumor microenvironment (e.g., ovarian cancer, breast cancer, kidney cancer and pancreatic cancer) (Mirzaei, Sarkar & Yong, 2017). Our results showed that Tregs were significantly increased in GBM and high WHO grade, which demonstrated Tregs may promote the development of gliomas (Banissi et al., 2009;Vandenberk & Van Gool, 2012). Previous study showed that tumor-infiltrating Tfh cells have protective roles by suppressing lymphoid tumor-promoting effects in breast cancer and colorectal cancer (Woroniecka et al., 2018). Our results also shown that Tfh cells was significantly decrease in gliomas. Previous study showed neutrophils and eosinophils are implicated in almost every stage of oncogenesis and display both anti-and pro-tumor properties (Curran, Evans & Bertics, 2011;Schneider, Kwan & Boockvar, 2018).

Myeloid cell lines:
Our results showed that neutrophils and eosinophils were significantly decreased in gliomas, which indicated neutrophils and eosinophils may have negative correlation with gliomas. Mast cells play an important role in the growth of tumors (Põlajeva et al., 2011). However, the contribution of mast cells in the microenvironment of solid malignancies remains controversial. Previous studies have illustrated that gliomas are associated with a profound accumulation of mast cells, and STAT5 play an important role on the recruitment of mast cells to gliomas (Põlajeva et al., 2014). Our results showed that resting mast cells were increased in gliomas, whereas activated mast cells were significantly decreased in gliomas. In addition, the infiltration of activated mast cells of GBM was prominently lower than AOA. Consistent with previous study, gliomas promoted mast cells recruitment, but gliomas may also block resting mast cell activation. And glioblstoma may show stronger ability of preventing mast cell activation than AOA. Dendritic cell vaccinations have emerged as newly strategies in the treatment of gliomas (Garg et al., 2016;Wen et al., 2019). Adjuvant dendritic cell-based immunotherapy for gliomas patients can induce long-term survival (Li et al., 2018). We found resting dendritic cells were significantly decreased in glioma, which was consistent with previous studies that dendritic cells recruitment correlated with favorable outcome. Previous study showed monocytes played an important role in cancer development and progression (Põlajeva et al., 2011). And different monocytes displayed both anti-and pro-tumor function (Prosniak et al., 2013). Our results showed monocytes were significant decreased in high grade glioma and GBM, which indicated stimulating monocytes recruitment may prevent progression of high grade gliomas. What's more, our results also showed monocytes were significantly lower in patients with radiotherapy than patients without radiotherapy, which may contribute to the resistance of gliomas radiotherapy.
Next, we further investigated tumor-infiltrating macrophages in gliomas. Based on their function, M0 macrophages can be differentiated into two categories: classically activated macrophages (M1 macrophages) and alternatively activated macrophages (M2 macrophages) by stimulation of various cytokines Jiang et al., 2019;Lailler et al., 2019). Previous studies show gliomas contain a plenty of macrophages, and tumor-associated macrophages correlate with progression and angiogenesis of gliomas by several ways, including: microRNA-macrophage feedback loop (Bao & Li, 2019;Liu et al., 2019), extracellular lipid loading (Offer et al., 2019) and cytokines like OPN, IRGM, IL-6 (Hernandez-SanMiguel et al., 2019;Wang et al., 2018;Xu et al., 2019). In addition, previous study successfully reeducated the pro-tumor M2 macrophages toward anti-tumor M1 macrophages by a dual-targeting biomimetic treatment strategy which provides a good method for the pharmacotherapy of gliomas . Our results showed that M0 and M2 were significantly increased in gliomas, while M1 were significant decreased. The results confirmed pro-tumor function of M0 and M2 as well as anti-tumor function of M1 in gliomas. And the infiltration of M0 macrophages showed correlated with decreased survival (P < 0.05). In addition, we found M0 macrophages were significantly lower in IDH mutant than IDH wildtype. This result may explain glioma patients with IDH-mutant always have favorable survival outcomes. After functional enrichment analysis, 3 gene sets includes 102 genes about angiogenesis were detected in high macrophage-infiltrating group, which was consistent with previous study.
Among 102 core genes, we found five stromal signatures correlated with high macrophages infiltration and angiogenesis indicated poor prognosis (P < 0.05). The five stromal signatures including HSPG2, FOXF1, KDR, COL3A1 and SRPX2. Sun et al. (2019) constructed the multicellular gene network between gliomas and macrophages, they found macrophage-related gene signature had good prognostic value for predicting resistance to targeted therapeutics and survival of glioma patients. We found the relationship between stromal signatures and macorphages signatures in gliomas. Based on precious studies, high expression of perlecan/HSPG2 in gliomas lead to tumor promotion through the transformation of brain extracellular matrix into tumor microenvironment (Hu et al., 2016). In addition, high expression of COL3A1 contributes to glioma cells proliferation and migration (Shin et al., 2017). And up-regulation of KDR promotes proangiogenic myeloid cells, result in low-grade to high-grade transition (Huang et al., 2017). Moreover, glioma with high tumor associated macrophages has increased expression of SRPX2 (Hung et al., 2016). We also found these five stromal signatures were correlated with clinical features like grades, histology, chemotherapy and radiotherapy. The expression of stromal signatures was significantly higher in patients with higher grade and worse pathological features like GBM. Then, we constructed a prognosis classifier based on HSPG2, FOXF1, KDR, COL3A1 and SRPX2. And the classifier showed strong predictive ability in both CGGA cohort and TCGA cohort. We confirmed the high infiltration of M2 macrophages in protein and histology level via CPTAC and the human protein atlas. Eventually, based on our investigation, up-regulation of five stromal signatures may lead to macrophage infiltration and angiogenesis in gliomas patients. And our classifier also provides strong predictive ability in prognosis. Our study may provide a novel immunotherapeutic strategy and a more clearly overview of tumor microenvironment in gliomas.
However, there still exited several limitations in our research. First, the calculation results from public database may show bias. Although we have verified the results in two independent cohorts, we should continue to deeper research. Second, the macrophage infiltration and angiogenesis function of five stromal signatures in gliomas need to further confirmation in vitro and in vivo.

CONCLUSIONS
Based on CGGA RNA-seq datasets, the landscape of 22 tumor-infiltrating immune cells in gliomas was analyzed by CIBERSORT. High infiltration of macrophages was correlated with poor outcome. After functional enrichment analysis, three gene sets includes 102 core genes about angiogenesis were detected in high macrophage-infiltrating group. Next, we found five stromal signatures indicated poor prognosis which including HSPG2, FOXF1, KDR, COL3A1, SRPX2. Five stromal signatures were adopted to construct a classifier. The classifier showed powerful predictive ability. Finally, we validated this classifier in TCGA and the result was consistent with CGGA. Our investigation of tumor microenvironment in gliomas may stimulate the new strategy in immunotherapy. The five stromal signatures correlated with poor prognosis also provide a strong predator of glioma patient outcomes.