Hypoxia coordinates the spatial landscape of myeloid cells within glioblastoma to affect survival

Abstract

Myeloid cells are highly prevalent in glioblastoma (GBM), existing in a spectrum of phenotypic and activation states. We now have limited knowledge of the tumor microenvironment (TME) determinants that influence the localization and the functions of the diverse myeloid cell populations in GBM. Here, we have utilized orthogonal imaging mass cytometry with single-cell and spatial transcriptomic approaches to identify and map the various myeloid populations in the human GBM tumor microenvironment (TME). Our results show that different myeloid populations have distinct and reproducible compartmentalization patterns in the GBM TME that is driven by tissue hypoxia, regional chemokine signaling, and varied homotypic and heterotypic cellular interactions. We subsequently identified specific tumor subregions in GBM, based on composition of identified myeloid cell populations, that were linked to patient survival. Our results provide insight into the spatial organization of myeloid cell subpopulations in GBM, and how this is predictive of clinical outcome.

Fig. 1.. Characterization of cell populations present in the GBM TME using IMC.

(A) Overview of the GBM patient samples obtained from Salford Royal Foundation Trust (STFT) and IMC workflow. (B) Representative IMC and H&E images from GBM cases taken from either the edge or core of the tumors, visualizing key microglial, macrophage, and neoplastic markers. (C) UMAPs visualizing the single-cell data acquired from the IMC workflow, demonstrating the distribution of the myeloid, neoplastic, and vascular markers over all cells. Each marker is normalized to the 99th percentile of its expression. (D) Heatmaps showing the mean marker expression for the populations identified in the single-cell IMC data using Leiden clustering. (E) UMAP showing the labeled cell populations resulting from Leiden clustering. (F) Comparison of the abundances of the different myeloid and nonmyeloid populations between the edge and core regions of the tumors. *P < 0.05, groups compared by multiple linear regression. Data shown as mean ± SEM. (G) Diffusion pseudotime and PAGA analysis demonstrating a pathway of myeloid differentiation starting at homeostatic microglia (Mg-Ho), through to proinflammatory activation microglia (TAM-Mg), proinflammatory macrophages (TAM-Mac), and ultimately to immunosuppressive myeloid cells (TAM-Supp).

Fig. 2.. Myeloid cells exhibit high homotypic and low heterotypic clustering behavior in GBM.

(A) Mapping of the myeloid populations identified from single-cell analyses to their locations in the TME in edge and core regions. (B) Diffusion pseudotime showing how differentiation away from a homeostatic microglial phenotype relates to myeloid cell positioning. (C) Spatial distribution analysis which shows whether populations are closer or further away from each other in the TME than would be expect by chance. (D) IMC image showing separation of microglia (P2RY12⁺) and macrophages (CD163⁺ and CD206⁺), alongside the position of the cells on the microglia-macrophage diffusion pseudotime axis, and how the cells are connected in the six–nearest neighbors (kNN) analyses. (E) Proportion of cellular interactions made by myeloid populations when each myeloid cell is connected to its six nearest neighbors. Hierarchical clustering of the interactions shows that cells with a similar phenotype also have similar proportions of interactions with other populations. (F) Clustering coefficient of the different populations in the edge and core of the tumor. (G) Correction of clustering coefficients for differences in cell abundance, in which a positive value suggests that cells are clustering at a greater rate than expected by chance. (H) Comparison of assortativity of myeloid populations in edge and core regions. Comparison made by Wilcoxon test with Benjamini-Hochberg correction (C and E), multiple linear regressions, corrected for multiple comparisons by Holm-Šídák (F and G), or Mann-Whitney U test (H). Box plots show range, interquartile range, and median of data (F to H).

Fig. 3.. Myeloid cell interactions with tumor, neuroglial, and vascular cells in GBM.

(A) Spatial density of different cell populations (calculated by Gaussian kernel density estimation) in the TME. (B) Myeloid cells mapped alongside tumor, neuroglial, and vascular cells in representative edge and core regions. (C) Quadratic entropy for the myeloid populations in the edge and core of the tumor. This quantifies how heterogeneous a cell is with respect to its interacting cell, with high values indicating that a cell interacts with several cells with different phenotypes. (D) Representative examples of myeloid populations mapped alongside nonmyeloid populations. (E) Difference between the observed rate of cell-to-cell interaction data (see fig. S3C) and the rate of interactions predicted by chance. This analysis corrects for the established differences in abundances of the myeloid and nonmyeloid cells between different regions. (F) Comparison between edge and core regions in the rate of infiltration of TAM-Supp cells into the parenchyma, defined here as being 10 μm away from nearest vascular cell. Comparisons made by linear mixed models (B, C, and E) or Wilcoxon tests (E) with Holm-Šídák corrections or Mann-Whitney U test (F). Box plots show range, interquartile range, and median of data (C, E, and F).

Fig. 4.. Association of hypoxia and fibrinogen with myeloid cell positioning in GBM.

Myeloid populations mapped alongside markers of (A) environmental hypoxia in GLUT1 and pERK1/2 staining, or (D) fibrinogen. Quantification of environmental GLUT1 (B), pERK1/2 (C), and fibrinogen (E) staining around each myeloid population in edge and core regions. A cell’s environment was defined as a 40-μm square centered on the cell. (F) Comparison of environmental stains in TAM-Supp cells that were either vascular associated, perivascular, or fully infiltrated into the tumor, in both edge and core regions. (G) Representative examples of TAM-Supp cells with different vascular associations in a normoxic edge region and a hypoxic core region. Comparisons made by linear mixed models with Holm-Šídák corrections (B, C, E, and F). Violin plots show range, interquartile range, and median of data (B, C, E, and F).

Fig. 5.. ST reveals hypoxia and chemokines as determinants of myeloid cell positioning in GBM.

(A) UMAP of the myeloid populations identified by Leiden clustering in the Ruiz-Moreno et al. scRNA-seq dataset. (B) Comparison between the transcriptomes of the myeloid populations and published gene lists for biological processes (MSigDB) (79), cell identities (PanglaoDB) (80), myeloid cell phenotype in GBM, and other conditions disease (, –78). Enrichment of the gene lists was calculated by overrepresentation analysis. (C) Alignment of the populations identified by IMC and scRNA-seq. (D) Schematic showing the strategy for deconvolving the Ravi et al. (37) ST datasets using the cell2location (82) and the transcriptomic signatures of the myeloid populations we identified in the Ruiz-Moreno et al. dataset (26). (E) Predicted abundances of the myeloid populations calculated using cell2location. Data shown as mean ± SEM. (F) Distribution of myeloid populations in ST spots within a representative deconvolved ST case. (G) Pearson’s R correlation between the different myeloid populations present in each spot. (H) Strategy for identifying explanatory factors that may control the positioning of myeloid cell populations. Using this strategy, myeloid cell abundances were correlated with the transcriptomic signatures of biological processes from the MSigDB database (I) or expression of chemokines (J). *P < 0.05, with correction for multiple Pearson’s R tests using the Benjamini-Hochberg procedure (G, I, and J). UV, ultraviolet; NF-κB, nuclear factor κB; TNFα, tumor necrosis factor–α; ROS, reactive oxygen species.

Fig. 6.. Spatial clustering of myeloid cells is associated with poor outcome in GBM.

(A) Strategy to identify environments in the deconvolved ST cases that have distinct patterns of myeloid cell distribution. (B) Identification of distinct patterns of myeloid cell distribution using k-means clustering, corresponding to five distinct myeloid cell environments. (C) Mapping of the five myeloid environments in example ST cases. (D) Proportion of interactions between spots from each myeloid environment, assuming that each spot is connected to its neighboring six spots. (E) Strategy to identify the transcriptomic changes arising from the nonmyeloid cells in each spot (see Materials and Methods). (F) Overrepresentation analysis of the signatures of biological process from the MSigDB database (79) in the remaining nonmyeloid genes in the different myeloid environments. (G) The top five enriched terms from the Gene Ontology (GO) database in the nonmyeloid genes in the different myeloid environments. (H) Strategy to deconvolve bulk RNA-seq GBM cases from the TCGA (The Cancer Genome Atlas) and CCGA (Chinese Glioma Genome Atlas) using the TAPE algorithm (–43), therefore allowing us to predict the proportions of myeloid environments in each case. (I) Proportion of myeloid environments predicted by the TAPE algorithm in the TCGA and CCGA cases. (J) Modeling the relationship between the abundance of each myeloid environment and GBM survival using Cox proportional hazards, with correction for multiple tests using Holm-Šídák. Hazard ratios are increase in risk of death per percentage point of myeloid environment abundance. Kaplan-Meier curves compare patients having the high (top 50%) or low estimates for the presence of that environment, and shaded areas are 95% confidence intervals.