Single-cell RNA-seq highlights intratumoral heterogeneity in primary glioblastoma

Abstract

Human cancers are complex ecosystems composed of cells with distinct phenotypes, genotypes, and epigenetic states, but current models do not adequately reflect tumor composition in patients. We used single-cell RNA sequencing (RNA-seq) to profile 430 cells from five primary glioblastomas, which we found to be inherently variable in their expression of diverse transcriptional programs related to oncogenic signaling, proliferation, complement/immune response, and hypoxia. We also observed a continuum of stemness-related expression states that enabled us to identify putative regulators of stemness in vivo. Finally, we show that established glioblastoma subtype classifiers are variably expressed across individual cells within a tumor and demonstrate the potential prognostic implications of such intratumoral heterogeneity. Thus, we reveal previously unappreciated heterogeneity in diverse regulatory programs central to glioblastoma biology, prognosis, and therapy.

Figure 1. Intratumoral glioblastoma heterogeneity quantified by single cell RNA-seq

(A) Workflow depicts rapid dissociation and isolation of glioblastoma cells from primary tumors for generating single cell and bulk RNA-seq profiles and deriving glioblastoma culture models. (B) Clustering of copy number variation (CNV) profiles inferred from RNA-seq data for all single cells and a normal brain sample. Clusters (dendrogram) primarily reflect tumor-specific CNV (colored bar coded as in panel D). Topmost cluster (red, arrow) contains the normal brain sample and 10 single cells, nine of which correlate with normal oligodendrocyte expression profiles and one with normal monocytes (‘Oligo’ and ‘Mono’, black and white heatmap). (C) Heatmap of CNV signal normalized against the ‘normal’ cluster defined in (B) shows CNV changes by chromosome (columns) for individual cells (rows). All cells outside the normal cluster exhibit chromosome 7 gain (red) and chromosome 10 loss (blue), which are characteristic of glioblastoma. (D) Multidimensional scaling illustrates the relative similarity between all 430 single cells and population controls. The distance between any two cells reflects the similarity of their expression profiles. Cells group by tumor (color code), but each tumor also contains outliers that are more similar to cells in other tumors. (E) RNA-seq read densities (vertical scale of 10) over surface receptor genes are depicted for individual cells (rows) from MGH30. Cell-to-cell variability suggests a mosaic pattern of receptor expression, in contrast to constitutively expressed GAPDH.

Figure 2. Unbiased analysis of intratumoral heterogeneity reveals coherent transcriptional modules

(A) Gene sets that vary coherently between cells in specific tumors or across the global dataset (colored boxes) were identified by principal component analysis or clustering (24). Hierarchical clustering of these gene sets across all cells (tree) reveals four meta-signatures related to hypoxia, complement/immune response, oligodendrocytes and cell cycle. (B) Heatmap shows expression of the cell cycle meta-signature, selected cell cycle gene sets and representative genes from the signature (rows) in individual glioblastoma cells (columns). Cells were grouped by tumor and ordered by meta-signature score. (C) Heatmap depicts hypoxia meta-signature as in (B).

Figure 3. Transcriptional signatures of a stem-like compartment in primary glioblastoma

(A) Stem-like (GSC) and differentiated (DGC) culture models were derived from patient tumor MGH26. GSCs grow as spheres (left, top), initiate tumors in xenotransplantation (right, top), and express the stem cell marker CD133 (right, bottom). (B) Heatmap depicts expression of genes (rows) from a stemness signature in differentiated models (DGC, left columns), stem-like models (GSC, right columns) derived from 3 tumors, and in 70 individual cells from MGH31 (middle). (C) Bar plot depicts the Pearson correlation coefficient (y-axis) between the stemness signature and selected transcriptional modules in each tumor (colored bars): cell cycle, transcriptional targets of POU3F2, SOX2, SALL2, OLIG2 (core TF) (43), NFI transcriptional targets (NFI) (42), and the proneural (PN), classical (CL), mesenchymal (MES) and neural (N) subtypes defined by the Cancer Genome Atlas (21). (D) Plot depicts stemness score (y-axis) computed from stemness signature gene expression in individual cells from each tumor (xaxis) ordered by score. Bar plots depict the overall variance (y-axis, standard deviation) in the stemness score (red) and the average variance of simulated control gene sets (blue), confirming the significance of the gradient.

Figure 4. Individual tumors contain a spectrum of glioblastoma subtypes and hybrid cellular states

(A) Heatmap depicts average expression of classifier genes for each subtype (rows) across all classifiable cells grouped by tumor (columns). PN: proneural, CL: classical, MES: mesenchymal, N: neural. Each tumor contains a dominant subtype, but also has cells that conform to alternate subtypes. (B) Hexagonal plots depict bootstrapped classifier scores for all cells in each tumor. Each data point corresponds to a single cell and is positioned along three axes according to its relative scores for the indicated subtypes (Supplemental Methods). Cells corresponding to each subtype are indicated by solid color, while hybrid cells are depicted by two colors. (C) Kaplan-Meier survival curves are shown for proneural tumors from the Cancer Genome Atlas (21). Intratumoral heterogeneity was estimated based on detected signal for alternative subtypes, and used to partition the tumors into a pure proneural group and three groups with the indicated additional subtype. Tumors with mesenchymal signal had significantly worse outcome than pure proneural (p<0.05). (D) Kaplan-Meier survival curves shown for proneural tumors partitioned based on the relative strength of alternative subtype signatures in aggregate (24). Tumors with high signal for alternative subtypes had significantly worse outcome than pure proneural (p<0.05).