Multiplexed imaging uncovers correlations between genotype and TME.

In our previous study, we showed that the single-cell mosaicism of genetic amplifications in GBM correlates with immune infiltration (5). To test whether we could identify features of the TME that drive the selection of specific amplifications, we conducted a multiplexed imaging study on formalin-fixed, paraffin-embedded (FFPE) samples from a cohort of 9 matched primary and recurrent GBM cases (Figure 1A and Supplemental Table 1; supplemental material available online with this article; https://doi.org/10.1172/jci.insight.179853DS1). The imaging was performed in 2 consecutive rounds on the same tissue, using modified CyCIF and STAR-FISH protocols (5–7) (see Methods for details). First-round imaging was aimed to acquire information about the TME state by assessing the presence of blood vessels (CD31), active infiltrating immune cells (CD45RO), and hypoxia (HIF1α) (Figure 1B). The next layer of imaging was focused on genetic features, including amplifications of EGFR, CDK4, and PDGFRA, and hotspot mutation in the TERT promoter (TERTp) (Figure 1C). The 3 amplifications were previously associated with distinct transcriptional states in GBM, namely NPC-like state with CDK4 amplification, AC-like state with EGFR amplification, and OPC-like state with PDGFRA amplification (2). Nuclear staining channels from the first and second rounds of imaging were overlaid to ensure no loss of cells occurred (loss of a single cell occurred in 9 imaged regions and between 2 and 6 cells in 3 imaged regions; total count of lost cells was 22 out of 20,205; Supplemental Figure 1). Image segmentation and quantification of both the nuclear and TME-related staining allowed us to classify each nucleus as belonging to a tumor cell with or without one or a combination of the 3 amplifications, a tumor cell with TERTp mutation, an immune cell, or an endothelial cell (EC) (Figure 1, B–E, and Supplemental Table 2). GBM cells can mimic ECs and pericytes (8, 9). Indeed, we identified 958 cells characterized by tumor-specific genetic markers and expressing the EC/pericyte marker CD31. The hypoxia marker, nuclear HIF1α, can be expressed both by tumor cells and immune cells in the TME and our analysis identifies both populations. Of note, cells harboring normal CDK4, EGFR, and PDGFRA copy number, no mutation, and no immune or endothelial markers, which may represent normal cells such as oligodendrocytes, neurons, or astrocytes, as well as tumor cells negative for our tumor markers, were excluded from quantitative analyses and only considered in spatial analyses for accurate measurement of regional cell density.

Immunogenotyping of primary and recurrent archival GBM samples. (A) Study outline. (B) Top panels: Representative images of immunofluorescent staining for markers of hypoxia (HIF1α), endothelial cells (CD31), and immune cells (CD45RO). Bottom panels: Image segmentation. (C) Top panels: Representative FISH images for CDK4, EGFR, and PDGFRA and STAR-FISH for TERT promoter (TERTp) mutation, respectively. Bottom panels: Image segmentation. Scale bars: 40 μm (B and C). (D) Left panels: Quantification of cell frequency based on phenotypes (top) and genotypes (bottom) in representative images in B and C. Right panels: Spatial distribution of cells classified into distinct phenotypes and genotypes corresponding to the genotype panels on the left. (E) Top panel: Quantification of cell frequency based on TERTp mutation status in representative images in B and C. Bottom panel: Spatial distribution of cells classified into TERTp WT and MUT.

First, we searched for changes in the frequencies of genotypes and phenotypes between matched primary and recurrent cases (Figure 2, A and B, and Supplemental Table 3). The frequency of cells with CDK4 amplification and cells with TERTp mutation was higher in recurrent samples (Figure 2, C and D; for CDK4 comparisons; 2-tailed paired Wilcoxon’s ranked test, P = 0.03; for TERTp mutant comparisons: 2-tailed paired Wilcoxon’s ranked test, P = 0.04). Among the phenotypes we identified, hypoxic EC mimicry was more prevalent in recurrent samples (Figure 2E; 2-tailed paired Wilcoxon’s ranked test, P = 0.02). To account for overall diversity of each imaged field of view, we calculated the Shannon index of diversity (see Methods), capturing the evenness of distribution of cells among the different genotypes or phenotypes in each tumor (Supplemental Figure 2). These indices were significantly different between primary and recurrent samples, when amplification or TERTp mutation status was taken into consideration (2-tailed paired Wilcoxon’s ranked test, P = 0.0023 and P value = 0.0025, respectively). Thus, the distribution of cells harboring TERTp mutation within each tumor is more heterogeneous than the variation in hypoxia, vasculature, or immune infiltration.

Differential cell frequencies and correlations between genotypes and phenotypes in primary and recurrent tumors. (A) Frequency of cells with distinct genotypes in primary and recurrent GBM for each of 5 imaged tumor regions per case (n = 9 matched cases). (B) Frequency of cells with distinct phenotypes in primary and recurrent GBM for each of 5 imaged tumor regions per case (n = 9 matched cases). (C) Frequency of cells with CDK4 amplification in matched primary and recurrent samples. Data points represent case average frequency after ROUT outlier removal. Dotted line: matched primary and recurrent cases. Violin plot shows mean and quartiles. P value of 2-tailed paired Wilcoxon’s ranked test is shown. (D) Frequency of cells with TERT promoter (TERTp) mutation in matched primary and recurrent samples. Data points represent case average frequency after ROUT outlier removal. Dotted line: matched primary and recurrent cases. Violin plot shows mean and quartiles. P value of 2-tailed paired Wilcoxon’s ranked test is shown. (E) Frequency of cells with TERTp mutation in matched primary and recurrent samples. Data points represent case average frequency after ROUT outlier removal. Dotted line: matched primary and recurrent cases. Violin plot shows mean and quartiles. P value of 2-tailed paired Wilcoxon’s ranked test is shown. (F) Pearson’s correlation between the frequencies of cells with distinct genotypes and phenotypes.

Next, we calculated Pearson’s correlation between the frequencies of cells with distinct genotypes and phenotypes (Figure 2F). The most notable difference between primary and recurrent tumor samples was the association between hypoxic tumor cells adopting an EC phenotype (hypoxic EC mimicry) and a genotype consisting of all 3 amplifications co-occurring in the same cell. This correlation was significant, but weak in primary tumor samples, with r = 0.3 and P = 0.041, and strengthening in recurrence, with r = 0.73 and P = 0.9 × 10–8, suggesting that co-amplification of EGFR, CDK4, and PDGFRA may allow more tumor cells to adopt this hypoxia-driven perivascular phenotype.

Previously, we showed that the relative frequency of co-amplification of EGFR and CDK4 at the single-cell level is linked with an immunosuppressive microenvironment in GBM (5). To test whether this association holds in recurrent tumors, we calculated the odds ratio (OR) for co-amplification of EGFR and CDK4 for each tumor. The classification of our tumors into ORhigh and ORlow revealed a significant difference in overall immune infiltration only in primary tumors, but not in the recurrent samples (Supplemental Figure 3). Thus, it is likely that these genetic drivers are more important in the early steps of establishing the immunosuppressive microenvironment.

To better understand the relationships between the genotypic and phenotypic diversity in our cohort, we next performed clustering analysis based on the frequencies of all cell types identified in each imaged region of interest. Clustering of individual images based on TERTp mutation frequencies revealed 4 classes of tumor areas. Similarly, amplification-based genotype clustering revealed 4 clusters of tumor areas, while phenotype frequencies divided the imaged tumor regions into 6 clusters (Figure 3, A and B). Most notable were the phenotype cluster 3, characterized by high frequency of immune cells, and cluster 4, with highest abundance of hypoxic cells (Supplemental Figure 4). Phenotype clusters 0 and 1 are both composed of intermediate frequencies of immune cells and ECs, but differ in overall cellularity (Supplemental Figure 4A). Interestingly, in addition to the highly hypoxic cluster 4, it was the more complex phenotypes of clusters 0 and 1 that were significantly associated with particular genotypes (Figure 3C). Tumor regions in phenotype cluster 0 are largely classified as genotype cluster 0 and mutation cluster 1, indicating a low frequency of all measured amplifications and depletion of TERTp-mutant cells despite overall high cellularity of these tissues (Figure 3, B and C, and Supplemental Figure 4, A–C). Tumor regions in phenotype cluster 1 are mainly comprised of genotype cluster 1 and mutation cluster 2, indicating an intermediate level of singly amplified EGFR and CDK4 cells, with a slightly higher frequency of TERTp-mutant cells compared with phenotype cluster 0 (Figure 3, B and C, and Supplemental Figure 4, A–C). Given the high cellularity of phenotype cluster 0, it is plausible that the majority of these tumor regions are enriched in cancer cells harboring genetic changes not captured by our probes. Cells with EGFR-only amplification were depleted from regions rich in hypoxic cells, which was not observed for EGFR/CDK4 co-amplified cells (Figure 3, D and E). This is in line with the previously reported angiogenic role of amplified EGFR signaling (10).

Genotype- and phenotype-based clustering of tumor microenvironments in primary and recurrent GBM. (A) Tumor region clustering based on cellular phenotypes, genotypes, and TERT promoter–mutant (TERTp-mutant) cells. Each point represents a tumor region. Numbers represent cluster identifier. (B) Connections between classification based on TERTp mutation, genotype, or phenotype clustering in primary and recurrent GBM. The width of each connection represents the number of tumor regions classified. (C) Contingency between the phenotypes and genotypes. Fisher’s exact test P values are: Cluster 0, P < 0.0001; Cluster 1, P < 0.0001; Cluster 4, P = 0.009. The color scale represents number of tumor regions classified. (D) Frequency of cells with EGFR amplification (left) and CDK4/EGFR co-amplification (right) across phenotype clusters. (E) Frequency of hypoxic (left) and immune cells (right) across phenotype clusters. The box-and-whisker plots in all bar graphs show the mean (midline) and 25th–75th (box) and 5th–95th (whiskers) percentiles. (F) Spatial distribution of cells with different genotypes and phenotypes. Columns in the heatmap represent spatial clusters determined based on XY coordinates of the cells.

We also noted that the connections between genotype- and phenotype-based classifications change between primary and recurrent tumors (Figure 3B). In recurrent tumors, genotype cluster 3 environments, with high frequencies of both EGFR-only amplified and EGFR/CDK4 co-amplified cells, lose their connection with hypoxic (phenotype cluster 4) and tumor-cell-rich nonvascularized immunodepleted environments (phenotype cluster 5). We instead observed a transition of genotype cluster 3 environments to tumor-cell-rich neighborhoods, with increased frequency of EC mimicry and intermediate levels of hypoxia and immune infiltration (phenotype cluster 2). Genotype cluster 2, rich in CDK4-amplified cells, associated with hypoxic and immune-rich environments in primary tumors switched to a tumor-enriched environment in recurrent tumors. These frequency-based analyses confirm that genomic amplification may constrict the evolvability of GBM cells and may impact the adaptation of the TME upon recurrence.

Since intercellular interactions depend on cellular proximity to their neighbors, we next performed clustering analysis considering the spatial localization of the cells analyzed in this study. We built a distance matrix and neighborhood feature vectors based on the Euclidian distances from each cell to its nearest neighbors for each tumor region analyzed in our study. Comparison of proximity-based clustering shows that spatial arrangement preferences for immune cells and vascular cells are comparable between primary and recurrent tumors (Figure 3F). However, CDK4 single-amplified cells have more genetically diverse neighbors in recurrent tumors than in primary tumors. Given that hypoxia is reported as contributing to the processes of EC mimicry and vasculogenic transformation, it was interesting to observe hypoxic EC-mimicking cells that were less abundant in primary tumors overall, possessing a higher diversity of neighbors in this setting. In recurrent tumors, hypoxic EC-mimicking cells increase in number but are primarily in proximity to other tumor cells and especially those with TERTp mutation. Hypoxia levels increasing over the course of disease may lead to an increased adaptation of malignant cells by adopting the EC phenotype. Initial stochasticity of this process could result in a more dispersed spatial arrangement of the EC-mimicking cells within the TME.

In summary, our results demonstrate that a significant shift in cellular composition from primary to recurrent GBM is reflected in single-cell genetic heterogeneity, but is also linked to tumor cells’ ability to adopt an EC/pericyte-like phenotype. This cellular adaptation is associated with local hypoxia and co-occurrence of all 3 oncogenic amplifications driving the distinct cellular states in GBM (2). It is plausible that the high level of aneuploidy does not allow these cells to thrive in a heterogeneous tumor, but provides a selective advantage under the pressure of chemotherapy and radiation.

Texture of ECM revealed by reflectance imaging associated with genetic diversity.

The differential composition of the ECM, with elevated levels of collagens, laminins, and hyaluronan, is another factor contributing to local microenvironment heterogeneity within GBM tissue (11). To image ECM structures without the need to add additional fluorophores to our staining panel, we took advantage of the differential refractive index of ECM components, which results in the scattering of light off the ECM structures generating a label-free contrast image by reflectance confocal microscopy (RCM) (12). Thus, we included RCM imaging to visualize the organization of the ECM in the first round of our imaging protocol. Across 90 imaged tumor areas, we identified 3 classes of ECM texture: “string-like,” “spotty,” and “mossy” (Figure 4A), with a majority of the fields of view (FOVs) containing a mixed texture (Figure 4B and Supplemental Table 3).

Reflectance imaging–based patterns of tissue organization predictive of cellular composition of the tumor. (A) Representative images of reflectance microscopy–based extracellular matrix textures: string-like, spotty, and mossy. Scale bars: 40 μm. (B) Representative images of mixed reflectance textures and their manual annotation. SL, string-like texture; SP, spotty texture; MO, mossy texture. Scale bars: 40 μm. (C) Influence of cell type on reflectance niche. Mean non-zero influences estimated from correctly classified quadrats (quadrat width and height = 199), weighted by the frequency at which each feature had a non-zero influence per class. Error bars indicate the 95% confidence intervals. n reps = 100, mean correctly classified quadrats per rep = 690.49, median ROC-AUC = 0.784.

Species distribution models (SDMs) are an ecological method to measure how environmental factors and species are spatially associated with a niche of interest and have previously been used to study the TME (13). We thus developed a neural network–based SDM to better understand the spatial relationship between cellular composition and the local ECM texture. The manual annotation of ECM texture was used to perform quadrat counting on each FOV’s cell segmentation data (see Methods for details). The network was then trained to predict a quadrat’s ECM texture given the abundance of each cell type in that quadrat. After training, the relationship between the abundance of each cell type and the ECM was quantified using the Integrated Gradients feature attribution method (14). This approach allows for the exclusion of incorrectly classified quadrats, meaning that the estimates of the relationship between cell types and the ECM are only based on correct classifications. After applying the Integrated Gradients method to 100 trained models, each cell type’s average attribution was calculated and multiplied by the frequency with which that cell type had an informative attribution (i.e., associated with a correct classification and had non-zero attribution; Supplemental Figure 5). The weighted non-zero attributions for each ECM texture show significantly different contributions from different cell types (Figure 4C). String-like texture is strongly associated with the presence of CDK4-amplified cells, EGFR-amplified cells, and TERTp-mutant cells, yet not associated with overall frequency of tumor cells (all cells with amplifications or mutations). In contrast, spotty ECM texture is linked to higher overall frequency of tumor cells and mildly related to presence of cells with co-amplified PDGFRA/EGFR. These results suggest that cells with PDGFRA/EGFR co-amplification may have distinct microenvironmental niche preferences compared with cells with other genotypes identified in this study. Thus, despite its heterogeneous origin, we show that reflectance can operate as a meaningful tool to classify TMEs.

Reflectance and spatial transcriptomics identify a cluster of perivascular immunosuppressive macrophages.

While classifying the reflectance images, we noted that the majority of the blood vessels in recurrent tumors are surrounded by a thick layer of ECM deposition (Figure 5, A–C). Vascular malformations, including vessel hyalinization, are frequent side effects of radiotherapy (15). Trichrome staining confirmed that these structures are highly enriched in collagen (Figure 5B). Interestingly, immunofluorescent staining revealed that the perivascular collagen rims are densely populated by immune cells, and that a vast majority of these cells are CD163+ (Figure 5A). CD163 is a marker of immunosuppressive polarization of tumor-associated macrophages, previously linked to GBM survival (16). The tight localization of these cells and protein around blood vessels possibly impede other immune players like CD8+ cells (Figure 5A) and hematogenous components from exiting the vasculature and could also have a supporting function for blood vessel structure.

Spatial profiling of perivascular regions in primary and recurrent GBM. (A) Representative images highlighting selected vascular regions of recurrent GBM samples. Reflectance imaging showing hyperintense signal around vessel lumen (left) and immunostaining for CD163 and CD8α of peripheral cells (right). Scale bars: 40 μm. (B) Trichrome staining of perivascular collagen in recurrent GBM of adjacent section of the tumor region imaged in panel A. Scale bar: 40 μm. (C) Quantification of collagen rim around blood vessels in FOVs containing vasculature for primary and recurrent samples. Color scale represents number of cases. (D) Semisupervised clustering generated from normalized gene expression of 52,588 cells (4 GBM tumors). (E) Top differentially expressed markers of cellular groups. Columns represent cell types and rows represent genes. Scaled expression data represented as z scores. (F) UMAP plot grouped by tissue of origin. Primary (top) or recurrent (bottom) tissues. (G) Cell type composition of primary and recurrent samples. (H) Representative images of linked IHC (top), cell type spatial plots (middle), and niche spatial plots (bottom) in primary and recurrent samples. IHC images represent matched tissue sample locations to spatial plots at serial section not more than 12 μm away. Scale bars: 120 μm. (I) Distribution of percentage of cell types present within TME niches. Yellow cells indicate highly represented outliers computed at α = 0.001 before normalization. (J) Niche composition of primary and recurrent tumor.

To test these hypotheses and elucidate the functional role of the perivascular immunosuppressive macrophages in recurrent GBM, we performed single-cell spatial transcriptomic characterization of 2 matched primary and recurrent human GBM samples using the CosMx platform (17). A total of 36 regions of interest containing blood vessels and surrounding tissue were selected (Supplemental Figure 6). CosMx enabled localization of 1,030 transcripts at single-cell resolution by performing multiple cycles of nucleic acid hybridization of fluorescent molecular barcodes on FFPE tissue. Our analysis was performed on 52,588 cells, with an average of 106 transcripts detected per cell (Supplemental Tables 4 and 5). Semisupervised Leiden clustering of both primary and recurrent samples identified 10 major cellular clusters, annotated based on established markers of cell type metaprograms (2) (Figure 5, D and E). Recurrent samples were enriched in mesenchymal-like tumor cells (“MES”), expressing high levels of NDRG1 and VEGFA that further distinguish them as hypoxic (Figure 5, F and G, and Supplemental Figure 7). We observed these MES cells arranging in spatial niches surrounding collagen-rimmed vessels (Figure 5H).

Of note, 2 different types of immune cells were associated with the blood vessels in recurrent tumor tissue: the perivascular macrophages (“PVM”), expressing CD68, CD14, CD74, and CD163, corresponding with CD163+ immunostaining, and the SPP1-, CCL7-, and CEACAM1-expressing cells (“Immune”) — most likely tumor-infiltrating lymphocytes, which would correspond to the spatial positioning of CD8a+ immunostaining (Figure 5, A, E, and H). While PVM cells do not appear to be the major source of collagen in the perivascular spaces (Supplemental Figure 7), they are in closest proximity to the fibroblast-like cell type that highly expresses COL1A1, COL3A1, and COL6A3. Notably, this fibroblast-like cell group also highly expresses DCN, encoding decorin, an ECM proteoglycan that binds collagen (Figure 5E).

Next, we performed niche analysis combining all the profiled tissue areas. Seven niches were identified, with variable abundance in primary and recurrent samples (Figure 5, H–J). The most interesting was the organization of the niches surrounding the blood vessels in recurrent tumors. In these areas, niche 1, which is a vascular niche composed of ECs and fibroblasts, was surrounded by niche 7, enriched with immune cells, glial cells, and fibroblasts (Figure 5I). PVM cells of recurrent tumors were primarily found in this zone. These areas were also encapsulated by niches rich in mesenchymal and OPC-like cells (Figure 5H). No such structured organization was observed around vasculature of primary tumor samples. Rather, niche 2, which is enriched with a different macrophage cell type (“MAC”) was observed to be more localized around vasculature and widely distributed in the TME of primary tumor samples (Figure 5, H–J). Surprisingly, niche 7, associated with perivascular regions in recurrence, was only identified in primary tumor areas possessing aggregated blood cells and PVM surrounded by glial cells, suggesting ruptured vessels with hemorrhage and glial scarring (18) (Figure 5H).

Since our initial unsupervised analysis identified 4 clusters of macrophages, 2 of which seemed to have the PVM phenotype, we asked whether this cluster could be further resolved to differentiate true perivascular macrophages from the ones linked to hemorrhagic areas in the tissue. Indeed, we found that these 4 clusters correspond to 4 spatially distinct locations (Figure 6, A–D; reanalysis of data presented in Figure 5). Cluster 5 macrophages are widely distributed in the tissue (Figure 6C). Clusters 8 and 9 are both associated with vasculature, but only cluster 9 cells are found next to vessels with collagenous deposition in recurrent samples (Figure 6C). Interestingly, cluster 9 cells express significantly higher levels of collagen-encoding COL1A1 and COL3A1, as well as fibronectin (FN1) and thus may contribute to the thickening of the perivascular ECM (Figure 6D). Cluster 10 macrophages are predominantly residing within the hemorrhagic areas (Figure 6C). The differences in perivascular niche composition between primary and recurrent tumor were also evident in cell-cell proximity analysis (Figure 6E). In primary tumor, ECs interact directly with fibroblasts and macrophage cluster 8, while in recurrent tumor there is enrichment of interactions between the ECs, fibroblasts, and macrophage clusters 8, 9, and 10. Thus, while transcript-based identification of distinct phenotypes of macrophages might be challenging with non–genome-wide spatial transcriptomics, the analysis of the localization of suspected cell types yields clear spatial differentiation between primary and recurrent GBM.

Perivascular niche macrophage cellular interactions in vessels with or without collagen rim. (A) UMAP plot highlighting cell type clusters 5, 8, 9, and 10 (higher resolution clustering of previous MAC and PVM cell type clusters). (B) Feature plots depicting macrophage/monocyte-derived cell type gene expression markers enriched in the areas of cluster 5, 8, 9, and 10. (C) Spatial plots of cell types from cluster 5, 8, 9, and 10 in primary and recurrent GBM. FOVs chosen possess vasculature validated through IHC and gene expression. Scale bar: 120 μm. (D) Top differentially expressed markers between clusters. Columns represent cell types and rows are genes. Scaled expression data represented as z scores. (E) Cell-cell interactions in primary and recurrent tumor tissue. Green lines show spatial proximity enrichment and red lines show depletion between pairs of cell types. Proximity enrichment derived by calculating the observed over the expected frequency of cell-cell proximity interactions. The expected frequency is the average frequency calculated from the spatial network simulations.