Subclone-specific microenvironmental impact and drug response in refractory multiple myeloma revealed by single-cell transcriptomics

Abstract

Virtually all patients with multiple myeloma become unresponsive to treatment over time. Relapsed/refractory multiple myeloma (RRMM) is accompanied by the clonal evolution of myeloma cells with heterogeneous genomic aberrations and profound changes of the bone marrow microenvironment (BME). However, the molecular mechanisms that drive drug resistance remain elusive. Here, we analyze the heterogeneous tumor cell population and its complex interaction network with the BME of 20 RRMM patients by single cell RNA-sequencing before/after treatment. Subclones with chromosome 1q-gain express a specific transcriptomic signature and frequently expand during treatment. Furthermore, RRMM cells shape an immune suppressive BME by upregulation of inflammatory cytokines and close interaction with the myeloid compartment. It is characterized by the accumulation of PD1⁺ γδ T-cells and tumor-associated macrophages as well as the depletion of hematopoietic progenitors. Thus, our study resolves transcriptional features of subclones in RRMM and mechanisms of microenvironmental reprogramming with implications for clinical decision-making.

Fig. 1. RRMM samples and scRNA-seq data set.

a Typical disease course of RRMM patients and sampling time-points of the study. b Workflow for scRNA-seq including CD138 sorting of myeloma cells. c Overview of scRNA-seq data analysis. Left: clustering and cell type identification was based on scRNA-seq of 212,404 cells from primary mononuclear bone marrow samples of RRMM patients (n = 20) as shown by an UMAP embedding colored by sample without batch-effect correction. Middle: myeloma cells subclones according to genetic aberrations were identified while changes of BME cells were evaluated against the HCA data set of healthy donors. Right: by integrating these data, subclone-specific interactions of myeloma cells with BME cells were revealed. d UMAP embedding as shown in panel (c) but colored according to expression of the plasma cell marker TNFRSF17. e Same as panel (d) but colored for plasma cells according to the annotation with SingleR.

Fig. 2. scRNA-seq analysis of RRMM tumor cells.

a UMAP embedding of RRMM PCs colored by patient. Dashed circle marks non-malignant plasma cells (nPCs). The pie chart inset shows the nPCs fraction colored according to patient. b Pearson correlation matrix of averaged gene expression levels per patient. Top, cytogenetic information; bottom, averaged gene expression levels of five MM-subtype-specific genes. c CNAs of exemplary patient sample. Top, heatmap of RRMM05 CNA signal normalized against nPCs derived from the HCA bone marrow reference data set. Horizontal lines divide subclones; bottom, coverage plot showing total copy number derived from whole-genome sequencing data of the same sample; * indicates agreement between both modalities in detecting major CNAs. d UMAP embedding of RRMM05 tumor cells colored by transcriptional cluster. e Same as panel (d) but colored according to CNA clone. f Same as panel (d) but with subclone 2 cells in black. g Bar plot of relative subclone abundances in RRMM05. h Heatmap of 348 differentially expressed genes between subclone 2 and 3. Example genes are listed and shown in bold letters if located on 1q. Subclone 3 has been downsampled to 1000 cells for visualization. Thresholds for differential expression using two-sided Wilcoxon rank-sum test were p-value < 0.05 (Bonferroni-adjusted) and logFC > 0.1. i Bar plot showing fraction of differentially upregulated genes in subclone 2 compared to subclone 3 by chromosomal location. j Scatter plot of +1q transcriptome signature score against fraction of +1q cells as determined by InferCNV. The Pearson’s correlation coefficient was R = 0.86 (p = 2.1 × 10⁻¹¹). Grouping of samples/patients into +1q “not-detected/rare” and “dominant” is indicated by vertical lines. k Violin plot of +1q signature scores for cells of the three different +1q groups “ND/rare” (n = 22,206 cells), “subclonal” (n = 13,345), and “dominant” (n = 45,946). The p-values from a Kruskal-Vallis test were <2 × 10⁻¹⁶. Box plot: center line, median; box limits, first and third quartile; whiskers, minimum/maximum values. l Genes of the +1q signature. The bar plot shows the fraction of subjects in which individual genes were upregulated in +1q clones when compared to the most similar clone without +1q.

Fig. 3. Treatment response of +1q subclones.

a Heatmap of RRMM13 CNA signal (pre- and post-treatment data combined). Horizontal lines divide subclones. b Bar plot of fraction of +1q cells pre- and post-treatment for RRMM13. c UMAP embedding of RRMM13 tumor cells before treatment. Coloring depicts transcriptional cluster (left), +1q cells (middle) and +1q signature score (right). d Violin plot and associated histogram of +1q signature score of RRMM13 tumor cells’ pretreatment. The positive predictive value (PPV) = 0.076 for the score of 0.2 indicated by the dashed line. e Line plot of subclone fraction derived from the CNA analysis per patient pre/post treatment. Red lines mark +1q subclones. Top, CNA stability score per patient; bottom, best treatment response.

Fig. 4. Analysis of cellular interaction of myeloma and BME cells.

a UMAP embedding of CD138^– BME cells colored by cell type (Supplementary Table 3). The inset provides a schematic overview for the location of major cell types. b Gene expression dot plot of major marker genes for individual cell types. c UMAP embedding as shown in panel (a) as point-density plot split in RRMM vs. healthy donors. d Bar plot of cell type fractions for RRMM (right) and healthy (left) donors individually. e Cellular interactions of myeloma tumor cells and BME cell types. Ligand–receptor expression was ordered according to the number of detected interactions. f Bar plot of the cumulative number of interactions detected in a given BME cell type in RRMM samples in comparison to healthy donors. Gray, all interactions; red, interactions gained in RRMM; blue, interactions lost in RRMM. Only selected interactions detected in ≥3 patients were included. g Gene expression dot plot of ligand–receptor expression. Left, expression in nPCs/myeloma tumor cells; right, expression in BME cell types. Only interactions that increased between myeloma tumor cells and immune cell types as compared to nPCs are shown. Interactions involving inflammatory cytokines are highlighted.

Fig. 5. T-cell heterogeneity in RRMM patient vs. healthy donor samples.

a UMAP embedding of subclustered T-cell populations colored by cell type of the combined RRMM/healthy data set. b UMAP embedding split and colored by RRMM/healthy status. c Heatmap showing averaged gene expression levels of T-cell population marker genes in the combined data set. d Gene expression dot plot showing exhaustion signature score levels in T-cell subpopulations in RRMM. e Volcano plot of differentially expressed genes using two-sided Wilcoxon rank-sum test in γδ T-cells (RRMM vs. healthy). Thresholds for differential expression were p-value < 0.05 (Bonferroni-adjusted) and logFC > 0.1. f Heatmap of clustered average gene expression and exhaustion signature score of γδ T-cells (RRMM vs. healthy). Selected genes are indicated to the right and signature expression levels of IFN response and ribosomal genes are shown at the bottom. Only samples with >60 profiled γδ T-cells were included. g Scatterplot of average effector score and average exhaustion score in γδ T-cells across patients. Regression line and 95% confidence interval are shown with data points colored according to iFISH status of t(11;14). The Pearson’s correlation coefficient was R = −0.63 (p = 0.0087). h Changes of gene expression in γδ T-cells and CD8⁺ cytotoxic T-cells upon treatment with IMiD. Top, violin plot of exhaustion signature score of donor BM1 (n = 369 cells) and patients RRMM01 (n = 65/425 cells) or RRMM07 (n = 130/95 cells) pre- and post-treatment. Pairwise Bonferroni-adjusted p-values from a two-sided Wilcoxon rank-sum test are shown. Box plot: center line, median; box limits, first and third quartile; whiskers, minimum/maximum values. Bottom, heatmap of average gene expression levels of genes of the exhaustion signature.

Fig. 6. CD16⁺ monocyte heterogeneity in RRMM.

a UMAP embedding of subclustered CD16⁺ monocytes colored by subtype. b UMAP embedding showing expression levels of major marker genes of subtypes indicated in panel (a). c Violin plot of selected marker genes of subtypes in panel (a) split into RRMM vs. healthy. d UMAP embedding as shown in panel (a) split and colored by donor status (RRMM/healthy). e Beeswarm plot for the comparison of CD16⁺ monocyte subtype fractions between RRMM (n = 19) and healthy (n = 8) individuals with Bonferroni-adjusted p-values from a two-sided Wilcoxon rank-sum test. Red center line indicates mean. f Heatmap showing selected differentially expressed genes between RRMM-enriched CD16⁺ monocyte subtypes. g Network plot of predicted cellular interactions between immune cell subsets in RRMM. Every cell type is connected to its 4 top interacting cell types based on the sum of interaction strengths. The node size corresponds to the number of connected cell types and the coloring corresponds to the interaction strength. h Gene expression dot plot of ligand/receptors of RRMM-enriched CD16⁺ monocyte subtypes (left) and associated interaction partners in the RRMM BME (right). Shown are selected interactions that are stronger in the IM/TAM3 subtypes as compared to TAM2 and TAM3.

Fig. 7. BME changes in +1q RRMM.

a Changes in cell-type composition in dependence of +1q. The bar plot in the center shows log₁₀ Bonferroni-adjusted p-values from a two-sided Wilcoxon rank-sum test of differences in cell type fractions between +1q ND/rare vs. subclonal/dominant groups. The box plots display the comparison of NK^dim cells (left) and TAM3 (right) subtype fractions between patients with +1q ND/rare (n = 7) vs. subclonal/dominant (n = 12). Box plot: center line, median; box limits, first and third quartile; whiskers, minimum/maximum values. b UMAP embedding of subclustered NK cells colored by subtype. c Gene expression dot plot for main marker genes of NK subtypes. d UMAP embedding showing GZMB and PRF1 expression levels (left) and point-density plot split in +1q ND/rare and +1q-subclonal/dominant groups (right). e Beeswarm plot for the comparison of the NK^dim effector cells fraction between patients with +1q ND/rare (n = 7), subclonal (n = 5), and dominant (n = 7). Pairwise Bonferroni-adjusted p-values from a two-sided Wilcoxon rank-sum test are indicated. Red center line: mean. f Scheme of transcriptional changes and altered cellular interactions in RRMM with +1q-specific changes colored in red.