Identification of Therapy-Induced Clonal Evolution and Resistance Pathways in Minimal Residual Clones in Multiple Myeloma through Single-Cell Sequencing

The advent of novel therapeutics, such as proteasome inhibitors and immunomodulatory drugs, has significantly enhanced survival outcomes for patients with multiple myeloma (MM) over the past two decades (1). Presently, the treatment objective for MM is transitioning from achieving complete response, traditionally assessed using conventional biochemical tests, to realizing undetectable minimal residual disease (MRD−) status, monitored via highly sensitive next-generation flow (NGF) or next-generation sequencing technologies (1, 2). A wealth of evidence substantiates the correlation between reduced disease progression or mortality in MM patients and the attainment of MRD− (3). Particularly for high-risk MM patients, achieving and maintaining MRD− status is crucial to counterbalance the negative prognostic implications of adverse cytogenetics (4).

However, recent studies exploring MRD dynamics have highlighted that only a subset of patients [those who have undergone autologous stem-cell transplant (5) or received continuous maintenance treatment (6)] experience a transition from detectable MRD (MRD+) to MRD−. This underscores the urgency of characterizing the biological properties of persistent MRD tumor cells. Such efforts can unveil new insights into the chemoresistance mechanisms of these cells and potentially disclose hitherto unknown vulnerabilities that can be exploited to eradicate persistent MRD tumor cells (4, 7).

A high degree of intraclonal heterogeneity is a key feature of MM. Patients with newly diagnosed MM (NDMM) invariably harbor multiple clones with distinct genetic backgrounds (8). Moreover, treatments exert varying impacts on intraclonal diversity, leading to the elimination of sensitive clones and the selection of resistant ones (7, 8). Our current understanding of the genetic and nongenetic biology of MRD tumor cells has been predominantly gleaned from bulk-level analysis (9–11). The advent of single-cell technologies, especially multi-omics single-cell technologies, has offered us a novel perspective on understanding the clonal evolution of MM (12). A recent single-cell study has revealed that high-risk major clones observed at relapse are already present at diagnosis as minor, undetectable subclones (13). Moreover, a subgroup of plasma cells (PC) exhibiting 1q alterations, TP53 mutations, and elevated expression of the PHF19 has been identified to be associated with the progression of MM (14). By integrating single-cell RNA sequencing (scRNA-seq) with various sequencing platforms, co-evolution of tumor cells and the tumor microenvironment (TME) has also been explored during the progression of MM (15).

However, two main limitations exist in previous MM clonal evolution studies. First, the focus has largely been on comparing clonal architecture at diagnosis and relapse using paired samples (16), with only a few studies investigating clonal evolution occurring early in MRD tumor cells. Second, clonal evolution has been primarily explored from a genetic perspective, including the acquisition of new genetic abnormalities or clonal competition and expansion through a Darwinian process (16, 17). Conversely, the role of transcriptional, metabolic, or epigenetic adaptation in clonal evolution has received comparatively limited attention. Despite the widespread acceptance that a tumor cell, or a group of tumor cells with a specific genetic alteration, will gain a competitive advantage to evade therapeutic pressure (18), a growing number of scRNA-seq studies are now focusing on the nongenetic determinants of malignant clonal fitness (19). As such, clonal evolution is progressively transitioning into a multidimensional concept, encompassing both genetic and nongenetic determinants that drive the therapeutic adaptation of clones. Despite this, the relationship between genetic and nongenetic resistance mechanisms in MM remains largely uncharted.

To fill these knowledge gaps, our study examined the transcriptional dynamics of PCs from diagnosis through to posttreatment. We collected and sequenced a total of 38 paired bone marrow (BM) samples. Our findings revealed three primary trajectories of transcriptional evolution posttreatment: clonal elimination in MRD− patients and clonal stabilization and clonal selection in MRD+ patients. Additionally, we observed disparate resistance pathways in residual malignant PCs (mPC). A comparison of the genetic and transcriptional dynamics unveiled a strong correlation between genetic and nongenetic factors driving the clonal evolution of MM. Lastly, we identified variations in cellular interactions between mPCs and the TME, with selective PCs exhibiting the strongest interactions with the TME.

Taken together, our study offered a deeper understanding of the clonal evolution and resistance mechanisms of MM. By examining these processes at a single-cell resolution on a large scale, we confirmed the link between genetic and nongenetic mechanisms driving clonal evolution, providing insights into potential therapeutic strategies aimed at disrupting the connections between MRD tumor cells and the TME.

BM samples were collected from patients with monoclonal gammopathy of undetermined significance (MGUS; n = 1), MM patients at diagnosis (n = 19), and a single healthy donor (HD). Paired BM aspirates were acquired after two to five cycles of induction therapy from 17 MM patients. The detailed clinical and demographic characteristics of the MGUS and MM patients are presented in Supplementary Table S1. All samples were obtained at the Institute of Hematology & Blood Diseases Hospital, Chinese Academy of Medical Sciences & Peking Union Medical College in Tianjin, China.

BM mononuclear cells (BMMC) were isolated from the BM aspirates using the density gradient centrifugation method with Ficoll-Paque PLUS lymphocyte isolation sterile solution (Cytiva, Sweden). For patients MM1 to MM8, CD138⁺ PCs were isolated using magnetic-activated cell sorting with CD138⁺ magnetic beads (Miltenyi Biotec, Paris, France). As indicated by several previous studies and ours, a ≥90% PC purity could be achieved after cell sorting (20–22). All samples were viably cryopreserved in dimethyl sulfoxide at a final concentration of 10%. This study was approved by the local institutional ethics committees under the oversight of the Institute of Hematology & Blood Diseases Hospital, Chinese Academy of Medical Science & Peking Union Medical College. In accordance with the guidelines of the Declaration of Helsinki, all patients and the HD provided informed consent prior to enrollment.

MRD was assessed at the time of paired posttreatment BM sample acquisition. Out of the 17 MM patients with paired posttreatment BM samples, 14 underwent MRD assessment at this defined time point. MRD was evaluated in accordance with the NGF method developed by EuroFlow (23, 24). The limit of detection achieved by NGF was calculated for each sample using the formula: (20/number of viable nucleated cells) × 100. A minimum of 2,000,000 nucleated cells were recorded, resulting in a sensitivity threshold ranging between 10⁻⁵ and 10⁻⁶. MRD detection was defined as when the percentage of phenotypically aberrant PCs was equal to or greater than the limit of detection achieved in the respective sample.

Frozen BM cells were rapidly thawed, washed, counted, and resuspended in PBS and 0.04% BSA to a final concentration of 1,000 cells per μL. In this study, samples were sequenced in two batches. Samples MM1 to MM8 constituted the first batch, whereas samples MM9 to MM19 comprised the second batch. The single-cell 3′ Library Kit V3 (10× Genomics, CA) was used for the cDNA synthesis of individual cells. These cells were then barcoded and ready for library construction following the manufacturer’s instructions. Upon preparation, the indexed cDNA libraries underwent sequencing on an Illumina NovaSeq 6000 platform using a paired-end mode. Novogene Co., Ltd. (Beijing, China) performed the sequencing according to the manufacturer’s guidelines, targeting an approximate sequencing depth of 20,000 reads per cell.

The processing of the raw sequencing data marked the initial stage, which involved alignment to the GRCh38 (UCSC) human reference genome utilizing the CellRanger pipeline (version 3.0.2, 10× Genomics, CA). Default parameters were employed for the alignment process. Postalignment, quality control procedures were conducted on the datasets using the Seurat package (version 3.1.1; ref. 25). The first step was the removal of doublet cells using DoubletFinder (version 2.0.3; ref. 26). Subsequently, low-quality cells were filtered out according to the following criteria: (i) numbers of detected unique molecular identifiers (UMI) and genes (log₁₀ scaled) that fell outside the median ± 3 median absolute deviation range and (ii) mitochondrial gene proportions below the median +3 median absolute deviation. Quality control was applied separately to CD138⁺ PCs and BMMCs, considering factors such as UMI count, gene count, and the mitochondrial ratio. After the exclusion of low-quality cells, the two filtered datasets were merged. Genes that were expressed in at least five cells were selected for downstream analysis.

After the elimination of low-quality single cells, the raw expression matrix was normalized using the NormalizeData function. The ScaleData function was then applied to mitigate variations in gene expression from cell to cell. Principal component analysis was used for linear dimensionality reduction, targeting the top 2,000 variable genes. To identify the appropriate number of principal components for subsequent clustering, we used the ElbowPlot function. Next, we applied the FindNeighbors function to calculate shared nearest neighbor graphs, and the Louvain algorithm was used, alongside Uniform Manifold Approximation and Projection for Dimension Reduction (UMAP) embedding in two-dimensional space, to cluster the single cells. For the clustered analysis of all PCs depicted in Fig. 1E, we did not perform any integration preprocessing. For all other cells, we applied the Harmony package (version 1.0; ref. 27) for integration prior to analysis, using the following parameters: Lambda = 1, theta = 1.

Cell types were determined by evaluating the average expression levels of established markers. These included SDC1, IGHG1, and MZB1 for PCs; CD79A, CD79B, and MS4A1 for B cells; S100A8, S100A9, and LYZ for myeloid cells; GNLY, NKG7, CD3D, and CD3E for NK cells and/or T cells; CD34 for hematopoietic stem and progenitor cells; and HBA1, HBA2, HBB, and PF4 for erythrocytes and/or megakaryocytes, as previously described (28–32).

Following this, PCs were isolated based on the expression of known marker genes such as SDC1, MZB1, and XBP1. A reclustering process was then implemented. The methods described above were used to cluster PCs into 22 subsets from both HD and patients. Malignant PCs from each patient, both pretreatment and posttreatment, were integrated for subsequent analysis.

In order to identify differential expressed genes (DEG), we used the FindMarkers or FindAllMarkers functions, using normalized data and setting parameters such as test.use = “wilcox” and logfc.threshold = log1.5. P values were adjusted using the Bonferroni correction method, which takes into account the total number of genes in the dataset. DEGs with adjusted P values below 0.05 were considered significant and chosen for subsequent analysis.

For the differential analysis depicted in Fig. 3, we used a pct threshold of 0.6 to identify DEGs across patients. For other differential analyses, a min.pct threshold of 0.3 was used. These DEGs were then uploaded to the Metascape (33) website (https://metascape.org), and default parameters were used to perform pathway enrichment analyses on the 50 hallmark pathways (34).

Copy number alteration (CNA) profiles in individual plasma/myeloma cells were identified using the InferCNV package (version 1.7.1), which calculates the average expression of adjacent genes over large genomic regions. PCs derived from HDs were used as the reference, and we profiled the CNA profile on the PCs of every patient individually.

To identify somatic large-scale copy number variants (CNV), the gene count matrix was used as the input with the following parameters: cluster_by_groups = False, cutoff = 0.1, analysis_mode = “subclusters,” and BayesMaxPNormal = 0.25. This allowed us to identify regions of the genome where there were significant increases or decreases in gene copy numbers, which can often be indicative of cancer-associated genetic alterations.

Developmental potential and cellular differentiation states in individual plasma/myeloma cells per patient were analyzed using the CytoTRACE package (version 0.3.3; ref. 35), using default parameters.

To determine the cell-cycle phase of each cell, we employed the CellCycleScoring function in the R package Seurat (version 3.1.1; ref. 25). This allowed us to assign a cell-cycle phase to each individual cell based on its gene expression profile.

The analysis of the proportion of Lambda+ and Kappa+ PCs was conducted as outlined by Andor and colleagues (36), albeit with some tailored modifications to the algorithm. To ascertain the expression of immunoglobulin Kappa (IGK) and Lambda (IGL) in PCs, IGK and IGL transcripts were quantified. For each individual cell, “m” was designated as the quantity of IGKL transcripts, and “n” as the highest quantity of transcripts from the three functional IGLC isotypes (IGLC2, IGLC3, and IGLC7). We calculated the absolute value of log₂(m/n) and retained cells that exceeded 1 for further examination. Cells were categorized as Kappa+ if m > n, and Lambda+ if m < n. Finally, the proportion of Lambda+ cells was calculated as L/(L + K), with “L” representing the count of Lambda+ cells and “K” the count of Kappa+ cells.

To explore the specific gene regulatory networks affiliated with PCs, we employed the pySCENIC package (version 0.11.2; ref. 37) in conjunction with the cisTarget database. The gene regulatory network of PCs was assembled using default parameters. To identify differentially activated regulons between cell groups, we utilized Wilcoxon rank-sum tests with “wilcox.test” and “p.adjust” for multiple hypothesis correction. Regulators meeting the following criteria were considered significant and selected: fold change >1 between resistant (selective, cycling-resistant group, or other resistant cells) and sensitive cells; adjusted P value < 0.05 with Bonferroni correction; and regulon frequency >0.5 in the resistant cells.

The gene expression from the Myeloma Research Foundation (MMRF) CoMMpass dataset was obtained from the publicly available database The Cancer Genome Atlas. The other three gene expression datasets along with corresponding clinical information of MM patients were obtained from the publicly available database Gene Expression Omnibus (GEO), under the accessing numbers GSE2658, GSE9782, and GSE57317. We then applied the gene set variation analysis package (version 1.48.2; ref. 38) to calculate the enrichment score of the cycling-resistant or selective gene signature for individual patients in each dataset. A univariate Cox proportional hazard regression model was implemented to determine the association of the cycling-resistant or selective gene score with overall survival (OS) in each cohort. HR from the univariate Cox regression analysis was used to determine the protective (HR < 1) and risky genes (HR > 1). In the Cox univariate analysis, the MMRF CoMMpass dataset was set as the test set, whereas the other three GEO datasets served as validation sets. The results of these Cox univariate regression analyses in each dataset were independent of each other.

We analyzed intercellular interactions between cell types using the CellPhoneDB package (version 2.1.4; ref. 39), leveraging ligand–receptor co-expression. The interaction score refers to the total average of the mean expression value of a single ligand–receptor partner in the corresponding interacting cell type. We merged the single-cell dataset of MM14 pretreatment and posttreatment and normalized the raw expression matrix using the NormalizeData function in Seurat. We then used the normalized counts as input for subsequent analysis, employing default parameters. In this study, we focused on the enriched cytokine/receptor interactions in selective PCs or resistant PCs and selected the cytokine/receptor interactions associated with highly significant (P value < 0.05) and highly expressed (exp_value > 0.1) cell–cell interaction pairs in selective PCs compared with sensitive PCs or resistant PCs compared with sensitive PCs.

All statistical analyses and graph generation were conducted using R (version 3.5.1, R Foundation, Vienna, Austria). All statistical tests employed in this study were two-sided. The identification of DEGs was performed using a two-sided Wilcoxon rank-sum test. The Bonferroni correction method was applied to control the false discovery rate and to correct P values for multiple testing. An adjusted P value less than 0.05 was considered statistically significant. The Kaplan–Meier method was utilized for survival data analysis, and differences in survival were evaluated using the log-rank test. HR and 95% confidence interval (CI) was calculated using univariable Cox regression models.

All patients included in this study were enrolled in the National Longitudinal Cohort of Hematological Diseases in China (ClinicalTrials.gov identifier: NCT04645199). In addition, the study complied with the Declaration of Helsinki and was approved by the Ethics Committee of the Institute of Hematology and Blood Diseases Hospital, Chinese Academy of Medical Science & Peking Union Medical College (Certificate: IIT2020023-EC-1).

The scRNA-seq data presented in this study have been deposited to the Genome Sequence Archive in the National Genomics Data Center, under BioProject ID: PRJCA018851 (https://ngdc.cncb.ac.cn/bioproject/browse/PRJCA018851), and Accession ID: HRA005218 (https://ngdc.cncb.ac.cn/gsa-human/browse/HRA005218). Please note that any data usage must fully comply with the Regulations on Management of Human Genetic Resources in China. The scRNA-seq matrix data from the BM of 8 HDs were obtained from the Human Cell Atlas (HCA) database (census of immune cells, https://data.humancellatlas.org/explore/projects/cc95ff89-2e68-4a08-a234-480eca21ce79). The scRNA-seq matrix data of PCs utilized for external validation of normal-like PCs (nlPC) can be accessed from the GEO database with the accession numbers GSE161195 and GSE117156. The bulk RNA-seq data of BM samples, along with clinical data from 858 NDMM samples, were procured from The Cancer Genome Atlas (https://tcga-data.nci.nih.gov/) Multiple MMRF CoMMpass dataset. The bulk RNA-seq data used for survival analysis can be found in the GEO database with the accession numbers GSE2658, GSE9782, and GSE57317. All other relevant data supporting the key findings of this study are available within the article and its Supplementary Information files or from the corresponding author at angang@ihcams.ac.cn upon reasonable request. All noncommercial reagents used in this article are available from the corresponding author at angang@ihcams.ac.cn upon reasonable request.

To comprehensively investigate the transcriptional dynamics of PCs throughout treatment, we leveraged scRNA-seq to dissect the transcriptional heterogeneity, subclonal structure, clonal evolution of transcriptional clusters and genetic subclones, and cellular interactions in MM at diagnosis and posttreatment. Our study sample set included 38 BM samples from patients with MGUS (n = 1), MM patients at diagnosis (n = 19), MM posttreatment (n = 17), and one HD (Fig. 1A; Supplementary Table S1).

We analyzed 161,403 CD138⁺ PCs and BMMCs derived from HD, MM pretreatment, and posttreatment samples, all of which passed stringent quality controls (Supplementary Fig. S1A and S1B; Supplementary Table S2). Additionally, we integrated our dataset with the scRNA-seq dataset from 8 HDs from the HCA (40). After normalizing the gene expression matrix, we identified six major cell subsets in the BM using graph-based clustering (Supplementary Fig. S1C and S1D). We observed a distinct separation of plasma/myeloma cells from immune cells, with notable enrichment of CD138⁺ cells in the MM TME (Fig. 1B–D). In total, we profiled 76,036 PCs (median = 902) and analyzed a number of immune cells per sample ranging from 839 to 10,137, with a median of 3,059 cells (Supplementary Table S3).

To delve deeper into the inter and intrapatient tumor heterogeneity in MM, we carried out unsupervised clustering for samples with a sufficient number of PCs (≥30 PCs). This process yielded 22 clusters of PCs. Interestingly, most PCs derived from HDs were clustered into cluster 1 (C1), a group that also contained PCs from various patients. Consequently, PCs in C1 were identified as HD/nlPCs (Fig. 1E). We observed significant transcriptional differences among patients, with the majority of PC clusters being composed exclusively of cells from a single patient. However, there were some exceptions: C15, for instance, contained PCs from 2 patients with distinct translocations [MM12 with t(6;14) and MM13 with t(4;14); Fig. 1F; Supplementary Table S1]. Furthermore, longitudinal samples from the same patients were grouped into separate clusters (MM5: C2 and C11; MM6: C5 and C19; Fig. 1F).

Despite the considerable interpatient tumor transcriptional heterogeneity, we noticed commonly overexpressed driver genes in mPCs across previously defined MM subgroups (41). The expression of translocation target genes, including CCND1, FGFR3 and NSD2, CCND2, MAFB, MAF, and CCND3, correlated with the translocations detected by interphase fluorescence in situ hybridization (iFISH; Fig. 1G; Supplementary Table S1). In addition, LAMP5, previously reported to be upregulated in MM patients with osteolytic lesions (42), was found to be overexpressed in 9 out of 18 patients (Fig. 1G).

In our present study, we identified a PC cluster (C1) that encompassed PCs from HDs and multiple MM patients. Specifically, 9.8% of PCs from MM2 pretreatment (MM2_pre; 339/3,449), 100% of PCs from MM9_pre, and 98.5% of PCs from MGUS were categorized into this subset (Supplementary Table S4). To further explore this, we initially characterized the clonotypic identity of PCs within this cluster. We analyzed the relative expression of the immunoglobulin κ and λ light chain variable regions (IGKV and IGLV, respectively), coupled with the immunoglobulin heavy chain variable region for every individual cell. Our findings revealed diverse immunoglobulin sequences from HD PCs. However, we noticed a specific immunoglobulin clonotype for some MM patients, including MM2_pre, MM3 posttreatment (MM3_post), and MM9_pre (Supplementary Fig. S2A–C). To further validate the monoclonal nlPCs in C1, we next determined the light chain restriction of PCs in each sample (see Supplementary Methods; ref. 36). For HD PCs, the frequency of immunoglobulin Lambda+ PCs ranged from 35.7% to 46.2%. For mPCs from each MM patient, the frequency of immunoglobulin Lambda+ PCs fell either below 2.3% or exceeded 85.5%. Conversely, nlPCs from MM2_pre, MM3_post, and MM9_pre exhibited an exclusive light chain assignment of 25.8%, 28.0%, and 1.4%, respectively (Supplementary Fig. S2D). The difference in light chain restriction between MM2_pre and MM9_pre suggested the potential admixture of normal cells within the nlPCs of MM2_pre. Collectively, this light chain restriction for MM9_pre supported our hypothesis that, despite clustering with HD PCs, at least a portion of the MM PCs within the seemingly healthy cluster were clonal PCs.

To further understand the transcriptional differences between nlPCs and mPCs, we subsequently performed unsupervised clustering on the nlPCs and mPCs from MM2_pre. Our results unveiled more transcriptional clusters than CNV clones. Moreover, nlPCs and mPCs were distinctly separated, whereas the transcriptome-based clustering displayed a high degree of consistency with CNV clones and cell groups (Fig. 2A). Further analysis revealed unique CNA profiles among three CNV subclones, with S1 and S3 demonstrating similar CNA profiles, and S2 exhibiting a seemingly normal CNA profile (Fig. 2B). In addition, cellular differentiation states inferred using CytoTRACE indicated a less differentiated state in mPCs and a more differentiated state in nlPCs (Fig. 2C). Further analysis, which included HD PCs, indicated an intermediate differentiated state of nlPCs (Supplementary Fig. S2E). DEG analysis in nlPCs and mPCs from MM2_pre compared with HD PCs revealed 55 shared upregulated DEGs. These genes included some previously described MM-specific targets, such as SLAMF7 (43) and MCL1 (Fig. 2D; Supplementary Table S5; ref. 44). Pathway enrichment analysis for the shared upregulated DEGs demonstrated significant enrichment for the NF-κB pathway (Fig. 2E). Conversely, compared with nlPCs, mPCs showed upregulated expression of endoplasmic reticulum and unfolded protein response (UPR) pathway gene PPIB, mitochondrial stress genes COX7C and COX6B1, and a previously reported MM-related resistance gene MIF (Fig. 2E; Supplementary Table S5). Furthermore, activation of the MYC signaling pathway was observed in mPCs (Fig. 2F). Selected MM-related driver genes, including ITGB7, LAMP5, KRAS, and BRAF, were expressed at the highest levels in mPCs, whereas nlPCs only exhibited low-level expression of these genes (Fig. 2G). Analysis of all nlPCs grouped in the C1 cluster demonstrated similar seemingly normal CNA profiles (Supplementary Fig. S2F), whereas the expression of MM-related genes was exclusively detected in patient-derived nlPCs (Fig. 2H). Further analysis of nlPCs from MM9_pre, MM3_post, and MGUS also demonstrated similar seemingly normal CNA profiles, whereas mPCs from MM3_post displayed various loci involved by CNV (Supplementary Fig. S3). Therefore, nlPCs might represent a transitional state of PCs, evolving from normal PCs to mPCs.

Finally, we validated our identification of nlPCs using two external scRNA-seq datasets (32, 46). Our analysis confirmed an exclusive heavy chain assignment in MM12 in the external datasets and validated the abnormal expression of MM-related driver genes in nlPCs from several patients (Supplementary Fig. S4). In conclusion, we identified a cluster of nlPCs present in several patients. This group of cells is most frequently seen when patients carry only minor clones of cytogenetic abnormalities (MM2_pre, MM3_post, and MM9_pre) or when patients do not exhibit cytogenetic abnormalities (MM12 in the external validation cohort).

Several distinguishing features were observed in this cluster: It lacked secondary genetic events as shown in the CNA profile and exhibited monoclonal immunoglobulin expression or light chain restriction, and despite clustering with HD PCs at the single-cell level, the nlPCs demonstrated significant transcriptional alterations. Our results suggested that nlPCs represented a transitional state of PCs, transitioning from normal PCs to mPCs. Given the multistep disease progression characteristic of MM, such nlPCs might represent a potential target for early intervention in MM and warrant further investigation.

To gain further insights into the transcriptional evolution and resistance pathways of residual mPCs, we performed unsupervised clustering on the mPCs from each patient separately, specifically focusing on those with sufficient posttreatment PCs (≥30 PCs; Supplementary Fig. S5). In accordance with the NGF MRD data (Supplementary Table S1), for patients who achieved MRD− status (sensitivity threshold at 10⁻⁵) posttreatment (MM4, MM7, and MM11), we observed only a very low percentage of residual mPCs among all PCs. Conversely, for MRD+ posttreatment patients, our results pointed to two distinct patterns: one where 50% of PCs in the BM posttreatment were mPCs and the other where almost all the PCs were residual mPCs (Supplementary Fig. S6A).

However, the percentage of mPCs within the total PC population does not fully reflect the transcriptional changes that occur in mPCs. Thus, we compared the dynamic shifts in mPC subclusters before and after treatment. In doing so, we identified two patterns of transcriptional evolution in 9 MRD+ patients. The transcriptional selection was observed in 4 patients (MM2, MM5, MM6, and MM14), with the major transcriptional cluster pretreatment becoming a minor cluster, whereas a minor cluster evolved into a predominant cluster posttreatment. In contrast, transcriptional stabilization was observed in the remaining 5 MRD+ patients (MM3, MM10, MM16, MM17, and MM19; Fig. 3A; Supplementary Fig. S5).

Thus, we categorized mPCs into three groups based on their response to treatment: sensitive PCs, such as those in clusters C0, C1, C3, and C5 in MM5, which were nearly eliminated following treatments; resistant PCs, which did not respond or only partially responded to treatment, such as those in clusters C0, C1, C2, and C4 in MM19; and selective PCs, which showed a growth advantage and became the predominant cluster posttreatment, including C2 in MM5, C1 in MM6, and C1 in MM14 (Supplementary Fig. S6B and S6C).

Although most of the resistant clusters remained in a state of cell-cycle arrest posttreatment, we identified a cell-cycle–related resistant cluster in 5 out of 7 patients who received bortezomib and lenalidomide in combination with dexamethasone treatment (Fig. 3B). To delve deeper into these cycling-resistant clusters, we carried out a DEG analysis between HD, sensitive, and cycling-resistant PCs. Interestingly, our results demonstrated that cycling-resistant PCs highly expressed several genes previously linked to myeloma drug resistance, such as CKS1B, STMN1, TUBA1B, and TYMS (47–50). We also noted high expression of PHF19, a gene recently implicated in the malignant progression of MM and PC leukemia (14).

The pathways enriched in cycling-resistant cells were associated with UPR, hypoxia, and cell-cycle/proliferation pathways. Notably, in line with recent reports that indicate an upregulation of antioxidant gene programs in cycling persister cells (51), we found that cycling-resistant cells also highly expressed MIF, a key antioxidant gene known to mediate MM cell resistance to proteasome inhibitors (Fig. 3C and D; Supplementary Table S6; ref. 45). However, our analysis of noncycling-resistant PCs highlighted only a few upregulated genes, such as S100A4, TXNDC5, and RPS10 (Supplementary Fig. S6D and S6E; Supplementary Table S7).

Compared with both HD PCs and sensitive PCs, the genes displaying upregulation in selective PCs were found to be significantly enriched in pathways such as NF-κB, as demonstrated by genes like FOSB, JUN, and EGR1. Furthermore, they were enriched in the antiapoptotic pathway, exemplified by genes such as MCL1. The gene PPP1R15A, crucial for controlling cell viability amid protein folding stress (52), was also observed to be upregulated in selective PCs. Notably, this gene has been previously associated with the progression from MGUS to MM (53). Furthermore, CD74, which is the primary surface receptor for MIF, was significantly upregulated in selective PCs (Fig. 3E and F; Supplementary Table S8). As previous studies have suggested that stimulation of surface CD74 contributes to NF-κB pathway activation (54), we thus hypothesized that CD74 played an important role in the interaction between selective cells and the TME.

Next, we investigated TFs that could regulate different drug resistance pathways in cycling-resistant and selective cells. Our analysis indicated that EZH2, E2F1, and E2F2, all of which are cell-cycle/proliferation-regulating TFs, were upregulated in cycling-resistant PCs (Fig. 3G). Conversely, the genes highly expressed in selective PCs were regulated by TFs such as FOS and JUN. It has been reported that the TF AP-1, formed by FOS and JUN, plays a crucial role in enhancing survival but decreasing the proliferation of MM cells (Fig. 3G; ref. 15).

Lastly, we compared the DEGs between selective PCs and cycling-resistant PCs (Supplementary Table S9). We found that selective PCs were enriched for pathways including NF-κB, apoptosis, and P53 pathways, whereas cycling-resistant PCs were enriched for pathways including MYC, E2F target, and UPR pathways (Fig. 3H).

We proceeded to analyze whether the molecular signatures that delineate cycling-resistant PCs (the cycling-resistant signature) and selective PCs (the selective signature) correlate with patient outcomes. We employed gene set variation analysis to calculate the expression of the cycling-resistant and selective signatures in 858 NDMM patients using the independent dataset from the Myeloma Research Foundation (MMRF) CoMMpass dataset. Survival analysis indicated that patients with a high cycling-resistant score had significantly shorter survival times compared with those with a lower score. The median OS was 39.0 months for the high-score group, compared with an unreached median for the low-score group, with an HR of 2.94 and 95% CI of 2.23–3.88 (P < 0.0001; Fig. 4A). However, it was observed that patients with a high selective score had a longer OS compared with those without a high cycling-resistant score, although both medians were not reached, and the HR was 0.56, with a 95% CI of 0.39–0.81 (P = 0.0015; Fig. 4B).

The adverse prognostic value of the cycling-resistant signature was further validated by utilizing three independent datasets (GSE2658, GSE9782, and GSE57317) consisting of 868 MM cases (Fig. 4C). Moreover, our results indicated a slightly longer OS for patients with a high selective score compared with those without a high selective score, although the difference did not reach statistical significance in the log-rank analysis (Fig. 4D).

To gain deeper insights into the interplay between genetic and nongenetic evolution, we applied inferCNV (55) to delineate the clonal architecture of genetic events in patients. A representative patient, MM5, exhibited branching evolution (16), with analysis of pretreatment and posttreatment mPCs identifying six CNV subclones among the eight transcriptional clusters (Fig. 5A). Consistent with iFISH results (Supplementary Table S1), we observed a common 13q deletion (13q-) in all six CNV subclones (Fig. 5B).

To investigate the clonal evolution pattern in MM5, we considered the clonal structure of each CNV subclone. Early genetic events were found in all or most of the CNV subclones, whereas later events occurred in only some of the subgroups. Additionally, ancestral malignant subclones displayed simpler clonal architectures, whereas descendant subclones exhibited more complexity. Specifically, in MM5, 13q- was identified as the earlier cytogenetic event, whereas 1p- was considered a later event. Moreover, we found that S2, S1, S3, S5, and S6 shared similar CNA profiles, with S1, S3, S5, and S6 acquiring some new cytogenetic abnormalities compared with S2. In contrast, the clonal architecture of S4 differed slightly from the S2 branch. Thus, the five main CNV subclones of MM5 at diagnosis could be divided into two branches (S4 and S1/S2/S3/S5/S6; Fig. 5C).

Further analysis revealed a high degree of concordance between transcriptional clusters and genetic subclones, such as C2 (transcriptional cluster) and S4 (genetic subclone), or between C4 (transcriptional cluster) and S6 (genetic subclone; Fig. 5D). This observation suggested a strong link between genetic and nongenetic mechanisms driving the clonal evolution of MM. Similarly, the transcriptional selection observed in MM5 posttreatment could also be interpreted from a genetic perspective. Notably, two groups of tiny subclones at diagnosis, namely, S4 and S6, became the dominant subclones posttreatment, whereas S1, S2, and S3 were virtually eliminated, and S5 was only partly preserved (Fig. 5E).

Viewed through the genetic lens of clonal evolution, S4 (corresponding to C2) emerged as the predominant major CNV subclone posttreatment, followed by S6 (corresponding to C4). Interestingly, MM5 exhibited a significantly shorter survival compared with other patients with MRD+ status, with an OS of only 10 months. We hypothesized that both the nongenetic resistance mechanisms of C2 and C4 contributed to the inferior outcome in MM5. To investigate this further, we inferred the cellular differentiation states using CytoTRACE, revealing a more differentiated state in C2 and a less differentiated state in C4. However, both C2 and C4 displayed higher differentiation compared with other transcriptional clusters (Fig. 5F).

Further DEG analysis between C2 and C4 clusters highlighted significant differences in gene expression patterns. C2 showed high expression of genes related to the “Hypoxia” pathway, such as PPP1R15A and ATF3, whereas C4 exhibited elevated expression of genes involved in “Oxidative phosphorylation,” including ATP6V0B, UQCRB, COX6C, and UQCRH. Additionally, the immunotherapeutic target TNFRSF17 (also known as BCMA) was highly expressed in the C4 cluster (Fig. 5G and H; Supplementary Table S10). These distinct patterns of metabolic adaptation in C2 (S4) and C4 (S6) clusters suggested their involvement in promoting prosurvival pathways posttreatment and potentially offering novel therapeutic opportunities to target residual malignant cells.

Similar genetic selection was observed in MM6 and MM16 (Fig. 5I; Supplementary Fig. S7), whereas MM10, MM16, MM17, and MM19 exhibited genetic stabilization (Fig. 5I; Supplementary Fig. S8). In summary, our findings demonstrated a close link between the genetic landscape and nongenetic mechanisms driving myeloma cell treatment resistance. Moreover, the distinct patterns of metabolic adaptation played crucial roles in driving the clonal selection of residual malignant cells.

To delve deeper into the TME remodeling after treatment, we analyzed 73,424 BMMCs derived from our HD, MM pretreatment, and posttreatment samples (MM09-MM17, MM19), all of which passed stringent quality controls, identifying 17 major cell types (Fig. 6A and B). Notably, we observed significant differences in cellular compositions between the two disease statuses. Consistent with previous studies, MM patients displayed trends toward increased proportions of NK cells, CD16⁺ monocytes, and CD8⁺ memory/effector cells pretreatment compared with HDs, whereas the proportions of B cells were decreased in MM pretreatment (28). Intriguingly, a substantial increase in the proportion of CD14⁺ monocytes was observed in MM posttreatment compared with MM at diagnosis (Fig. 6C).

To understand the functional implications of these changes, we conducted a DEG analysis on CD14⁺ monocytes posttreatment. The results revealed significantly elevated expression of chemokines, such as CCL3, CCL4, and CXCL2, by CD14⁺ monocytes (Fig. 6D; Supplementary Table S11). This heightened expression of chemokines suggested an enhanced recruitment effect of CD14⁺ monocytes on T cells. Conversely, genes highly expressed by CD16⁺ monocyte pretreatment were primarily related to the “Interferon alpha response” pathway, including IFI6, IFL44L, and IFITM1 (Fig. 6E; Supplementary Table S12), indicating that antimyeloma treatment partially mitigated the inflammation in the TME (56).

Taken together, our results indicated that the treatment not only remodeled the composition of the MM TME but also altered the transcriptional status of immune cells, reflecting a complex interplay between the immune system and myeloma progression.

Considering the potential immunosuppressive effects of the TME on providing shelter for MRD cells from antimyeloma treatment (56), we conducted a detailed investigation of the cellular interactions between TME immune cells and residual mPCs using the accumulated ligand–receptor interaction database CellPhoneDB (39). Consistent with previous studies (29), our findings highlighted the strongest interactions between mPCs and myeloid cells, particularly monocytes and dendritic cells (DC; Fig. 6F and G). Further examination revealed even stronger interactions between selective mPCs and TME immune cells, compared with cycling-resistant and/or sensitive mPCs, mediated by APP/CD74, COPA/CD74, FAM3C/HLA-C, MIF/CD74, and SELL/SELPLG interactions (Fig. 6H; Supplementary Fig. S9A). Notably, the most differential interactions we identified involved the upregulated expression of CD74 on the surface of selective mPCs and its corresponding receptors in the TME, suggesting that CD74 overexpression in selective cells might contribute to their immune escape.

Moreover, we found enhanced interactions between selective mPCs and NK cells, mediated by HLA-E/CD94:NKG2A interactions, which represent inhibitory signals to NK cells (Supplementary Fig. S9B). Additionally, HLA-F/LILRB1 interactions, representing an inhibitory signaling axis in the TME (57), were enriched between selective mPCs and myeloid cells (Fig. 6I). These findings raised intriguing questions about the potential involvement of these differential cellular interactions in the observed patterns of transcriptional evolution and distinct resistance pathways among the residual mPCs. Further studies are warranted to elucidate the functional implications of these interactions in the immune escape and survival of selective mPCs.

This scRNA-seq study of paired diagnostic and posttreatment samples represented the most extensive single-cell dataset to date. It elucidated transcriptional dynamics and uncovered resistance pathways within MRD tumor cells in patients with MM. There were five main findings in this study. First, a population of nlPCs was identified in MM. Second, three main trajectories of transcriptional evolution posttreatment were observed. These were clonal elimination in MRD− patients, as well as clonal stabilization and clonal selection in MRD+ patients. Third, a close link between genetic and nongenetic mechanisms behind the clonal evolution of MM was confirmed, which further supported the multidimensional background of clonal evolution. Fourthly, transcriptional and metabolic adaptations were observed as important activated prosurvival pathways in MRD tumor cells. Finally, differential cellular interactions between mPCs and TME cells were observed, showing that the TME might provide shelter for MRD tumor cells.

Genetic studies employing iFISH and whole-exome sequencing have demonstrated that MM is not a single disease entity but also a spectrum of molecular conditions, each with distinct clinical outcomes (58). Moreover, gene expression profiling has revealed various subgroups within these conditions, defined by the cytogenetic features of PCs and the aberrant expression of a D-group cyclin (41). Previous phenotypic studies have identified several subclones within patients with MM, each associated with unique chemoresistant profiles, cytogenetic characteristics, and clonogenic potentials (7). Recently, an MGUS-like phenotype has been described in MM patients based on the prevalence of BM PC and clonal PC within the BM PC compartment (59). These findings underscore the remarkable intertumor heterogeneity of MM and suggest the existence of a cell population evolving from normal to mPCs identifiable at the early stages of MM pathogenesis. In alignment with previous scRNA-seq studies (29, 32, 46), our study identified a PC cluster composed of PCs from HDs and multiple MM patients. Interestingly, even clonal PCs from MM patients co-clustered with HD PCs. We designated such clonal PCs as nlPCs. Further analyses revealed several distinguishing features of nlPCs, such as the absence of secondary genetic events as shown in their CNA profiles, the expression of monoclonal immunoglobulins or light chain restriction, and extensive transcriptional alterations.

Therapy-induced clonal evolution plays a crucial role in MM drug resistance (58). As treatments exert selective pressure, subclones either carry or acquire specific genetic alterations, conferring a competitive fitness that enables these PCs to evade therapeutic pressure. Utilizing quantitative multigene fluorescence in situ hybridization, we have charted the phylogeny of tumor subclones and identified four patterns of evolutionary architecture in 129 longitudinal samples from 57 MM patients. We found that patients exhibiting clonal stabilization had superior OS, whereas those with linear evolution had the worst outcomes (16). Increasingly, studies are focusing on clonal evolution in MRD tumor cells (6, 7). Previous research indicates that standard-risk patients are more susceptible to clonal selection events, whereas high-risk patients more often acquire mutations (4). A recent single-cell study has revealed that genetic abnormalities observed at relapse are already present at diagnosis within subclonal populations (13). Based on these findings, we investigated the transcriptional dynamics between disease diagnosis and posttreatment. Our results revealed three main trajectories of transcriptional evolution posttreatment: clonal elimination in MRD− patients, as well as clonal stabilization and clonal selection in MRD+ patients. Furthermore, by delineating the clonal architecture of genetic events per patient using inferCNV (55), we found a high correlation between genetic and nongenetic mechanisms underpinning the clonal evolution of MM. We proposed further investigations to elucidate the clinical significance of nlPCs.

There is growing evidence suggesting that drug resistance in cancer cannot be solely attributed to simple genetic causes (19). The concept that nongenetic mechanisms contribute to therapeutic resistance is increasingly gaining recognition (19). In the context of MM, subclones with identical genetic backgrounds can result in multiple transcriptional clusters (29), which points to the nongenetic mechanisms that influence transcriptional status. A recent scRNA-seq study has reported upregulation of the UPR, endoplasmic reticulum, and mitochondrial stress pathways in primary refractory MM patients (32). Aligning with these findings, we observed a metabolic shift to fatty acid oxidation in resistant PCs in a state of cell-cycle arrest, whereas selective PCs demonstrated a preference for the NF-κB pathway in our study. Through this mechanism of transcriptional selection, MM cells could rapidly adapt to therapeutic pressure and attain a selective advantage within the therapeutic landscape.

MM progression involves aberrant interactions between PCs and the TME (60). Comprehensive evidence suggests that both hematopoietic and nonhematopoietic cells, including regulatory T cells (61), regulatory B cells (62), myeloid-derived suppressor cells (63, 64), and stromal cells (56), contribute to PC proliferation. Our study identified a potential supportive role of the TME for MRD tumor cells. The most pronounced interaction between MRD tumor cells and the TME was observed in selective PCs. Furthermore, the upregulated expression of CD74 by selective PCs and the potential contribution of CD74-related PC-myeloid cell interactions to their immune evasion were noted. Given that myeloid cells or monocyte-derived cells, including osteoclasts (65) and DCs (66), comprise the primary cellular component of the microenvironment that promotes immunosuppression, disrupting the interactions between selective PCs and myeloid cells by blocking CD74 could be a potential therapeutic approach to eradicate MRD tumor cells.

This study also highlights several remaining challenges. First, because only one MGUS sample was sequenced in our study, the interpretation of the results of our study needs to take into account the limited number of MGUS patients in our cohort. Nevertheless, the identification of nlPCs added substantial weight to the concept of intertumor heterogeneity in MM. Moreover, the inclusion of steroids in the induction regimen may lead to transcriptional and metabolic adaptation in residual mPCs. Further experimental studies are necessary to confirm whether the transcriptional and metabolic adaptation in mPCs is in response to a specific drug or represents a general phenomenon.

In conclusion, single-cell profiling of MRD tumor cells could provide novel insights into the transcriptional evolution and molecular mechanisms underlying therapeutic resistance. This knowledge might help identify new vulnerabilities that can be targeted, thereby contributing to improved long-term disease control in MM.

L. Qiu reports grants from Pfizer, nonfinancial support and other support from BeiGene, and grants, nonfinancial support, and other support from Johnson & Johnson outside the submitted work. No disclosures were reported by the other authors.

J. Cui: Conceptualization, data curation, formal analysis, validation, investigation, visualization, writing–original draft, writing–review and editing. X. Li: Conceptualization, data curation, software, formal analysis, investigation, visualization, writing–original draft, writing–review and editing. S. Deng: Conceptualization, funding acquisition, investigation, writing–original draft, project administration, writing–review and editing. C. Du: Resources, data curation. H. Fan: Resources, data curation. W. Yan: Resources, data curation. J. Xu: Resources, data curation. X. Li: Writing–review and editing. T. Yu: Resources, data curation. S. Zhang: Resources, data curation. R. Lv: Resources, data curation. W. Sui: Resources, data curation. M. Hao: Resources, data curation. X. Du: Writing–review and editing. Y. Xu: Data curation, supervision. S. Yi: Data curation, supervision. D. Zou: Data curation, supervision. T. Cheng: Data curation, supervision. L. Qiu: Conceptualization, resources, data curation, supervision, funding acquisition, writing–original draft, writing–review and editing. X. Gao: Conceptualization, resources, data curation, supervision, funding acquisition, writing–original draft, writing–review and editing. G. An: Conceptualization, resources, data curation, supervision, funding acquisition, writing–original draft, project administration, writing–review and editing.

We thank Prof. Nikhil C. Munshi for his insightful suggestions for this study. We thank all the patients and families for their participation in this research study. This study was supported by the National Natural Science Foundation (grants 82270218, U22A20291 to G. An), Youth Fund of the National Natural Science Foundation of China (grant 81900214, to S. Deng), the National Natural Science Foundation (grant 32370722 to X. Gao), the CAMS Innovation Fund for Medical Sciences (CIFMS; grant 2022-I2M-1-022 to L. Qiu, grant 2021-I2M-C&T-B-079 to G. An, and grant 2021-I2M-1-041), and the International Cooperation Projects of National Natural Science Foundation (grant 81920108006, to L. Qiu). J. Cui is supported by a grant from International Myeloma Society.