Spatial-reprogramming derived GPNMB+ macrophages interact with COL6A3+ fibroblasts to enhance vascular fibrosis in glioblastoma

Abstract

Background: Neoadjuvant therapy plays an important role in the treatment of glioblastoma (GBM), but a considerable proportion of patients remain unresponsive to the combination of immune checkpoint blockade (ICB) and antiangiogenic therapy. Understanding the mechanisms underlying resistance to this treatment and developing novel therapeutic strategies are crucial.

Methods: We integrate extensive single-cell and spatial transcriptomic data to dissect the cellular composition and spatial organization of the GBM tumor microenvironment and validate our findings through experiments such as multiplex immunohistochemistry and atomic force microscopy. We applied 101 machine learning algorithms to evaluate the prognostic and immunological value of COL6A3⁺ tumor-associated fibroblasts (TAFs) and GPNMB⁺ monocyte-derived macrophages (MDMs) in multiple GBM cohorts and immunotherapy cohorts.

Results: We constructed a stromal cell atlas in GBM and identified a distinct subset of COL6A3⁺ TAFs with functional characteristics of matrix fibroblasts. We found that COL6A3⁺ TAFs are significantly enriched in non-responders to neoadjuvant combination therapy. These fibroblasts drive the spatial-reprogramming of anti-tumorigenic MDMs into a pro-tumorigenic phenotype. In turn, these reprogrammed immunosuppressive GPNMB⁺ MDMs promote vascular fibrosis mediated by COL6A3⁺ TAFs through the GPNMB-ITGB5 interaction.

Conclusions: Our findings highlight the critical role of COL6A3⁺ TAFs in regulating MDM function and spatial distribution, as well as their contribution to fibrotic tumor vasculature formation. Additionally, we propose targeting COL6A3⁺ TAFs with cilengitide as a potential therapeutic strategy.

Keywords: Anti-angiogenic agents; Cancer-associated fibroblasts; Immune checkpoint blockade; Single-cell RNA sequencing; Tumor-associated fibroblasts; Tumor-associated macrophages.

Fig. 1

A significant reduction of COL6A3⁺ TAFs in neoadjuvant responders. A Uniform manifold approximation and projection (UMAP) plot showing different subtypes of stromal cells (n = 13,143 cells). B Dot plot showing the top 3 highly expressed genes for each subtype of stromal cells. C UMAP plot showing the relative abundance of subtypes based on the sample origin. D Pie charts displaying the frequencies of cell subtypes in different patient groups. E Gene network diagram showing functional analysis of COL6A3⁺ TAF characterized genes. F Representative immunohistochemical staining of pathological slides and box plots showing increasing COL6A3 expression and average optical density (AOD) in responder (n = 3) and nonresponder (n = 3) samples. G Representative immunohistochemical images and AOD quantification of COL6A3 expression in peritumoral regions (normal brain, n = 8), pilocytic astrocytoma (WHO grade I, n = 8), astrocytoma (WHO grade II/III, n = 8), oligodendroglioma (WHO grade II/III, n = 8), and glioblastoma (WHO grade IV, n = 8)

Fig. 2

Assessment of the degree of hypoxia and T cell activation function in two MDM subpopulations with opposite prognostic characteristics. A UMAP plot showing different subtypes of MDMs (n = 40,184 cells). B Dot plot showing the top 3 highly expressed genes for each subtype of MDMs. C Forest plot showing the results of multivariate cox regression analysis of overall survival in patients in the TCGA-GBM/LGG cohort (n = 494), with error bars showing 95% confidence intervals. D Kaplan–Meier plot showing that a lower proportion of ICAM1⁺MDMs and a higher proportion of GPNMB⁺MDMs based on GSVA scores grouping is associated with shorter overall survival. E Raincloud plots displaying the MPI scores of different MDM subgroups. F Experimental design of isolating patient-derived ICAM1⁺ MDMs and GPNMB⁺ MDMs. G Representative FCM plots showing the percentages of patient-derived ICAM1⁺ MDM and GPNMB⁺ MDM. H Representative mIHC images of DAPI, CD68, GPNMB, ICAM1, and HIF-1α. Scale bar = 50 μm. The left dashed box shows ICAM1⁺ MDM staining, and the right box shows GPNMB⁺ MDM staining. I Ridge plot showing the expression level of HIF-1α. J Experimental design of the hPBMCs-derived T cells and primary MDM co-culture system. K Representative FCM scatterplots showing the T cell cytotoxic activity level. L Representative FCM scatterplots showing the T cell exhaustion level. M Representative FCM scatterplots showing the T cell proliferation level

Fig. 3

Spatial transcriptomics reveals the co-localization of cell subpopulations in different pathological locations. A Box plots showing the estimated proportions of each cell cluster in different histological regions of GBM, obtained through deconvolution of bulk RNA-seq data from the Ivy-GBM cohort (n = 270). B H&E staining (left) and surface plots (right) showing the distribution of different cell subtypes in newly diagnosed GBM patients and neoadjuvant non-responders. C Representative mIHC images of GLUT1, GPNMB, CD68, and ICAM1 in the pseudopalisading necrosis region and α-SMA, CD68, ICAM1, and GPNMB in the peri-vascular region. Box plots (right) displaying the proportions of ICAM1⁺ MDMs and GPNMB⁺ MDMs among all CD68⁺ MDMs in randomly selected high-power fields (HPF, n = 6)

Fig. 4

Pseudotime trajectory analysis reveals the differentiation pathways of monocytes/MDMs. A UMAP plot showing the monocyte/MDM trajectory derived from Slingshot. B Trajectory plots showing the monocyte/MDM trajectory derived from Monocle2, colored based on cell subpopulations (left), cell state (top right), and pseudotime progression (bottom right). C Pseudotime projection of transcriptional changes in key genes during differentiation along the trajectory. D Heatmap (left) showing the differential gene expression of cells along the differentiation trajectory. On the right, the top five significantly enriched hallmark biological terms at each differentiation stage are displayed. E Chord diagram showing the interactions between seven important cell subtypes in the neoadjuvant therapy responders (left) and non-responders (right). F Heatmap (left) showing the top-ranked ligands that regulate MDMs by COL6A3⁺ TAFs according to Nichenet. Dot plot (middle) displaying the expression levels of these ligands in different cell subpopulations. Heatmap (right) showing the downstream genes that these ligands activate

Fig. 5

COL6A3⁺ TAFs mediate the spatial-reprogramming of ICAM1⁺ MDMs into GPNMB⁺ MDMs through TGFβ3 and CSF1. A Experimental design of the process to obtain TAF-conditioned medium (TAF-CM) and tumor cell-conditioned medium (TCM). B Immunofluorescence showing the expression of PDGFRA and COL6A3 in tumor cells and COL6A3⁺ TAFs. C Box plots showing the concentrations of TGF-β3 and CSF1 in TAF-CM and TCM. D Surface plot displaying the spatial analysis trajectory. Yellow dots represent spatial locations included in the gene gradient analysis. Continuous line plot illustrating the changes in ICAM1⁺ MDM and GPNMB⁺ MDM signature genes from the MVP region through the CT region to the PAN region. E Experimental design of sorting of CD44⁺ MES-like tumor cells (above). Representative images showing the expression levels of CD44, CD24, PDGFRA, and EGFR in the sorted cells (below). F Line charts showing the temporal changes in the expression of ICAM1 and GPNMB measured via flow cytometry. G Line charts showing the temporal expression changes in the signature genes of ICAM1⁺MDM and GPNMB⁺MDM measured by qRT-PCR. H Box plot displaying the expression levels of ICAM1⁺ MDM and GPNMB⁺ MDM signature genes under different treatment conditions measured by qRT‒PCR. I Representative FCM scatterplots showing the expression levels of ICAM1 and GPNMB under different treatments. J Experimental design of establishing patient-derived GBO cultures. K Representative mIHC images of COL6A3, ICAM1, and GPNMB (left) and the box plot quantifying the fluorescence intensity (right). The fluorescence intensity of GPNMB was divided by the fluorescence intensity of ICAM1 to assess the extent of MDM reprogramming

Fig. 6

COL6A3⁺ TAFs enhance vascular fibrosis and lead to T cell exclusion through GPNMB/ITGB5 interplay. A Pearson correlation analysis showing the associations between GPNMB and fibrosis-related genes in the CGGA-GBM/LGG cohort (n = 322). B Scatter plot showing the differentially expressed genes between COL6A3⁺ TAFs and other stromal cells. C Dot plot showing the expression levels of ITGB5 across different stromal cell subpopulations. D Venn diagram showing ITGB5 as a key interacting gene of GPNMB. The data sources include characteristic genes of COL6A3⁺ TAFs from our scRNA-seq dataset, GPNMB-binding proteins from the STRING database, and the ECM receptor interaction gene set from MSigDB. E Multiplex immunohistochemistry revealing the co-localization of GPNMB and ITGB5 in non-responder. G The line graph quantitatively depicting the fluorescence intensity of GPNMB (yellow) and ITGB5 (red) within the analyzed region. H Immunofluorescence staining showing the expression levels of PDGFRA and COL6A3 in different treatment groups (left), and the magnified images were quantified (right). The relative fluorescence intensity was calculated by dividing the total fluorescence intensity of each group by the number of nuclei and normalizing it to that of the control group. The treatment group with only sGPNMB added was the comparison group for statistical significance. I Representative atomic force microscopy images showing the cell morphology (above) and three-dimensional distribution of the cell height (below). J Box plots quantifying the Young’s modulus and roughness of the cell surface. K Representative mIHC images of CD31, COL6A3, and CD3 in newly diagnosed GBM patient (left), neoadjuvant combination therapy non-responder (middle), and neoadjuvant combination therapy responder (right)

Fig. 7

High infiltration of COL6A3⁺ TAFs and GPNMB⁺ MDMs is correlated with immunotherapy resistance. A Schematic diagram showing the computational framework for establishing the cell subtype infiltration gene set. B Through a tenfold cross-validation framework, a total of 101 combinations of machine learning algorithms were used for the cell subpopulation infiltration gene set. TCGA-GBM/LGG (n = 691) was used as the training dataset; CGGA-GBM/LGG (n = 693/325), GLASS (n = 371), and GSE108474 (n = 487) were used as validation datasets. The C-index for each model was calculated based on its application across all datasets. C Kaplan–Meier survival curve shows the overall survival of patients in the training cohort, grouped according to the risk scores calculated by the model. D Kaplan–Meier survival curve shows the overall survival of patients in the validation cohort, grouped according to the risk scores calculated by the model. E Kaplan–Meier survival curves demonstrating the associations between cell infiltration, which is based on GSVA scores, and overall survival of patients across different immunotherapy cohorts (left). Box plots combined with violin plots showing the relationship between cell subtype infiltration scores and treatment response in patients with different therapeutic outcomes (middle), as well as the correlation between infiltration scores and TIDE scores (right)

Fig. 8

TME landscape of non-responding GBM patients to combination therapy. COL6A3⁺ TAFs and MDMs form a positive feedback loop. Peri-vascular COL6A3⁺ TAFs induce functional and spatial changes in ICAM1⁺ MDMs by secreting TGFβ3 and CSF1. GPNMB⁺ MDMs, which are located in pseudo-palisading necrotic regions, promote COL6A3⁺ TAF-mediated collagen deposition through GPNMB-ITGB5 ligand-receptor interaction, leading to vascular fibrosis and T cell exclusion. The abundance of COL6A3⁺ TAFs can determine the efficacy of combination therapy