Identification of resistance pathways and therapeutic targets in relapsed multiple myeloma patients through single-cell sequencing

Abstract

Multiple myeloma (MM) is a neoplastic plasma-cell disorder characterized by clonal proliferation of malignant plasma cells. Despite extensive research, disease heterogeneity within and between treatment-resistant patients is poorly characterized. In the present study, we conduct a prospective, multicenter, single-arm clinical trial (NCT04065789), combined with longitudinal single-cell RNA-sequencing (scRNA-seq) to study the molecular dynamics of MM resistance mechanisms. Newly diagnosed MM patients (41), who either failed to respond or experienced early relapse after a bortezomib-containing induction regimen, were enrolled to evaluate the safety and efficacy of a daratumumab, carfilzomib, lenalidomide and dexamethasone combination. The primary clinical endpoint was safety and tolerability. Secondary endpoints included overall response rate, progression-free survival and overall survival. Treatment was safe and well tolerated; deep and durable responses were achieved. In prespecified exploratory analyses, comparison of 41 primary refractory and early relapsed patients, with 11 healthy subjects and 15 newly diagnosed MM patients, revealed new MM molecular pathways of resistance, including hypoxia tolerance, protein folding and mitochondria respiration, which generalized to larger clinical cohorts (CoMMpass). We found peptidylprolyl isomerase A (PPIA), a central enzyme in the protein-folding response pathway, as a potential new target for resistant MM. CRISPR-Cas9 deletion of PPIA or inhibition of PPIA with a small molecule inhibitor (ciclosporin) significantly sensitizes MM tumor cells to proteasome inhibitors. Together, our study defines a roadmap for integrating scRNA-seq in clinical trials, identifies a signature of highly resistant MM patients and discovers PPIA as a potent therapeutic target for these tumors.

Extended Data Fig. 1. Sorting strategy of plasma cells and overview single cell data collection from the patients.

a-c, Flow cytometry plots showing sorting strategy (CD38⁺CD138⁺) for plasma cells after doublet exclusion from 3 representative patients. Plots were generated using FlowJo software. (Methods). d, Schematic diagram showing the statists of patient enrollment to KYDAR clinical trial and their current status.

Extended Data Fig. 2. Overview of quality control and data analysis for the single bone marrow plasma cells.

a, Dot plots showing number of reads, number of UMIs, and percentage of cells analyzed per batch of 384 cells (that were pooled for library construction) for all CD38+CD138+ bone marrow single cells from 66 participants (11 control donors, 15 NDMM patients, and 40 PRMM patients). (Methods). b, Heat map showing clustering analysis of 3,862 ‘contamination’ cells that pass quality control but do not express key plasma genes. Shown are representative genes of non-PC. c, Heat map showing clustering analysis of 51,297 PC sorted from all the 66 participants featuring normalized single cell expression levels of a selected set of most variable genes. Clustering is performed using 2,038 genes as features. d, Heat maps showing clustering analysis of PC from patient KYDAR 21 (left) together with the same number of PC from healthy control group (right) using 1,393 gene features.

Extended Data Fig. 3. The robustness of differential expression z-score.

a, comparison of standard log2 fold change vs. our z-score for patient KYDAR 24 clone 1 relative to healthy control. Size of the dots corresponds to the mean expression (UMIs) in the healthy control samples. b, Scatter plot showing the -log₁₀ p-value from the two-sided Mann-Whitney test and the z-score. c, scatter plot comparing our z-score to variance/mean showing no clear bias with high z-score and cell-cell variability. d, z-score scatter plots comparing different sub-sampling of cells (100, 400 and 800 cells). e, z-score scatter plots comparing different sub-sampling of UMIs per cell (300, 400 and 500 cells). f, The distribution of the z-score from patient KYDAR 24.

Extended Data Fig. 4. Differential expression and functional enrichment analyses comparing NDMM and PRMM patients.

a, Volcano plot showing a selected set of differential genes between NDMM and PRMM patients. X axis showing the difference between the average z-score of NDMM patients and PRMM patients; Y axis showing the -log₁₀ p-value from a bootstrap test using t-statistic (Methods). b, Heat map showing the correlation map of all the differential genes identified between the NDMM and PRMM patients. Genes are clustered using hierarchical clustering. c, Functional enrichment of the 3 gene modules using Metascape. d, Graphical illustration of non-responder overall PI resistant pathway.

Extended Data Fig. 5. Clinical response to DARA-KRD treatment.

a, Distribution of international myeloma working group (IMWG) response grades in KYDAR trial. b, Waterfall plot shows the best overall responses according to the IMWG in MM. c,d, Kaplan-Meir curve showing the progression free survival (panel c) and overall survival (panel d) of PRMM patients treated with DARA-KRD.

Extended Data Fig. 6. Differential expression analysis of DARA-KRD resistant patients.

a, Volcano plot showing a selected set of differential genes between PRMM responder group and PRMM non-responder group. X axis showing the difference between the average z-score of PRMM responder group and PRMM non-responder group; Y axis showing the log₁₀ p-value from a bootstrap test using t-statistic (Methods). b, Scatter plot on the left showing the overlap between gene module 1 in Fig. 2a and resistance signature 1 in Fig. 4a. X axis showing the difference between the average z-score of NDMM patients and PRMM patients; Y axis showing the difference between the average z-score of PRMM responder group and PRMM non-responder group. Scatter plot on the right showing the same as the left one, but for gene module 3 and resistance signature 2. c, Venn diagram showing the overlap between the PRMM patient groups classified in Fig. 2a and PRMM responder and non-responder groups.

Extended Data Fig. 7. Longitudinal data collection of PRMM patients.

a, Representative flow cytometry plots showing sorting strategy (CD38+CD138+) for post treatment patients’ plasma cells after doublet exclusion, same as in Extended Data Fig. 1a, but for the samples from patient KYDAR 10 after 4 cycles DARA-KDR treatment. b, Same as in panel a, but for patient KYDAR 10 after 10 cycles DARA-KDR treatment. c, Heat map depicting the relative frequency of the immunoglobulin sequences for each clone (including 11 control donor and 15 newly diagnosed patients); left, bottom, immunoglobulin heavy chain constant region (IGHC); middle, immunoglobulin light chain constant region (IGKC and IGLC); right, immunoglobulin light chain variable region (IGKV and IGLV).

Extended Data Fig. 8. Longitudinal single cell analysis of MM patients pre- and post-treatment.

a, Upper part showing 2D projection of malignant PC of patient KYDAR 21, including healthy PC as reference, highlighting all cells (left), cells collected from baseline (middle left), cells collected from Cycle 4 (middle right), cells collected from Cycle 10 (right). The color of the cells represents healthy PC, and malignant clone 1, 2, and 3. Lower part showing 2D projection of expression of a selected set of differential genes over the metacell model. b, Heatmap showing clone composition in each sample from the patient KYDAR 21 and expression of selected differential genes in each clone. (c) and (d) showing the same as panel (a) and (b), but for patients KYDAR 34 with corresponding samples and clones. (e) and (f) showing the same as panel a and b, but for patients KYDAR 30 with corresponding samples and clones. g, The dynamics of S100A11 (left) and S100A10 (right) expression in individual patients during DARA-KRD treatment.

Extended Data Fig. 9. CsA potentiates Carfilzomib and induced cell death in myeloma cell lines.

a, Graphical illustration of PRMM non-responder resistant pathway, highlighting PPIA. b, Bar plots showing normalized RNAseq PPIA gene transcription in MM cell lines, highlighted in red lines − RPMI-8226 and U266 cell lines (Methods). c, Barplot showing PPIA mRNA levels measured by qPCR in RPMI-8226 PPIA+/+ and PPIA-/-. d, Western blot of PPIA protein levels in RPMI-8226 PPIA+/+ and PPIA-/- cell lines. Shown are cropped western blots of PPIA protein levels in RPMI-8226 PPIA+/+ and PPIA-/- cell lines (tubulin is used as control). Full image can be found in raw data section. Western Blot experiment was conducted once with 3 successful replicates per cell line group. e, Dose-response curves to determine the half maximal inhibitory concentrations (IC₅₀) of CsA in RPMI-8226 and U266 cell lines (Methods). The dose-response curves and IC₅₀ values were determined by a non-linear regression fit. Results are shown as the mean±SEM (Standard Error of the Mean) of 3 independent experiments. f, Immunofluorescence staining of RPMI 8226 (upper) and U266 (lower) cell lines post treatment with Carfilzomib or in combination with CsA for 48h. DAPI (Blue) staining for control live cells, Propidium Iodide (Red) staining for late apoptosis (death), FITC (Green) staining for early apoptosis, and merged gate showing all staining. Representative images from two independent experiments are shown. g, Dose-response curves to determine the half maximal inhibitory concentrations (IC₅₀) of CsA in RPMI-8226 PPIA+/+ and PPIA-/- cell lines (Methods). Results are shown as the mean±SEM (Standard Error of the Mean) of 3 independent experiments.

Extended Data Fig. 10. Cyclosporin A, known inhibitor for PPIA shown synergistic effect with the proteasome inhibitor carfilzomib.

a, Representative flow cytometry plots showing sorting strategy for single RPMI-8226 cells. b, 2D projection showing metacell analysis of scRNA-seq of RPMI-8226 cells during CsA, CFZ, and CFZ+CsA treatment at 4 h and 8 h. Different cell states are indicated by colors. c, 2D projection showing the distribution of the cells (down sampled to 299 cells) in each sample over the metacell map in panel b. d, heatmap showing z-score of differentially expressed genes in each cell state. e, gene ontology enrichment of the 115 up regulated genes and 78 down regulated genes. f, Heatmap depicting the z-score of resistance signature 1 and 2 (Fig. 4a) genes in 15 NDMM, 21 PRMM responder, and 7 PRMM non-responder patients, together with ex-vivo cells of one highly resistant PRMM patient Z01. g, barplot showing scores of resistance signature 1 in 15 NDMM, 21 PRMM responder, and 7 PRMM non-responder patients, together with one resistant PRMM patient Z01. h, barplot showing Z-scores of PPIA expression in the same NDMM and PRMM patients in panel g, together with the resistant PRMM patient Z01.

Fig. 1. Single-cell atlas of PCs derived from NDMM and PRMM patients.

a, The 2D projection showing expression profiles of 51,297 QC-positive single BM PCs derived from 60 patients: 11 controls, 15 NDMM and 34 PRMM patients. b, The 2D projection of a selected set of key PC genes and main MM driver genes over the metacell model. c, Barplot depicting percentage of seemingly healthy PCs within total collected PCs among 15 individual NDMM (N) and 34 PRMM (K) patients, with the total number of collected PCs from each patient in the brackets of the labels on the x axis. d, Boxplots showing normalized expression (Methods) of canonical PC genes and common MM driver genes in 11 controls (H), and 15 NDMM (N) and 34 PRMM (K) individual patients. The center line represents the median, the box limits the upper and lower quartiles, and the whiskers 1.5× the interquartile range. e, Barplot depicting malignant PC clonality in individual NDMM and PRMM patients.

Fig. 2. New MM gene modules define subsets of PRMM patients.

a, Heatmap depicting the z-score of 66 differential genes with an FDA-adjusted P <0.05 (statistical test, see Methods) between the malignant PCs of NDMM and PRMM patients. b, Dot plots showing average expression of a selected set of highly differential genes between NDMM and PRMM patients, in four different patient groups: healthy, NDMM patients, PRMM group 1 and PRMM group 2. The center line represents the mean, the box limits the 95% confidence intervals, and the points all the data. c, The 2D projection of the same patient groups in b over the metacell model. Individual patients are located in the mean coordinates of their corresponding cells. d, The 2D projection of gene module scores of the three gene modules in a over the metacell model. Individual cells are represented by small dots, whereas metacells are represented by bigger dots proportional to the cell number of each metacell.

Fig. 3. Gene signatures of PRMM patients can predict clinical response to DARA treatment.

a, KM plot and analysis of PFS comparing patients who achieved PR and better with patients who did not achieve PR (logrank test, two-sided P = 0.0039, HR = 5.7, 95% CI = 1.5−21.7). b, KM plot and analysis as in a for OS (logrank test, two-sided P = 0.068, HR = 3.2, 95% CI = 0.9−11.9). c, KM plot and analysis of PFS for patients with double high-risk cytogenetic aberration on FISH analysis compared with single or no high-risk aberration (logrank test, two-sided P = 0.00038, HR = 5.1, 95% CI = 1.9−14.1). d, KM plot and analysis as in c for OS (logrank test, two-sided P = 0.033, HR = 2.9, 95% CI = 1.1−7.9). e, KM plot and analysis of PFS for patients with high module-1 genes score compared with low module-1 score (logrank test, two-sided P = 0.00012, HR = 6.4, 95% CI = 2.2−18.6). f, KM plot and analysis as in e for OS (logrank test, two-sided P = 0.0036, HR = 2.9, 95% CI = 1.02−8.2). g, KM plot and analysis comparing PFS of patients with no double hit and low module-1 score (00, blue), patients with double hit and low module-1 score (10, red), patients with no double hit and high module-1 score (01, green) and patients with double hit and high module-1 score (11, orange) (logrank test, two-sided P = 0.000004). h, KM plot and analysis as in g for OS (logrank test, two-sided P = 0.022). i, Histogram of the log₂(average module-1 gene expression) in the MMRF’s CoMMpass data for 783 NDMM patients (upper panel), 69 patients after first line of treatment (middle panel) and 56 patients after at least two lines of treatment (lower panel). j, KM plot and analysis for PFS comparing NDMM patients in the CoMMpass data with high module-1 score (red) and patients with low module-1 gene score (blue) (logrank test, two-sided P = 8.12 × 10⁻⁷, HR = 2.9, 95% CI = 1.9−4.5). k, KM plot and analysis as in j for OS (logrank test, two-sided P = 4.57 × 10⁻¹⁷, HR = 3.9, 95% CI = 2.22−6.87). l, Multivariate KM plot and analysis for PFS of CoMMpass data with module-1 high and FISH double-hit genetic risk (logrank test, two-sided P = 3.72 × 10⁻¹⁵). m, KM plot and analysis as in l for OS (logrank test, two-sided P = 1.8 × 10⁻³¹).

Fig. 4. Molecular pathways of broadly resistant MM patients.

a, Heatmap depicting the z-scores of 133 differential genes in NDMM and PRMM patients. The gene list is obtained by comparing the malignant PCs of the PRMM responder and nonresponder patient groups with an FDA-adjusted P < 0.05 (for statistical test, see Methods). b, Dot plots showing expression of a selected set of highly differential genes between PRMM responder and nonresponder patients. The center line represents the mean, the box limits 95% CIs and the points all the data. c, GO enrichment analysis of resistance signature 1 on the left and signature 2 on the right. d, The 2D projection of resistance signature scores over the metacell model. Individual cells are represented by small dots, whereas metacells are represented by bigger dots proportional to the cell number of each metacell. e, Boxplot showing the response score of individual cells from 24 individual patients using a shallow neural network classifier predicting the response to treatment after 4 months. The center line represents the median, the box limits the upper and lower quartiles, and the whiskers 1.5× the interquartile range. PD, progressive disease.

Fig. 5. Longitudinal single-cell analysis of MM patients pre- and post-treatment.

a, Graphical illustration of clonal evolution and dynamics of malignant PCs into three main trajectories during DARA-KDR treatment: sensitive, resistant and selection. b, Barplots depicting longitudinal malignant PC clonality in pre-treatment (baseline) and post-treatment (cycle 4, 3 months; cycle 10, 9 months) of individual PRMM patients. Healthy PCs are shown in green, whereas different clones of malignant PCs are shown in different shades of gray. Patient outcome at 4 months is indicated on the right of the patient labels. c, Top: 2D projection of malignant PCs of patient KYDAR 24, including healthy PCs as a reference, highlighting all cells (left), cells collected from baseline (middle) and cells collected from cycle 4 (right). The color of the cells represents healthy PCs, and malignant clones 1 and 2. Bottom: 2D projection of expression of a selected set of differential genes over the metacell model. d, Heatmap showing expression of selected differential genes between the two malignant PC clones of patient KYDAR 24, with healthy PCs as a reference. e, The dynamics of CD38 expression in individual patients during DARA−KRD treatment. f, Boxplot showing the response scores of individual cells of each malignant PC clone in 24 individual patients using a shallow neural network classifier predicting the response to treatment after 4 months. The center line represents the median, the box limits the upper and lower quartiles, and the whiskers 1.5× the interquartile range.

Fig. 6. CsA, known inhibitor for PPIA shows synergistic effect with the proteasome inhibitor CFZ.

a, Dose-response curves for RPMI-8226 PPIA^+/+ and RPMI-8226 PPIA^−/− cell lines to determine sensitivity to the first-line proteasome inhibitor, CFZ. Cells were incubated with CFZ at serial dilution, starting at doses of 125 nM for 48 h before cell viability was determined using an MTS proliferation assay. Results are shown as the mean ± s.e.m. of three independent experiments. IC₅₀ values are 34 nM for RPMI-8226 PP/A^+/+ and 14 nM for RPMI-8226 PPIA^−/−. P values of the statistical comparison between IC₅₀ values are indicated (P< 0.001, extra sum-of-squares F test). b, Dose-response curves for RPMI-8226 and U266 cell lines. Cells were incubated with CFZ with or without CsA (1.5 μg) for 48 h before cell viability was determined using an MTS proliferation assay. Results are shown as the mean ± s.e.m. of three independent experiments. IC₅₀ values for CFZ versus combination therapy are 28 nM and 3.8 nM in RPMI-8226 cell line (left), and 31 nM and 5.6 nM in U266 cell line (right), respectively. Results are shown as the mean ± s.e.m. of three independent experiments. P values of the statistical comparison between IC₅₀ values are indicated (P < 0.001, extra sum-of-squares F test). c, FACS apoptosis assay of RPMI-8226 and U266 cell lines using annexin and DAPI staining for early and late apoptosis, respectively. d, Immunofluorescence analysis and quantification of proliferating cells. Statistical significance was determined using exact Wilcoxon’s rank-sum test (Mann-Whitney U-test, *P = 0.0286, two sided; n = 4 per group; results are shown as the mean ± s.e.m.). e, The 2D projection showing metacell analysis of scRNA-seq of RPMI-8226 cells during CsA, CFZ and CFZ + CsA treatment at 4 h and 8 h. Different cell states are indicated by colors. f, Barplot showing the cell state composition in each culture at 4 h and 8 h. g, Barplot showing normalized expression (Methods) of selected genes in each cell state. The results are shown as the mean±s.d. h, FACS apoptosis assay using annexin and DAPI staining for early and late apoptosis, respectively, on ex vivo cultured, patient malignant PCs, treated with either CFZ or CFZ + CsA. i, Barplot showing the quantification of live, early apoptotic and late apoptotic cells. Statistical significance was determined using exact Wilcoxon’s rank-sum test (Mann-Whitney U-test); *P = 0.0286, two sided; n = 4 replicates from patient Z01).