Pericyte mediates the infiltration, migration, and polarization of macrophages by CD163/MCAM axis in glioblastoma.

Microenvironment cells (MCs) play a critical role in tumor proliferation, progression, and metastasis. However, it has not been adequately addressed whether MCs could be used as a reliable prognostic marker in glioblastoma (GBM). In the current study, the cell pair (CP) score was constructed in 1137 GBM samples based on the cell pair algorithm and Gaussian finite mixture model (GMM) and was verified in 73 GBM samples from the Xiangya cohort. CP score predicted GBM patients' survival and response to anti-PD-1 treatment with high sensitivity. Macrophage markers CD68 and CD163 were co-expressed with pericyte markers MCAM and MG2. Pericyte could mediate the infiltration, migration, and M2 type polarization of macrophages by MCAM. The CP score was a valuable tool for predicting survival outcomes and guiding immunotherapy for GBM patients. Cell pair pericyte/macrophage and gene pair CD163/MCAM were biologically significant in the tumor microenvironment of GBM.

Introduction

The World Health Organization (WHO) classification defines grade I and II glioma as low-grade glioma (LGG). It defines grade III and IV glioma as high-grade glioma (HGG)(Louis et al., 2016), among which glioblastoma (GBM) has been recognized as the most devastating primary brain tumor with an extremely high mortality rate. Typically, the 10-year survival rate of patients with LGG is 47%, with a median survival time of 11.6 years (Smoll et al., 2012), while the median survival time of patients with GBM is less than 15 months (Yang et al., 2022). Despite the surgical resection with adjuvant chemoradiotherapy, the prognosis of GBM patients remains poor. So far, biomarkers including isocitrate dehydrogenase (IDH), 1p19q, and O-6-methylguanine-DNA methyltransferase (MGMT) have been used for precise classification of GBM patients to promote the clinical management and fulfill individualized treatment (Neftel et al., 2019; Wang et al., 2017).

Microenvironment cells (MCs), including T cells, mast cells, tumor-associated macrophages (TAMs), cancer-associated fibroblasts (CAFs), and natural killer (NK) cells, could elicit a robust immune response against tumor (Hiraoka, 2010; Zhang et al., 2020, 2021b). MCs play a central role in regulating the immunosurveillance of cancer (Beatty and Gladney, 2015) and creating a permissive microenvironment that accelerates tumor progression(Schreiber et al., 2011). Notably, MCs have been proposed to be either the mediator of immunotherapy (Peng et al., 2020) or the immunotherapeutic target(Sabado et al., 2017). Besides, as the rapid development of bioinformatics provides insight into cancer research based on large-scale analysis, numerous studies have established the microenvironment cell-based risk signatures in various cancer types. However, the comprehensive role of MCs in the tumor microenvironment (TME) of GBM lacks in-depth understanding (Chen and Hambardzumyan, 2018). Therefore, developing an MC-based signature can be promising to help determine the prognostic value of MCs in GBM and improve the efficacy of immunotherapeutic approaches. However, although algorithms such as xCell, CIBERSORT, and TIMER have been developed to quantify the expression level of MCs based on bulk sequencing datasets to facilitate the research on MCs, these methods are restricted by the different reference genomes that might lead to discrepant research results from other studies. Some computational deconvolution methods might even require validation of the results using single-cell genomic and/or imaging methods. Given that the fraction of each MC in the TME is within a relatively stable range due to intra-sample heterogeneity or inter-sample heterogeneity, exploring the ratio of different MCs could potentially optimize the quantification of MCs in the research of TME.

In this study, 65 immune cell types were collected to construct a risk signature (Zhang et al., 2021c). These cell types were quantified in GBM cohorts to estimate their prognostic value. A cell pair (CP) score was constructed based on the relative abundance of identified immune cell types. As a result, a high CP score predicted worse overall survival in GBM patients. Besides, the CP score was profoundly associated with various tumorigenic and immunogenic factors and could sensitively predict response to anti-PD-1 immunotherapy. Meanwhile, cell pair macrophage/pericyte and gene pair CD163/MCAM was expected to be the potential prognostic markers and therapeutic targets for GBM. Pericyte was further found to mediate the infiltration, migration, and polarization of macrophages in the TME of GBM.

Results

Construction of CP score and its prognostic value

The overall study design was shown in Figure S1. 26 prognostic immune cell types were identified by Univariate Cox regression analysis in TCGA GBM samples and paired with 65 integrated immune cell types collected from previously published studies. Each cell pair was assigned 1 or 0 as the score according to the relative expression level. After calculating the 2-year area under the curve (AUC) of all cell pairs, 13 immune cell pairs were identified with the highest 2-year AUC. After performing GMM, the CP score based on 6 immune cell pairs finally stood out with the highest AUC (Figure 1A). The HR of the 13 cell pairs with the highest 2-year AUC value was shown in Figure 1B. AUC of the CP score models sorted by the GMM classifier in all the 8191 formulas was shown in Figure 1C. CP score predicted worse survival in LGG samples, GBM samples, and pan-glioma samples from TCGA (logrank test, p < 0.001; Figures 1D–1F, respectively). Moreover, CP score was a hazardous factor in GBM samples and pan-glioma samples from the Xiangya cohort (logrank test, p < 0.001; Figures 1G and 1H, respectively). The ROC analyses with the 2-year, 3-year, 4-year, 5-year AUC of 0.703, 0.738, 0.767, 0.797 confirmed that CP score served as a prognostic marker in predicting the survival status of GBM patients from TCGA (Figure 1I). The CP score was further constructed and verified in meta-cohort (Figure S2A) and each of the six GBM datasets (Figure S2B). Notably, the CP score significantly stratified GBM patients’ survival and predicted worse outcomes. The Univariate Cox regression analyses confirmed that CP score was a hazardous factor in all six GBM datasets and meta-cohort (Figure S2C). The Univariate Cox regression analyses and Multivariate Cox regression analyses further proved that CP scores and age were hazardous factors in TCGA and CGGA. At the same time, IDH mutation and 1p19q codeletion were favorable factors (Figure S2D).

Figure 1

Construction of CP score

(A) Flow diagram of the cell pair algorithm.

(B) Forest plot depicting the 13 cell pairs.

(C) The pattern of the logistic regression model correlated with the AUC values and was identified by a Gaussian mixture. There are nine clusters of 8191 combinations. Kaplan–Meier curves for two CP score groups in (D). LGG samples, (E). GBM samples, and F. glioma samples in the TCGA dataset. Logrank test, p< 0.001.

(G) Kaplan–Meier curves for two CP score groups in GBM samples from the Xiangya cohort. Logrank test, p< 0.001.

(H) Kaplan–Meier curves for two CP score groups in glioma samples from the Xiangya cohort. Logrank test, p< 0.001.

(I) ROC curve measuring the sensitivity of CP score in predicting GBM patient’s 2-year, 3-year, 4-year, and 5-year survival in TCGA dataset. The area under the ROC curve was 0.703, 0.738, 0.767, and 0.797, respectively.

Genomic features of the CP score groups

Given the prognostic value of the CP score, somatic mutation analysis and copy number variation (CNV) analysis were performed to characterize the genomic features of the CP score in TCGA. A global CNV profile exhibited insignificant differences regarding mutation frequency in two CP score groups (Figure S3A and Table S2). Specifically, in the high CP score group, the high CP score group, and the mutation rates of PTEN, TP53, EGFR, and TTN were 33%, 27%, 26%, and 26%, respectively, in the high CP score group (Figure S3B). On the contrary, mutation rates of TP53, TTN, and EGFR were 43%, 35%, and 22%, respectively, in the low CP score group (Figure S3C). Besides, the CP score showed high sensitivity in predicting the mutation status of IDH with an AUC of 0.778 (Figure S3D).

Immune regulatory mechanisms related to CP score

Metabolism has been proposed to influence the activation and quiescence of immune cell (Pearce and Pearce, 2013). A high CP score was significantly associated with glutathione metabolism, kynurenine metabolism, and prostanoid biosynthesis in the Xiangya cohort and meta-cohort (Figures 2A and S4A, respectively). The cancer immunity cycle has been proposed to comprehensively reflect the functions of several chemokines and immunomodulators (Chen and Mellman, 2013; Xu et al., 2018). Of note, most of the steps in the cancer immunity cycle were upregulated in the high CP score group, including the release of cell antigens (Step 1), tumor antigen presentation (Step 2), recruitment of immune cells (Step 4) (CD8 T cell, dendritic cell, Macrophage, myeloid-derived suppressor cell (MDSC), monocyte, neutrophil, NK cell, Th1 cell, Th17 cell, Th22 cell), and infiltration of immune cells into tumors (Step 5) in Xiangya cohort and meta-cohort (Figures 2B and S4B, respectively). The enhanced activities of the multiple steps in the cancer immunity cycle may subsequently increase the infiltration levels of MCs in the high CP score group. The RNA modification N6-methyladenosine (m6A) has been proved to significantly regulate immune recognition, immune activation, and immune cell fate decisions (Shulman and Stern-Ginossar, 2020). Of interest, the high CP score group had a higher level of CBLL1, ELAVL1, WTAP, RBM15B, YTHDC2, YTHDF1, YTHDF2, and YTHDF3 in meta-cohort (Figure S4E), while the expression differences of m6A regulators were not significant in Xiangya cohort (Figure S4C). Besides, a high CP score group was associated with inflammatory activities in the Xiangya cohort and meta-cohort (Figures S4D and S4F, respectively).

Figure 2

Immunogenic and tumorigenic characteristics of CP score in TCGA

(A) Heatmap illustrates the expression pattern of metabolic signatures in the CP score in the Xiangya cohort.

(B) Differences in the various steps of the cancer immunity cycle between high- and low-CP score groups in Xiangya cohort.

(C) GEP score in high- and low-CP score groups.

(D) HRD in high- and low-CP score groups.

(E)Number of segments in high- and low-CP score groups.

(F) TGF-beta response in high- and low-CP score groups.

(G) Leukocyte fraction in high- and low-CP score groups.

(H) Stromal fraction in high- and low-CP score groups.

(I) IFNG score in high- and low-CP score groups.

(J) IFNG.GS score in high- and low-CP score groups.

(K) ISG.RS scores in high- and low-CP score groups.

(L) TCR Shannon in high- and low-CP score groups.

(M) TCR richness in high- and low-CP score groups.

(N) APM score in high- and low-CP score groups.

The immune escape mechanisms included four significant aspects: tumor immunogenicity, antigen presentation capacity, microenvironment cell, and regulation of immune checkpoint (Schumacher and Schreiber, 2015). A series of tumorigenic and immunogenic factors were first evaluated. The high CP score group exhibited a higher T cell-inflamed gene expression profile (GEP), indicating higher response rates of anti-PD-1 therapy (Figure 2C) (Ayers et al., 2017). The high CP score group also displayed lower homologous recombination deficiency (HRD), an indicator of cell death (Figure 2D). Interestingly, a low CP score group correlated with more segments (Figure 2E). Stroma signatures, including leukocyte fraction, stromal fraction, TGF-β response, interferon-gamma (IFNG), IFNG.GS, and ISG.RS-related signatures, were higher in the high CP score group (Figures 2F–2K, respectively). In terms of antigen presentation capacity, the high CP score group exhibited higher levels of T cell receptor (TCR) Shannon, TCR Richness, and higher antigen processing and presenting machinery (APM) score (Figures 2L–2N, respectively).

The immune infiltration characteristics of CP score groups were also assessed. As a result, the high CP score group correlated with a higher ESTIMATE score, immune score, and stromal score in the Xiangya cohort (Figure 3A). Based on three various algorithms, the high CP score group was significantly associated with immune suppressive cells, including regulatory T cells (Tregs), TAMs, CAFs, and T helper 2 cells (Th2) in the Xiangya cohort (Figure 3A). The association between CP score and seven types of immune checkpoint molecules was next explored (Schreiber et al., 2011; Wang et al., 2020). The high CP score group had significantly higher expression of the majority of immune checkpoint molecules, including ICOSLG, PDCD1, CTLA4, and CD40, and might potentially evade immune response via these classical immune checkpoint molecules in the Xiangya cohort (Figure 3B).

Figure 3

The predictive value of CP score in immunotherapy

(A) Heatmap illustrates the expression pattern of immune infiltrating cells in CP score.

(B) Heatmap illustrates the expression pattern of seven types of immunomodulators in the CP score.

(C) Expression differences of CP score in patients with or without anti-PD-1 responses. p = 0.012.

(D) Kaplan–Meier curves for two CP score groups in the GSE78220 dataset. Logrank test, p< 0.1412.

(E) CP score in groups with different anti–PD-1 clinical response statuses (CR/PR and SD/PD). Statistical significance was based on the Wilcoxon test (Wilcoxon, p =0.036).

(F) Kaplan–Meier curves for two CP score groups in the IMvigor210 dataset. Logrank test, p = 0.00174.

(G) CP score in groups with different anti–PD-1 clinical response statuses (CR, PR, SD, PD). Statistical significance was based on the Kruskal–Wallis test (Kruskal-Wallis, p =0.013).

The GO results of GSVA confirmed that tumorigenic pathways, including regulation of ERBB signaling pathway, Toll-like receptor signaling pathway, NF-kB transcription factor activity, glial cell activation, and immunogenic pathways, including regulation of macrophage, chemokine production, mast cell activation, were more activated in high CP score group in Xiangya cohort (Figure S5A). The GO results of GSVA were also validated in TCGA and meta-cohort (Table S3). Additionally, the high CP score group was significantly associated with PD-1 therapeutic effect, T cell signaling, Hypoxia signaling, exosome signaling, and signalings of immunosuppressive cells in the Xiangya cohort (Figure S5B).

CP score predicted immunotherapeutic responses

Immunotherapy has revolutionized cancer treatment. Therefore, the predictive value of CP score in immunotherapeutic response was also explored. In the cohort exploring the anti-PD-1 immunotherapy response in GBM patients, patients responding to anti-PD-1 immunotherapy were less likely to present a high CP score (Figure 3C). A high CP score in the melanoma dataset, GSE78220, predicted a worse survival outcome (Figure 3D). Likewise, patients with stable and progressive diseases were likely to present high CP scores (Figure 3E). CP score was also constructed in the IMvigor210 cohort (urothelial carcinoma dataset). A high CP score predicted a worse survival outcome (Figure 3F). Patients with stable and progressive diseases were likely to present high CP scores (Figure 3G).

Functional annotation of cell pair macrophage/pericyte

Subsequently, we identified macrophage/pericyte as the most prognosis-related cell pair based on the 2-year AUC value. GBM patients were separated into two cell groups based on the relative abundance of macrophage and pericyte, in which macrophage > pericyte was termed cell group M and macrophage < pericyte was termed cell group P. The GO results of GSVA confirmed that immunogenic pathways were more activated in the cell group M in TCGA (Figure S6A). And also, cell group M was significantly associated with multiple metabolic pathways (Figure S6B) and inflammatory activities (Figure S6C) in TCGA. Based on six various algorithms, the two cell groups were significantly differentially associated with immune suppressive cells, including Tregs, TAMs, and CAFs (Figure 4A). In addition, cell group M had significantly higher expression of most immune modulator molecules, including TGF-β, CD27, CD80, and CD40 (Figure S7A). The expression differences of immune checkpoint molecules in cell groups were absent from somatic mutation and CNV but were strongly associated with methylation (Figure S7B). Further, cell group M predicted worse survival in TCGA (Figure 4C). Moreover, the univariate and multivariate cox regression analysis proved that cell pair macrophage/pericyte was a valuable prognostic marker in GBM patients (Figure S8B).

Figure 4

Prognostic value of gene pair CD163/MCAM

(A) Estimation of the correlation between immune cells and cell pair macrophage/pericyte in different algorithms in TCGA.

(B) Estimation of the correlation between immune cells and gene pair CD163/MCAM in different algorithms in TCGA.

(C) Kaplan–Meier curves for two cell groups in TCGA. Logrank test, p< 0.001.

(D) Kaplan–Meier curves for two gene groups in TCGA. Logrank test, p = 0.00115.

(E) Kaplan–Meier curves for two gene groups based on sequencing data of Xiangya cohort. Logrank test, p = 0.04908.

(F) Kaplan–Meier curves for two gene groups based on IHC staining of Xiangya cohort. Logrank test, p = 0.02014.

(G) IHC staining for CD163 and MCAM in four representative samples of the Xiangya cohort (scale bar: 100um for 10X and 25um for 40X).

Functional annotation of gene pair CD163/MCAM

After pairing the specific markers of macrophage and pericyte, CD163/MCAM was identified as the most prognostic related gene pair based on the 2-year AUC value. GBM patients were then separated into two gene groups based on the relative abundance of CD163 and MCAM, in which CD163 > MCAM was termed a high group, and CD163 < MCAM was termed a low group. The GO results of GSVA confirmed that immunogenic pathways were more activated in the high group in TCGA (Figure S9A). The high group was significantly associated with multiple metabolic pathways in TCGA (Figure S9B). Besides, the high group was related to inflammatory activities in TCGA (Figure S9C). Based on six various algorithms, the two gene groups were significantly differentially associated with immune suppressive cells, including Tregs, TAMs, and CAFs (Figure 4B). The high group positively correlated with most immune modulator molecules, including CD80, TGF-β, and CD40 (Figure S10A). The expression difference of immune checkpoint molecules in cell groups was absent from somatic mutation and CNV but was strongly associated with methylation (Figure S10B).

Furthermore, the high group predicted worse survival in TCGA (Figure 4D). The prognostic value of gene pair CD163/MCAM was also explored in the five included datasets. Notably, the high group was associated with significantly decreased survival in four datasets (Figure S8A). Furthermore, the univariate and multivariate cox regression analysis proved that gene pair CD163/MCAM was a valuable prognostic marker in GBM patients (Figure S8C).

Validation of gene pair CD163/MCAM in Xiangya cohort

The high group was associated with decreased survival in the sequencing data of 73 GBM samples from the Xiangya cohort (Figure 4E). Besides, IHC staining was performed on the 45 GBM samples from the Xiangya cohort. Based on the H-scores of CD163 and MCAM in IHC staining results, 45 GBM samples were divided into a high group (CD163>MCAM) and a low group (CD163

Besides, the high group was significantly associated with PD-1 therapeutic effect, T cell signaling, hypoxia signaling, exosome signaling, and signalings of immunosuppressive cells in the Xiangya cohort (Figure S11A). Based on six various algorithms, the high group was significantly associated with immune suppressive cells such as TAMs and CAFs (Figure S11B). Additionally, the high group had higher classical immune modulator molecule expression, including TIM3, TIGIT, PD-L2, PD-L1, PD-1, LAG3, IDO1, and CTLA-4 (Figure S11C). The enhanced activities of the multiple steps in the cancer immunity cycle were also observed in the high group (Figure S11D).

Pericytes mediate the infiltration, migration, and polarization of macrophages

To further elucidate the biological significance of gene pair CD163/MCAM and cell pair macrophage/pericyte, multiplex immunofluorescence staining of macrophage markers CD68, CD163, and pericyte markers MCAM, NG2 was performed in GBM samples. Notably, macrophage markers CD68 and CD163 were almost coexpressed with pericyte markers MCAM and NG2 in the tumor region and perivascular area (Figures 5A and S12). Moreover, the expression of M2 macrophage marker CD163 decreased along with the reduced expression of pericyte markers MCAM and NG2, which suggested that pericytes might biologically regulate the activity of macrophages. NG2/MCAM recognized pericytes, and cells at different distances from pericytes were quantified (Figures 5B and 5E). The number of DAPI⁺ cells, CD68⁺ cells, CD163⁺ cells, CD68⁺CD163⁺ cells, NG2⁺ cells, MCAM⁺ cells, and NG2⁺MCAM⁺ cells was calculated (Figures 5C and 5F). Notably, increasing active pericytes were potentially surrounded by M2 macrophages (Figures 5D and 5G). A co-cultured system between macrophage and pericyte was designed (Figure 6B). Based on the ssGSEA algorithm, the ratio of macrophages and pericytes in the TME of GBM was approximately 1:1 (Figure 6A). As microglia share thea similar property with macrophages in the CNS, the HMC3 cell line was used. So, HMC3 cells were co-cultured with primary pericyte cells at the ratio of 1:1 in the subsequent experiments. The qPCR results showed that siRNA-1 and siRNA-2 could significantly suppress the expression of MCAM (Figure 6C). The western blotting results further proved the knock-down efficiency of this two siRNA (Figure 6D). In the co-cultured system for transwell assay, HMC3 cells had significantly decreased ability in migration in two siRNA groups compared with the control group (Figure 6E). In the co-cultured system for multiplex immunofluorescence staining, HMC3 cells were more likely to polarize into M1 macrophages in two siRNA groups (Figure 6F). Correspondingly, HMC3 cells were less likely to polarize into M2 macrophages in two siRNA groups (Figure 6G). The above results suggested that pericytes could potentially mediate the M2-type polarization of macrophages in the TME of GBM.

Figure 5

Interaction between macrophages and pericytes

(A)Multiplex immunofluorescence staining of macrophage markers CD68, CD163, and pericyte markers MCAM, NG2 in GBM samples.

(B) StrataQuest software recognizes pericytes by NG2/MCAM. Quantify the number of cells at different distances from pericytes: 0–25μm (green area), 25–50μm (brown-yellow area).

(C) The scatterplot shows the number of DAPI⁺ cells, CD68⁺ cells, CD163⁺ cells, CD68⁺CD163⁺ cells, NG2⁺ cells, MCAM⁺ cells, NG2⁺MCAM⁺ cells.

(D) The scatterplot shows the number of CD68⁺CD163⁺ cells at different locations from NG2⁺MCAM⁺ cells: 0–25μm, 25–50μm.

(E) StrataQuest software recognizes pericytes by NG2/MCAM. Quantify the number of cells at different distances from pericytes: 0–25μm (green area), 25–50μm (brown-yellow area).

(F) The scatterplot shows the number of DAPI⁺ cells, CD68⁺ cells, CD163⁺ cells, CD68⁺CD163⁺ cells, NG2⁺ cells, MCAM⁺ cells, and NG2⁺MCAM⁺ cells.

(G) The scatterplot shows the number of CD68⁺CD163⁺ cells at different locations from NG2⁺MCAM⁺ cells: 0–25μm, 25–50μm.

Figure 6

Pericytes mediate the migration and polarization of macrophages

(A) The relative expression of macrophages and pericytes in the TME of GBM based on the ssGSEA algorithm.

(B) The flow diagram of the co-culture system between macrophages and pericytes.

(C) The relative RNA expression of MCAM in the NC and three siRNA target groups.

(D) The protein expression of MCAM in the NC group and three siRNA target groups.

(E) Transwell assay for the co-cultured pericytes. Statistical analysis of the migrated pericytes in different siRNA groups.

(F)Multiplex immunofluorescence staining of M1 macrophage markers CD68, CD11c in the cocultured system between macrophages and pericytes. Statistical analysis of the positive M1 macrophages in different siRNA groups (scale bar: 25um).

(G)Multiplex immunofluorescence staining of M2 macrophage markers CD68, CD163 in the cocultured system between macrophages and pericytes (scale bar: 25um). Statistical analysis of the positive M2 macrophages in different siRNA groups. Data are represented as mean ±SD.

Discussion

MCs regulate both the immunosurveillance and immune escape of cancer cells. The prognostic value of MCs has been reported in various cancer types. However, the overall survival benefits of multiple MCs in GBM have not been fully explored. Moreover, considering the non-uniform reference genomes and immune cell signatures, previous immune cell-derived prognostic models are limited in the cross-validation of different transcriptomic datasets. To resolve this issue, we introduced the concept of cell pairs for constructing a prognostic immune signature. We explored the probability of using the relative abundance of immune cells to calculate the CP score, which extensively eliminated the need for data normalization and increased the accuracy of developing a signature.

In this study, CP score was established in GBM samples from TCGA and could significantly stratify the OS of GBM patients. In the Xiangya cohort, a high CP score was associated with inferior survival in GBM and glioma patients. Likewise, ha high CP score predicted decreased survival in CGGA and GEO datasets. CP score was a hazardous marker in all included GBM samples. Furthermore, the annotated genomic features of the CP score revealed that the IDH1 missense mutations more frequently occurred in the low CP score group (14%), in line with the fact that IDH mutations confer favorable survival outcomes in glioma patients(Yan et al., 2009). PTEN and EGFR, two of the most overrepresented mutated genes in the high CP score group (33 and 26%, respectively), have been frequently activated in GBM and confer terrible survival outcomes in glioma patients(Benitez et al., 2017; Cancer Genome Atlas Research Network, 2008). Surprisingly, TP53 (43%) and TTN (35%) were more enriched in the low CP score group. CP score could also sensitively predict the IDH mutation status and may thus be a potent predictor for the oncogenic process.

The immune characteristics of two CP score groups were summarized. The cancer immunity cycle was highly active in the high CP score group, indicating an increased level of MCs. Notably, immunosuppressive cells, including Tregs, MDSCs, TAMs, and CAFs, were more actively and abundantly observed in the high CP score group, indicating the potential immune evasion in the TME. Moreover, GBM samples with high CP scores presented higher tumor immunogenicity, antigen presentation capacity, and stroma signatures than GBM samples with low CP scores, all of which have been used as a reliable marker in predicting the immune escape of cancer (Tauriello et al., 2018). The significant correlation between CP score and classical immune checkpoint molecules such as ICOSLG, PDCD1, CTLA4, and CD40 also suggested that CP score could effectively predict immune checkpoint blockade (ICB) therapy response (Zhang et al., 2021a). The above findings suggested an orientation for including CP score as the indicator of immunosuppression. Immunotherapy, especially ICB, has tremendously revolutionized tumor treatment. Based on a recent clinical trial exploring the anti-PD-1 response rates in GBM, patients with higher CP scores exhibited lower response rates. Notably, high CP scores also predict worse survival outcomes and worse immunotherapy response rates in two widely used immunotherapy cohorts, the IMvigor210 cohort, and the GSE78220 dataset.

Subsequently, macrophage/pericyte was identified as the most prognostic related cell pair. Based on the specific markers of macrophage and pericyte, CD163/MCAM was further identified as the most prognostic corresponding gene pair. The immunogenic and tumorigenic features of macrophage/pericyte and CD163/MCAM revealed that macrophage and CD163 were more involved in promoting tumor progression and suppressing the immune response. Macrophage and CD163 correlated with more immune checkpoint molecules, and the difference in macrophage/pericyte and CD163/MCAM were not caused by mutation but were closely connected to methylation. The prognostic value of macrophage/pericyte and CD163/MCAM were also proved in both sequencing data and IHC staining of the Xiangya cohort.

The biological function of the interaction between macrophage and pericyte has been proved during tissue repair and inflammation (Minutti et al., 2019; Shibahara et al., 2020). However, whether the macrophage/pericyte axis plays a role in the TME of GBM remains unknown. The remarkable biological significance of cell pair macrophage/pericyte and gene pair CD163/MCAM in our bioinformatic analyses was further validated in our in vitro experiments. Surprisingly, macrophage markers CD68, CD163, and pericyte markers MCAM, NG2 was highly coexpressed in GBM samples. Pericytes also mediate the migration and M2-type polarization of macrophages. Taken together, pericyte was proposed to correlate with macrophage and regulate its biological activity in GBM. When the balance between pericyte and macrophage was broken, cell pair macrophage/pericyte and gene pair CD163/MCAM became prognostic biomarkers.

We comprehensively collected the immune cell types in the TME and introduced a cell pair algorithm for developing a robust immune signature in GBM. The immune signature could help identify GBM patients with better immunotherapy responses. Moreover, macrophage/pericyte and CD163/MCAM significantly affected GBM patients’ survival. Notably, pericytes could mediate the infiltration, migration, and M2 type polarization of macrophages in GBM.

Limitations of the study

The complex regulatory network leading to the prognostic cell pair macrophage/pericyte and gene pair CD163/MCAM remains further explored. An in vivo mouse model is needed to verify the biological significance of cell pair macrophage/pericyte and gene pair CD163/MCAM.

STAR★Methods

Key resources table

Resource availability

Lead contact

Further information could be directly consulted to the lead contact, Zhixiong Liu (zhixiongliu@csu.edu.cn).

Materials availability

This study did not generate new materials or reagents.

Experimental model and subject details

Cell lines

HMC3 cells were purchased from iCell and cultured in 1640 complete medium (Sigma, R8758-500ML) with 10% FBS at 37°C with 5% CO2. HUM-iCell-n011 were purchased from iCell and cultured in a specific medium (iCell, Primed-icell-015) at 37°C with 5% CO2.

Method details

Datasets collecting and preprocessing

The publicly available GBM cohorts were collected from The Cancer Genome Atlas (TCGA; https://xenabrowser.net/), Chinese Glioma Genome Atlas (CGGA; http://www.cgga.org.cn/), and Gene Expression Omnibus (GEO; https://www.ncbi.nlm.nih.gov/geo/). A total of 523 GBM patient samples were from TCGA. A total of 1137 GBM patient samples from 6 cohorts were defined as meta-cohort. The information on the platforms and patient samples are provided in Table S1.

Raw data from the microarray dataset generated by Agilent and Affymetrix was downloaded from the GEO database. Gene expression profiles generated by Illumina and corresponding clinical information were obtained from TCGA and CGGA. Our previous finding provided a detailed procedure for processing the raw data(Zhang et al., 2021d).

Datasets collecting of the Xiangya cohort

Glioma tissues were collected and written informed consent was obtained from all patients. The included glioma tissues were approved by the Ethics Committee of Xiangya Hospital, Central South University.

Immune cell signature collection

Immune cell signatures were integrated from publicly available resources(Zhang et al., 2021c). The gene sets of the immune cell types from different literature were combined and considered reliable.

Development of a reliable risk signature in GBM

Univariate Cox analysis was performed to screen out prognosis-associated immune cell types with prognostic values in the GBM dataset, TCGAGBM (523 samples). Prognosis-associated immune cell types (Ci) were paired with all 65 microenvironment cell types (Cj). For a cell pair started with Ci, Ci and Cj, Score_ij = 1 (exp_Ci – exp_Cj > 0) and Score_ij = 0 (exp_Ci – exp_Cj < 0). 2 year- Area Under the Curve (AUC) was adopted to estimate the performance of each Score_ij and find out the cell pair with statistically significant prognosis and highest 2 year-AUC (Harrell et al., 1996). For each Ci, Score_ij was identified with the highest 2 year-AUC. The identified cell pairs with the highest 2 year-AUC were further sorted with the hazard ratio (HR)>1, and duplicate cell pairs were removed. Subsequently, classification was conducted with cell pair model-based hierarchical agglomerative clustering based on the Gaussian finite mixture model (GMM). Then, the CP score was calculated with these selected Score_ij:

CP score = Σ Score_ij

GBM patients were grouped into high and low CP score groups based on the cutoff value of the CP score according to the optimal cut-off value determined by the survminer package.

Genomic alterations in CP score

Somatic mutations and copy number variations (CNVs) data corresponding to the GBM samples were downloaded from TCGA. Somatic mutations were visualized using R package maftools. CNVs associated with the two CP score groups were determined using GISTIC 2.0 analysis.

Functional annotation of CP score

The gene signatures for 115 metabolism-relevant signaling pathways were from published work(Rosario et al., 2018). Seven immune checkpoint molecules were from a previous study(Thorsson et al., 2018). The diverse immune modulators were collected(Thorsson et al., 2018). The xCell algorithm(Aran et al., 2017), TIMER algorithm(Zhang et al., 2017), EPIC algorithm(Racle and Gfeller, 2020), MCPcounter algorithm(Becht et al., 2016), quanTlseq algorithm(Moll et al., 2014), and CIBERSORT algorithm(Newman et al., 2015) were used for identifying microenvironment cells in the tumor microenvironment of GBM.

Prediction of CP score in immunotherapy response

GBM samples receiving anti-PD1 immunotherapy in the PRJNA482620 dataset were collected for evaluating the predictive value of CP score(Zhao et al., 2019). The IMvigor210 cohort (urothelial carcinoma cohort) and the GSE78220 (melanoma dataset) were further used to predict the immunotherapy response (Hugo et al., 2017; Wang et al., 2020). The processing of raw data from both datasets was based on the DEseq2 R package. The expression value of the raw matrix was transformed to the TPM value. CP score was calculated in these two cohorts, respectively.

Identification of cell pair macrophage/pericyte and gene pair CD163/MCAM

Based on the 2y-AUC, the most prognosis-related cell pair was explored. Functional annotation, including biological processes, metabolic pathways, inflammatory signature, and immune infiltration, was performed for the identified cell pair macrophage/pericyte.

CD31, NG2, PDGFR beta, MCAM, and Nestin were used as the markers for pericyte, while CD11b, CD68, CD163, CD14, and CD16 were used as the markers for macrophage. Markers from macrophages and markers from pericyte were then paired. The most prognosis-related gene pair was also explored based on 2y-AUC. Functional annotation, including biological processes, metabolic pathways, inflammatory signature, and immune infiltration, was performed for the identified gene pair CD163/MCAM.

Transcriptomic sequencing of Xiangya cohort

RNAstore-fixed tumor tissues from 105 glioma patients were used for RNA sequencing(Zhang et al., 2021d). The detailed procedures were as follows. RNA was sheared followed by sequencing library preparation using the NEBNext Ultra RNA Library Prep Kit. After targeting region capture by biotin-labeled probes, the captured libraries were sequenced on an Illumina Hiseq platform to generate 125/150 bp paired-end reads. In-house perlscripts were used to process raw data (raw reads). Then, reads containing adapter and ploy-N, and low-quality reads were removed to obtain clean data (clean reads). Reference genome and gene model annotation files were obtained from the genome website. The index of the reference genome was built using Hisat2 v2.0.5 and paired-end clean reads were aligned to the reference genome. FeatureCounts v1.5.0-p3 was then used to count the read numbers mapped to each gene. The TPM value of each gene was calculated on the basis of the gene length and reads count.

Immunohistochemistry

Patients (n = 45) undergoing the surgical removal of GBM in Xiangya Hospital, Central South University, were the sources of tissues. Tissues were then fixed by formalin and embedded in paraffin for subsequent obtaining of slices (4 μm). The sections were dewaxed in xylene 3 times for 10 minutes each. The dewaxed sections were then dehydrated in 100% ethanol, 95% ethanol, 90% ethanol, 85% ethanol, 80% ethanol, and 70% ethanol for 5 minutes each time. The sections were then placed in a citric acid buffer and boiled for antigen retrieval, and 3% H2O2 was adopted as a blocker of HRP activity from endogenesis. 5% BSA was used for section blocking. Rabbit polyclonal anti-CD163 (Rabbit, 1:1000, 16646-1-AP, Proteintech, China) and anti-MCAM antibody (Rabbit, 1:200, 17564-1-AP, Proteintech, China) were the primary antibodies, while the HRP-labeled Goat Anti-Rabbit IgG (PV9000, ZSGB-BIO, China) was the secondary antibody. Sections with the primary antibody were incubated at 4 degrees Celsius overnight. Substrate (3, 3′-diaminobenzidine, DAB) with a mixture of solutions 1 and 2 at 1 drop/1 mL was used for checking the signal. Hematoxylin was used in slices-staining. An optical microscope was finally used for observation after staining. As for the intensity score, four intensity levels, negative, weak, moderate, and strong, were assigned with 0, 1, 2, and 3, respectively. As for extent score (proportion of stained cells), 10%, 10–25%, 25–50%, 50–75%, >75% were assigned with 0, 1, 2, 3, and 4, respectively. The H-score was calculated as extent∗intensity and had the range of 0–12.

Multiplex immunofluorescence staining

The tissue microarray was purchased from the Wuhan Bioyeargene Biotechnology company (NCT805, Wuhan, China) and the ethics was approved. Paraffin sections of the tissue microarray of 80 cores from 35 GBM samples and 5 normal samples (duplicated cores per case) were deparaffinized. After antigen retrieval, sections were blocked with 3% H2O2 and 2% BSA. Different primary antibodies, NG2 (Rabbit, 1:3000, 55027-1-AP, Proteintech, China), MCAM (Rabbit, 1:100, 17564-1-AP, Proteintech, China), CD68 (Rabbit, 1:3000, GB113150, Servicebio, China), CD163 (Rabbit, 1:3000, 16646-1-AP, Proteintech, China), were sequentially applied, followed by horseradish peroxidase-conjugated secondary antibody incubation (1:500, GB23301, GB23303, Servicebio, China) and tyramide signal amplification (TSA) (FITC-TSA, CY3-TSA, 594-TSA and 647-TSA (Servicebio, China)). After labeling the human antigen, the nucleus was coated with 4′,6-Diamidino2-phenylindole dihydrochloride (DAPI), and an antifade mounting medium was applied. Multispectral images were obtained by scanning stained slides using the Pannoramic Scanner (3D HISTECH, Hungary). DAPI emits blue light through UV excitation wavelength 330–380 nm and emission wavelength 420 nm in fluorescence analysis. CY3 emits red light through excitation wavelength 510–560 nm and emission wavelength of 590 nm. 594 emits fuchsia light through excitation wavelength 594 nm and emission wavelength 615 nm. FITC emits green light through excitation wavelength 465-495nm and emission wavelength 515–555nm. 647 emits pink light through excitation wavelength 608–648nm and emission wavelength 672–712nm. Multispectral images were analyzed using case viewer (C.V 2.3, C.V 2.0) and pan luminosity viewer (P.v 1.15.3) image analysis software, and positive cells were quantified at the single-cell level. The spatial analysis of the stained cells was performed using the StrataFAXS software.

RT-qPCR assay

The primers of β-actin (F ACCCTGAAGTACCCCATCGAG; R AGCACAGCCTGGATAGCAAC) and MCAM (F CACCGTCCCTGTTTTCTACCC; R TCCCCTTCCTTCAGCATTCCCA) were designed using the primer 5.0. Total RNAs were extracted and reversely transcribed into cDNA by HiScript Q RT SuperMix for RT-qPCR. The expression levels of β-actin and MCAM were quantified using 2-ΔΔCT.

Western blotting assay

The western blotting assay assessed the expression level of MCAM and β-actin. Anti-MCAM (Rabbit, 1:2000, 17564-1-AP, Proteintech, China) and anti-β-actin (Mouse, 1:5000, 66009-1-Ig, Proteintech, China) were used as the primary antibody. HRP goat anti-mouse IgG (Mouse, 1:5000, SA00001-1, Proteintech, China) and HRP goat anti-rabbit IgG (Rabbit, 1:6000, SA00001-2, Proteintech, China) were used as the secondary antibody. ECL development was used for visualization.

Coculture of HMC3 and human microvascular pericyte for transwell assay

Human microvascular pericyte (HUM-iCell-n011) and human microglia cell line (HMC3) were purchased from iCell (http://www.icellbioscience.com). At the density of 2 × 105/mL, human microvascular pericytes were added to the 6-well plate and transfected with si-NC and si-MCAM. The treated human microvascular pericytes were then digested and resuspended using 10% DMEM and were added to the lower chamber. After human microvascular pericytes were transfected for 48h, HMC3 cells were also digested and resuspended at 1 × 106/mL. 100ul HMC3 cells were added to the upper chamber. At the ratio of 1:1, human microvascular pericytes and HMC3 cells were cocultured for 48h. After being washed with phosphate buffer saline (PBS) twice, the upper chamber was fixed using acetone and methyl alcohol at the ratio of 1:1 for 20 min. The upper chamber was then stained with 0.5% crystal violet for 5 min for photographing.

Coculture of HMC3 and human microvascular pericyte for multiplex immunofluorescence staining

The transfected human microvascular pericytes were added to the upper chamber (1 × 105 each well), while HMC3 cells were added to the lower chamber and placed on the cell slides (1 × 105 each fall). After coculturing for 24h, the cell slides were washed twice with PBS and fixed with 4% paraformaldehyde for 30 min. 0.3% triton was added to the cell slides for 30 min under 37°C. The cell slides were blocked with 5% BSA for 60 min. After being washed with PBS three times, the cell slides were added with primary antibodies CD68 (1: 50, Mouse, 14-0688-82, ThermoFisher, America), CD11c (1: 50, Rabbit, 17342-1-AP, Proteintech, China), and CD68 (1: 50, Mouse, 14-0688-82, ThermoFisher, America), CD163 (1: 50, Rabbit, 16646-1-AP, Proteintech, China), respectively. Anti-mouse (1:200, A11008, ThermoFisher, America) and anti-rabbit (1:200, A11005, ThermoFisher, America) IgG secondary antibodies were then added. The cell nucleus was stained with DAPI for 10 min. The cell slides were finally blocked with buffer glycerin and observed with a microscope.

Quantification and statistical analysis

Survivorship curves were generated using the R package survminer. The univariate and multivariate Cox regression analyses determined the clinical significance of prognostic factors. Correlation coefficients were calculated by Pearson correlation to analyze. The receiver operating characteristic (ROC) analysis was visualized using the R package pROC. The R package maftools were used for depicting the mutation landscape of TCGA via OncoPrint(Gu et al., 2016). For normally distributed variables, significant quantitative differences between and among groups were determined by a two-tailed t-test or one-way ANOVA, respectively. For nonnormally distributed variables, significant quantitative differences between and among groups were determined by a Wilcoxon test or a Kruskal–Wallis test, respectively. All statistical analyses were performed on R project 3.6.3 and GraphPad Prism 8.4.3. p < 0.05 was considered statistically significant.

REAGENT or RESOURCE	SOURCE	IDENTIFIER
Antibodies
Rabbit polyclonal anti-CD163	Proteintech	16646-1-AP
Rabbit polyclonal anti-MCAM antibody	Proteintech	17564-1-AP
HRP-labeled Goat Anti-Rabbit IgG	ZSGB-BIO	PV9000
Rabbit polyclonal anti- NG2 antibody	Proteintech	55027-1-AP
Rabbit polyclonal anti- CD68 antibody	Servicebio	GB113150
Rabbit polyclonal anti- CD163 antibody	Proteintech	16646-1-AP
horseradish peroxidase-conjugated secondary antibody	Servicebio	GB23301
Mouse polyclonal anti-β-actin antibody	Proteintech	66009-1-Ig
HRP goat anti-mouse IgG	Proteintech	SA00001-1
HRP goat anti-rabbit IgG	Proteintech	SA00001-2
Mouse polyclonal anti- CD68 antibody	ThermoFisher	14-0688-82
Rabbit polyclonal anti- CD11c antibody	Proteintech	17342-1-AP
Anti-mouse IgG secondary antibody	ThermoFisher	A11008
Anti-rabbit IgG secondary antibody	ThermoFisher	A11005
Deposited data
TCGAGBM	UCSC Xena platform	https://xena.ucsc.edu/
CGGA311	The Chinese Glioma Genome Atlas	http://www.cgga.org.cn
CGGA668	The Chinese Glioma Genome Atlas	http://www.cgga.org.cn
GSE4412	Gene Expression Omnibus	https://www.ncbi.nlm.nih.gov/geo
GSE83300	Gene Expression Omnibus	https://www.ncbi.nlm.nih.gov/geo
GSE108474	Gene Expression Omnibus	https://www.ncbi.nlm.nih.gov/geo
Experimental models: Cell lines
HMC3	iCell	http://www.icellbioscience.com
HUM-iCell-n011	iCell	http://www.icellbioscience.com
Oligonucleotides
β-actin F     ACCCTGAAGTACCCCATCGAG	This paper	N/A
β-actin R    AGCACAGCCTGGATAGCAAC	This paper	N/A
MCAM F     CACCGTCCCTGTTTTCTACCC	This paper	N/A
MCAM R    TCCCCTTCCTTCAGCATTCCCA	This paper	N/A
Software and algorithms
R (version 3.6.3)	Downloaded from the R Project for Statistical Computing	https://www.r-project.org/
GraphPad Prism (version 8.4.3)	Downloaded from the GraphPad Prism for Statistical Computing	https://www.graphpad-prism.cn