Longitudinal single-cell analysis of a myeloma mouse model identifies subclonal molecular programs associated with progression

Abstract

Molecular programs that underlie precursor progression in multiple myeloma are incompletely understood. Here, we report a disease spectrum-spanning, single-cell analysis of the Vκ*MYC myeloma mouse model. Using samples obtained from mice with serologically undetectable disease, we identify malignant cells as early as 30 weeks of age and show that these tumours contain subclonal copy number variations that persist throughout progression. We detect intratumoural heterogeneity driven by transcriptional variability during active disease and show that subclonal expression programs are enriched at different times throughout early disease. We then show how one subclonal program related to GCN2 stress response is progressively activated during progression in myeloma patients. Finally, we use chemical and genetic perturbation of GCN2 in vitro to support this pathway as a therapeutic target in myeloma. These findings therefore present a model of precursor progression in Vκ*MYC mice, nominate an adaptive mechanism important for myeloma survival, and highlight the need for single-cell analyses to understand the biological underpinnings of disease progression.

Fig. 1. A single-cell transcriptional map of malignant cells from progressing Vκ*MYC mice.

a Schematic of Vκ*MYC mouse cohort and experimental workflow for collection of single cells. Graphics created in part using BioRender.com. b Disease burden across cohort as determined by M-protein measurements from SPEP. Statistical comparison of multiple groups was performed using a Wilcoxon rank-sum test (two-sided) corrected for multiple testing (Benjamini−Hochberg). Boxplots represent the distribution of each measurement within defined groups, where the central rectangle spans the interquartile range, the central line represents the median, and “whiskers” above and below the box show the value 1.5× the interquartile range. Only P values for statistically significant comparisons are listed. Data points represent measurements from biologically-independent animals (Control (n = 3), early-MM (n = 5), int-MM (n = 3), and active-MM (n = 7)). c tSNE visualization of 17,504 B lineage cells, coloured by sample ID. d tSNE visualization of B lineage cells, coloured by plasma cell gene signature score (Hay SB et al.). e tSNE visualization of B lineage cells, coloured by the relative expression of indicated cell-type-specific genes. f tSNE visualization of plasma cells identified by criteria displayed in (d) and (e) coloured by the relative expression of Vκ*MYC transgene (tMYC). g Bar plot showing the distribution of normal vs. malignant plasma cells across samples No malignant cells were detected in EMM2 and thus it was removed from downstream malignant cell analyses. Source data are provided in SourceData_Fig. 1.xlsx. EMM: early-MM, IMM: intermediate-MM, AMM: active-MM.

Fig. 2. Core versus disease-stage specific gene expression programs in malignant cells from Vκ*MYC mice.

a Heatmap of differentially expressed genes shared by all malignant cells in Vκ*MYC mice compared to normal plasma cells (FDR < 0.05). Heatmap is split vertically to show normal plasma cells (nPC) versus malignant plasma cells (mPC), the latter of which is further split by disease stage group. The upper and lower panels of the heatmap separate upregulated and downregulated genes, respectively. A subset of 100 randomly selected cells per disease stage group are shown and data represent scaled expression values (any values outside a range of −2 to 2 were clipped). b Top 20 positively/negatively enriched terms from MSigDB gene set enrichment analysis (H, C2, C6, FDR < 0.05) computed using core upregulated/downregulated genes identified by DE analysis in (a). c–e Disease stage-specific genes that are significantly differentially expressed between disease stage groups. Coloured dots represent the mean expression of disease stage samples for each gene, with error bars depicting the standard error of the mean. Statistical comparisons were performed using a two-sided t-test with subsequent correction for multiple testing (Bonferroni). Grey data points represent mean expression of respective genes in cells from each biologically-independent animal (Cont1 = 44 cells, Cont2 = 72 cells, Cont3 = 148 cells, EMM1 = 45 cells, EMM4 = 52 cells, EMM5 = 71 cells, IMM1 = 206 cells, IMM2 = 88 cells, IMM3 = 149 cells, AMM1 = 2,003 cells, AMM2 = 830 cells, AMM3 = 1,379 cells, AMM4 = 822 cells, AMM5 = 302 cells, AMM6 = 323 cells, AMM7 = 310 cells). Genes are grouped according to the pattern of expression throughout progression. Source data are provided in SourceData_Fig. 2.xlsx.

Fig. 3. Drivers of intratumoural heterogeneity in the malignant cell compartment of Vκ*MYC mice.

a Heatmap of genome-wide copy number variations (CNVs) inferred from scRNA-seq data of malignant plasma cells as determined using InferCNV. Columns represent genome position across chromosomes. Rows represent CNVs averaged by intra-sample CNV subpopulation, which were identified using the ward.D2 hierarchical clustering/random forest method implemented by analysis_mode = ‘subclusters’ in inferCNV. CNV-driven subpopulation sizes ranged from 1 to 979 cells (median 44). The height of each CNV subpopulation is proportionate to its fractional composition within a given tumour. b UMAP visualization of malignant cells from each active-MM mouse coloured by transcriptional cluster. Gene expression-driven cluster sizes ranged from 10 to 474 cells (median 79). c UMAP visualization of malignant cells from each active-MM mouse coloured by CNV subpopulation. d Bar plot showing the distribution of CNV subpopulations (fill) across transcriptional clusters (x-axis). Results for (b–d) are organized for each active-MM mouse in columns, with subject names and the number of cells/transcriptional clusters listed above. Source data are provided in SourceData_Fig. 3.xlsx. EMM: early-MM, IMM: intermediate-MM, AMM: active-MM.

Fig. 4. Molecular programs driving intratumoural heterogeneity in Vκ*MYC mice with active-MM.

a Heatmap of Jaccard Index between significantly enriched Reactome terms across 41 intra-tumour malignant cell clusters. Groupings of clusters with increased similarity (Similarity Programs) were determined according to the complete linkage method for hierarchical clustering and are labelled below the heatmap. Columns are annotated with information related to sample identity, cell cycle phase, and Mki67/Top2a expression. b Map of Reactome terms with significant enrichment in malignant cell clusters from Similarity Program A (ISR-GCN2). The full hierarchy of each Reactome pathway is shown for context but only significantly enriched shared pathways are highlighted in orange. c Gene set scoring for indicated Reactome signatures calculated using Seurat’s AddModuleScore across disease groups in Vκ*MYC data (Control = 264 cells, early-MM = 168 cells, int-MM = 443 cells, active-MM = 5,969 cells). Boxplots within violin plots represent the distribution of each measurement within defined groups, where the central rectangle spans the interquartile range, the central line represents the median, and “whiskers” above and below the box show the value 1.5× the interquartile range. Statistical comparison of multiple groups (normal PCs vs. each Vκ*MYC disease group) was performed using the Wilcoxon rank-sum test (two-sided) corrected for multiple testing (Benjamini−Hochberg). Source data are provided in SourceData_Fig. 4.xlsx.

Fig. 5. Proposed model for patterns of tumour heterogeneity during myeloma disease progression in Vκ*MYC mouse model.

Graphics created in part using BioRender.com.

Fig. 6. ISR-GCN2 pathway is progressively activated in myeloma patients and supports myeloma cell survival.

a ISR-GCN2 gene signature scoring in publicly-available microarray patient data from Chng et al. calculated by taking the mean of scaled expression values for genes from the ISR-GCN2 gene set. b ISR-GCN2 gene signature scoring calculated using Seurat’s AddModuleScore in publicly-available scRNA-seq patient data from Ledergor et al. across disease groups. Statistical comparisons of multiple groups (Healthy vs. each disease group) in (a) and (b) were performed using the Wilcoxon rank-sum test (two-sided) corrected for multiple testing (Benjamini−Hochberg). Boxplots in (a) and (b) represent the distribution of each measurement within defined groups, where the central rectangle spans the interquartile range, the central line represents the median, and “whiskers” above and below the box show the value 1.5× the interquartile range. c MTT analysis of GCN2iB treatment in HMCLs (1 μM, 72 h). Data represent the mean of biological replicates (n = 2 for MY5, JJN3, XG6, AMO1, and RPMI8226; n = 3 for XG7, OPM2; n = 4 for MM1S, U266) with error bars representing standard deviation shown for samples with more than two biological replicates. d Pearson correlation (cor.test, two-sided) between the ISR-GCN2 gene signature score (x-axis) and viability relative to DMSO (1 μM GCN2iB, y-axis). The linear regression line is plotted in black with confidence interval shaded grey. Each dot represents one HMCL. e Flow cytometric analysis of apoptosis in U226 cells after 48 h treatment with 5 μM GCN2iB. f Flow cytometric analysis of apoptosis in Vκ12598 tumour cells (CD138⁺/B220⁻) after 48 h treatment with 5 μM GCN2iB. Scatter plots in (e) and (f) are representative of multiple independent experiments (see Supplementary Fig. 8c, d). Quadrants corresponding to viable tumour cells are highlighted in red. g Bar plots depicting cell survival, as determined by trypan blue assay, in GCN2 knockout HMCLs. Bar heights represent mean relative viability from two independent experiments (as shown by individual data points). Source data are provided in SourceData_Fig. 6.xlsx.