Mesenchymal Stem Cells and Fibroblasts Contribute to Microvascular Proliferation in Glioblastoma and are Correlated with Immunosuppression and Poor Outcome Available to Purchase

Microvascular proliferation (MVP) is an aggressive histopathologic phenotype that is a defining hallmark of glioblastoma (GBM), the most common and deadly primary adult brain tumor. MVP can independently upgrade lower-grade gliomas (LGG) into World Health Organization (WHO) grade 4 gliomas and is detectable in a variety of solid tumors, including melanoma and breast, endometrial, and prostate cancers, in which it also correlates with poor prognosis (1–3). Even though MVP is detectable by histopathology and has substantial clinical relevance, little is known about the biology underlying the regulation of MVP or its contribution to the tumor microenvironment, especially at the single-cell level.

Classically, MVP was thought to occur because of aberrant endothelial cell hyperplasia (4, 5). However, multiple groups have reported that glioma stem cells may also contribute to MVP by directly differentiating into endothelial cells and pericytes, a phenomenon known as vascular mimicry (6–9). Based on our previous work (10–15) and the work of others (16–23), we suspected the involvement of another cell type in MVP, mesenchymal stem cells (MSC). Beginning in the 19th century, MSCs were characterized as multipotent stem cells from extracranial locations like the bone marrow (BM) that could differentiate into cells of mesenchymal heritage such as osteocytes and chondrocytes (24–27). We have not only found that MSCs reside in the perivascular niche in mouse brains (13, 14) but also demonstrated that intravascularly administered human BM-derived MSCs are capable of homing to glioma stem cell xenografts (28). Furthermore, MSCs may be the progenitors of fibroblasts under noncancerous conditions (13, 16–22), presenting the possibility that glioma cancer-associated fibroblasts (CAF; ref. 17) may originate from MSCs. Thus, we hypothesized that MVP is composed of a rich network of cell types stemming from MSCs and CAFs, offering therapeutic targets for GBM. Although single-cell technology has refined our understanding of certain cellular subsets in the GBM microenvironment, the substrates of MVP have yet to be characterized at this resolution.

Although MVP is classified as an abnormality of vasculature, it is unknown whether this phenomenon contributes directly to the angiogenic needs of the tumor. Previous studies have failed to find substantive evidence linking MVP to microvessel density, the standard for assessing neovasculature in tumors (1). In fact, the presence or absence of MVP was a superior predictor of poor survival compared with microvessel density. This suggested that MVP has alternative contributions to benefit tumor survival. Therefore, we also hypothesized that the cellular substrates of MVP, namely MSCs and CAFs, may also contribute to an immunosuppressive tumor microenvironment.

To begin to deconvolute the biology of MVP in the context of glioma and enrich for MSCs, we performed single-cell RNA sequencing (scRNA-seq) on paired CD45⁺CD105^± and CD45⁻CD105⁺ cells from 16 primary glioma patient samples with and without MVP. We chose endoglin (CD105) as an enriching marker for MSCs because it is one of a few canonical markers of human BM MSCs (alongside CD73, CD90, and PDGFRα; refs. 16, 29), we have previously found CD105⁺ cells in the mouse glioma perivascular niche (10), and because CD105 is associated with abnormal angiogenesis (30). Furthermore, CD105 is expressed by endothelial cells (31) and other cell types we suspected are associated with the perivascular niche such as fibroblasts (32) and myeloid cells (31). In line with our hypothesis, we found an appreciable enrichment of MSCs and CAFs in the CD45⁻CD105⁺ cell fraction, with specific MSC subsets predominantly from MVP samples. Multiplex immunofluorescence (IF) confirmed the presence of MSCs in the glioma perivascular niche. Furthermore, we identified platelet-derived growth factor receptor β (PDGFRB) as a putative driver gene of MVP. Signal transduction network analyses between perivascular stromal cells (PVSC) and immune cell components along with spatial transcriptomics implicated increased TGF-β pathway activity and a higher frequency of immunosuppressive immune cell subsets under MVP conditions.

Together, this study (i) refines the paradigm of MVP by highlighting the role of PDGFRβ⁺ MSCs that may differentiate into other PDGFRβ⁺ PVSCs, including cancer-associated fibroblasts (CAFs), pericytes, fibromyocytes, and smooth muscle cells, (ii) implicates a potential role whereby MVP promotes an immunosuppressive microenvironment in GBM, and (iii) provides a single-cell atlas in which both human glioma PVSCs and immune cells were enriched for in the context of MVP (17, 33, 34).

This study aimed to elucidate which cells were involved in MVP at the single-cell level through analysis of 16 fresh human glioma specimens using scRNA-seq complemented by multiplex IF and digital spatial profiling (DSP). FACS was used to enrich for CD45⁻CD105⁺ and CD45⁺ cellular populations. Bioinformatic analyses included top differential gene expression, RNA velocity, and de novo signal transduction network analyses. We then confirmed the prognostic relevance of MVP-associated PVSC gene signatures across multiple solid tumors using publicly available databases.

Sixteen patient tumor samples were obtained under MD Anderson Cancer Center Institutional Review Board–approved protocol 2012 to 0441 between July 2021 to July 2022. All patients signed informed consents for participation in protocol 2012 to 0441 before surgery or sample collection. Samples were collected consecutively without randomization. The clinical demographics, including age and sex, of individual patients are shown in Table 1. There was no attrition of subjects.

To verify whether the samples had MVP, we not only reviewed the clinical report generated by the board-certified neuropathologist on call but also the scientific specimen submitted for assay with Dr. Gregory Fuller (expert neuropathologist). The diagnosis was based on the 2021 WHO Classification of Tumors of the central nervous system (CNS; 35). Four GBMs and one IDH-mutant WHO grade 4 astrocytoma had MVP (MDA3671, MDA3753, MDA3841, MDA3807, and MDA3944, respectively), whereas three GBMs and one IDH-mutant WHO grade 4 astrocytoma did not (MDA3707, MDA3692, MDA3822, and MDA3868, respectively). None of the WHO grade 2 and 3 gliomas had MVP. WHO grade 2 and 3 gliomas were included in the analysis because all have the potential to dedifferentiate and develop MVP. All tumors were chemoradiation-naïve, with the majority exposed to no/low doses of corticosteroids to promote an accurate reflection of the biological, especially immunologic, processes investigated in this dataset.

Fresh tumor tissue was dissociated using the gentleMACS Dissociator (Miltenyi Biotec) according to the manufacturer’s instructions. Percoll density gradient centrifugation at 512 g for 20 minutes at 18°C was then used to deplete the myelin layer, and the resulting single-cell suspension was counted in Cellometer vision CBA (Nexcelom Biosciences). Cells were resuspended in freezing media containing 90% AB serum (purchased from Sigma, Cat. # H4522) and 10% DMSO and stored in liquid nitrogen until analysis.

For scRNA-seq analysis, cells were thawed, washed, and sorted for viable CD45⁺and CD105⁺ cells using the BD FACSAria III SORP sorter (BD Biosciences) and a surface staining cocktail of fluorescently conjugated antibodies comprised of CD45 Alexa Fluor 532 (clone HI30, eBioscience, 58-0459-42; RRID: AB_11218673), live/dead discrimination viability dye Pacific Orange (Invitrogen, L34968), and CD105 (Invitrogen, SN6, allophycocyanin). Cells from each sample were counted to load 15,000 cells per sample on the 10× chromium chip (Chromium platform, 10× Genomics) with a target of 10,000 cells per sample for the downstream analysis. As experiments used precious patient samples, we did not perform technical replicates. Single-cell mRNA libraries were built using Chromium Next GEM Single Cell 3’ Library Construction V3 Kit. Libraries were sequenced using NovaSeq 6000 Reagent Kit (flow cell type – S2-100, run format – read 1: 28 bp, read 2: 91 bp paired-end sequencing, and 8 i7 index).

The CellRanger v6.0.0 software (10× Genomics) was used to process the sequencing reads. The “cellranger count” pipeline was used to align the reads to the human transcriptome GRCh38-2020-A genome and compute the count matrix. The generated cell-by-gene unique molecular identifier (UMI) count matrix was analyzed using the Seurat R package v.4.1.0 (RRID: SCR_007322). We kept only the cells expressing at least 200 genes. The genes were also filtered by the maximum of 8,000 expressed genes and of 15% mitochondrial genes. The UMI counts were then normalized for each cell by the total expression and log-transformed. We used Seurat’s default method to identify highly variable genes and scale data for regressing out variation from UMIs and mitochondrial genes. Cell-cycle regression was also performed to regress out the differences between the G2M and S phase scores. The scaled data with variable genes were used to perform principal component analysis. The top 30 principal components were chosen for further analysis, including clustering to identify cell populations. Uniform Manifold Approximation and Projection values were calculated in the Seurat R package using the top 30 principal components and min_dist = 0.8. Harmony was used to perform batch-effect correction and to integrate samples from individual patients. Seurat was used to identify cluster-specific marker genes and visualization with dot, violin, and feature plots.

To annotate the cell type for each cluster, we identified marker genes using the Wilcoxon rank-sum test by comparing one cluster with the others. Marker genes were upregulated significantly by at least a log₂ fold change value of 0.25 in each cluster compared with all other cells. Next, we checked canonical markers in each cluster to determine the cell type. Finally, we performed Gene Ontology (GO) enrichment analysis. Cell-type annotation revealed that cluster 9 was enriched with long noncoding RNAs MALAT1 and NEAT1. In addition, cluster 9 was negative for canonical markers of immune cells and the nFeature density was low (mean = 707.8 units). Contributions to this cluster came from 15/16 patients. Thus, this cluster consisted of stressed cells and was excluded from further analysis.

To estimate RNA velocity, we first generated loom files containing the spliced and unspliced RNA from the original CellRanger output BAM files and the GRCh38 (version refdata-gex-GRCh38-2020-A, 10X Genomics) human genome reference using velocyto (36). Dynamical modeling, latent time, and phase analyses and visualizations were all performed using scVelo (37).

To identify enriched molecular pathways based on differentially expressed genes (DEG), an overrepresentation analysis was performed on DEGs from each cluster using clusterProfiler V4.2.2 (RRID: SCR_016884). Gene sets from GO biological processes were used.

Signal transduction networks were generated as previously described (38). In brief, cross-talk scores were comprised of preferential expression measures, reflecting gene expression across all cell types, and non–self-talk scores, reflecting the Shannon entropy for a ligand and receptor pair. The product of the min–max normalized preferential expression measures and non–self-talk score values yields a cross-talk score for each ligand–receptor pair and for each cell-type pair. The signaling network between two cell types was generated by treating it as if it were a prize-collecting Steiner forest problem. The “edge cost” represents the probability of a functional interaction between two genes and the “node prize” represents cell-specific gene activity.

The network, nodes, and edges were imported into Cytoscape (RRID: SCR_003032). Nodes width and height were locked. Different shapes were selected to identify cell type. Size of the node was selected to represent gene expression. Minimum values were set at 0 (size of 20.0) and maximum values at 1.5 (size of 80.0). Node fill color represented the node prize, which was defined on both the cell type gene expression specificity and network distance to the ligand/receptor gene. Minimum values were set at 0 and 0.35. Edge thickness was inversely proportional to edge cost, which is the probability of a functional interaction between two genes. Minimums were set at 0 (edge width of 40) and maximums at 1.0 (edge width of 1). Cross-talk edges were denoted with magenta, indicating a ligand–receptor interaction. The compound spring embedder layout was applied.

We used CellChat to analyze cell–cell communication networks. Processed scRNA-seq data, normalized and log-transformed, served as input. Gene expression data were mapped to the CellChat database, which annotates signaling pathways based on ligand–receptor interactions. CellChat inferred intercellular communication by identifying significant interactions and calculating communication probabilities using a built-in statistical framework. Identified signaling networks were visualized through chord diagrams to reveal dominant pathways and their contributions. Source and target cells were defined based on the biological question and as detailed in CellChat vignettes.

We used hallmark gene sets from the Molecular Signatures Database to analyze biological processes and pathways. Gene expression matrices were prepared by normalizing raw counts and filtering low-expressed genes. AUCell, an R package, was used to quantify the activity of hallmark gene sets within individual cells. Cell clusters were downsampled for computational efficiency in AUCell. Gene set activity scores were computed by ranking each cell’s expression data and evaluating the enrichment of hallmark genes using the AUC. These scores enabled comparative analyses of pathway activity across cells. Hallmark gene sets were cross referenced with CellChat signaling pathways to infer the ligands, receptors, and cell types involved in biological processes.

Approximately 30% of the volume of tissue specimens used to generate scRNA-seq data was first split into smaller pieces to generate formalin-fixed paraffin-embedded tissue sections which were stained with hematoxylin and eosin (H&E) following standard methods and mounted with Cytoseal 60 (Thermo Fisher Scientific). H&E-stained sections were imaged using an Olympus BX53 microscope with Olympus CellSens software. For IF, serial sections were deparaffinized in xylene and rehydrated in graded ethanols. Slides were submerged in citrate-based antigen retrieval buffer (pH 6.0), heated for 25 minutes in a vegetable steamer, and then cooled to room temperature. Tissues were blocked with 5% BSA in PBS for 1 hour at room temperature and then incubated overnight at 4°C with anti-GFP (Abcam, ab13970, diluted 1:1,000) and anti-CD31 (R&D Systems, AF3628, diluted 1:200) primary antibodies diluted in blocking buffer. Slides were washed and incubated for one hour with Alexa Fluor 488 donkey anti-chicken and Alexa Fluor 594 donkey anti-goat secondary antibodies (Invitrogen/Thermo Fisher Scientific), each diluted 1:500 in blocking buffer. Slides were washed and mounted with Vectashield antifade mounting medium with 4',6-diamidino-2-phenylindole (DAPI, Vector Laboratories), and coverslips were sealed with clear, acrylic nail polish. IF was imaged using an Olympus FV1000 confocal microscope with Olympus CellSens imaging software.

Moribund mice were euthanized by CO₂ inhalation and tissues were fixed by cardiac perfusion of 4% paraformaldehyde followed by 10% neutral buffered formalin. Fixed tissues were processed and paraffin-embedded by standard methods. Formalin-fixed paraffin-embedded tissues were sectioned and stained with modified Harris hematoxylin (Epredia) and eosin-Y (Thermo Fisher Scientific) following standard methods and mounted using Cytoseal XYL (Thermo Fisher Scientific). H&E-stained sections were imaged using an Olympus BX53 microscope with Olympus CellSens software. For IF, serial sections were deparaffinized in xylene and rehydrated in graded ethanols. Slides were submerged in IHC Antigen Retrieval Solution (low pH; Invitrogen/Thermo Fisher Scientific), heated for 30 minutes in a vegetable steamer and then cooled to room temperature. Tissues were blocked with 10% normal donkey serum in PBS for 1 hour at room temperature and then incubated overnight at 4°C with anti-GFP (Abcam, ab13970, diluted 1:1,000) and anti-CD31 (R&D Systems, AF3628, diluted 1:200) primary antibodies diluted in blocking buffer. Slides were washed and incubated for one hour with Alexa Fluor 488 donkey anti-chicken and Alexa Fluor 594 donkey anti-goat secondary antibodies (Invitrogen/Thermo Fisher Scientific), each diluted 1:500 in blocking buffer. Slides were washed and incubated for 5 minutes with Hoechst 33,342 nuclear stain and mounted with Prolong Diamond antifade mounting medium. IF was imaged using a Lunaphore COMET microscope with HORIZON viewer software (Lunaphore Technologies).

For multiplex IF staining, we followed the Opal protocol staining method (39). Cells were stained for multiplex IF using an RX-BOND (Leica) automated multiplex slide stainer. Tissue sections were stained with antibodies for the following markers: CD31 (Abcam, 647-JC/70A, 1:100 dilution; RRID: AB_2890260), CD105 (Cell Signaling Technology, 3A9, 1:25 dilution; RRID: 3A9), Olig-2 (Abcam, EPR2673, 1:100 dilution), Iba1/AIF1 (Sigma-Aldrich, 20A12.1, 1:100 dilution), CD73 (Cell Signaling Technology, D7F9A, 1:50 dilution), CD90 (Abcam, EPR133, 1:50 dilution), CD34 (Abcam, EP373Y, 1:200 dilution), TAGLN (Thermo Fisher Scientific, GT336, 1:50 dilution), DCN (Proteintech, Poly 144667-1-AP, 1:25 dilution), and PDGFRβ (Cell Signaling Technology, 28E1, 1:25 dilution), with subsequent visualization using Akoya Opal fluorophores (650, 480, 570, 540, 690, 620, and 520, respectively), and nuclei were visualized with DAPI (1:2000 dilution). All sections were coverslipped using Vectashield Hardset895 mounting medium. Slides were scanned using a Vectra/Polaris slide scanner (PerkinElmer). For multispectral analysis, each of the individually stained sections was used to establish the spectral library of the fluorophores using Inform (Akoya). Images were acquired at ×20 magnification and unmixed using the acquired spectral libraries using Inform software (Akoya).

For multispectral IF analysis, each of the individually stained sections was used to establish the spectral library of the fluorophores using Inform software (Akoya); this library was then used to unmix all multiplexed images acquired at 20X magnification. Furthermore, quantitative analysis was performed using HALO software (version 3.6.4). A pathologist (S. Jindal) annotated margins around normal blood vessels (n = 60), and in cases with MVP, proliferative blood vessels were annotated (n = 60). Fluorescence multiplex analysis algorithm, HighPlex FL v4.04 (Indica Labs), was used to determine the number of markers coexpressing in cells of associated blood vessels. The data were plotted using Prism v.8.0 (GraphPad Software; RRID: SCR_002798). The infiltration analysis algorithm from the HALO Spatial Analysis Module was used to determine the distance of marker-expressing cells from blood vessels. Statistical significance for nonparametric data was calculated using a Wilcoxon rank-sum or Mann–Whitney U test, and P value <0.05 was considered statistically significant. Statistical significance for parametric data was calculated using an unpaired t-test, and P value <0.05 was considered statistically significant.

DSP (NanoString) experiments were performed according to the manufacturer's protocol and as previously described (40). Briefly, slides were stained and imaged to visualize morphology markers CD45 (2B11+PD7/26, Novus, Cat. # NBP2-34528AF532, 10 µg/mL), PDGFRβ (D-6, Santa Cruz Biotechnology, Cat. # sc-374573AF594, 1:20), CD31 (JC/70A, Abcam, Cat. # ab215912-100ul, 1:50), and syto13. Images at 20X magnification were assembled to yield high-resolution regions of interest in which appropriate masks were applied for selection of areas of interest. Areas of interest were illuminated, and released tags were collected into 96-well plates as previously described. Sequential sections were used for GeoMx Protein and RNA profiling. For protein detection, GeoMx Immune Profile Core, immuno-oncology drug target, immune activation status, immune cell typing, and pan tumor modules were used to identify 54 biomarkers. Detection of mRNA was performed for morphology marker–defined compartments according to the NanoString GeoMX RNA assay protocol described previously (41) using the Whole Transcriptome Atlas probe reagent. Library preparation was also performed according to the manufacturer's protocol. The Illumina NovaSeq 6000 S1-100 platform was used to sequence the resulting libraries. The raw counts were converted to digital count conversion files using NanoString’s GeoMx NGS pipeline v.2.0. Digital count conversion files were uploaded onto the GeoMx DSP software v2.5 for further analysis. Quality control and data processing were performed, followed by Q3 normalization.

GEPIA2 analysis Gene Expression Profiling Interactive Analysis (GEPIA, http://gepia2.cancer-pku.cn, version 2) is an open-access online tool that allows the analysis of RNA-seq data of 9,736 tumors and 8,587 normal samples from The Cancer Genome Atlas (TCGA) and the Genotype-Tissue Expression programs, respectively (42). The RNA-seq datasets GEPIA2 uses is based on the UCSC Xena project and computed by a standard pipeline. In this study, it was used to compare RNA expression between cancer subtypes and normal tissue as well as perform a survival analysis associated with the PVSC-associated MVP gene signature (a median cutoff of 50% was used).

Tg(Pdgfrb-EGFP)JN169Gsat/Mmucd mice were purchased from the Mutant Mouse Resource and Research Centers (MMRRC) at the University of California, Davis (stock number MMRRC:031796-UCD). An EGFP reporter gene derived from Aequorea victoria, followed by a polyadenylation sequence, was inserted into a BAC clone (BAC address for this line is RP23-10L16) at the initiating ATG codon of the first coding exon of the PDGFRB gene so that EGFP expression was driven by the regulatory sequences of the PDGFRB promoter. The modified BAC DNA was purified and injected into fertilized mouse eggs, surviving embryos were implanted into foster mothers, and pups were screened for the transgene. The resulting founder mice were bred to establish the line, and hemizygous progeny were mated to Crl:CD1(ICR) mice (The Jackson Laboratory, MPD ID: 866) each generation thereafter. These mice were then crossed with Tg(NES-TVA)J12Ech/J (a gift from Dr. Ganesh Rao, Baylor College of Medicine) mice. Gliomas were then induced via injection of an RCAS-PDGFB-HA virus into the striatum of 1-day-old pups. Mice without an EGFP tag were generated in the same manner except without manipulation of the EGFP gene. Mice were sacrificed when moribund. All animal manipulations were performed in the veterinary facilities at MD Anderson Cancer Center in accordance with institutional, state, and federal laws and ethics guidelines with a protocol explicitly approved by the Institutional Animal Care and Use Committee of MD Anderson Cancer Center.

We injected mice expressing the ntva-transgene from a modified nestin promoter with RCAS-PDGFB-HA (43) viruses to generate an immunocompetent model of PDGFB/PDGFRβ-driven MVP. To confirm our results, we also generated a reporter mouse model wherein we tagged PDGFRβ with constitutively expressed eGFP, crossed these mice with ntva^+/+ mice, and then injected RCAS-PDGFB-HA (43) viruses.

R v4.1.2 and GraphPad Prism software v9 were used for the statistical analyses. The tests performed have been indicated in the individual figure legends.

The data generated in this study are publicly available in the NCBI (accession number PRJNA1229526) and the code is on Github at https://github.com/ccpoon-ucalgary/PVSCs.

To characterize MVP at the single-cell level, CD105 and CD45 were used to enrich for perivascular MSCs, CAFs and immune cells, respectively, in paired glioma samples (Supplementary Fig. S1). Independent preparation of scRNA-seq libraries from CD105⁺ and CD45⁺ cells enhanced the resolution of the analyses (Fig. 1A). This study analyzed 16 treatment-naïve human glioma specimens. These included five WHO grade 4 gliomas with MVP and four without (MVP is a sufficient, but not necessary, criterion to confer a WHO grade 4 diagnosis (Table 1; ref. 35). Two WHO grade 3 and five WHO grade 2 specimens without MVP were also included, as all gliomas can dedifferentiate and acquire MVP. A board-certified pathologist (G.N. Fuller) classified the specimens. Most tumors were corticosteroid-naïve or received low doses to minimize confounding effects on immune responses.

Uniform Manifold Approximation and Projection and k-means clustering of 76,141 cells after quality control processes identified 27 cellular clusters (Fig. 1B). The contribution of individual tumors is shown in Fig. 1C. These clusters were annotated based on known cell-type marker expression (Fig. 1B; Supplementary Table S1; Supplementary Figs. S2 and S3), with transcriptionally similar clusters grouped together, such as clusters 5, 16, and 23. Clusters 5, 16, and 23 exhibited a signature consistent with the perivascular cells described in a scRNA-seq study of patients with temporal lobe epilepsy (44), but the major markers for fibroblasts (MGP and DCN), pericytes (HIGD1B and NDUFA4L2), fibromyocytes (IGFBP5), and smooth muscle cells (TAGLN and MYL9) were shared between cells, suggesting transitional cell states within glioma. Thus, these cells were designated as “fibroblastic cells.” To support this designation, we also noted that CAF markers, ACTA2, PDGFRB, and VIM, identified in GBM-associated fibroblasts reported by Jain and colleagues (17), were expressed by clusters 5, 16, and 23 (Supplementary Table S1; Supplementary Fig. S4). Notably, clusters 5, 16, 23, and 22 were negative for endothelial cell markers PECAM1 and CLDN5. Cluster 22 expressed pericyte markers SLC38A11, CARMN, and ATP1A2 (44), as well as PDGFRB. However, this cluster lacked expression of other fibroblastic markers, suggesting it represented a more differentiated pericyte cell state than the fibroblastic cells. A separate cluster of smooth muscle cells (MT2A⁺MT1X⁺MT1M⁺ (44); cluster 18) was also identified which was negative for both PDGFRB and fibroblastic markers.

Clusters 8 and 13 expressed the traditional MSC markers CD73 (NT5E), CD90 (THY1; ref. 29), and PDGFRA (16). CD105, another traditional MSC marker, was not significantly upregulated in clusters 8 and 13 because of the initial CD105 enrichment. Modules of the triple-positive CD73⁺CD90⁺CD105⁺ and PDGFRA⁺ MSCs corresponded to clusters 8 and 13 (Supplementary Fig. S4), further confirming the identification of these clusters as MSCs. An additional designation of astrocytic and fibroblastic subtypes was assigned to clusters 8 and 13, respectively, because these clusters expressed markers for astrocytes and radial glial cells (such as GFAP, S100B, and FABP7) or fibroblastic cells (such as DCN and TAGLN). MSC clusters 8 and 13 were bridged by pericyte cluster 22, suggesting a potential shared transcriptional ontogeny. In addition, fibroblastic MSCs clustered with fibroblastic cells, suggesting that fibroblastic MSCs, fibroblastic cells, and pericytes were transcriptionally similar to cells of the perivascular niche.

We next wanted to determine whether fibroblastic MSCs and fibroblastic cells resided in a perivascular location. However, this first required identification of a marker with high sensitivity and specificity for PVSCs. We performed a differential gene expression analysis of our scRNA-seq dataset between PVSCs and all other cells in the dataset. Thus, the top 50 differentially expressed candidate PVSC markers were identified in an unbiased manner (Fig. 1D). The most upregulated markers were TAGLN, a marker reported to be highly upregulated in smooth muscle cells but also all other perivascular cell classes (fibroblasts, fibromyocytes, and pericytes), RGS5, a pericyte marker, and DCN, a fibroblast marker that is also found in all perivascular cells (44, 45). Although TAGLN, RGS5, and DCN were highly sensitive markers, they were not specific. However, PDGFRB was highly sensitive and specific for PVSCs, being expressed in only 1.7% of all other cell types in the dataset (Supplementary Table S2).

Thus, antibodies to PDGFRβ in addition to the classic MSC markers CD73, CD90, and CD105 were chosen to further characterize the anatomic location of PVSCs using multiplex IF in gliomas with (n = 3) and without (n = 4) MVP (Fig. 1E; Supplementary Fig. S5). All gliomas with MVP were WHO grade 4 GBMs. A WHO grade 4 astrocytoma, a WHO grade 3 oligodendroglioma, and two WHO grade 2 oligodendrogliomas represented gliomas without MVP. To obtain a more holistic view of the perivascular niche, antibodies to markers of oligodendrocytes and tumor cells (OLIG2), microglia and macrophages (AIF1), and endothelial cells (CD31) were also selected. Gliomas with MVP had an expanded PDGFRβ⁺ PVSC cell morphology compared with gliomas without MVP. The PDGFRβ⁺ cells were distinct from the CD31⁺ endothelial cells comprising the blood vessels. To verify whether PVSC staining was distinct from endothelial cells, we also investigated other markers upregulated in our scRNA-seq dataset associated with PVSCs (TAGLN/transgelin, DCN/decorin) and endothelial cells (CD34/CD34; Supplementary Fig. S6). Only 1.37% ± 0.69% of all PDGFRβ⁺ cells colocalized with CD34 (P = 0.002), unlike the similar levels of colocalization observed between TAGLN (45.07% ± 24.0%) and DCN (28.91% ± 21.6%; P = 0.24). Classic MVP morphology seemed to be attributed to PVSCs (PDGFRβ⁺ cells). PDGFRβ⁺ cells were present but not expanded morphologically in gliomas without MVP (Supplementary Fig. S5). Two populations of MSCs were again identified, converging with our scRNA-seq findings (PDGFRβ⁻CD73⁺CD90⁺CD105⁺ cells corresponded with astrocytic MSCs (Supplementary Fig. S5), and PDGFRβ⁺CD73⁺CD90⁺CD105⁺ cells corresponded with fibroblastic MSCs (Fig. 1E)]. Both astrocytic and fibroblastic MSCs were identified at a low frequency as expected of stem cell populations.

We manually annotated 20 vessels per specimen to quantitatively evaluate fibroblastic and astrocytic MSCs in tumors with and without MVP and in nonproliferative, “normal-adjacent” areas of tumors with MVP to address intratumoral heterogeneity (Supplementary Fig. S7). Only PDGFRβ⁺ MSCs were significantly elevated (P = 0.009, unpaired t-test with the Welch correction) in gliomas with MVP versus without, indicating that fibroblastic MSCs are associated with the perivascular compartment and are associated with MVP, not just MSCs of any subtype (Fig. 1F). No differences were found between control and normal-adjacent perivascular microenvironments, although the trend supported a possible upregulation of fibroblastic MSCs in the normal-adjacent regions. Furthermore, CD105⁺/PDGFRβ⁺ cells were enriched in MVP areas (P = 0.0073; Supplementary Fig. S8), indicating that in addition to morphologic expansion, they were also present in increased numbers. To verify whether CD105⁺ and CD105⁺/PDGFRβ⁺ cells were found in a perivascular location, a distance analysis was performed to identify the quantity of cells within or beyond 50 µm increments from a blood vessel (Supplementary Fig. S9). The majority of CD105⁺ and CD105⁺/PDGFRβ⁺ cells were found in the immediate proximity of blood vessels, validating the use of these markers to enrich for perivascular cells. In summary, it seemed PVSCs played a major role in MVP morphology.

To understand the functional role of PDGFRB as a marker of glioma-associated PVSCs, MSCs (clusters 8 and 13) and adjacent clusters (clusters 3, 5, 12, 16, 18, 22, and 23) were reclustered to identify potentially associated cell states (Fig. 2A). Similar to our previous clustering, endothelial cells formed a distinct cluster (Fig. 2A), indicating transcriptional separation from PVSCs. Again, MSCs exhibited transcriptional clustering as either astrocytic (Fig. 2A, cluster 2; NT5E+PDGFRA+GFAP+FABP7+) or fibroblastic (Fig. 2A, cluster 4; THY1+PDGFRA+IGFBP5+APOD+DCN+PDGFRB+) subtypes (Fig. 1B). Now, however, we were also able to capture a cluster of less-differentiated MSCs (Fig. 2A, cluster 5) expressing markers associated with stemness, including SOX2, CD44, and PDPN. To estimate cellular trajectories, we performed an RNA velocity analysis with dynamic modeling (Fig. 2B) and found that cluster 5 demonstrated bidirectionality toward clusters 2 and 4. The fibroblastic cells (clusters 0, 4, 6, 8, and 10) clustered together and exhibited multidirectionality between clusters, supporting a transitional nature. Furthermore, latent time analysis [which represents the cell’s internal clock and position in an underlying biological process (37)] supported cluster 5 as a more nascent cell population (Fig. 2C), with cluster 6 (fibroblastic, pericytic, and fibromyocyte cells) representing a more differentiated cell population (Fig. 2A). Furthermore, RNA velocity analysis (Fig. 2D and E) revealed PDGFRB was a putative driver gene that was induced transcriptionally during the fibroblastic MSC stage (Fig. 2A, cluster 4) through the intermediate fibroblastic, pericytic, and smooth muscle cell stages (Fig. 2A, clusters 0, 8, and 10) until it was downregulated but still present in the fibroblastic, pericytic, and fibromyocyte cell stages (Fig. 2A, cluster 6). PDGFRB is expressed throughout PVSC differentiation (Fig. 2F). PDGFRB expression is highest immediately following divergence from less-differentiated MSCs (Fig. 2A, Cluster 5). Although downregulated, PDGFRB persists in terminal fibroblastic, pericytic, and fibromyocyte cell states (Fig. 2A, Cluster 6), supporting its role as a PVSC marker in human glioma. A schematic of the hierarchy of differentiation between cell states is proposed in Fig. 2G, although it is important to note that dedifferentiation may also occur. Taken together, our data support the involvement of PDGFRB/PDGFRβ in MVP.

To validate our findings that PDGFRB can be a marker of PVSCs associated with MVP in vivo, we reverse-translated our findings in human tissue and generated a reporter mouse model of MVP (Supplementary Fig. S10). The morphology of MVP was recapitulated by the eGFP⁺ (PDGFRβ⁺) cells adjacent to the CD31⁺ cells (Supplementary Fig. S10), akin to what we observed in the human condition (Fig. 1E). In mice, MVP was not seen in the contralateral, non–tumor-containing hemisphere despite the presence of eGFP⁺ cells (Supplementary Fig. S10F), suggesting that microenvironmental factors also contribute to MVP. Taken together, our data supports the involvement of the PDGFRB/PDGFRβ pathway in MVP.

We next investigated whether cellular composition differed between gliomas with and without MVP [33,295 cells associated with MVP (43.7% of all cells), n = 5 gliomas versus 42,846 cells not associated with MVP (56.2% of all cells), n = 11 gliomas; Fig. 3A and B]. Greater than 75.0% of fibroblastic MSCs, brain-resident perivascular macrophages, smooth muscle cells, and HLA-DRA-low (HLA-DRA^lo) suppressive monocytes were found in gliomas with MVP. Conversely, greater than 75.0% of all intermediate monocytes, plasmacytoid dendritic cells, and IFN-stimulated gene–expressing effector T cells were found in gliomas without MVP. The enrichment of fibroblastic MSCs, smooth muscle cells, HLA-DRA^lo suppressive monocytes, and brain-resident perivascular macrophages coupled with depletion of IFN-stimulated gene–expressing effector T cells, intermediate monocytes, and plasmacytoid dendritic cells in gliomas with MVP suggests an association between MVP and immunosuppressive modulation. Brain-resident perivascular macrophages expressed CD105, further endorsing its use to enrich for diverse cell types within the glioma perivascular niche, including myeloid cells, endothelial cells, and PVSCs.

To investigate functional differences between PVSCs (clusters 5, 13, 16, 22, and 23) between gliomas with and without MVP, we performed GO analysis to identify the top differentially enriched biological processes (Fig. 3C). PVSCs associated with MVP demonstrate enriched pathways for extracellular matrix (ECM) organization, extracellular structure, and external encapsulating structure organization. This enrichment aligns with the substantial morphologic changes observed in PVSCs of gliomas with MVP compared with those without MVP. Other important processes observed preferentially in MVP-associated PVSCs included the positive regulation of cell adhesion, apoptotic signaling pathways, and cellular oxidant detoxification, all of which are known to promote cancer invasion and propagation. One of the top enriched biological processes was the negative regulation of immune system processes, concordant with our observations of cellular subset differences predisposing to immunosuppression (as discussed in Fig. 3A and B).

TGF-β pathway genes constitute a key component of the GO definition of “negative regulation of immune system processes.” To corroborate whether PVSCs are associated with immunosuppression and share similarities with CAFs outside the CNS, the gene signature of glioma PVSCs was compared with that of ovarian cancer fibroblasts, known to promote immunosuppression via TGF-β pathway upregulation (28, 46). We found that the ovarian cancer fibroblast signature aligned with glioma PVSCs and was more highly upregulated in glioma PVSCs with MVP compared with those without MVP (Supplementary Fig. S11). This finding supports the association of glioma PVSCs, including CAFs, with immunosuppression in MVP conditions. Glioma PVSCs are transcriptionally similar to CAFs in non-CNS cancers, suggesting the generalizability of these results.

To further investigate the immunosuppressive networks active in PVSCs, CytoTalk was used to construct unbiased, de novo signal transduction networks between PVSCs and all immune cell subsets. More TGF-β pathway ligand–receptor interactions were observed in the MVP microenvironment than in the non-MVP microenvironment (Supplementary Tables S3 and S4). Additionally, in the MVP cohort of gliomas, the brain-resident perivascular macrophages and HLA-DRA^lo suppressive monocytes were shown to possess receptors involved with the TGF-β axis (TGFBR1/TGFBR2/ENG), whereas the PVSCs possessed the corresponding ligands (TGFB2/TGFB3; Fig. 3D and E; Supplementary Tables S3 and S4, respectively). Moreover, the type of ligand–receptor interactions in the TGF-β pathway varied in gliomas with and without MVP. In gliomas with MVP, PVSCs communicated with immune cells via TGFB2/TGFB3 and TGFBR1/TGFBR2/ENG, respectively. In gliomas without MVP, these interactions occurred through SDC2 on the immune cells and TGFB1 from the PVSCs, suggesting that when attempting to inhibit the TGF-β pathway, certain ligands and receptors will be higher-value targets and that the TGF-β pathway plays different roles in MVP versus non-MVP microenvironments (Supplementary Tables S3 and S4). Overall, upregulated TGF-β pathway activity and increased frequencies of brain-resident perivascular macrophages and HLA-DRA^lo suppressive monocytes are associated with PVSCs in the MVP microenvironment.

Fatty acid metabolism (47, 48) and glycolysis (49) are metabolic states associated with immunosuppression. Thus, we investigated whether PVSCs were associated with altering immune metabolic states using the fatty acid metabolism and glycolysis hallmark gene sets (50). Expressions of both gene sets were higher in HLA-DRA^lo suppressive monocytes and brain-resident perivascular macrophages compared with other immune cells (Supplementary Figs. S12 and S13). Each gene set was cross-referenced with signaling pathways in CellChat, a tool for delineating intercellular communication networks (51). The two communication networks with the greatest gene overlap with the fatty acid metabolism pathway were the bone morphogenetic protein (BMP) and thrombospondin pathways, whereas collagen and EGF network genes overlapped with the glycolytic pathway (Supplementary Figs. S12 and S13). To determine which ligand–receptor pairs were involved in these pathways, we identified the associated PVSC cell types and quantified the strength of these associations using CellChat. Of the PVSCs, fibroblastic MSCs drove the BMP pathway in fatty acid metabolism via BMP8A, affecting not only HLA-DRA^lo suppressive monocytes and brain-resident perivascular macrophages but also brain-associated macrophages, microglia, and intermediate monocytes through multiple receptors (Supplementary Fig. S12C). Similarly, fibroblastic MSCs had a sizeable effect on HLA-DRA^lo suppressive monocytes and NKT cells via THBS1 and CD47 (Supplementary Fig. S12D). The collagen pathway had considerable overlap with the glycolysis hallmark gene set, with fibroblastic MSCs and fibroblastic cells acting as sources of collagen for various integrins expressed by immune cells (Supplementary Fig. S13C). The EGF pathway was driven mostly by smooth muscle cells expressing TGFA and EGF interacting with EGFR on myeloid and T cells (Supplementary Fig. S13D). There was no overlap observed between oxidative phosphorylation hallmark gene sets and CellChat pathways upregulated in PVSCs, suggesting that PVSCs are not involved in the oxidative phosphorylation of immune cells. In sum, PVSCs may play a role in influencing innate and adaptive immune cell types toward fatty acid metabolism and glycolysis, metabolic states associated with immunosuppression.

To compare the immunosuppressive contribution of other vascular/perivascular cell types such as astrocytic MSCs and endothelial cells with PVSCs, we investigated CellChat networks upregulated in our dataset known to be associated with immune evasion (TGFb, IL10, VEGF, and MIF), checkpoint inhibition (VISTA), and immune cell recruitment (CXCL, CCL, and CX3C) in GBM (Supplementary Fig. S14). The expression of individual checkpoint molecules is displayed in Supplementary Fig. S15. The majority of the communication networks investigated revealed substantial contributions from PVSCs, with fibroblastic MSCs often driving the contribution. The EGF, VISTA, and CXCL pathways were the exception in which smooth muscle cells and endothelial cells contributed more significantly. Overall, PVSCs participate in the majority of vascular/perivascular immunosuppressive activities investigated in this dataset. Furthermore, besides the CXCL pathway, fibroblastic MSCs may play a role in recruiting myeloid cell and T-cell populations through CCL26-, CCL2-, and CX3CL1-mediated signaling (Supplementary Fig. S14G and S14H), offering a potential mechanism whereby immunosuppressive populations are maintained.

To support the involvement of HLA-DRA^lo suppressive monocytes and brain-resident perivascular macrophages in the perivascular niche, we performed DSP. We used antibodies to identify cells positive for CD45, PDGFRβ, and CD31 to optimally delineate areas of CD45 positivity around blood vessels (Fig. 4A). Three GBMs represented samples with MVP. Controls included a GBM, a WHO grade 3 oligodendroglioma, and a WHO grade 2 oligodendroglioma. When examining cell types, the top differentially upregulated protein markers included CD14 and CD163, supporting HLA-DRA^lo suppressive monocytes (CD14⁺) and brain-resident perivascular macrophages (CD163⁺) were upregulated in the perivascular niche in gliomas with MVP (Fig. 4B), corresponding with our scRNA-seq dataset (Fig. 3A and B). Similar to our transcriptomic dataset, NK cells (CD56) were significantly downregulated.

When examining inflammatory markers, we observed that the immunologic activity of CD45⁺ cells was starkly higher in the MVP-associated perivascular regions compared with the regions without MVP (Fig. 4B). Although the immunosuppressive markers ARG1 and B7-H3 (CD276) were some of the top differentially expressed proteins, immunostimulatory markers STING (STING1), 4-1BB (TNFRSF9), ICOS, OX40L (TNFSF4), and GZMB were also upregulated. We examined the level of expression of these markers in our scRNA-seq dataset. We found that STING1, TNFSF4, and ARG1 were not detectable or barely detectable. T-cell subsets predominantly expressed TNFRSF9, ICOS, and GZMB, indicating that adaptive, rather than innate, immune populations contributed to these immunostimulatory markers in the perivascular niche (Supplementary Fig. S16). CD276 was expressed by both brain-associated macrophages and PVSCs. Overall, the decreased quantity of NK cells, the increased presence of CD14⁺ and CD163⁺ myeloid cells, and the expression of CD276 by macrophages and PVSCs supports a myeloid- and PVSC-driven immunosuppressive phenotype in the perivascular niche of gliomas with MVP.

DSP was also used to obtain a transcriptomic pathway analysis of the perivascular region. Upregulation of PDGF signaling (P = 0.0008), TGF-β receptor complex signaling (P = 0.0008), ECM organization (P = 0.0008), and collagen formation pathways (P = 0.0008) were observed in the perivascular region of gliomas with MVP (Fig. 4C; Supplementary Table S5). Further examination of pathways involving PDGFRβ revealed an upregulation of MAPK family signaling cascades (P = 0.0008), diseases of signal transduction by growth factor receptors and second messengers (P = 0.0008), signaling by receptor tyrosine kinases (P = 0.0001), MAPK1/MAPK3 signaling (P = 0.004), RAF/MAP kinase cascade (P = 0.005), and FLT3 Signaling (P = 0.01; Supplementary Table S5). When examining which cells provided the ligands to drive PDGFB signaling, CellChat analysis revealed a substantial contribution from myeloid populations (Supplementary Fig. S17). These findings are consistent with the results of our previous assays and further support involvement of the PDGFR and TGF-β pathways in MVP, along with an intimate association with immune activity.

To assess the clinical relevance of MVP-associated PVSCs, we investigated whether hallmark gene sets involved with tumor progression were enriched. The angiogenesis, epithelial–mesenchymal transition, and hypoxia gene sets were all significantly upregulated in MVP-associated PVSCs compared with non–MVP-associated PVSCs (P ≤ 2.2 × 10⁻¹⁶; Supplementary Fig. S18). We also investigated the contribution of PVSC subsets to communication networks involved with tumor growth and vascular integrity, such as NOTCH, FGF, ANGPT, and WNT (Supplementary Fig. S19). The NOTCH, FGF, and ANGPT networks had large contributions from fibroblastic MSCs and fibroblastic cells, followed by a small contribution from pericytes. The WNT pathway, however, saw a larger contribution from smooth muscle cells, indicating that PVSC subtypes play different roles in tumor growth. In addition, a gene signature was developed to facilitate investigation of the clinical relevance of MVP-associated PVSCs in TCGA datasets. The top DEGs in PVSCs with MVP compared with those without MVP were identified to derive a gene signature (Fig. 5A). GEPIA2 (42) was used to analyze the expression levels of the top five DEGs—COL3A1, PCOLCE, TIMP1, NNMT, and COL1A1—in the TCGA-GBM, TCGA-LGG, and Genotype-Tissue Expression (ref. 52) datasets. These genes demonstrated the highest expression in GBM, followed by LGG, and had the lowest expression in normal tissues (Fig. 5B). In the TCGA cohorts of LGG, there was a significant difference in overall (P = 3.9 × 10⁻⁹; Fig. 5C) and disease-free survival (P = 0.00029; Fig. 5D), suggesting that this signature could identify LGGs that are at risk of developing MVP and progressing to a higher grade. In the TCGA cohorts of GBM, there were no differences in disease-free (P = 0.22; Supplementary Fig. S20A) and overall survival (OS; P = 0.42; Supplementary Fig. S20B), as expected. These GBM specimens were classified using older nonmolecular WHO guidelines, and likely all had MVP. Investigating GBM classified using newer guidelines in a published cohort (53), our signature correlated with disease-free survival (P = 0.045; Supplementary Fig. S20C) and demonstrated a trend toward OS (P = 0.052; Supplementary Fig. S20D). Expression of COL3A1, PCOLCE, TIMP1, NNMT, and COL1A1 individually correlated with decreased OS in LGG (Supplementary Fig. S21), supporting the robustness of the combined gene signature. This signature also predicted OS in TCGA datasets of clear cell renal cell carcinoma (P = 0.0095), squamous cell carcinoma of the lung (P = 0.043), muscle-invasive bladder cancer (P = 0.044), mesothelioma (P = 0.033), and uveal melanoma (P = 0.0029; Fig. 5E–I). These findings suggest the broader applicability of this signature to non-CNS malignancies and that the MVP process is not limited to gliomas. Overall, COL3A1, PCOLCE, TIMP1, NNMT, and COL1A1 represent an ECM molecule–rich signature that has the potential to be targeted or used as a biomarker in glioma and other cancers.

Although MVP is a crucial clinical marker for prognosis, it remains poorly characterized at the single-cell level, and its impact on the glioma microenvironment is not well understood. In this study, we captured MSCs in both low- and high-grade human gliomas for the first time. Our analysis implicates MSCs and their progeny, including CAFs (collectively termed “PVSCs”), in the pathogenesis of MVP.

Our study found that MSCs, CAFs, pericytes, fibromyocytes, and smooth muscle cells (i.e., PDGFRB⁺ PVSCs) all appeared to exist in transitional cell states, whereas similar transitional states are not seen in nonmalignant, human brain–perivascular niches (44). Establishing this connection between the cell states of PVSCs at the single-cell level has many implications, because MSCs originate from the periphery, particularly the BM (16, 20, 21, 54), and can differentiate into various cell types. Furthermore, we used unbiased approaches to identify PDGFRB as a candidate driver gene for PVSCs that can become involved in MVP. Our study raises the possibility that MVP may be driven, at least partly, by MSCs attracted from the periphery that can populate the perivascular niche with other PVSCs. The implication that PVSCs are derived from MSCs suggests that MSCs may represent a common pathway that can be targeted to disrupt MVP. As such, systemic therapeutics could target MSCs in the periphery without needing to cross the blood–brain barrier. Antiangiogenics targeting MVP and endothelial cells, such as bevacizumab, have shown limited success in multiple clinical trials (55–57). Our findings suggest that an added focus on PVSCs may be warranted to develop more effective treatments. Although we recognize a major limitation of the nature of clinical specimen-based work renders it correlative and lacking in functional validation, our goal was to propose PVSCs (which includes MSCs and CAFs) as an overlooked cellular population in glioma that may have important implications for MVP genesis and create a platform for future hypothesis generation.

CAFs were only recently identified in the GBM microenvironment by Jain and colleagues, but trajectory analyses and ontogeny could not be addressed (17). Evaluating CAFs in GBM has been challenging due to inconsistent markers, extensive tissue manipulation required for isolation, selection biases introduced by current isolation techniques for late-stage α-SMA (ACTA2)-expressing CAFs, and cell scarcity that may be a product of these other challenges (16, 17, 58–65). This scarcity hampers the isolation of sufficient CAF quantities for comprehensive studies such as trajectory and ontogeny analyses (17). In our work, by enriching for the MSC marker CD105 and leveraging the well-defined states of glioma progression (WHO grades 2–4), we captured enough CAFs and transcriptionally related cells to begin dissecting their ontogeny and potential roles in glioma. Our dataset demonstrates that isolating CD105⁺ cells is a viable strategy for capturing not only PVSCs in glioma but also other vascular/perivascular niche cell types, including endothelial cells and brain-resident perivascular macrophages, for scRNA-seq. Unexamined CD105^- PVSCs may influence the perivascular niche and currently have an unknown role in MVP, meriting additional research.

Prior to the advent of scRNA-seq, mouse studies suggested BM precursors, potentially “pericyte progenitor cells" (PPC), migrated to non-CNS tumors to become CAFs (18, 66–68). However, the incorporation of PPCs into CNS tumors has been controversial as no incorporation was observed in orthotopic GL261 (67) and neural stem cell/U87(68) glioma models, whereas only low frequencies of PPCs were seen in Rag1-deficient mice with tumors derived from transformed astrocytes (69). Our dataset nominates MSCs/PVSCs as the potential human correlate to this PPC population, as PVSCs and PPCs share the marker, PDGFRβ (66, 69), and MSCs can originate from the BM (16, 20, 21, 54). Our human dataset suggests that CAFs and related cells can migrate from the periphery into the glioma perivascular space despite the controversial evidence from previous mouse models of glioma.

MSCs are multipotent stromal cells that, like CAFs, exhibit diverse protumorigenic and immunomodulatory properties (70–73). We identified MSC subtypes, including an astrocytic and fibroblastic population, that may colonize the perivascular niche. An even less-differentiated MSC population that potentially gives rise to the astrocytic and fibroblastic subtypes was also discerned. Exploring markers of this less-differentiated MSC population could inform future lineage tracing studies to identify the cell of origin for glioma PVSCs. RNA velocity analysis identified distinct PVSC subtypes. PVSC subtypes were also demonstrated in our RNA velocity analysis. Our analyses revealed that fibroblastic MSCs and fibroblastic cells contributed heavily to immunosuppression, immune cell metabolism, and tumor progression. However, smooth muscle cells showed strong involvement in EGF and CXCL pathways. These findings underscore the distinct roles of PVSC subtypes in cellular processes and highlight the need to develop cell type–specific therapies.

Additionally, we sought to investigate how PVSCs affected the microenvironment. PVSCs undergoing MVP were associated with immunosuppressive changes. For instance, more than 75% of the HLA-DRA^lo suppressive monocytes and fewer than 25% of IFN-stimulated gene–expressing effector T cells were found in MVP-associated gliomas. This suggests MVP is linked to the upregulation of immunosuppressive and downregulation of proinflammatory immune cell subsets. This is despite the inclusion of WHO grade 4 gliomas in the MVP-free control group and therefore not a result of comparing WHO grade 4 gliomas with LGGs. We used unbiased de novo signal transduction analyses to identify TGF-β pathway interactions between PVSCs and HLA-DRA^lo suppressive monocytes/perivascular macrophages. Ligands were attributed to PVSCs, whereas receptors were associated with myeloid cells, suggesting PVSCs instigated these interactions. Our protein-level DSP suggested enrichment of these specific myeloid subsets in the perivascular niche. However, not all immune interactions in the perivascular niche were prominently immunosuppressive, a fact more easily uncovered in our dataset as steroid use was limited. Although ARG1 and B7-H3, an emerging checkpoint inhibitor molecule, were highly differentially expressed, there was also upregulation of immunostimulatory markers STING, 4-1BB, ICOS, OX40L and GZMB. Cross-referencing with our scRNA-seq dataset showed that most proinflammatory markers were expressed by T and NKT cells, not myeloid cells. This indicates opposing immunosuppressive and immunostimulatory forces ascribed to myeloid and adaptive immune cell types, respectively, actively vying for control in the perivascular niche. In addition to being a key passageway for immune cell trafficking into and out of the CNS, the glioma perivascular niche is thus poised as an important battleground for treatments that can help sway the microenvironment toward more proinflammatory, antitumor phenotypes. Our study expands upon findings from previous literature (74, 75), proposing specific myeloid cell subsets in the perivascular niche are associated with MVP in human gliomas.

Beyond glioma, the generalizability of our results is suggested not only by the shared gene signature between TGF-β ovarian CAFs (46) and MVP-associated PVSCs but correlation of our MVP-derived signature (COL3A1, PCOLCE, TIMP1, NNMT, and COL1A1) with decreased overall and disease-free recurrence in LGG and five other solid-tumor types in TCGA. It is expected that OS was not significantly different in TCGA GBMs with low versus high MVP signatures, as essentially all these tumors have MVP as per the older WHO definition at the time of collection (76, 77). Recent work on paired primary and recurrent GBM specimens showed an increase in ECM-associated genes at recurrence that was attributed to pericytes (78). As per our earlier discussion, we suspect these “pericytes” identified by PDGFRβ are actually PVSCs. In that study, COL1A1 was also identified as an ECM molecule associated with a more mesenchymal, aggressive phenotype, previously linked to “oncostreams” with “fibroblast-like morphologic changes” (79). Our study demonstrates a convergence of knowledge and offers a missing puzzle piece, suggesting that COL1A1 upregulation is attributed to PVSCs and that the ECM molecules COL3A1 and PCOLCE may be other druggable targets relevant to glioma progression. Additional targets identified in our study, TIMP1 and NNMT, are involved with ECM remodeling (80, 81) and known to be expressed by non-CNS fibroblasts and CAFs (82, 83). A recent study found that the presence of glioma-associated fibroblasts alone did not seem to be directly linked to changes in OS, as the presence of any one of five glioma-associated fibroblast markers did not alter outcome in the TCGA cohort (17). Instead, our data support that it is specifically the signature of PVSCs associated with MVP that is connected with prognosis. Each gene in the MVP-associated signature independently correlated with decreased OS in LGG, emphasizing the robustness of this signature. This signature may help stratify patients who could be more likely to progress and benefit from earlier, more aggressive interventions. In addition to using this signature as a biomarker, the identification of PDGFRβ and CD73/CD90/CD105 as markers of glioma stromal activity suggests these markers could be valuable in clinical practice, particularly in the context of antiangiogenic therapies. The finding that PVSCs, especially fibroblastic MSCs and fibroblastic cells, may contribute to immunosuppression indicates that these markers could predict therapeutic response, especially in immunotherapy-resistant gliomas. Gliomas without MVP that exhibit low expression of these markers, as well as the COL3A1, PCOLCE, TIMP1, NNMT, and COL1A1 signature, may represent optimal candidates within a therapeutic window for aggressive immunotherapy before MVP develops. This study proposes that targeting PDGFRβ+ stroma, using commercially available inhibitors for example, could enhance concurrent immunotherapies such as immune checkpoint blockade, especially if future functional studies confirm the role of PVSCs in immune evasion. Similarly, combining previously unsuccessful antiangiogenic therapies such as bevacizumab with PDGFRβ inhibitors may offer a new strategy to target angiogenesis and tumor progression.

As alluded to earlier, one major limitation of our study is its largely correlative nature, as it is heavily based on clinical specimens. Future functional validation will be crucial for establishing the direct role of PVSCs in immunomodulation. Using the mouse models proposed in this study, it will also be important to elucidate the mechanisms underlying PDGFRβ+ stroma-associated immunosuppression and tumor progression. However, given the lack of fibroblast/CAF atlases derived from human brain/glioma tissue (33, 34) and the paucity of databases focusing on perivascular biology, our data help fill a gap in the current literature. Especially because fibroblasts/CAFs and MSCs exhibit organ-specific transcriptomic profiles (62, 84, 85), it is imperative to characterize them in anatomically distinct regions such as the brain in states of normalcy and pathology. To our knowledge, this dataset is one of the largest collections of PVSCs in high- and low-grade human gliomas. By excluding chemoradiation-treated samples and mitigating steroid effects, we attempted to accurately represent the “natural” glioma stromal and immune microenvironment. Lastly, to definitively prove PVSCs are not dedifferentiated endothelial cells or have not fused with endothelial cells will require further study. Our transcriptomic data showed that endothelial cell markers consistently clustered separately from PVSCs, and our multiplex IF data suggest that if dedifferentiated or fused endothelial cells did play a role, their role is very minor.

In conclusion, our work provides a high-dimensional perivascular cell atlas and a single-cell resolution link between MSCs, CAFs, fibromyocytes, pericytes, and smooth muscle cells that helps to start unraveling the confusion and complexity surrounding the identification and ontogeny of PVSCs in human glioma. Our data suggest that MSCs of the fibroblastic subtype are found in the perivascular niche, and that these MSCs may be able to colonize this niche with other PVSC subtypes that exist in overlapping cell states. Because our dataset enriched for both PVSCs and the immune milieu, it is also poised to help elucidate the intricate interactions between PVSCs and immune cells. PVSCs are potentially important targets and biomarkers not only because of their role in MVP but because they may be associated with immunosuppression and ECM changes that promote progression (Supplementary Fig. S22). Overall, our work supports focusing on PVSCs to reveal mechanisms to combat glioma.

F.F. Lang reports other support from DNAtrix outside the submitted work; in addition, F.F. Lang has a patent for Delta-24-RGD pending. P. Sharma is a Scientific Advisory Board Member for Achelois, Adaptive Biotechnologies, Affini-T, Akoya Biosciences, Apricity, Asher Bio, BioAtla LLC, BioNTech, Candel Therapeutics, Catalio, C-Reveal Therapeutics, Dragonfly Therapeutics, Earli Inc, Enable Medicine, Glympse, Henlius/Hengenix, Hummingbird, ImaginAb, InterVenn Biosciences, JSL Health, LAVA Therapeutics, Lytix Biopharma, Marker Therapeutics, Matrisome, Oncolytics, Osteologic, PBM Capital, Phenomic AI, Polaris Pharma, Soley Therapeutics, Sporos, Spotlight, Time Bioventures, Trained Therapeutix Discovery, Two Bear Capital, and Xilis Inc, and the work in this article is not related to these companies. No disclosures were reported by the other authors.

C.C. Poon: Conceptualization, data curation, formal analysis, validation, investigation, visualization, methodology, writing–original draft, writing–review and editing. S.M. Herbrich: Formal analysis, validation, methodology, writing–review and editing. Y. Chen: Formal analysis, validation. A. Hossain: Investigation, methodology. G.N. Fuller: Data curation, validation, methodology. S. Jindal: Data curation, formal analysis, validation, investigation, visualization, methodology. S. Basu: Data curation, formal analysis, validation, investigation, methodology. D. Ledbetter: Investigation. M. Macaluso: Formal analysis, investigation, visualization, methodology. L.M. Phillips: Validation. J. Gumin: Data curation. Z. He: Data curation, formal analysis, investigation, methodology. B.C. Parker Kerrigan: Resources, data curation. S.K. Singh: Validation. P. Singh: Data curation, formal analysis, validation, investigation, methodology. M.F. Zaman: Investigation. D. Ng Tang: Investigation, project administration. S. Goswami: Conceptualization, validation, investigation. F.F. Lang: Conceptualization, resources, data curation, supervision, funding acquisition, validation, investigation, methodology, writing–original draft, writing–review and editing. P. Sharma: Conceptualization, resources, data curation, formal analysis, supervision, funding acquisition, validation, investigation, visualization, methodology, writing–original draft, writing–review and editing.

We acknowledge Liangwen Xiong, Baoxiang Guan, Lihong Long, and Jing Yang for technical assistance. We acknowledge the CATALYST-working group at MD Anderson Cancer Center for human glioma samples. C.C. Poon reports grants from the Parker Institute for Cancer Immunotherapy and T.C. and Jeanette Hsu Foundation, the NCI, the Broach Foundation for Brain Cancer Research, the Elias Family Fund, the TLC2 Foundation, and the Ergon Foundation Award during the conduct of the study. P. Sharma is a member of the James P. Allison Institute. This work was supported by the Parker Institute for Cancer Immunotherapy and T.C. and Jeanette Hsu Foundation (P. Sharma), the NCI (1R01CA214749, 1R01CA247970, and 2P50CA127001; F.F. Lang), the Broach Foundation for Brain Cancer Research (F.F. Lang), the Elias Family Fund (F.F. Lang), the TLC2 Foundation (F.F. Lang), R37 (1 R37 CA279192-01), the M.D. Anderson Physician Scientist Award (S. Goswami), and the Ergon Foundation Award (S.M. Herbrich). CATALYST is supported by the MD Anderson GBM Moon Shot.