Tumor Microenvironment Is Critical for the Maintenance of Cellular States Found in Primary Glioblastomas

Abstract

Glioblastoma (GBM), an incurable tumor, remains difficult to model and more importantly to treat due to its genetic/epigenetic heterogeneity and plasticity across cellular states. The ability of current tumor models to recapitulate the cellular states found in primary tumors remains unexplored. To address this issue, we compared single-cell RNA sequencing of tumor cells from 5 patients across four patient-specific glioblastoma stem cell (GSC)-derived model types, including glioma spheres, tumor organoids, glioblastoma cerebral organoids (GLICO), and patient-derived xenografts. We find that GSCs within the GLICO model are enriched for a neural progenitor-like cell subpopulation and recapitulate the cellular states and their plasticity found in the corresponding primary parental tumors. These data demonstrate how the contribution of a neuroanatomically accurate human microenvironment is critical and sufficient for recapitulating the cellular states found in human primary GBMs, a principle that may likely apply to other tumor models. SIGNIFICANCE: It has been unclear how well different patient-derived GBM models are able to recreate the full heterogeneity of primary tumors. Here, we provide a complete transcriptomic characterization of the major model types. We show that the microenvironment is crucial for recapitulating GSC cellular states, highlighting the importance of tumor-host cell interactions.See related commentary by Luo and Weiss, p. 907.This article is highlighted in the In This Issue feature, p. 890.

Figure 1.. Glioblastoma cerebral organoids better recapitulate GBM bulk patient tumors than glioma spheres or tumor organoids.

a. Schematic overview of model generation. b. Representative histology of the glioma cerebral organoid and xenograft models. Scale bar: 50um. c. Representative histology of the patient sample (for all H&E images, see Supplementary Figure 1). Scale bar: 50um. d. Combined principal component analysis of each model type with the bulk sample of the patient tumor. e. Unsupervised hierarchical clustering (Ward’s method, Euclidean distance) of the expression signatures of all genes in the bulk patient tumor with the average expression values of each available model type. f. Distributions of the Spearman correlation of each model type with the bulk patient tumor using the expression signatures of all genes in the bulk sample (**** indicates GLICO is significantly greater than the other three models using one-way ANOVA test with Tukey’s multiple comparisons, P < 0.0001. Minimum dispersion cutoff of 0.5. Using highly variable genes also results in GLICO with a significantly higher correlation, P < 0.0001). g. Hierarchical clustering (Ward’s method, Cosine distance) of the average transcriptomes of high variation GSCs for each model type and published primary scRNA-seq (14).

Figure 2.. Glioblastoma cerebral organoids are enriched for a proneural subpopulation present in primary GBM.

a. UMAP projection of batch-corrected cells across all samples and models colored by Louvain clustering assignments, and additional UMAP projections of well-known genes. b. Expression of top marker genes defining Louvain clusters as determined by Wilcoxon test of each gene between each cluster versus all other clusters. c. Heatmap and hierarchical clustering of Louvain clusters by relevant gene set z-scores. d. Pie charts depicting percentage of cells in each cluster by model. Legend keys are functional group name followed by cluster numbers in that group. e. Percentage of cells belonging to each model type by GBM subtype. Subtypes were computed by correlation to TCGA reference data of Classical, Mesenchymal, and Proneural subtypes. Error bars represent one standard deviation.

Figure 3.. Glioblastoma cerebral organoids are enriched for NPC-like and OPC-like cellular states.

a. Stemness, oligodendrocyte, and astrocyte scores based on curated lineage gene sets(5). b. Expression of marker genes in each of the three subpopulations. c. Box plots of cell type scores by model type. Significance was determined using one-way ANOVA (*P < 0.05; ***P < 0.001). d. UMAP projection colored by cell assignments using Neftel gene sets(7). Cells were assigned to the cell type that had the maximum score within that cell above a minimum evidence threshold of 0.3. e. Relative meta-module scores for cellular states plot(7). Axes are the log1p transformation of the absolute difference between cell type scores where each quadrant represents an OPC, AC, NPC, or MES-like cell state. Colors indicate density of cycling cells. f. Normalized confusion matrix between our Louvain cluster cell assignments and cellular state assignments(7) as in Figure 3D. g. UMAP projection of batch-corrected cells with addition of eight primary GBM samples(14) colored by (left) unsupervised Louvain clustering assignments and (right) condition, subset to GLICO, XE, and Primary tumor cells. Arrow indicates proneural cluster. h. Bar chart depicting percentage of cells in the NPC/OPC cluster (C3 in Figure 3G) by model (*P < 2.2 × 10⁻¹⁶, hypergeometric enrichment test). Error bars represent one standard deviation.

Figure 4.. Glioblastoma cerebral organoids express Notch pathway members and GBM invasiveness markers.

a. Volcano plot of GLICO model vs. 2D and TO models. Genes with an absolute log fold change > 0.6 are labeled. b. Gene set enrichment analysis of genes differentially expressed in GLICO from Figure 4A. Showing gene sets for proneural signature (P < 2.2×10⁻¹⁶, FDR < 2.2×10⁻¹⁶) and spheroid vs. adherent growth (P < 2.2×10⁻¹⁶, FDR=0.08). c. Expression of SOX4, BCAN, and DLL3 across models (* indicate p < 0.05 compared to all other models). d. Volcano plot of subsequent pathway analysis of MAST differential expression results in Figure 4A. Pathways with an activation Z-score > 0 are considered activated in GLICO and an activation Z-score < 0 are considered activated in 2D and TO models. Size of points proportional to the number of molecules supporting the annotation. e. Immunofluorescence staining of BCAN, DLL3, and SOX4 in GSCs isolated from the GLICO model (Scale bar 100um; inset: Scale bar 10um). f. UMAP projection of cell type assignment and cell velocity streams as determined by RNA Velocity(30) in individual samples. g. Differential expression between RNA Velocity “root” cells and all other cells.

Figure 5.. Cerebral organoid influences GBM cellular states.

a. Heatmap depicting genes differentially expressed between GLICO and 2D, showing log-normalized expression of genes in GLICO, 2D, and the re-plating GLICO into 2D coculture experiment (CC). Genes are clustered using hierarchical clustering. Representative genes are labeled. b. Principal component analysis of patient samples cultured in GLICO and then re-plated in 2D for 30 days (CC), as compared to GLICO and 2D samples, colored by condition, cell cycle assignment, and NPC or AC-like cell type score using Neftel gene sets(7). c. Representative fluorescence images of GFP+ GSCs (1.) in GLICO, (2.) 7 days after and (3.) 30 days after the re-plating experiment from GLICO to 2D. d. Schematic of re-plating experiment results.