Single-cell RNA sequencing reveals tumor heterogeneity, microenvironment, and drug-resistance mechanisms of recurrent glioblastoma

Abstract

Glioblastomas are highly heterogeneous brain tumors. Despite the availability of standard treatment for glioblastoma multiforme (GBM), i.e., Stupp protocol, which involves surgical resection followed by radiotherapy and chemotherapy, glioblastoma remains refractory to treatment and recurrence is inevitable. Moreover, the biology of recurrent glioblastoma remains unclear. Increasing evidence has shown that intratumoral heterogeneity and the tumor microenvironment contribute to therapeutic resistance. However, the interaction between intracellular heterogeneity and drug resistance in recurrent GBMs remains controversial. The aim of this study was to map the transcriptome landscape of cancer cells and the tumor heterogeneity and tumor microenvironment in recurrent and drug-resistant GBMs at a single-cell resolution and further explore the mechanism of drug resistance of GBMs. We analyzed six tumor tissue samples from three patients with primary GBM and three patients with recurrent GBM in which recurrence and drug resistance developed after treatment with the standard Stupp protocol using single-cell RNA sequencing. Using unbiased clustering, nine major cell clusters were identified. Upregulation of the expression of stemness-related and cell-cycle-related genes was observed in recurrent GBM cells. Compared with the initial GBM tissues, recurrent GBM tissues showed a decreased proportion of microglia, consistent with previous reports. Finally, vascular endothelial growth factor A expression and the blood-brain barrier permeability were high, and the O⁶ -methylguanine DNA methyltransferase-related signaling pathway was activated in recurrent GBM. Our results delineate the single-cell map of recurrent glioblastoma, tumor heterogeneity, tumor microenvironment, and drug-resistance mechanisms, providing new insights into treatment strategies for recurrent glioblastomas.

Keywords: cancer microenvironment; drug resistance; recurrent glioblastoma; scRNA-seq; tumor heterogeneity.

FIGURE 1

Landscape of cell composition in primary and recurrent glioblastoma. (A) Workflow of sample collection, sequencing, and analysis for single‐cell RNA‐seq data. Three samples of recurrent glioblastoma sequenced by Singleron platform. (B) Uniform Manifold Approximation and Projection (UMAP) plot showing the major clusters of all cells in a glioma. Each dot of the UMAP plot represents a single cell. Cells are color coded by their associated cell clusters. (C) UMAP plot showing the illness state of individual tumor cells in glioblastoma patients. Primary glioblastoma patients were collected from CGGA, and recurrent glioblastoma patients were collected from an independent cohort from the Sun Yat‐sen University Cancer Center. (D) UMAP plot showing the major lineages in glioblastoma (GBM). Each dot of the UMAP plot represents a single cell. Cells are color coded by their associated cell lineages. (E) UMAP plots showing the origins of individual tumor cells in GBM patients. (F) The CNV pattern of each cell in immune cells and tumor cells. Red means amplification and blue indicates deletion. The intensity of the red and blue color components correlates with the gain and loss values.

FIGURE 2

Cellular landscape of the tumor cell in primary and recurrent glioblastoma patients. (A) UMAP plot showing the illness state of individual tumor cells in glioblastoma patients. (B) UMAP plot showing the illness state of individual tumor cells in glioblastoma patients. Primary patients were collected from CGGA, and relapse patients were collected from an independent cohort from the Sun Yat‐sen University Cancer Center. (C) UMAP plot showing the major clusters of tumor cells in glioblastoma patients. Each dot of the UMAP plot represents a single cell. Cells are color coded by their associated cell clusters. (D) Normalized Shannon diversity index (NSDI) of major clusters that indicates the evenness of major clusters in each patient. (E) Heatmap showing the illness state preferences of major clusters estimated by the Ro/e score that denotes the ratio of observed to expected cell number.

FIGURE 3

The frequency and evenness of each cell type in samples. (A) UMAP plot showing the major lineages of immune cells in a glioblastoma. Each dot of the UMAP plot represents a single cell. Cells are color coded by their associated cell lineages. (B) Dot plot showing the 16 signature gene expressions across the nine cellular clusters. The size of dots represents the proportion of cells expressing the particular marker, and the spectrum of color indicates the mean expression levels of the markers (log1p transformed). (C) Bar plot showing the fraction of major lineages of tumor cells in recurrent glioblastoma. (D) Bar plot showing the fraction of major lineages of tumor cells in primary glioblastoma. (E) Heatmap showing the illness state preferences of major tumor cell lineages estimated by the Ro/e score that denotes the ratio of observed to expected cell number. (F) Heatmap showing logarithmic interaction scores between all cell subsets.

FIGURE 4

Developmental trajectory and identification of proliferation‐related tumor clusters in recurrent GBM. (A) The left UMAP plot shows the developmental trajectories (black lines). The right UMAP plot shows a pseudo‐time score from dark blue to yellow, indicating developmental states, respectively. Tumor cell subpopulations inferred by monocle3. (B) Gene set enrichment analysis in cluster 12, ordered by generation, enrichment of stem‐like cells, and tumor invasiveness in recurrent GBM. (C) UMAP plot showing the expression of proliferation gene markers in tumor subpopulations. (D) Bar plot showing the cell‐cycle distribution in different illness states. (E) Box plot showing GSVA of cluster 12 using a CGGA cohort. (F) Kaplan–Meier curves show probability of OS among patients with high radial glia signature (top 50%) versus low radial glia signature (bottom 50%) in patients from the CGGA cohort. (G) Kaplan–Meier curves show probability of OS among patients with high radial glia signature (top 50%) versus low radial glia signature (bottom 50%) in patients from TCGA‐GBM LGG cohort.

FIGURE 5

Identification of treatment resistance‐related tumor clusters in glioblastoma. (A) Heat map of differentially expressed genes (based on two‐sided Wilcoxon test) in primary and recurrent GBM. (B) Box plot showing the expression of VEGFA in tumor subpopulations. Each comparison is made using a two‐sided Wilcoxon test. (C) UMAP plot showing the expression of VEGFA in tumor subpopulations. (D) Gene set variation analysis (GSVA) of BBB signature using single‐cell data. Each comparison is made using a two‐sided Wilcoxon test. (E) Violin plots showing the expression of MGMT in tumor subpopulations. Each comparison is made using a two‐sided Wilcoxon test. (F) GSVA of MGMT signature using single‐cell data. Each comparison is made using a two‐sided Wilcoxon test. (G) Kaplan–Meier curves showing the probability of OS among patients with high radial glia signature (top 50%) versus low radial glia signature (bottom 50%) in patients from the CGGA cohort. (H) Kaplan–Meier curves showing the probability of OS among patients with high MGMT signature (top 50%) versus low MGMT signature (bottom 50%) in patients from the CGGA cohort.

FIGURE 6

The changes of microglial cell number and VEGF protein expression in primary and recurrent GBM. (A) Immunohistochemistry images of VEGF protein expression in GBM and r‐GBM. (B) Immunofluorescence images of TMEM119 ⁺ cells and DAPI ⁺ in GBM and r‐GBM. (C) Quantitation of VEGF ⁺ cell relative to total DAPI ⁺ cells and of TMEM119 ⁺ cell relative to total DAPI ⁺ cells in GBM and r‐GBM (n = 6). Two‐tailed Student's t‐test: ***p < 0.001.