Spatially resolved transcriptomic profiles reveal unique defining molecular features of infiltrative 5ALA-metabolizing cells associated with glioblastoma recurrence

Abstract

Background: Spatiotemporal heterogeneity originating from genomic and transcriptional variation was found to contribute to subtype switching in isocitrate dehydrogenase-1 wild-type glioblastoma (GBM) prior to and upon recurrence. Fluorescence-guided neurosurgical resection utilizing 5-aminolevulinic acid (5ALA) enables intraoperative visualization of infiltrative tumors outside the magnetic resonance imaging contrast-enhanced regions. The cell population and functional status of tumor responsible for enhancing 5ALA-metabolism to fluorescence-active PpIX remain elusive. The close spatial proximity of 5ALA-metabolizing (5ALA +) cells to residual disease remaining post-surgery renders 5ALA + biology an early a priori proxy of GBM recurrence, which is poorly understood.

Methods: We performed spatially resolved bulk RNA profiling (SPRP) analysis of unsorted Core, Rim, Invasive margin tissue, and FACS-isolated 5ALA + /5ALA - cells from the invasive margin across IDH-wt GBM patients (N = 10) coupled with histological, radiographic, and two-photon excitation fluorescence microscopic analyses. Deconvolution of SPRP followed by functional analyses was performed using CIBEROSRTx and UCell enrichment algorithms, respectively. We further investigated the spatial architecture of 5ALA + enriched regions by analyzing spatial transcriptomics from an independent IDH-wt GBM cohort (N = 16). Lastly, we performed survival analysis using Cox Proportinal-Hazards model on large GBM cohorts.

Results: SPRP analysis integrated with single-cell and spatial transcriptomics uncovered that the GBM molecular subtype heterogeneity is likely to manifest regionally in a cell-type-specific manner. Infiltrative 5ALA + cell population(s) harboring transcriptionally concordant GBM and myeloid cells with mesenchymal subtype, -active wound response, and glycolytic metabolic signature, was shown to reside within the invasive margin spatially distinct from the tumor core. The spatial co-localization of the infiltrating MES GBM and myeloid cells within the 5ALA + region indicates PpIX fluorescence can effectively be utilized to resect the immune reactive zone beyond the tumor core. Finally, 5ALA + gene signatures were associated with poor survival and recurrence in GBM, signifying that the transition from primary to recurrent GBM is not discrete but rather a continuum whereby primary infiltrative 5ALA + remnant tumor cells more closely resemble the eventual recurrent GBM.

Conclusions: Elucidating the unique molecular and cellular features of the 5ALA + population within tumor invasive margin opens up unique possibilities to develop more effective treatments to delay or block GBM recurrence, and warrants commencement of such treatments as early as possible post-surgical resection of the primary neoplasm.

Keywords: 5ALA; Glioblastoma; Glycolysis; Mesenchymal subtype; Myeloid; Recurrence; Spatial transcriptomics; Wound response.

Fig. 1

Differential gene expression analysis in spatially distinct GBM regions. Schematic figure delineating the steps for tissue collection from distinct GBM regions (Core, Rim, and Invasive margin) followed by 5ALA-based FACS isolation of Invasive margin cells into 5ALA + and 5ALA − subpopulations (A). The spatially resolved RNA profile (SPRP) from each unsorted region and sorted cells were interrogated by gene set enrichment analysis (GSEA), Deconvolution algorithm, single-cell gene signature scoring (SC sig-scoring), Exon–intron split analysis (EISA), network inference, and complemented by spatially resolved transcriptomics (stRNA-seq) (A). Representative hematoxylin and eosin-stained images of 5ALA + GBM infiltrative margin, scanned using NanoZoomer, magnification × 40, scale bar 200 μm for patient 28 and 100 μm for patient 37. Top (patient 28): Dashed line boundary demarcates cellular tumor-filling white matter. Individual tumor cells are observed within the cerebral neocortex (black circles), with neurons normally distributed (normally formed cerebral neocortex) (green circles). Microvascular proliferation (dark blue arrow) beyond the margin of the cellular tumor. Representative normal blood vessel (red arrow). Bottom (patient 37): Individual tumor cells are observed within the cerebral neocortex (black circles), with a normal distribution of white matter blood vessels (red circles). Deep white matter diffusely infiltrated by tumor cells at low density (green arrows). Other cells observed within this deep white matter region are oligodendrocytes with macrophages and astrocytes (B). The Inv margin representing the non-enhancing on T1 with gadolinium region located outside the MRI contrast region (C). The two-photon excitation fluorescence image demonstrates a distribution of PpIX across radiologically defined spatially distinct tumor core (Core) (D, Upper panel) and Invasive margin (Inv) (D, Lower panel) regions. Volcano plot representing differential gene expression between 5ALA + and Tumor Core (E), Rim (F), Invasive margin (G), and 5ALA − cells (H). Heatmap showing the normalized enrichment scores (NES) representing significantly enriched hallmarks (padj < 0.05) in a specific GBM region (I). GSEA was performed between a particular region and all other regions. The color represents the value of NES where yellow and black indicate the highest (NES = 3.5) and lowest (NES = 0) NES values, respectively. Only significant NES values are shown (adj. p-values < 0.05). Ki67 immunohistochemistry (IHC) of tumor core—superficial medial (J), anterior medial (K), tumor rim (L), and Invasive margin (M) to estimate the fraction of proliferating cells in spatially distinct regions of GBM. The scale bar indicates 25 µm. CD31-IHC represents the tumor vascularity identifying the number of vascular structures in the Core (N), Rim (O), and Invasive margin (P)

Fig. 2

Enrichment of GBM subtypes, molecular and metabolic gene signatures in distinct GBM regions. The normalized enrichment scores (NES) of the significantly enriched (padj < 0.05) GBM subtypes are shown in distinct GBM regions and 5ALA sorted cells (A). Retrieval of gene sets specifying different GBM subtypes (Verhaak et al.) was followed by gene set enrichment analysis (GSEA). The color code indicates the differential NES values (yellow and black represent higher and lower NES, respectively). GSEA plot shows that the GBM mesenchymal subtype is significantly enriched (NES: 2.1, padj = 2.3 × 10⁻⁶) in 5ALA + cells (B). NeuN immunohistochemistry (IHC) of Core—(C), Rim (D), and Invasive margin (E) to estimate the proportion of NeuN positive (neuronal) cells (arrows). The scale bar indicates 25 µm. Differential z-scored normalized expression Log₂(TPM + 1) of the significantly regulated (Limma, padj < 0.05) leading edge genes of GBM subtypes (Verhaak et al.)—Neural (F), and Mesenchymal (G)—are shown as a heatmap. Heatmap illustrating the normalized enrichment scores (NES) representing enriched cellular and metabolic states (padj < 0.05) in distinct GBM intratumor regions and 5ALA sorted cells (H). The gene signatures of cellular states (Developmental and Inflammatory wound response) and metabolic states (Glycolytic—GPM, and Mitochondrial—MTC) were retrieved and subjected to GSEA. GSEA plots represent the enrichment of GPM (I), and Inflammatory wound response (J) in 5ALA + cells compared to 5ALA − cells and Invasive margin, respectively. Heatmaps showing the differential expression of significantly regulated (Limma, padj < 0.05) leading edge genes of GPM (K) and Inflammatory wound response (L) in Core, Rim, Invasive margin, 5ALA − and 5ALA + cells. Stacked bar plot representing the transcriptional program estimates across 10 GBM samples (M). Each transcriptional program (Developmental and Inflammatory wound response) was divided into three expression-based categories based on their gene expression pattern (High, Intermediate, and Low). Each color indicates a specific expression-based category of a transcriptional program. The log₁₀ ratios of High Inflammatory wound response and High Developmental transcriptional program in the unsorted Core region and 5ALA + cells across 10 GBM patients are shown (N). Box plots represent the median expression in Developmental and Inflammatory wound response genes enriched in cells with the high expression of different transcriptional programs (Developmental-High and Invasive margin-High) across Core and 5ALA + cell populations (O). P-values calculated from Student’s t tests are shown

Fig. 3

Enrichment of 5ALA + gene signatures in myeloid-like aneuploid, diploid myeloid, and non-myeloid malignant cells. tSNE plot representation of the sing cells (N = 261,092) annotated by GBmap as malignant and myeloid cell types (A). The color code represents the proposed cell annotation from GBmap. Copy number variation (CNV) categories (aneuploid and diploid) based on were mapped as different colors (aneuploid: Orange and diploid: Green) (B). tSNE plot with Louvain clustering of myeloid cells from GBmap dataset (N = 134,405, as annotated by Gbmap) (C). tSNE representation of the myeloid-like aneuploid (N = 5171) (left) and diploid myeloid (N = 127,483) (right) cells. The color code represents the 5ALA + UCell enrichment score calculated based on the expression of the 5ALA + -specific genes (N = 251) (D). UCell scores of different gene signatures (5ALA + , GPM, Inf. wound, Inf. response, MES, and TNFα) across aneuploid myeloid (AM), diploid myeloid (DM), and non-myeloid-like malignant (malignant) cells (E). CIBERSORTx-derived signature matrix based on three cell annotations (DM, MLA, and Malignant) (F). Color code represents the z-scored expression. Selected differentially expressed genes are shown. Estimated fractions of different cellular states (MLA, DM, and Malignant) across unsorted Core and sorted 5ALA + and 5ALA − cells (G). Louvain clustering based on the single-cell RNA-seq data from GBmap (H). The color code represents 18 different clusters (cluster 0 to cluster 11). tSNE plots of GBmapdataset are shown where the color code represents the single-cell wise UCell scores for different gene signatures—MES1-like (I), MES2-like (J), Inf. wound response (K), GPM (L), MTC (M), and 5ALA + (N)

Fig. 4

5ALA + cells represent CD44 expressing mesenchymal cells. mRNA-based stemness index (mRNAsi) values across distinct GBM regions (Core, Rim, Invasive margin, 5ALA − and 5ALA +) for each patient are represented in a heatmap (A). Row (samples) and columns (GBM regions) are clustered by using a correlation algorithm. Comparison of the mRNAsi values according to brain regions (left) and TCGA-GBM samples (right) are shown as bar diagrams (B). The TCGA-GBM samples were pre-stratified according to GBM subtypes as described by Verhaak et al. [10]. Kruskal–Wallis tests showed a significantly higher mRNAsi in 5ALA + cells compared to Core and Rim, with p-values shown. Pearson correlation coefficient values between the mRNAsi and mRNA expression of selected genes are shown (C). The genes that showed a significant (p-value < 0.05) positive correlation with mRNAsi in 5ALA + cells for each patient were selected. tSNE clustering plots based on scRNA-seq dataset (Couturier et al.) are shown (D–G). The color code represents the single-cell wise UCell scores for different cell-signatures—glial progenitor cells (GP) (D), oligo-lineage cells (OLC) (E), mixed population including truncated radial glial cells, and cancer mesenchymal cells (F), and 5ALA + cells (G). tSNE clustering plots (Couturier et al.) with the color code representing the single-cell wise gene expression values for different marker genes—CD44 (H), AQP4 (I), FAM107A (J), and SOX9 (K), GLI3 (L), and TIMP1 (M) are shown. Heatmap showing the z-scored log2 TPM expression of selected marker genes (GIL3, TIMP1, FAM107A, SOX9, AQP4, and CD44) based on scRNA-seq dataset (Couturier et al.) across different cell types (5ALA + , GP, OLC, and mixed population) (N). qPCR validation results showing Log₂ gene relative gene expression of CD44 gene in spatially resolved RNA profiles across Core, Invasive margin, and 5ALA + cells (O). P-values are shown as calculated by paired T-test

Fig. 5

Identification of rare infiltrative 5ALA + cell cluster within GBM invasive margin. GSEA score of 5ALA + gene signature (A), MES (B), and Inf. wound response (C) are shown for patient sample UKF#334, with enriched (red) and random (blue) spots. Spatial localization of the 5ALA + gene signature is marked (dashed box), where red (5ALA-INV) and green (5ALA-CT) boxes represent 5ALA + enriched spots. GSEA between 5ALA-CT and 5ALA-INV is shown as a bubble plot (D). Enrichment scores (− log₁₀ FDR) are color-coded (yellow—high enrichment; black—low enrichment). H&E-stained tissue section of patient UKF#313 showing the necrotic core, pseudopalisading cells, and cellular tumor (CT) regions (E). Inferred CNV analysis of chromosomes 7 (F) and 10 (G) presented by a spatial surface plot where gain and loss of chromosomes 7 and 10 respectively are color-coded. GSEA of 5ALA + gene signature illustrated by a surface plot where GSEA score is depicted by color code and 5ALA + enriched spots marked (dashed box) (H). H&E-stained tissue section of patient UKF#269 showing CT and adjacent non-tumor (NT) regions (I). Inferred CNV analysis of chromosomes 7 (J) and 10 (K) is shown by a surface plot. GSEA of 5ALA + gene signature is shown by a surface plot where 5ALA + enriched spots are marked (dashed box) (M). H&E-stained tissue section of patient UKF#275 showing CT with surrounding pseudopalisading cells (M). The scaled spatial program score indicates the expression of the reactive immune program (N). GSEA score representing spatial arrangements of 5ALA + gene signature (O) and GPM (P) enriched spots (dashed box). Spatial locations of 5ALA + enriched spots are represented by the GSEA score where red indicates high enrichment of 5ALA + gene signature (Q). The 5ALA-INV and 5ALA-CT spots are indicated by red and green boxes, respectively. GSEA scores representing GPM (R), MTC (S), and hypoxia (T). Spatially weighted correlation analysis and spatial overlap of 5ALA + gene signature with established transcriptional signatures (U). Gene sets with high cross-correlation with 5ALA + gene signature are marked (dashed box) (U). CNV analysis of Chr7 (left) and Chr10 (right), across four samples (UKF#275, UKF#243, UKF#251, and UKF#334) (V). The color code indicates the gain or loss of chromosomes. Each box represents the average CNV value from 50 selected spots across 5ALA and NT regions. The stars indicate the significance of the p-value < 0.001

Fig. 6

Association of 5ALA + gene signatures with recurrence in GBM patients. Stacked bar plots representing the CIBERSORTx-derived average cell-type estimates from Netfel, Richards, and Garofano et al. datasets across primary and recurrent GBM patients from TCGA (A). Heatmap showing the − log10 adjusted p-value of the enriched gene signatures across recurrent and primary GBM tumors (B). The conditions representing the comparisons between recurrent and primary GBM tumors from TCGA and CCGA are given in columns (TCGA—Recurrent vs. Primary and CCGA—Recurrent vs. Primary). GBM data from TCGA has been analyzed in a paired and unpaired manner. Columns are divided into upregulated (Red) or downregulated (Blue) segments based on the regulation of genes between Recurrent vs. Primary samples. Each row represents the different gene signatures. Scatter plots representing the correlation between 5ALA + gene signature scores and overall survival (months) in recurrent and primary patients from GLASS (C and D) cohort. Spearman correlation coefficient (R) and p-values are shown. Forest plot representing the hazard ratio of different factors in recurrent GBM patients (E). tSNE coupled with Louvain clustering of the primary and recurrent GBM cells from two patients (F). 5ALA + UCell score mapping onto the tSNE (G). Comparison of 5ALA + UCell score between primary and recurrent cells (H). (Bar represents the mean). Schematic diagram showing the 5ALA + cells with MES subtype and distinct transcriptional programs such as wound response signatures are associated with recurrence of GBM tumor (I)