Single-cell RNA sequencing reveals the impact of chromosomal instability on glioblastoma cancer stem cells

Abstract

Background: Intra-tumor heterogeneity stems from genetic, epigenetic, functional, and environmental differences among tumor cells. A major source of genetic heterogeneity comes from DNA sequence differences and/or whole chromosome and focal copy number variations (CNVs). Whole chromosome CNVs are caused by chromosomal instability (CIN) that is defined by a persistently high rate of chromosome mis-segregation. Accordingly, CIN causes constantly changing karyotypes that result in extensive cell-to-cell genetic heterogeneity. How the genetic heterogeneity caused by CIN influences gene expression in individual cells remains unknown.

Methods: We performed single-cell RNA sequencing on a chromosomally unstable glioblastoma cancer stem cell (CSC) line and a control normal, diploid neural stem cell (NSC) line to investigate the impact of CNV due to CIN on gene expression. From the gene expression data, we computationally inferred large-scale CNVs in single cells. Also, we performed copy number adjusted differential gene expression analysis between NSCs and glioblastoma CSCs to identify copy number dependent and independent differentially expressed genes.

Results: Here, we demonstrate that gene expression across large genomic regions scales proportionally to whole chromosome copy number in chromosomally unstable CSCs. Also, we show that the differential expression of most genes between normal NSCs and glioblastoma CSCs is largely accounted for by copy number alterations. However, we identify 269 genes whose differential expression in glioblastoma CSCs relative to normal NSCs is independent of copy number. Moreover, a gene signature derived from the subset of genes that are differential expressed independent of copy number in glioblastoma CSCs correlates with tumor grade and is prognostic for patient survival.

Conclusions: These results demonstrate that CIN is directly responsible for gene expression changes and contributes to both genetic and transcriptional heterogeneity among glioblastoma CSCs. These results also demonstrate that the expression of some genes is buffered against changes in copy number, thus preserving some consistency in gene expression levels from cell-to-cell despite the continuous change in karyotype driven by CIN. Importantly, a gene signature derived from the subset of genes whose expression is buffered against copy number alterations correlates with tumor grade and is prognostic for patient survival that could facilitate patient diagnosis and treatment.

Keywords: Aneuploidy; CIN; CNV; CSCs; Cancer stem cells; Chromosomal instability; Copy number variation; Glioblastoma; Heterogeneity.

Fig. 1

CB660 NSCs and GliNS2 CSCs have distinct gene expression profiles. a, Heatmap showing the normalized gene expression of the most variably expressed genes used for hierarchical clustering analysis. The top dendrogram shows that CB660 NSCs and GliNS2 CSCs cluster into two distinct populations. In the row above the heatmap, blue indicates CB660 NSCs and red indicates GliNS2 CSCs. b, Principal component analysis (PCA) and cell cycle phase analysis of CB660 NSCs and GliNS2 CSCs. The graph shows the separation of CB660 NSCs (circles) and GliNS2 CSCs (triangles) into distinct groups. The color of the circles or triangles corresponds to the predicted cell cycle phase of each cell. c, Pie plots showing the fraction of CB660 NSCs and GliNS2 CSCs in each phase of the cell cycle. P > 0.05; Chi-square test comparing the cell cycle profiles of CB660 NSCs and GliNS2 CSCs

Fig. 2

Gene expression scales with chromosome copy number. a, Heatmap of the estimated copy number (ECN) of all chromosomes (columns) in single CB660 NSCs and GliNS2 CSCs (rows). On the scale, ECN = 0 indicates diploid gene expression levels. The column adjacent to the heatmap shows the cell cycle phase of each cell as determined in Fig. 1b with the color of the bar corresponding to the predicted cell cycle phase. b, Quantification of chromosomal instability in CB660 NSCs and GliNS2 CSCs. Bar, median; box 25th to 75th percentile; whiskers, minimum and maximum. P < 2E-16; Mann-Whitney U test comparing CB660 NSCs and GliNS2 CSCs. c and d, Heatmaps showing normalized gene expression in CB660 NSCs and GliNS2 CSCs for chromosome 7 (C) and chromosome 13 (D)

Fig. 3

Identification of copy number dependent and independent differentially expressed genes between CB660 NSCs and GliNS2 CSCs. a, Volcano plot showing both unadjusted and copy number adjusted differentially expressed genes between CB660 NSCs and GliNS2 CSCs. The dashed line shows the statistical significance cut-off (P.adjust< 0.05, Mann-Whitney U test and Bonferroni adjustment) used for differential gene expression analysis. Dark red points indicate genes that remain significantly differentially expressed after copy number adjustment while light red points indicate genes that are not significantly differentially expressed after copy number adjustment. b, Venn diagram showing the number of overlapping differentially expressed unadjusted and copy number adjusted genes. c, Scatter plot showing the correlation between the negative log10 transformed p-value of unadjusted and copy number adjusted differential gene expression analysis. Spearman correlation coefficient = 0.77. d, Scatter plot showing the correlation between c-MYC expression and c-MYC copy number in CB660 NSCs and GliNS2 CSCs. Spearman correlation coefficient = 0.54. e, Bar graphs quantifying c-MYC expression in CB660 NSCs and GliNS2 CSCs before (left graph) and after copy number adjustment (right graph). Bar, median; box 25th to 75th percentile; whiskers, minimum and maximum. P.adjust> 0.05 after copy number adjustment, Mann-Whitney U test and Bonferroni adjustment. f, Scatter plot showing the correlation between SLC23A2 expression and SLC23A2 copy number in CB660 NSCs and GliNS2 CSCs. Spearman correlation coefficient = 0.51. g, Bar graphs quantifying SLC23A2 expression in CB660 NSCs and GliNS2 CSCs before (left graph) and after copy number adjustment (right graph). Bar, median; box 25th to 75th percentile; whiskers, minimum and maximum. P.adjust = 4E-3 after copy number adjustment, Mann-Whitney U test and Bonferroni adjustment. h, Scatter plot showing the correlation between TFAP2C expression and TFAP2C copy number in CB660 NSCs and GliNS2 CSCs. Spearman correlation coefficient = − 0.02. i, Bar graphs quantifying TFAP2C expression in CB660 NSCs and GliNS2 CSCs before (left graph) and after copy number adjustment (right graph). Bar, median; box 25th to 75th percentile; whiskers, minimum and maximum. P.adjust = 0.02 after copy number adjustment, Mann-Whitney U test and Bonferroni adjustment

Fig. 4

A copy number independent gene signature is prognostic for tumor grade and patient survival. a, Heatmap showing the relative gene expression of the copy number independent differentially expressed genes that are involved in pathways that negatively regulate cell proliferation in GliNS2 CSCs compared to CB660 NSCs. b, Graph showing the copy number independent (CI) gene signature score for glioma samples in data set GSE1993 stratified by histological grade. Grade I/II = pilocytic astrocytoma (n = 2) and diffuse astrocytoma (n = 5), Grade III = anaplastic astrocytoma (n = 19), and Grade IV = glioblastoma (n = 39). Bar, median; box 25th to 75th percentile; whiskers, minimum and maximum. P = 1E-3, ANOVA analysis. c, Kaplan-Meier plots showing that the CI gene signature score is prognostic for patient survival in four independent data sets. The median CI score was used as the cutoff to dichotomize patients into CI-Hi and CI-Lo groups with the number of patients in each group indicated in parentheses. Hazard ratios were calculated using a Cox regression model, and p-values were calculated by using log-rank tests to determine statistical differences between survival curves. Below each Kaplan-Meier plot is a table showing the number of patients at risk over time