Prognosis and Immunotherapy Significances of a Cancer-Associated Fibroblasts-Related Gene Signature in Gliomas

As a cold tumor, malignant glioma has strong immunosuppression and immune escape characteristics. The tumor microenvironment (TME) provides the “soil” for the survival of malignant tumors, and cancer-associated fibroblasts (CAFs) are the architects of matrix remodeling in TME. Therefore, CAFs have potent regulatory effects on the recruitment and functional differentiation of immune cells, whereby they synthesize and secrete numerous collagens, cytokines, chemokines, and other soluble factors whose interaction with tumor cells creates an immunosuppressive TME. This consequently facilitates the immune escape of tumor cells. Targeting CAFs would improve the TME and enhance the efficacy of immunotherapy. Thus, regulation of CAFs and CAFs-related genes holds promise as effective immunotherapies for gliomas. Here, by analyzing the Chinese Glioma Genome Atlas and the Cancer Genome Atlas database, the proportion of CAFs in the tumor was revealed to be associated with clinical and immune characteristics of gliomas. Moreover, a risk model based on the expression of CAFs-related six-gene for the assessment of glioma patients was constructed using the least absolute shrinkage and selection operator and the results showed that a high-risk group had a higher expression of the CAFs-related six-genes and lower overall survival rates compared with those in the low-risk group. Additionally, patients in the high-risk group exhibited older age, high tumor grade, isocitrate dehydrogenase wildtype, 1p/19q non-codeletion, O-6-methylguanine-DNA methyltransferase promoter unmethylation and poor prognosis. The high-risk subtype had a high proportion CAFs in the TME of glioma, and a high expression of immune checkpoint genes. Analysis of the Submap algorithm indicated that the high-risk patients could show potent response to anti-PD-1 therapy. The established risk prediction model based on the expression of six CAFs-related genes has application prospects as an independent prognostic indicator and a predictor of the response of patients to immunotherapy.

Introduction

Glioma is the most common form of a primary malignant brain tumor in adults. In particular, glioblastoma (GBM) is one of the most lethal and highly aggressive cancer, associated with low survival, less than 1 year (Lapointe et al., 2018). Although current evidence indicates surgical resection as the main treatment approach for gliomas, postoperative radiotherapy, and chemotherapy can also be administered according to the specific condition of patients. However, most gliomas are difficult to completely be removed by surgery without affecting normal brain functions due to the growth properties and its special local anatomy characteristics. Some patients also do not respond well to radiotherapy and chemotherapy, and consequently, develop rapid recurrence after standard therapy (Tan et al., 2020). Despite recent progress in the development of new drugs, researchers are facing challenges in developing therapeutics for gliomas due to biological properties, including the blood-brain barrier, tumor specificity of gliomas, and immune environments (Aldape et al., 2019).

The tumor microenvironment (TME), as the “soil” for tumor survival, is crucial for tumor survival and is closely associated with the malignant behavior of tumor cells equivalent to “seed.” Compelling evidence shows that tumor cells interact with the extracellular matrix (ECM), immune cells, chemokines, and cytokines to create a favorable microenvironment for the proliferation and metastasis of tumors (Maman and Witz, 2018). The microenvironment of different tumors is diverse. For instance, in glioma, tumors can cooperate with peritumoral cells through chemokines and cytokines, direct contact, extracellular vesicles, nanotubes and microtubules, to promote tumor proliferation, brain invasion, angiogenesis and immunosuppression. This consequently creates a microenvironment conducive to the growth of malignant tumors (Broekman et al., 2018). Therefore, an in-depth understanding of the interactions between tumor cells and peritumoral cells may provide a new perspective to managing gliomas.

Interestingly, cancer-associated fibroblasts (CAFs) are an important component of stromal cells in TME and an “architect” of matrix remodeling, which are closely related to the prognosis of solid tumors (Cox, 2021). In TME, CAFs play a key role in the induction of Epithelial-mesenchymal Transition (EMT), and the maintenance of the pool of cancer stem cells and drug resistance by interacting with tumor cells and immune cells and releasing a variety of soluble factors (Su et al., 2018; Erin et al., 2020). Moreover, CAFs regulate tumor immunity and promote immune escape and resistance to cancer immunotherapy (Liu et al., 2019a). Therefore, targeting CAFs would not only inhibit the “seed” of cancer but also transform the “soil” into a microenvironment that inhibits tumor growth, and consequently transform the “enemy” that promotes tumor progression into a “friend” that inhibits tumor growth or metastasis (Chen and Song, 2019). Also, CAFs are activated to different degrees at different stages of tumor development. Studies have shown that different cytokines secreted by CAFs can play a pro-cancer or anti-cancer role (Liu et al., 2019b; Wang et al., 2021).

It is also noteworthy that because glioma is a “cold” tumor, immunotherapy has poor efficacy against the malignancy (Jackson et al., 2019). As such, to improve the efficacy of immunotherapy, transforming the “cold” environment into a “hot” one without causing neurotoxicity is imperative. Additionally, novel immunotherapeutic approaches, including oncolytic virus and adoptive T cell therapy, may exploit the T cell response to overcome the “cold” state of glioma, such as GBM (Buerki et al., 2018). Intriguingly, CAFs are co-expressed with the Fibroblast Activation Protein α (FAPα) and the Platelet-derived Growth Factor Receptor (PDGFR), which are the main components of stromal cells in GBM. Oncolytic adenoviruses have also been shown to specifically target GBM cells and CAFs (Li et al., 2020).

Cancer-associated fibroblasts are highly heterogeneous in terms of tissue origin, phenotype and function (Louault et al., 2020), but their function in the glioma microenvironment is yet to be fully elucidated. In the present study, we analyzed the sequencing data of glioma cohorts from the Chinese Glioma Genome Atlas (CGGA) database by the Estimate Proportion of Immune and Cancer cells (EPIC) algorithm (Racle et al., 2017; Liu et al., 2020), which is a common method to analyze the cell types based on the gene expression profile of the different type of cells, to quantify the expression of CAFs, then further analyze the relationship among CAFs and clinical features, tumor purity, immune score, stromal score, ESTIMATE score and stemness score of gliomas. Differentially expressed genes related to prognosis in the high and low expression subtypes of CAFs were also screened. Moreover, the risk scores of gliomas were assessed according to the expression of six genes identified as independent prognostic factors via the Least Absolute Shrinkage and Selection Operator (LASSO) regression analysis and validated in the Cancer Genome Atlas (TCGA) database. The established risk signature containing six CAFs-related-genes strongly correlates with the clinical and immune characteristics of glioma, including immune cells, immune checkpoints, and immunotherapy.

Materials and Methods

Data Collection

The RNA sequencing datasets [mRNAseq_693, mRNAseq_325 and mRNA sequencing (non-glioma as a control)] and the corresponding clinical and molecular information, including sex, age, grade, IDH status, 1p/19q status, MGMT promoter, and survival status information, were retrieved from the CGGA database¹ (Zhao et al., 2021). The data were categorized into the training cohort [including 625 cases of low-grade gliomas (LGG) and 388 cases of GBM] and the testing cohort (comprising transcriptome data from 698 cases of gliomas in TCGA database from the Genomic Data Commons Data Portal (GDC).² Additionally, the FPKM data was converted into TPM for subsequent analysis. The mRNAseq_693 and mRNAseq_325 data were merged into one metadata set, and batch effects were removed using the combat function in the SVA R package.

Quantification of Different Cell Types in the Tumor Microenvironment

Estimate Proportion of Immune and Cancer cells (Racle et al., 2017) is an effective algorithm used to simultaneously estimate the proportion of cancer and immune cell types according to the gene expression in tumors based on a unique set of RNA-seq reference gene expression profiles described previously (Tirosh et al., 2016). It allows for accurate prediction of the proportion of cancer and non-malignant cell types even in the absence of a priori information about cancer cells. In this study, the reference profiles from tumor-infiltrating cells were used as the parameter, whereas the TPM data of CGGA glioma were used as input data. The proportion of CAF was estimated using the R package “EPIC” to explore the changes of matrix components in gliomas. CIBERSORT is a tool that deconvolutes the expression matrix of human immune cell subtypes based on the principle of linear support vector regression. Researchers also apply this tool to determine immune cell composition, which comprises 22 immune cell subtypes, based on the specific gene expression data of the cells (Newman et al., 2015). We analyzed gene expression data with standard annotation using CIBERSORT source code in relative mode. The algorithm was run using the LM22 signature and 100 permutations. For each sample, the final CIBERSORT output estimates were normalized to sum to one such that it could be interpreted directly as cell fractions for comparison across different immune cell types and datasets. Immune infiltration analysis based on a single-sample gene set enrichment analysis (ssGSEA) score can also be employed to explore the degree of immune infiltration of gliomas. It defines an enrichment score representing the absolute degree of enrichment of the gene set in each sample in each dataset. ssGSEA-based evaluation of the level of immune infiltration in a sample according to the expression levels of immune cell-specific marker genes demonstrated the immune infiltration landscape of gliomas (Bindea et al., 2013). The ESTIMATE package was used to estimate tumor immune score and tumor purity (Yoshihara et al., 2013). The stemness index of glioma was calculated according to the expression of tumor stem cell genes (Miranda et al., 2019). Submap³ algorithms (Hoshida et al., 2007) were further utilized to predict the clinical response to immune checkpoint blocking therapy for PD-1 and CTLA4 in the low-and high-risk score groups. P-value < 0.05 denoted statistical significance.

Identification of Risk Genes and Calculation of Risk Score

Data were grouped into high subtype and low subtype according to the median value of CAFs. This was followed by the analysis of the differences between the mRNAseq_325 and mRNAseq_693 CGGA data. Differentially expressed genes (DEGs) were obtained according to a fold change (logFC) > 1.5 and P < 0.05. Genes with P < 0.001 and genes associated with prognosis were screened through univariate Cox regression analysis. A total of 329 genes were identified from mRNAseq_325 and 116 from the mRNAseq_693 CGGA database. Crossing the two sets of genes yielded 104 genes. Furthermore, multivariate Cox regression analysis was performed to identify 11 genes related to glioma overall survival (OS) (P < 0.001). Subsequently, LASSO analysis was employed, to minimize the risk of overfitting a prognostic model and construct a risk model (Tibshirani, 1997). The LASSO was utilized for variable selection and shrinkage via the “glmnet” package in R software (Simon et al., 2011). The independent variable in the regression was the normalized expression matrix of candidates for prognostic factors. The response variable was the OS time and state of the patients in the CGGA cohort. The penalty parameter (λ) of the model was constructed through tenfold cross-validations, followed by the minimum criteria (i.e., the value of λ corresponding to the lowest partial likelihood deviance). The risk score of the patients was calculated according to the normalized expression level of the prognostic gene signature and their corresponding regression coefficients according to the formula: risk score = e^{sum (each gene}′^{s expression × corresponding coefficient)}. Finally, patients were divided into high and low risk groups according to the median risk score.

Establishment of the Nomogram

The nomogram incorporated age, grade, 1p/19q status based on the CGGA cohort. The prognostic risk score model was established via the “RMS” package in R. The consistency between predicted survival rate and actual survival rate using time-dependent calibration curves, and verified in the TCGA cohort. The concordance index (C index) was calculated to evaluate the effectiveness of the model in prognosis prediction. The C index ranged between 0.5 and 1.0; notably, a higher index denoted the better the performance of the model in predicting survival rate.

Functional Enrichment Analysis

Spearman’s correlation analysis was used to reveal the genes associated with the risk score; here, the correlation coefficient ≥ 0.4 was selected. Genes related to risk scores were sorted according to the calculated correlation coefficient. Next, Gene Set Enrichment Analysis (GSEA) was performed with the “clusterProfiler” R package (Yu et al., 2012), using “h.all.v7.0.entrez.gmt” as a reference gene set. P values were adjusted using the Benjamini and Hochberg methods. P-value < 0.05 implied statistically significant differences. Lastly, the results of the first five enrichment analysis were visualized using “enrichplot” and “ggplot2” (Wickham, 2016) R software package.

Statistical Analysis

Patients in the CGGA training and TCGA validation cohorts were categorized into the high- and low-risk subtypes according to the median risk score. Wilcoxon rank-sum test was applied to compare the high- and low-risk subtypes. Differences among three or more subtypes of patients were tested using the K-W test. Kaplan-Meier analysis and log-rank test were employed to analyze survival rates between low-risk and high-risk subtypes. Univariate and multivariate Cox regression analyses identified the independent factors associated with the OS of glioma. All statistical analyses were conducted in R software (version 4.0.3), and P < 0.05 denoted statistical significance (^∗P < 0.05, ^∗∗P < 0.01, ^∗∗∗P < 0.001, ^****P < 0.0001).

Results

Cancer-Associated Fibroblasts Are Closely Related to the Clinical Progression of Gliomas

A detailed flow chart of this analysis is shown in Figure 1. Using EPIC, the proportion of CAFs was calculated in different glioma samples in the CGGA cohort. The mRNAseq_693 and mRNAseq_325 datasets in the CGGA cohort were merged and batch discrepancies were eliminated. Subsequently, the proportion of CAFs was calculated in gliomas in the training cohort. Patients were then stratified into high and low CAFs subtypes according to the median of CAFs. The prognosis of the high-CAFs subtype was worse than the low-CAFs subtype (log-rank, P < 0.0001; Figure 2A). According to the WHO (2016) grading guidelines, patients with high CAFs showed poor prognosis in Oligodendroglioma with IDH mutation and 1p/19q co-deletion, Astrocytoma with IDH mutant, Astrocytoma with IDH wildtype, GBM with IDH mutant, and GBM with IDH wildtype (Supplementary Figure 1). In addition, both univariate and multivariate Cox regression analyses verified that the proportion of CAFs is an independent risk predictor for gliomas (Figures 2B,C). On the other hand, the high proportion of CAF is mainly enriched in glioma patients with higher age (P < 0.0001), high WHO grade (P < 0.001), IDH wildtype status (P < 0.001), 1p/19q non-codeletion status (P < 0.001), and MGMT promoter un-methylated status (P < 0.01). However, CAFs did not differ between genders (Figures 2D–I). To eliminate the effect of data consolidation in mRNAseq_693 and mRNAseq_325 datasets, the proportion of CAFs in both datasets was calculated, respectively, which yielded a similar result to the merged dataset (Supplementary Figure 2). Results suggest that the proportion of CAFs is of promise as an independent predictor for the prognosis and progression of gliomas.

FIGURE 1

Figure 1. Flowchart of the study.

FIGURE 2

Figure 2. Predictive values of cancer-associated fibroblasts (CAFs) in the survival of glioma patients and its relationship with pathological characteristics. (A) Comparison of overall survival between low and high CAFs subtypes in the Chinese Glioma Genome Atlas (CGGA) database. (B) Univariate Cox regression analysis. Forest plot of the association between CAFs and glioma survival. (C) Multivariate Cox regression analysis showed that the CAFs was an independent predictor of gliomas. (D–N) Distribution of the CAFs in patients stratified according to age (D), gender (E), WHO grade (F), isocitrate dehydrogenase (IDH) status (G), 1p/19q status (H), O-6-methylguanine-DNA methyltransferase (MGMT) promoter status (I), Tumor purity (J), Immune score (K), Stromal score (L), ESTIMATE score (M), and Stemness score (N). ****P < 0.0001; ***P < 0.001; **P < 0.01; ns, no significant.

A High Proportion of Cancer-Associated Fibroblasts Is Associated With Immune Landscape and Stemness of Gliomas

Cancer-associated fibroblasts as the architect of matrix remodeling in TME may affect cell components of TME and promote progression of tumor malignancy. In this study, the ESTIMATE and GSVA packages in R were used to estimate tumor purity, immune score, stromal score, ESTIMATE score, and stemness score of gliomas in the CGGA cohort. A high proportion of CAFs was revealed to be associated with the glioma patients with low tumor purity, high stromal score, high ESTIMATE score, and high stemness score (P < 0.0001; Figures 2J–N). At the same time, the prognosis of glioma patients characterized by low tumor purity, high immune score, high stromal score, high ESTIMATE score, and high stemness score was poor (P < 0.0001, Supplementary Figures 3A–E). The findings suggest the association of low tumor purity, high immune score, high stromal score, high ESTIMATE score, and high stemness score with a high proportion of CAFs and the OS of glioma patients.

Construction of Prognostic Gene Signatures of Cancer-Associated Fibroblasts Which Are Related to the Status Gliomas

After calculating the proportion of CAFs in the CGGA mRNAseq_693 and mRNAseq_325 cohorts, the training cohorts were stratified into high and low subtypes according to the median of CAFs. Subsequently, univariate Cox regression analysis was conducted on both datasets of DEGs to screen for genes related to prognosis (with P < 0.001). Through a crossover analysis of the two DEGs sets, 104 genes were obtained. Eleven genes related to the total survival of glioma patients were obtained via multivariate Cox regression analysis (P < 0.001). The risk of overfitting was minimized using the LASSO regression algorithm. The risk score was calculated according to the expression level and regression coefficient of six genes (ABCC3, CTHRC1, MSR1, PDLIM1, TNFRSF12A, and CHI3L2) (Figures 3A,B). The risk score = (ABCC3 × 0.00530703547540297) + (PDLIM1 × 0.0015 8225241997429) + (CHI3L2 × 0.0000527762285183538) + (MSR1 × 0.00118065457455651) + (CTHRC1 × 0.004799633 48819953) + (TNFRSF12A × 0.00212454022154963). The formula was also used to calculated risk scores for the glioma patients in the TCGA validation cohorts. Of note, the six genes associated with the risk of glioma were highly expressed in the high-risk group (Figure 3C). There were significant differences in risk scores among patients with different ages in terms of first diagnosis (P < 0.001), WHO grade (P < 0.0001), IDH (P < 0.0001), 1p/19q (P < 0.0001) and MGMT promoter (P < 0.001). These six genes were also mainly highly expressed in high-grade glioma and IDH wildtype groups (Supplementary Figures 4A–L). These results strongly demonstrated that the risk model, based on the six genes is closely associated with the clinical progression of glioma, therefore, could be employed as a risk prediction model for glioma.

FIGURE 3

Figure 3. Identification of a CAFs-related six-gene risk signature for overall survival risk characteristics by least absolute shrinkage and selection operator (LASSO) regression analysis in CGGA cohort. (A) Cross-validation for tuning parameter (lambda) screening in the LASSO regression model. (B) LASSO coefficient spectrum of six genes in gliomas. (C) The heatmap shows the association between risk and clinic pathological characteristics of the six-gene risk signature. LASSO, least absolute shrinkage and selection operator. ****P < 0.0001; ***P < 0.001; ns, no significant.

Survival Analysis and Pathological Features Between High- and Low-Risk Patients

Kaplan-Meier analysis showed that the OS of the high-risk subtype was worse compared to that of the low-risk subtype (P < 0.0001; Figure 4A). Additionally, the risk score was related to a prognostic value in gliomas, stratified according to the WHO guidelines (2016) for the grading of tumors. The high-risk subtype was related to the low OS of the patients with Oligodendroglioma with IDH mutant and 1p/19q co-deletion, Astrocytoma with IDH mutant, Astrocytoma with IDH wildtype, GBM with IDH mutant, and GBM with IDH wildtype (Supplementary Figure 5). Next, we performed a multivariate Cox regression analysis of the risk score and clinical-pathological features of glioma patients. Results demonstrated that risk score is an independent risk factor to predict the OS in patients with glioma (HR = 1.731, 95% CI = 1.341–2.237, P < 0.001; Figure 4B). The heatmap showed that high glioma mortality is related to an increased risk score (Figure 4C). Furthermore, through ROC curve analysis, the accuracy of risk score as a prognostic factor for glioma was validated. Results showed that the risk score could predict the OS of the CGGA cohorts (5-year, AUC = 0.789). Of note, higher CAFs (AUC = 0.7994), tumor histology (AUC = 0.7822), IDH status (AUC = 0.8066), 1p/19q status (AUC = 0.8163) and MGMT promotor status (AUC = 0.5641) was shown in the high-risk group than in the lower group (Figures 4D–I). Moreover, glioma data from TCGA was utilized to validate the risk score. The LASSO regression analysis was performed on the TCGA data to calculate the patients’ risk scores using similar regression coefficients. Subsequently, K-M survival analysis of the TCGA data was performed to assess the risk model. Results demonstrated lower OS in the high-risk subtype than that in the low-risk subtype (Supplementary Figure 6A). Univariate Cox regression analysis was conducted to explore the prognostic value of the risk score. Results showed a significant correlation of the risk score with OS in TCGA LGG-GBM (HR = 4.255, 95% CI = 3.538–5.117, P < 0.001, Supplementary Figure 6B). Moreover, multivariate Cox regression analysis demonstrated that the risk score was an independent prognostic indicator (HR = 1.888, 95% CI = 1.438–2.4, P < 0.001, Supplementary Figure 6C). ROC curve analysis revealed that the risk model had a strong predictive value for the OS of glioma patients (Supplementary Figure 6D). Furthermore, six CAFs-related genes were highly expressed in TCGA LGG-GBM cohorts, particularly in the high-risk subtype (Supplementary Figure 6E), consistent with the results of CGGA cohorts. Taken together, the results provided evidence that the high-risk subtype is associated with low OS rate and some clinical-pathological features, therefore, has application prospects as a risk model based on the specific expression of the six genes.

FIGURE 4

Figure 4. Prognostic significance of risk signature derived risk score in the CGGA cohort. (A) Kaplan-Meier analysis of CGGA glioma patients was stratified by median risk. (B) Multivariate Cox regression analysis showed that the risk signature was an independent predictor of gliomas. (C) Risk score distribution, patient survival time, and glioma status. The black dotted line is the optimal threshold for classifying patients into low-and high-risk subtypes. (D–I) A high-risk score is associated with a lower survival rate for gliomas. ROC curves showed the predictive efficiency of risk characteristics, overall survival (D), CAFs (E), Histology (F), IDH status (G), 1p/19q status (H), and MGMT promoter methylation status (I).

The Risk Model for Individual Prognostic Prediction

Nomogram is a powerful tool used to quantitatively determine individual risk in the clinical setting by integrating various risk factors. Herein, using a six CAFs-related genes signature, a nomogram was constructed based on age, grade, 1p/19q gene deletion status, and risk score to predict the probability of 3- and 5-year OS. Meanwhile, the calculated C index was 0.760. Each factor was graded according to its contribution to OS (Figure 5A). The calibration curve showed a consistent actual survival rate with the predicted survival rate (Figures 5B,C). The accuracy of the prognostic prediction model was verified in the TCGA cohort, and the calculated C index was 0.844. The correction chart demonstrated that the 3- and 5-year OS corroborated with the predicted values (Figures 5D,E). These data demonstrate that the risk model is relatively accurate and can improve the ability to estimate individual prognosis.

FIGURE 5

Figure 5. The nomogram can predict the prognosis of gliomas. (A) A nomogram of the gliomas cohort (training set) was used to predict overall survival. (B,C) A calibration plot was used to predict the 3-year (B) and 5-year survival (C) in the CGGA cohorts (training set). Calibration plots for 3-year (D) and 5-year survival (E) in the Cancer Genome Atlas (TCGA) cohort (testing set). The X- and Y-axis represent the Nomogram prediction and actual survival, respectively. The solid line represents the nomogram of the forecast and the vertical line represents the 95% confidence interval.

The Immune Landscape Between Different Risk Subtypes

To determine the potential biological processes of risk score, Spearman’s correlation analysis was performed to identify the genes related to the risk score. Genes with Spearman correlation coefficient ≥ 0.4 were used for GSEA analysis. Results revealed that genes associated with risk score were mainly associated with epithelial-mesenchymal transformation (P < 0.0001, q < 0.0001), hypoxia (P < 0.0001, q < 0.0001), inflammation (P < 0.0001, q < 0.0001), interferon-gamma response (P < 0.0001, q < 0.0001), and NFkB-mediated TNFα signal transduction (P < 0.0001, q < 0.0001; Figure 6A). Furthermore, the ESTIMATE and GSVA packages in R were employed to estimate tumor purity, immune score, stromal score, ESTIMATE score, stemness score, and the proportion of immune cells in the CGGA cohort. Results showed that the high-risk subtype was associated with a high proportion of CAFs, high stemness score, high stromal score, high immune score, the high ESTIMATE score, and low tumor purity (P < 0.0001; Figure 6B). In addition, aDC, pDC, iDC, T helper cells, macrophage, Th2 cells, Treg cells, and B cells were rich in high-risk subtype, whereas neutrophil, Th1 cells were rich in low-risk subtypes (Figure 6B). To further verify that immunosuppressive cells were mainly enriched in the high-risk group, a histogram was generated and the proportion of immune cells in the CGGA data was calculated via the CIBERSORT algorithm. Results revealed that immunosuppressive cells, including Tregs and M2 macrophages, were mainly concentrated in the high-risk subtype (Supplementary Figures 7A,B). Intriguingly, the proportion of immune cells calculated by ssGSEA and the CIBERSORT algorithm showed that CD4⁺, CD8⁺ T cells, and NK cells were also enriched in the high-risk subtype (Supplementary Figure 7C).

FIGURE 6

Figure 6. Gene set enrichment analysis (GSEA) and immune cell infiltration analysis. (A) The GSEA of genes positively correlated with risk score. (B) Heatmap showing immune cells infiltration by ssGSEA algorithm between low- and high-risk groups in the CGGA cohort. NES, normalized enrichment score. ****P < 0.0001.

The Immune Checkpoint Landscape Between Different Risk Subtypes

Previously, T cells and NK cells exhaustion have been demonstrated to potentially facilitate cancer cells to escape host immunity, which is linked to poor prognosis (Gonzalez-Gugel et al., 2016; Zarour, 2016). In the present work, most of these exhaustion markers [TIGIT (Manieri et al., 2017), CEACAM1 (Huang et al., 2015), CTLA4 (Krummel and Allison, 1996; Postow et al., 2015), LAG3 (Ruffo et al., 2019), PD-1 (Ishida et al., 1992; Postow et al., 2015), PD-L1 (Postow et al., 2015), and TIM3 (Huang et al., 2015)] were highly expressed in the high-risk subtype (Figure 7A), indicating an elevated level of immune exhaustion in the tumors of high-risk glioma patients. It is well known that CAFs can promote the transformation of M1 macrophages to M2 and induce the secretion of related cytokines to promote tumor invasion, angiogenesis, and change the immune landscape of tumors (Farhood et al., 2019; An et al., 2020). Therefore, we explored the expression of chemokines and cytokines secreted by CAFs in the high-risk subtype and the low-risk subtype. M2 macrophage chemokines (IL10, IL13, CSF1, TGFB1, TGFB2, and TGFB3) were found to be highly positively correlated, whereas M1 macrophage chemokines (HMGB1, TNF) were weakly positively correlated with risk score (Zhang et al., 2021; Figures 7B,C). Additionally, analysis of the expression of chemokines and cytokines secreted by CAFs in the high-risk subtype and the low-risk subtype demonstrated that most of the cytokines secreted by CAFs were highly expressed in the high-risk subtype (Figure 7D). Patients with a high-risk score exhibited a type I interferon response, type II interferon response, and the activation of pro-inflammatory functions (Supplementary Figure 7C). These observations provided evidence that patients with a high-risk score are eligible for and respond to immunotherapy. Furthermore, subclass mapping was applied to compare the expression profile of the high- and low-risk subtypes with another published dataset containing 47 patients with melanoma that responded to immunotherapies (Roh et al., 2017). The high-risk subtype showed potent responses to anti–PD-1 therapy (Bonferroni correction P = 0.01; Figure 7E).

FIGURE 7

Figure 7. The role of risk signature in immune microenvironment and immunotherapy. (A) Expression of immune cell exhaustion marker genes in the high- and low-risk subtypes. (B,C) Correlation between representative macrophage chemokines and risk score, CGGA dataset (B), TCGA dataset (C). M1 macrophage chemokines (TNF, HMGB1), M2 macrophage chemokines (IL10, IL13, CSF1, TGFB1, TGFB2, and TGFB3). (D) Cytokines secreted by CAFs high expression in the high-risk subtype. (E) Predicting response to immunotherapy (anti-PD-1 and anti-CTLA4) in the high-and low-risk subtype based on the Submap algorithms.

Discussion

The present study, through analysis of the CGGA and TCGA databases of glioma, revealed that CAFs are associated with clinical and immunological characteristics of gliomas. More evidence shows that CAFs are an independent prognostic factor and CAFs are enriched in patients with a poor survival rate, older age, high tumor grade, IDH wild-type, 1p/19q non-codeletion, and MGMT promoter unmethylation (Figures 2A-I). This work has also demonstrated the association of a high proportion of CAFs in glioma patients with low tumor purity and high immune score, high stromal score, ESTIMATE score, and high stemness score (Figures 2J-N). By intersecting the two DEGs sets and conducting a further screen through multivariate Cox regression analysis, the remaining prognostic genes were used for LASSO regression analysis. Assessment of the risk score was based on the expression of six CAFs-related genes (ABCC3, CTHRC1, MSR1, PDLIM1, TNFRSF12A, and CHI3L2) (Figure 3). Similar to the proportion of CAFs, the risk score was correlated with a survival rate, clinical characteristics, tumor purity, immune score, stromal score, ESTIMATE score, and stemness score. These data strongly demonstrate that the risk score holds great promise as an independent prognostic factor for glioma (Figure 4). Furthermore, independent prognostic factors were used to construct a nomogram to directly predict individual OS based on six CAFs-related-gene signature risk score (Figure 5). It was revealed that the main enriched pathways in the tumor, including epithelial-mesenchymal transformation (Shintani et al., 2016), hypoxia (Lappano et al., 2020), inflammatory Response (Ershaid et al., 2019), interferon-gamma response (Broad et al., 2021), and NFkB-mediated TNFα signal transduction (Katanov et al., 2015) were closely related to the function of CAFs. These observations indicate that risk scores can fully reflect the function of CAFs (Figure 6A). To explore the role of the risk score in the TME, ssGSEA, and CIBERSORT were employed to evaluate immune cell infiltration. Percentage analysis of immune cells demonstrated that gliomas were mainly enriched with macrophages and few T cells. Also, significant differences in the abundance of immune cells between the high-and low-risk subtypes were noted. Th1 cells were mainly enriched in the low-risk subtype while immunosuppressive cells such as Th2 cells, Tregs, and M2 macrophages were mainly enriched in the high-risk category. Furthermore, CD4⁺ T cells, CD8⁺ T cells, and NK cells were highly enriched in the high-risk subtype (Figure 6B and Supplementary Figure 7). In addition, high-risk patients were shown to express high levels of markers of immune cell exhaustion (Figure 7A). These data demonstrate that although immune-activated cells are highly expressed in the high-risk group, it is possible that they are in a state of functional inhibition.

High expression of immune-related cytokines and chemokines in high-risk subtypes provide evidence on the relationship between CAFs and tumor-associated macrophages (TAMs) (Figures 7B–D). Studies show that recruitment of TAMs into the glioma environment may induce immunosuppression and tumor promotion effects (Hambardzumyan et al., 2016). Moreover, CAFs may promote the polarization of macrophages to M2 macrophages and exert an immunosuppressive effect in TME, and more, CAFs can secret cytokines and chemokines to interact with various immune cells in the tumor environment, contributing to the malignant transformation of tumors and resistance to treatment (Linares et al., 2020). Therefore, secretion factors from CAFs are promising indicators for tumor diagnosis and prognosis, and as drug targets.

Previous evidence indicates that CCL2 in the glioma microenvironment promotes the recruitment of Tregs and myeloid-derived suppressor cells (MDSCs) (Chang et al., 2016). On the other hand, CD70⁺ CAFs are independent markers for poor prognosis in invasive colorectal cancer, they potentially increase the infiltration of Tregs and promote the immune escape of tumor cells (Jacobs et al., 2018). Immunosuppression mediated by high expression of PD-L1 in gliomas has been demonstrated to be potentially associated with the infiltration of TAMs and M2 polarization (Zhu et al., 2020). Studies have also shown that the differentiation and survival of macrophages are dependent on the colony-stimulating Factor-1 (CSF-1), and IL6 and the granulocyte-macrophage colony-stimulating factor (GM-CSF) secreted by CAFs, which promote the differentiation of monocytes into M2-like TAMs (Cho et al., 2018). Consistently, Pyonteck et al. (2013) found that targeting TAMs with the CSF-1R inhibitor significantly improves the survival rate of GBM mice. These pieces of evidence indicate that CAFs and their secreted factors orchestrate with the immune cells in the TME to promote glioma progression.

Increasing evidence has shown that T cell dysfunction contributes to tumor immune escape in patients with gliomas (Mirzaei et al., 2017; Woroniecka et al., 2018). In the presence of adenosine, activation of CD8⁺ T cells effectuates a decrease in the expression of IFN-γ and tumor necrosis factor-α (TNF-α), thereby inhibiting anti-tumor response (Takenaka et al., 2019). Moreover, CD8⁺ T cells during tumorigenesis mediate the immune editing of immunogenic tumor clones, contributing to immune escape in murine glioma (Kane et al., 2020). The present study, based on the six-gene signature of CAFs, reported a higher expression of immune checkpoint molecules in high-risk subtype tumors. It was notable that the functions of CAFs rely on immune checkpoint activation, which induces the loss of antigen specificity in CD8⁺ T cells and blocks the activity of T cells (Lakins et al., 2018). Furthermore, evidence indicates that to inhibit CD8⁺ T cells, CAFs can secrete high levels of IL-6, which remodel the immunosuppressive TME (Kato et al., 2018). Previous researchers also demonstrated that CAFs can recruit and balance CD4⁺ effector T cell subsets (Th1 and Th2), promoting the recruitment and differentiation of Tregs, and facilitate the transformation of macrophages to the M2 phenotype through frequent interactions with TAMs. Additionally, CAFs can inhibit the ability of CD8⁺ T cells to kill tumor cells. To achieve this, CAFs reduce T cell infiltration, blocking both the cytotoxic activity of T cells and T-cell communication in the TME (An et al., 2020; Baker et al., 2021). Taken together, these findings support the view that CAFs potentially inhibit CD8⁺ T function, promote the formation of MDSC, and establish an immunosuppressive TME to facilitate the immune escape of tumor cells.

Additionally, the present work revealed the activation of the inflammation-promoting function, type I IFN response, and type II IFN response in the high-risk subtype (Figure 6B and Supplementary Figure 7C). The expression levels of immune checkpoint-related genes were also higher in the high-risk subtype (Figure 7A), demonstrating that immunosuppressed patients in the high-risk category may respond to immune checkpoint blockers. By predicting the response and resistance of different risk subtypes to immune checkpoint blockers therapy, it was intriguing that patients with high-risk scores showed higher responses to anti-PD-1 therapy (Figure 7E). As such, patients in the high-risk subtype are likely to be more responsive to immunotherapy. These results demonstrate that the risk score based on the expression of the six CAFs-related genes is promising as a novel and reliable method for evaluating the prognosis and clinical response to immunotherapy of glioma patients.

More evidence shows that CAFs-induced inhibition of immunosuppression may further enhance the response of tumors to immunotherapy. Current evidence shows several therapeutic strategies that target CAFs, including (i) targeting cytokines and chemokines (such as TGF-β or IL6) through immunotherapy, directly consuming CAFs (such as FAP-DNA vaccine) via cell surface labeling; (ii) targeting CAFs through the elimination of endothelial progenitors using bevacizumab, normalizing activated CAFs (e.g., using VDR ligand calcipotriol); (iii) targeting CAFs-derived extracellular matrix proteins or their related signal transduction in animal models (Chen and Song, 2019; Liu et al., 2019a). Moreover, selectively targeting CAFs using nanomedicine has been revealed to enhance the infiltration of cytotoxic T cells and inhibits tumor growth (Zhen et al., 2017). More importantly, the synergistic effect of blocking immune checkpoint molecules and targeting CAFs can be achieved by remodeling the immunosuppressive microenvironment and achieving an immunotherapeutic response (Feig et al., 2013). Normalization of CAFs may eliminate the tumor-promoting effect and increase the sensitivity of treatment (Öhlund et al., 2014; Vennin et al., 2019). Collectively, CAFs secrete cytokines in the TME during glioma progression, interact with immune cells, mediate the formation of the immunosuppressive microenvironment and induce the transformation to the malignant phenotype. Therefore, an improved understanding of the interaction of CAFs with anti-tumor immunity is crucial in establishing effective immunotherapy. In this regard, approaches such as spatial transcriptomics, single-cell RNA sequencing (Moncada et al., 2020), and organoids (Xu et al., 2021) can be employed to comprehensively understand the Spatio-temporal dynamics of CAFs as they interact with tumor and immune cells and their role in TME of gliomas.

To the best of our knowledge, this is the first study to explore the clinical features of CAFs in gliomas and establish a prognostic signature, based on CAFs for predicting the survival outcome of glioma patients and immunotherapy efficacy. The results are based on public data sets analysis, therefore, more exploration is warranted on heterogeneity in different patient groups, including intertumoral or intratumoral heterogeneity. Also, the biological functions of six CAFs-related genes as a prognostic signature should be explored deeply. Notably, the roles of CAFs vary in different patients given the complex and diverse immune microenvironment of gliomas. Thus, the precise effects of CAFs and their interaction with tumor cells and immune cells deserve further clarification.

Conclusion

Increased CAFs infiltration in gliomas is significantly correlated with older age, high tumor grade, IDH status, 1p/19q status and MGMT promoter status, tumor purity, immune score, stromal score, ESTIMATE score, stemness score, and patient prognosis. The established risk prediction model, based on the expression of six CAFs-related genes, has application prospects as an independent prognostic indicator. The risk model holds great promise in predict prognosis and immunotherapy response in glioma patients.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

Author Contributions

HX and KY conceived and designed the study and drafted the manuscript. ZC, SZ, and GH performed data analysis and wrote the manuscript. JT, WH, and W-QG revised the manuscript. All authors reviewed the manuscript.

Funding

This study was supported by the National Natural Science Foundation of China (Grant Nos. 31571399 and 82060456) and Finance Science and Technology Project of Hainan Province (Grant No. ZDYF2019129).

Conflict of Interest

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

Publisher’s Note

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

Acknowledgments

The authors gratefully acknowledge contributions from the CGGA and TCGA network.

Supplementary Material

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

Footnotes

1. ^ http://www.cgga.org.cn/
2. ^ https://portal.gdc.cancer.gov/
3. ^ https://cloud.genepattern.org/gp

References

Aldape, K., Brindle, K. M., Chesler, L., Chopra, R., Gajjar, A., Gilbert, M. R., et al. (2019). Challenges to curing primary brain tumours. Nat. Rev. Clin. Oncol. 16, 509–520. doi: 10.1038/s41571-019-0177-5

PubMed Abstract | CrossRef Full Text | Google Scholar

An, Y., Liu, F., Chen, Y., and Yang, Q. (2020). Crosstalk between cancer-associated fibroblasts and immune cells in cancer. J. Cell Mol. Med. 24, 13–24. doi: 10.1111/jcmm.14745

PubMed Abstract | CrossRef Full Text | Google Scholar

Baker, A. T., Abuwarwar, M. H., Poly, L., Wilkins, S., and Fletcher, A. L. (2021). Cancer-associated fibroblasts and T cells: from mechanisms to outcomes. J. Immunol. 206, 310–320. doi: 10.4049/jimmunol.2001203

PubMed Abstract | CrossRef Full Text | Google Scholar

Bindea, G., Mlecnik, B., Tosolini, M., Kirilovsky, A., Waldner, M., Obenauf, A. C., et al. (2013). Spatiotemporal dynamics of intratumoral immune cells reveal the immune landscape in human cancer. Immunity 39, 782–795. doi: 10.1016/j.immuni.2013.10.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Broad, R. V., Jones, S. J., Teske, M. C., Wastall, L. M., Hanby, A. M., Thorne, J. L., et al. (2021). Inhibition of interferon-signalling halts cancer-associated fibroblast-dependent protection of breast cancer cells from chemotherapy. Br. J. Cancer 124, 1110–1120. doi: 10.1038/s41416-020-01226-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Broekman, M. L., Maas, S. L. N., Abels, E. R., Mempel, T. R., Krichevsky, A. M., and Breakefield, X. O. (2018). Multidimensional communication in the microenvirons of glioblastoma. Nat. Rev. Neurol. 14, 482–495. doi: 10.1038/s41582-018-0025-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Buerki, R. A., Chheda, Z. S., and Okada, H. (2018). Immunotherapy of primary brain tumors: facts and hopes. Clin. Cancer Res. 24, 5198–5205. doi: 10.1158/1078-0432.CCR-17-2769

PubMed Abstract | CrossRef Full Text | Google Scholar

Chang, A. L., Miska, J., Wainwright, D. A., Dey, M., Rivetta, C. V., Yu, D., et al. (2016). CCL2 Produced by the Glioma Microenvironment is essential for the recruitment of regulatory T Cells and myeloid-derived suppressor cells. Cancer Res. 76, 5671–5682. doi: 10.1158/0008-5472.CAN-16-0144

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, X., and Song, E. (2019). Turning foes to friends: targeting cancer-associated fibroblasts. Nat. Rev. Drug Discov. 18, 99–115. doi: 10.1038/s41573-018-0004-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Cho, H., Seo, Y., Loke, K. M., Kim, S. W., Oh, S. M., Kim, J. H., et al. (2018). Cancer-stimulated CAFs enhance monocyte differentiation and protumoral TAM activation via IL6 and GM-CSF secretion. Clin. Cancer Res. 24, 5407–5421. doi: 10.1158/1078-0432.CCR-18-0125

PubMed Abstract | CrossRef Full Text | Google Scholar

Cox, T. R. (2021). The matrix in cancer. Nat. Rev. Cancer Online ahead of print. doi: 10.1038/s41568-020-00329-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Erin, N., Grahovac, J., Brozovic, A., and Efferth, T. (2020). Tumor microenvironment and epithelial mesenchymal transition as targets to overcome tumor multidrug resistance. Drug Resist. Updat. 53:100715. doi: 10.1016/j.drup.2020.100715

PubMed Abstract | CrossRef Full Text | Google Scholar

Ershaid, N., Sharon, Y., Doron, H., Raz, Y., Shani, O., Cohen, N., et al. (2019). NLRP3 inflammasome in fibroblasts links tissue damage with inflammation in breast cancer progression and metastasis. Nat. Commun. 10:4375. doi: 10.1038/s41467-019-12370-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Farhood, B., Najafi, M., and Mortezaee, K. (2019). Cancer-associated fibroblasts: secretions, interactions, and therapy. J. Cell Biochem. 120, 2791–2800. doi: 10.1002/jcb.27703

PubMed Abstract | CrossRef Full Text | Google Scholar

Feig, C., Jones, J. O., Kraman, M., Wells, R. J., Deonarine, A., Chan, D. S., et al. (2013). Targeting CXCL12 from FAP-expressing carcinoma-associated fibroblasts synergizes with anti-PD-L1 immunotherapy in pancreatic cancer. Proc. Natl. Acad. Sci. U.S.A. 110, 20212–20217. doi: 10.1073/pnas.1320318110

PubMed Abstract | CrossRef Full Text | Google Scholar

Gonzalez-Gugel, E., Saxena, M., and Bhardwaj, N. (2016). Modulation of innate immunity in the tumor microenvironment. Cancer Immunol. Immunother. 65, 1261–1268. doi: 10.1007/s00262-016-1859-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Hambardzumyan, D., Gutmann, D. H., and Kettenmann, H. (2016). The role of microglia and macrophages in glioma maintenance and progression. Nat. Neurosci. 19, 20–27. doi: 10.1038/nn.4185

PubMed Abstract | CrossRef Full Text | Google Scholar

Hoshida, Y., Brunet, J. P., Tamayo, P., Golub, T. R., and Mesirov, J. P. (2007). Subclass mapping: identifying common subtypes in independent disease data sets. PLoS One 2:e1195. doi: 10.1371/journal.pone.0001195

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, Y. H., Zhu, C., Kondo, Y., Anderson, A. C., Gandhi, A., Russell, A., et al. (2015). CEACAM1 regulates TIM-3-mediated tolerance and exhaustion. Nature 517, 386–390. doi: 10.1038/nature13848

PubMed Abstract | CrossRef Full Text | Google Scholar

Ishida, Y., Agata, Y., Shibahara, K., and Honjo, T. (1992). Induced expression of PD-1, a novel member of the immunoglobulin gene superfamily, upon programmed cell death. EMBO J. 11, 3887–3895. doi: 10.1002/j.1460-2075.1992.tb05481.x

CrossRef Full Text | Google Scholar

Jackson, C. M., Choi, J., and Lim, M. (2019). Mechanisms of immunotherapy resistance: lessons from glioblastoma. Nat. Immunol. 20, 1100–1109. doi: 10.1038/s41590-019-0433-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Jacobs, J., Deschoolmeester, V., Zwaenepoel, K., Flieswasser, T., Deben, C., Van den Bossche, J., et al. (2018). Unveiling a CD70-positive subset of cancer-associated fibroblasts marked by pro-migratory activity and thriving regulatory T cell accumulation. Oncoimmunology 7:e1440167. doi: 10.1080/2162402X.2018.1440167

PubMed Abstract | CrossRef Full Text | Google Scholar

Kane, J. R., Zhao, J., Tsujiuchi, T., Laffleur, B., Arrieta, V. A., Mahajan, A., et al. (2020). CD8(+) T-cell-Mediated immunoediting influences genomic evolution and immune evasion in murine gliomas. Clin. Cancer Res. 26, 4390–4401. doi: 10.1158/1078-0432.CCR-19-3104

PubMed Abstract | CrossRef Full Text | Google Scholar

Katanov, C., Lerrer, S., Liubomirski, Y., Leider-Trejo, L., Meshel, T., Bar, J., et al. (2015). Regulation of the inflammatory profile of stromal cells in human breast cancer: prominent roles for TNF-alpha and the NF-kappaB pathway. Stem Cell Res. Ther. 6:87. doi: 10.1186/s13287-015-0080-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Kato, T., Noma, K., Ohara, T., Kashima, H., Katsura, Y., Sato, H., et al. (2018). Cancer-associated fibroblasts affect intratumoral CD8(+) and FoxP3(+) T Cells Via IL6 in the tumor microenvironment. Clin. Cancer Res. 24, 4820–4833. doi: 10.1158/1078-0432.CCR-18-0205

PubMed Abstract | CrossRef Full Text | Google Scholar

Krummel, M. F., and Allison, J. P. (1996). CTLA-4 engagement inhibits IL-2 accumulation and cell cycle progression upon activation of resting T cells. J. Exp. Med. 183, 2533–2540. doi: 10.1084/jem.183.6.2533

PubMed Abstract | CrossRef Full Text | Google Scholar

Lakins, M. A., Ghorani, E., Munir, H., Martins, C. P., and Shields, J. D. (2018). Cancer-associated fibroblasts induce antigen-specific deletion of CD8 (+) T Cells to protect tumour cells. Nat. Commun. 9:948. doi: 10.1038/s41467-018-03347-0

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

Lappano, R., Talia, M., Cirillo, F., Rigiracciolo, D. C., Scordamaglia, D., Guzzi, R., et al. (2020). The IL1beta-IL1R signaling is involved in the stimulatory effects triggered by hypoxia in breast cancer cells and cancer-associated fibroblasts (CAFs). J. Exp. Clin. Cancer Res. 39:153. doi: 10.1186/s13046-020-01667-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, M., Li, G., Kiyokawa, J., Tirmizi, Z., Richardson, L. G., Ning, J., et al. (2020). Characterization and oncolytic virus targeting of FAP-expressing tumor-associated pericytes in glioblastoma. Acta Neuropathol. Commun. 8:221. doi: 10.1186/s40478-020-01096-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Linares, J., Marin-Jimenez, J. A., Badia-Ramentol, J., and Calon, A. (2020). Determinants and functions of CAFs secretome during cancer progression and therapy. Front. Cell Dev. Biol. 8:621070. doi: 10.3389/fcell.2020.621070

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, B., Chen, X., Zhan, Y., Wu, B., and Pan, S. (2020). Identification of a gene signature for renal cell carcinoma-associated fibroblasts mediating cancer progression and affecting prognosis. Front. Cell Dev. Biol. 8:604627. doi: 10.2139/ssrn.3688868

CrossRef Full Text | Google Scholar

Liu, T., Han, C., Wang, S., Fang, P., Ma, Z., Xu, L., et al. (2019a). Cancer-associated fibroblasts: an emerging target of anti-cancer immunotherapy. J. Hematol. Oncol. 12:86. doi: 10.1186/s13045-019-0770-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, T., Zhou, L., Li, D., Andl, T., and Zhang, Y. (2019b). Cancer-associated fibroblasts build and secure the tumor microenvironment. Front. Cell Dev. Biol. 7:60. doi: 10.3389/fcell.2019.00060

PubMed Abstract | CrossRef Full Text | Google Scholar

Louault, K., Li, R. R., and DeClerck, Y. A. (2020). Cancer-associated fibroblasts: understanding their heterogeneity. Cancers 12:3108. doi: 10.3390/cancers12113108

PubMed Abstract | CrossRef Full Text | Google Scholar

Maman, S., and Witz, I. P. (2018). A history of exploring cancer in context. Nat. Rev. Cancer 18, 359–376. doi: 10.1038/s41568-018-0006-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Manieri, N. A., Chiang, E. Y., and Grogan, J. L. (2017). TIGIT: a key inhibitor of the cancer immunity cycle. Trends Immunol. 38, 20–28. doi: 10.1016/j.it.2016.10.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Miranda, A., Hamilton, P. T., Zhang, A. W., Pattnaik, S., Becht, E., Mezheyeuski, A., et al. (2019). Cancer stemness, intratumoral heterogeneity, and immune response across cancers. Proc. Natl. Acad. Sci. U.S.A. 116, 9020–9029. doi: 10.1073/pnas.1818210116

PubMed Abstract | CrossRef Full Text | Google Scholar

Mirzaei, R., Sarkar, S., and Yong, V. W. (2017). T cell exhaustion in glioblastoma: intricacies of immune checkpoints. Trends Immunol. 38, 104–115. doi: 10.1016/j.it.2016.11.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Moncada, R., Barkley, D., Wagner, F., Chiodin, M., Devlin, J. C., Baron, M., et al. (2020). Integrating microarray-based spatial transcriptomics and single-cell RNA-seq reveals tissue architecture in pancreatic ductal adenocarcinomas. Nat. Biotechnol. 38, 333–342. doi: 10.1038/s41587-019-0392-8

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Öhlund, D., Elyada, E., and Tuveson, D. (2014). Fibroblast heterogeneity in the cancer wound. J. Exp. Med. 211, 1503–1523. doi: 10.1084/jem.20140692

PubMed Abstract | CrossRef Full Text | Google Scholar

Postow, M. A., Callahan, M. K., and Wolchok, J. D. (2015). Immune checkpoint blockade in cancer therapy. J. Clin. Oncol. 33, 1974–1982. doi: 10.1200/JCO.2014.59.4358

PubMed Abstract | CrossRef Full Text | Google Scholar

Pyonteck, S. M., Akkari, L., Schuhmacher, A. J., Bowman, R. L., Sevenich, L., Quail, D. F., et al. (2013). CSF-1R inhibition alters macrophage polarization and blocks glioma progression. Nat. Med. 19, 1264–127237. doi: 10.1038/nm.3337

PubMed Abstract | CrossRef Full Text | Google Scholar

Racle, J., de Jonge, K., Baumgaertner, P., Speiser, D. E., and Gfeller, D. (2017). Simultaneous enumeration of cancer and immune cell types from bulk tumor gene expression data. Elife 6:e26476. doi: 10.7554/eLife.26476.049

CrossRef Full Text | Google Scholar

Roh, W., Chen, P. L., Reuben, A., Spencer, C. N., Prieto, P. A., Miller, J. P., et al. (2017). Integrated molecular analysis of tumor biopsies on sequential CTLA-4 and PD-1 blockade reveals markers of response and resistance. Sci. Transl. Med. 9:eaah3560. doi: 10.1126/scitranslmed.aah3560

PubMed Abstract | CrossRef Full Text | Google Scholar

Ruffo, E., Wu, R. C., Bruno, T. C., Workman, C. J., and Vignali, D. A. A. (2019). Lymphocyte-activation gene 3 (LAG3): the next immune checkpoint receptor. Semin. Immunol. 4:101305. doi: 10.1016/j.smim.2019.101305

PubMed Abstract | CrossRef Full Text | Google Scholar

Shintani, Y., Fujiwara, A., Kimura, T., Kawamura, T., Funaki, S., Minami, M., et al. (2016). IL-6 secreted from cancer-associated fibroblasts mediates chemoresistance in NSCLC by increasing epithelial-mesenchymal transition signaling. J. Thorac. Oncol. 11, 1482–1492. doi: 10.1016/j.jtho.2016.05.025

PubMed Abstract | CrossRef Full Text | Google Scholar

Simon, N., Friedman, J., Hastie, T., and Tibshirani, R. (2011). Regularization paths for Cox’s proportional hazards model via coordinate descent. J. Stat. Softw. 39, 1–13. doi: 10.18637/jss.v039.i05

PubMed Abstract | CrossRef Full Text | Google Scholar

Su, S., Chen, J., Yao, H., Liu, J., Yu, S., Lao, L., et al. (2018). CD10(+)GPR77(+) cancer-associated fibroblasts promote cancer formation and chemoresistance by sustaining cancer stemness. Cell 172, 841–856.e816. doi: 10.1016/j.cell.2018.01.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Takenaka, M. C., Gabriely, G., Rothhammer, V., Mascanfroni, I. D., Wheeler, M. A., Chao, C. C., et al. (2019). Control of tumor-associated macrophages and T cells in glioblastoma via AHR and CD39. Nat. Neurosci. 22, 729–740. doi: 10.1038/s41593-019-0370-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Tan, A. C., Ashley, D. M., López, G. Y., Malinzak, M., Friedman, H. S., and Khasraw, M. (2020). Management of glioblastoma: state of the art and future directions. CA Cancer J. Clin. 70, 299–312. doi: 10.3322/caac.21613

PubMed Abstract | CrossRef Full Text | Google Scholar

Tibshirani, R. (1997). The lasso method for variable selection in the Cox model. Stat. Med. 16, 385–395. doi: 10.1002/(SICI)1097-0258(19970228)16:4<385::AID-SIM380>3.0.CO;2-3

CrossRef Full Text | Google Scholar

Tirosh, I., Izar, B., Prakadan, S. M., Wadsworth, M. H., Treacy, D., Trombetta, J. J., et al. (2016). Dissecting the multicellular ecosystem of metastatic melanoma by single-cell RNA-seq. Science 352, 189–196. doi: 10.1126/science.aad0501

PubMed Abstract | CrossRef Full Text | Google Scholar

Vennin, C., Mélénec, P., Rouet, R., Nobis, M., Cazet, A. S., Murphy, K. J., et al. (2019). CAF hierarchy driven by pancreatic cancer cell p53-status creates a pro-metastatic and chemoresistant environment via perlecan. Nat. Commun. 10:3637. doi: 10.1038/s41467-019-10968-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, Z., Yang, Q., Tan, Y., Tang, Y., Ye, J., Yuan, B., et al. (2021). Cancer-associated fibroblasts suppress cancer development: the other side of the coin. Front. Cell Dev. Biol. 9:613534. doi: 10.3389/fcell.2021.613534

PubMed Abstract | CrossRef Full Text | Google Scholar

Wickham, H. (2016). ggplot2: Elegant Graphics for Data Analysis. New York: Springer-Verlag. doi: 10.1007/978-3-319-24277-4

CrossRef Full Text | Google Scholar

Woroniecka, K. I., Rhodin, K. E., Chongsathidkiet, P., Keith, K. A., and Fecci, P. E. (2018). T-cell dysfunction in glioblastoma: applying a new framework. Clin. Cancer Res. 24, 3792–3802. doi: 10.1158/1078-0432.CCR-18-0047

PubMed Abstract | CrossRef Full Text | Google Scholar

Xu, R., Zhou, X., Wang, S., and Trinkle, C. (2021). Tumor organoid models in precision medicine and investigating cancer-stromal interactions. Pharmacol. Ther. 218:107668. doi: 10.1016/j.pharmthera.2020.107668

PubMed Abstract | CrossRef Full Text | Google Scholar

Yoshihara, K., Shahmoradgoli, M., Martínez, E., Vegesna, R., Kim, H., Torres-Garcia, W., et al. (2013). Inferring tumour purity and stromal and immune cell admixture from expression data. Nat. Commun. 4:2612. doi: 10.1038/ncomms3612

PubMed Abstract | CrossRef Full Text | Google Scholar

Yu, G., Wang, L. G., Han, Y., and He, Q. Y. (2012). clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS 16, 284–287. doi: 10.1089/omi.2011.0118

PubMed Abstract | CrossRef Full Text | Google Scholar

Zarour, H. M. (2016). Reversing T-cell dysfunction and exhaustion in cancer. Clin. Cancer Res. 22, 1856–1864. doi: 10.1158/1078-0432.CCR-15-1849

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, Y., Liu, L., Qiu, Q., Zhou, Q., Ding, J., Lu, Y., et al. (2021). Alternative polyadenylation: methods, mechanism, function, and role in cancer. J. Exp. Clin. Cancer Res. 40:51. doi: 10.1186/s13046-021-01852-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, Z., Zhang, K. N., Wang, Q., Li, G., Zeng, F., Zhang, Y., et al. (2021). Chinese Glioma Genome Atlas (CGGA): a comprehensive resource with functional genomic data from Chinese Gliomas. Genomics Proteomics Bioinform. 19, 1–12. doi: 10.1016/j.gpb.2020.10.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhen, Z., Tang, W., Wang, M., Zhou, S., Wang, H., Wu, Z., et al. (2017). Protein nanocage mediated fibroblast-activation protein targeted photoimmunotherapy to enhance cytotoxic T Cell infiltration and tumor control. Nano Lett. 17, 862–869. doi: 10.1021/acs.nanolett.6b04150

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhu, Z., Zhang, H., Chen, B., Liu, X., Zhang, S., and Zong, Z. (2020). PD-L1-mediated immunosuppression in glioblastoma is associated with the infiltration and M2-polarization of tumor-associated macrophages. Front. Immunol. 11:588552. doi: 10.3389/fimmu.2020.588552

PubMed Abstract | CrossRef Full Text | Google Scholar