Spatial transcriptomics reveals niche-specific enrichment and vulnerabilities of radial glial stem-like cells in malignant gliomas

Abstract

Diffuse midline glioma-H3K27M mutant (DMG) and glioblastoma (GBM) are the most lethal brain tumors that primarily occur in pediatric and adult patients, respectively. Both tumors exhibit significant heterogeneity, shaped by distinct genetic/epigenetic drivers, transcriptional programs including RNA splicing, and microenvironmental cues in glioma niches. However, the spatial organization of cellular states and niche-specific regulatory programs remain to be investigated. Here, we perform a spatial profiling of DMG and GBM combining short- and long-read spatial transcriptomics, and single-cell transcriptomic datasets. We identify clinically relevant transcriptional programs, RNA isoform diversity, and multi-cellular ecosystems across different glioma niches. We find that while the tumor core enriches for oligodendrocyte precursor-like cells, radial glial stem-like (RG-like) cells are enriched in the neuron-rich invasive niche in both DMG and GBM. Further, we identify niche-specific regulatory programs for RG-like cells, and functionally confirm that FAM20C mediates invasive growth of RG-like cells in a neuron-rich microenvironment in a human neural stem cell derived orthotopic DMG model. Together, our results provide a blueprint for understanding the spatial architecture and niche-specific vulnerabilities of DMG and GBM.

Fig. 1. Overview of workflow and sample information.

a Illustration of the overall workflow. Examples images for SpotClean-adjusted gene expression (OLIG2) and BANKSY-based spatial clustering for DMG1 are shown. b UMAP of all Harmony-integrated spatial transcriptomic spots (n = 26, 460) from all samples (n = 11). Colors represent different samples. Clusters of peritumor spots are highlighted by blue dashed lines. c Overview of the workflow for CNV-based estimation of tumor cell content. GBM4 is shown as an example. Colors indicate the predicted tumor cell content score. d UMAP of all Harmony-integrated spatial transcriptomic spots from all samples (n = 26,460 spots). Colors indicate the predicted tumor cell content score. The same clusters of peritumor spots in b are highlighted by blue dashed lines. e The colocalization of CNV and H3K27M mutation in spots from DMG samples (n = 5). The number of spots from each sample is indicated.

Fig. 2. Niche-specific spatial transcriptional programs and their relationship to subclonal architecture.

a Heatmap of horizontally integrated spatial transcriptional signatures/modules expressed in all malignant spots (n = 19,767). The sample identify for each spot is shown by different colors in the top bar. b Heatmap of GWR-based, spatially weighted correlation of the 48 spatial cluster signatures, marked by different colors on the left. Hierarchical clustering confirmed 4 clustered programs that correspond to the four modules in a, marked by different colors on the top. c The percentages of cluster signatures per tumor type contributing to each program are shown. d Representative examples of the correlation between the module gene expression and the histopathology of glioma niches across all samples (n = 10). e Heatmap of spatial weighted correlation between our modules and Ravi et al. modules across all malignant spots (n = 19,767). Correlation and P values were determined by non-parametric Spearman’s correlation test, Bonferroni-corrected for multiple comparisons. Source data are provided as a Source Data file, and individual P values are included. f Dot plot indicating subclonal distribution across different glioma niches marked by transcriptional modules from all subclones (n = 24). Dot size and color reflect the percentage.

Fig. 3. RNA isoform diversity across glioma niches revealed by long-read spatial transcriptomics.

a Workflow to identify niche-specific isoforms and splicing junctions (SJs), predict regulatory splicing factors (SFs), and analyze their prognostic value for patient survival. b Gene Ontology (GO) enrichment of the host genes of niche-specific isoform across samples (n = 10). P values were determined by two tailed hypergeometric test. Significantly enriched (Benjamini–Hochberg-adjusted P < 0.05) biological processes are shown. c–e Differential isoform enrichment of CHI3L1 in different niches. c Individual transcripts of CHI3L1 isoforms 201 and 205 in high-quality malignant spots from all glioma samples (n = 19,071 spots) are represented by transcript tracks, colored by their niche identities. Transcript structures based on Ensembl annotation are shown at the bottom, with colored regions representing exons and gray regions representing introns. d Tile plots showing scaled mean expression of CHI3L1-201 and CHI3L1-205 across niches to compare their niche-specific enrichment patterns. e As an example, spatial plots of PSI values (colors) in DMG3 quantifying the relative expression of CHI3L1-201 and CHI3L1-205 are shown. f Venn diagram showing the overlap of splice junctions detected in our long reads sequencing dataset and the TCGA-GBM dataset. Differential isoform enrichment of TOMM6-201 and TOMM6-202 in different niches of GBM samples (n = 6727 spots), shown by tile plots (g) and individual transcripts in each niche (h). i Kaplan–Meier survival curves comparing the overall survival in the TCGA-GBM cohort^,, stratified by high (n = 62) vs. low expression (n = 102) of TOMM6-202 specific splice junction. Curve comparison P value was determined by two tailed Log-rank (Mantel–Cox) test. j Predicted splicing regulation of TOMM6 by splicing factor YBX3. Peaks indicate eCLIP and shRNA-KD RNA-seq read density in HepG2 cell line from ENCODE. Biological replicates have similar results. Gray box shows the alternatively spliced region of TOMM6-202. k Kaplan–Meier survival curves comparing the overall survival between YBX3 high (n = 17) versus low groups (n = 149) in the TCGA-GBM cohort. Curve comparison P value was determined by two tailed Log-rank (Mantel–Cox) test.

Fig. 4. Profiling ecosystems in glioma niches.

a The average of the deconvoluted percentage of individual cell types in spots from each glioma niche (n = 19,767). Integrated GBM, DIPG, CTX, and CB public scRNA-seq datasets were used as a reference in the deconvolution analysis. b Dot plot of deconvoluted percentage of neurons in individual malignant spots from different glioma niches, quantified in all tumor samples or by tumor type (n = 10,342, 1182, 6290, and 1951 spots; n = 3917, 211, 1796, and 774 spots; n = 3713, 48, 2533, and 347 spots; n = 2712, 923, 1961, and 830 spots, respectively). Boxes indicate quartiles, horizontal bar indicates median, and whiskers indicate range, up to 1.5-fold inter-quartile range. P values, Kruskal–Wallis test for between group differences with Holm’s correction for multiple comparisons. c Dot plot of deconvoluted percentage of RG-like cells in individual malignant spots from different glioma niches, quantified in all tumor samples or by tumor type (n = 10,342, 1182, 6290, and 1951 spots; n = 3917, 211, 1796, and 774 spots; n = 3713, 48, 2533, and 347 spots; n = 2712, 923, 1961, and 830 spots, respectively). Boxes indicate quartiles, horizontal bar indicates median, and whiskers indicate range, up to 1.5-fold inter-quartile range. P values, Kruskal–Wallis test for between group differences with Holm’s correction for multiple comparisons. d Re-analysis of Filbin et al. DIPG scRNA-seq dataset presented by UMAP scatter plot. Cell identities for each cluster are consistent with the original study. The red-dashed circle highlights the malignant cell populations. OC_N, normal oligodendrocytes. e Violin plots of the relative RG score expression in individual cell types in d. n = 206, 87, 1437, 480, 49, 94, and 96 cells, respectively. Two tailed Wilcox test for between group differences for comparisons. f Violin plots of the relative RG-Vs-AC score expression in AC1 (n = 206 cells) and AC2 (n = 87 cells). Two tailed Wilcox test for between group differences for comparisons. g Re-analysis of Neftel et al. GBM scRNA-seq dataset in scatter plot showing the expression of RG score. The cellular states for each single cell were obtained from the annotations by Neftel et al.

Fig. 5. Identification of regulatory programs of RG-like cells in a neuron-rich invasive niche.

a Selected spots from the tumor core (Tumor), leading edge (LE), and tumor-infiltrated GCL (GCL_TI), and normal spots without CNV from GCL (GCL_N) of DMG1 are shown, distinguished by color. b Spatial plot of the RG score in selected spots. b′ The RG score in each region (n = 256, 74, 71, and 200 spots, respectively) is quantified and compared. Boxes indicate quartiles, horizontal bar indicates median, and whiskers indicate range, up to 1.5-fold inter-quartile range. P values, Kruskal–Wallis test for between group differences with Holm’s correction for multiple comparisons. c Spatial plot of the Neuron score in selected spots. c′ The Neuron score in each region (n = 256, 74, 71, and 200 spots, respectively) is quantified and compared. Boxes indicate quartiles, horizontal bar indicates median, and whiskers indicate range, up to 1.5-fold inter-quartile range. P values, two tailed Wilcox test for between group differences for comparisons. d The spatial trajectory from Tumor to GCL_N inferred by SPATA2. Highlighted spots were selected for analysis (upper left). Heatmap of genes (row) dynamically expressed in selected spots (column, marked by different colors representing regions) along the spatial trajectory is shown on the right, clustered into three modules. Genes upregulated in GCL_TI exhibiting an ‘one-peak’ pattern is highlighted in red. The expression of FAM20C along the spatial trajectory is shown as an example (bottom left). The shaded area represents the 95% confidence interval. e Venn chart showing the overlap between genes upregulated in GCL_TI (d) and the RG signature genes. The 73 overlapping genes are termed GCL_TI RG signature. f Gene Ontology enrichment analysis of the GCL_TI RG signature genes. P value was determined by two tailed hypergeometric test. Top enriched (Benjamini–Hochberg-adjusted P < 0.05) biological processes are shown. g Kaplan–Meier curves comparing the overall survival between GCL_TI RG signature high (n = 136) versus low (n = 30) groups in the TCGA-GBM cohort. Curve comparison p value was determined by two tailed Log-rank (Mantel–Cox) test.

Fig. 6. FAM20C mediates invasive growth of RG-like cells in a neuron-rich microenvironment.

a Spatial expression pattern of FAM20C in selected spots of DMG1. b The two tailed Paired Samples Wilcoxon Signed Rank Test comparing the expression of FAM20C between paired tumor core and invasive niches from different samples (n = 10). Boxes indicate quartiles, horizontal bar indicates median, and whiskers indicate range, up to 1.5-fold inter-quartile range. c Kaplan–Meier curves comparing the overall survival between FAM20C high (n = 32) versus low (n = 134) groups in the TCGA-GBM cohort. Curve comparison P value was determined by two tailed Log-rank (Mantel–Cox) test. d Left: Illustration of transwell assay with HPT-hNSCs in growth factor deprived medium in the top well and 10% fetal bovine serum in the bottom well. Middle: representative images of crest violet staining for migrated LacZ, F1, and F2 HPT-hNSCs. Right: Quantification of the OD values of invading cells for each group. n = 3 for each group. P values were determined by Student t-test. Error bars indicate mean ± SEM. e Left: Illustration of transwell assay with HPT-hNSCs in the top well and hNSC-derived neurons (day 40) in the bottom well in growth factor deprived medium. Middle: Representative images of crest violet staining for migrated LacZ (n = 3), F1 (n = 3), and F2 (n = 4) HPT-hNSCs. Right: Quantification of the OD values of invading cells for each group. P values were determined by Student t-test. Error bars indicate mean ± SEM. f Left: representative low-magnification images of whole-mount brain sections from LacZ and F1 group (n = 10 for each group). The LacZ and F1 mice shown in this panel both reached the endstage at around 80 days post HPT-hNSC xenograft. White dashed lines mark the tumor inside the brain parenchyma. Red-dashed lines mark the enlarged lateral ventricles (LV). Right: quantification of mCherry⁺ tumor area inside (n = 3 for LACZ group and n = 5 for F1 group) and outside (n = 3 for each group) the brain parenchyma for each group. P values were determined by Student t-test. Error bars indicate mean ± SEM. g Representative images of immunofluorescence (IF) co-labeling of FAM20C/Human Nuclear Antigen (hNu) in tumor areas inside the brains of LacZ and F1 mice (n = 10 for each group). The dashed line marks the border between the tumor core (T) and tumor-infiltrated brain area (TI). Arrows point to FAM20C/hNu co-labeled tumor cells in TI, arrowheads point to adjacent hNu-negative neurons with small, rounded nuclei. h Left: Representative images of IF co-labeling of mCherry/Ki67 in tumor areas inside and outside the brains of LacZ and F1 mice (n = 10 for each group). V, blood vessel. The dashed line marks the border of the brain stem. Right: the ratio of Ki67⁺mCherry⁺ cells among total mCherry⁺ cells in tumor areas inside (n = 3 for LacZ group, n = 5 for F1 group) and outside (n = 3 for LacZ group, n = 5 for F1 group) the brains of LacZ and F1 mice. P values were determined by Student t-test. Error bars indicate mean ± SEM. Source data for d-h are provided as a Source Data file. Scale bars, 50 μm in d, e, g, h, 1 mm in f.