Malignant clonal evolution drives multiple myeloma cellular ecological diversity and microenvironment reprogramming
Molecular Cancer volume 21, Article number: 182 (2022)
Multiple myeloma (MM) is a heterogeneous disease with different patterns of clonal evolution and a complex tumor microenvironment, representing a challenge for clinicians and pathologists to understand and dissect the contribution and impact of polyclonality on tumor progression.
In this study, we established a global cell ecological landscape of the bone marrow (BM) from MM patients, combining single-cell RNA sequencing and single-molecule long-read genome sequencing data.
The malignant mutation event was localized to the tumor cell clusters with shared mutation of ANK1 and IFITM2 in all malignant subpopulations of all MM patients. Therefore, these two variants occur in the early stage of malignant clonal origin to mediate the malignant transformation of proplasmacytes or plasmacytes to MM cells. Tumor cell stemness index score and pseudo-sequential clonal evolution analysis can be used to divide the evolution model of MM into two clonal origins: types I and IX. Notably, clonal evolution and the tumor microenvironment showed an interactive relationship, in which the evolution process is not only selected by but also reacts to the microenvironment; thus, vesicle secretion enriches immune cells with malignant-labeled mRNA for depletion. Interestingly, microenvironmental modification exhibited significant heterogeneity among patients.
This characterization of the malignant clonal evolution pattern of MM at the single-cell level provides a theoretical basis and scientific evidence for a personalized precision therapy strategy and further development of a potential new adjuvant strategy combining epigenetic agent and immune checkpoint blockade.
Multiple myeloma (MM) is the most common type of plasma cell malignancy (93%), characterized by the abnormal proliferation of terminally differentiated clonal plasma cells in the bone marrow (BM), and accompanied by chromosomal instability and cytogenetic abnormalities [1,2,3]. The progress of MM is a multi-level and multi-stage dynamic process involving a wide range of clonal genetic variation and molecular evolutionary dynamics, which drive the heterogeneity of the tumor cell ecological composition and spatial clonal structure, mediating treatment resistance and recurrence [4,5,6]. Efforts to link the progress and prognosis of MM with gene resistance cloning have motivated further study into identifying the drivers of genetic variation in functionally heterogeneous cloning . Nonetheless, existing studies on the clonal evolution of MM are mainly based on bulk-sequencing or computer simulation data at the tissue level of a high-volume BM aspirate [8,9,10,11]. To date, based on single-cell sequencing, Amit et al. identified resistance pathways and therapeutic targets in relapsed MM , and Lohr et al. revealed metabolic reprogramming as a resistance mechanism in BRAF-mutated MM . However, few studies reported the clonal evolution trajectory of MM malignant cells at the single-cell level. Moreover, such analyses ignore the intratumor cell heterogeneity that is not only an important result of tumor subclonal evolution at the cellular level but also the micro-scale basis of individual heterogeneity, mediating the diverse prognostic outcomes of patients [14,15,16].
Existing studies have only focused on myeloma evolution with respect to primary variation and acquired secondary variation in longitudinal samples, without exploring the interaction mechanism of variation and selection [17, 18]. Variation and selection play equally important complementary and promotive roles in establishing tumor evolutionary patterns. Additionally, there remain many challenges in transforming the observational studies of rich heterogeneous cell-state phenomenology from single-cell analysis to a mechanical understanding of cancer evolution dynamics and causal mechanism patterns . For example, defining a common feature to reflect the evolutionary process of a tumor would help to obtain more recognizable tumor biological characteristics toward developing better intervention strategies.
Therefore, in this study, we combined single-molecule long-read genome sequencing with single-cell RNA sequencing (scRNA-seq) to identify genomic instability events and the single-cell ecological landscape in MM patients. Bioinformatics approaches were applied for the direct reconstruction and in-depth global description of the cloning pattern of MM malignant cells in natural and drug-driven states to determine the change and selection of the microenvironment during the evolutionary process. The main findings were then verified in large-scale clinical MM samples, offering significant robustness in the potential to obtain a general framework of tumor cell evolution, replacing the current coarse-grained, step-wise, and deterministic cell and tissue models .
Materials and methods
Among the seven patients, three had newly diagnosed MM (NDMM), three had relapsed and/or refractory MM (RRMM), and one had CD20-positive RRMM (CD20+ RRMM). Among the three RRMM patients, RRMM1 was resistant to thalidomide and melphalan; RRMM2 received bortezomib, pegylated liposomal doxorubicin, cyclophosphamide, and dexamethasone but failed to respond; and RRMM3 received lenalidomide, bortezomib, pegylated liposomal doxorubicin, cyclophosphamide, and dexamethasone but relapsed after all regimens. Additionally, a patient was previously diagnosed with primary central nervous system lymphoma without BM invasion and had been in complete remission for 1 year without any treatment at the time of BM sampling and was considered as a control.
BM aspirates were collected from seven MM patients and one control donor who agreed to a multiplex library and sequencing protocol that covered all study procedures, which were approved by the Ethics Review Committee of Beijing Tongren Hospital and Sun Yat-sen University Cancer Center. All patients provided written informed consent. The clinical data of all patients are shown in Supplementary Table 1. All sequencing was performed at Biomarker Technologies Corporation (Beijing, China).
Single-cell transcriptome profiling
Sample preparation and cDNA library construction were performed with the 10× Genomics Single Cell 3’v3.1 kit according to the manufacturer’s instructions. Based on microfluidic technology, individual cells and reagents required for the reaction were wrapped in GEM droplets with a bead (containing a cell barcode) on the microchip, and the droplets containing the cells were collected. The cells were lysed to release mRNA, which binds to the cell barcode primer on the bead to complete the reverse transcription reaction. The GEMs were broken, and cDNA was recovered and amplified by polymerase chain reaction to construct the cDNA library. The cDNA product and library concentration were detected using a Qubit 4.0 fluorescence quantification instrument, and the insert size was detected using a Qseq400 Bioanalyzer to ensure a single peak type, no spurious peak, no junction, and no primer dimer. Finally, the sample library was sequenced using the NovaSeq 6000 instrument on the Illumina platform. After identifying the cassava base, the original image file was converted into a sequence file and stored in FASTQ format.
Single-molecule long-read sequencing and Nanopore sequencing
The sequencing procedure was performed according to standard Oxford Nanopore Technologies protocols . First, high-molecular-weight genomic DNA was extracted from the BM aspirate using HiPure Tissue & Blood DNA Kit (D3018–03, Magen), and Nanodrop, Qubit, and 0.35% agarose gel electrophoresis were used for purity, concentration, and integrity quality inspection. Next, 2 μg of high-quality nucleic acid was prepared and fragmented with a G-TUBE tube, thereby breaking the genomic DNA to an average of approximately 8-kb fragments. The genomic DNA ligation reaction kit (SQK-LSK109) was used to construct the library. NEBNext FFPE DNA Repair Mix and NEBNext Ultra II End Repair/dA-Tailing Module were used to repair the damage of nucleic acid fragments, end repair, and to add A bases to the end of the fragments, thereby purifying the magnetic beads. After adding the barcode sequences and ligating the sequencing junction using Amplification Free Barcode Expansion Kit 1–12 (EXP-NBD104, Nanopore), the beads were purified to complete the library construction. Finally, the library was quantified and sequenced based on the Qubit fluorometer. For the original electrical signal obtained by sequencing, Guppy software  was used for neural network base calling (https://github.com/rrwick/Basecalling-comparison) to obtain the original sequence file, which was stored in FASTQ format.
Cell-surface labeling was performed using flow cytometry to identify and verify MM malignant plasma cells. To detect the expression of cell antigens, cell size, and intracellular particle content, a 100 μL sample of fresh BM with EDTA was collected and mixed with fluorescein-labeled monoclonal antibodies (BD and Beckman Coulter Company) to the cell-surface markers including CD45, CD38, CD56, CD138, CD20, and CD19. Then, the mixture was incubated for 15 min at room temperature and protected from light. Red blood cell lysate (2 mL) was added, and incubated for 10 min at room temperature in the dark, followed by the addition of 2 mL phosphate-buffered saline (PBS), and the sample was centrifuged at 352 g for 5 min. After removing the supernatant, 3 mL PBS was added to wash and discard the supernatant, and then 300 μL PBS was added for machine detection. Finally, the results were analyzed on the FACSCanto II flow cytometer (BD) with BD FACSDiva software.
Sequence alignment and gene quantification
Cell Ranger 5.0.1 software was used for sequence comparison and quantification of the sequencing data (per official recommendation of 10× Genomics), and the sequence read was mapped and aligned to the reads of the reference human genome reference (Hg38) using STAR software . All unique gene names of the transcripts were recorded, the cells were labeled by the barcode, and the transcripts were labeled by unique molecular identifiers (UMIs) to quantify the number of cells and genes after comparison with the reference. All reads that mapped to the same gene and had the same UMI sequence were folded and different UMIs corresponding to the same gene were quantified, which produced a digital matrix for cell gene expression quantification. For all downstream analyses, we selected cells that have at least 1000 UMIs (indicating the number of captured transcripts) mapped to at least 200 unique genes and ensured that each gene is expressed in more than three cells. We excluded cells with poor viability and quality by removing more than 10% of the cells whose gene counts reflected mitochondrial genes or ribosomal RNA.
Construction of a single-cell atlas
The IntegrateData function in the R Seurat package  was used to merge single-cell data, and cell clustering analysis was performed according to default parameters (http://satijalab.org/seurat/). Principal component analysis and t-distributed stochastic neighbor embedding (t-SNE) methods were used for dimensionality reduction and visualization of the clustering results, and the results were projected to a two-dimensional image, which was defined as a single-cell atlas. The “FindAllMarkers” function in the Seurat package was performed to identify the specific genes expressed in each cell cluster, with P < 0.05 considered statistically significant.
Identification of cell types
For control donor BM cells, we used the SingleR package in R  to annotate the cell types. SingleR assigns cell identities to single-cell transcriptomes by comparison with the reference datasets of pure cell types for microarray or RNA-seq sequencing. Here, we used the previously defined single-cell transcriptome expression profile as a reference system . We identified and extracted malignant plasma cells based on clinical and laboratory features of immunophenotypes in MM (CD38+CD56+CD138+CD19−CD20−). For non-classical CD20+ RRMM patients, a phenotype of CD38+CD56+CD138+CD19−CD20+ was used as the screening condition for malignant plasma cells. The extracted malignant plasma cells were verified using flow cytometry. Clusters were annotated based on the expression of known marker genes (Supplementary Table 2).
For label scoring of the cell cycle and cell characteristics, we adopted non-parametric and unsupervised scoring assumptions based on the expression patterns of characteristic genes, similar to a previously described gene set scoring strategy . For scoring of the cell cycle and cell characteristics, an unlabeled sample gene expression matrix was used as input, which included the scRNA-seq count expression profile and log2-standardized chip expression profile of large-scale clinical MM patients or a bulk-seq count expression profile. First, the algorithm performed a non-parametric nuclear density estimation test on the overall gene expression profile. Second, based on the kernel density estimation results, the samples were sorted according to their expression levels. Next, the samples were ranked for expression levels based on the results of the nuclear density estimation, and the cell cycle and the rank statistic for each cell characteristic were calculated, similar to the Kolmogorov–Smirnov test. Finally, the cell cycle and enrichment scores of each cell feature were obtained and the output was provided as a data matrix corresponding to each sample.
The stemness score of the single-cell atlas was obtained using the TCGAanalyze_Stemness function in the R package TCGAbiolinks , which evaluates the degree of carcinogenic differentiation by extracting a series of marker genes to quantify the characteristics of stem cells from the transcriptional expression and epigenetic patterns of non-transformed pluripotent stem cells and their differentiated offspring using a publicly available molecular atlas . Subsequently, one-class logistic regression machine learning algorithm (OCLR) was used for multi-platform analysis of these transcriptomic, methylomic, and transcription factor (TF) binding sites to obtain two independent indices of stem cell characteristics: DNA methylation-based stemness index (mDNAsi), which reflects epigenetic characteristics, and gene expression-based stemness index (mRNAsi), which reflects gene expression patterns. We applied a stemness score on gene expression patterns (mRNAsi) to the single-cell atlas of MM malignant subclones and control donor BMs to identify the evolutionary patterns of malignant origin and intratumor molecular heterogeneity in MM.
The stemness score can determine the origin of the clonal evolution of MM malignant plasma cells, and pseudotime analysis can infer the trajectory of its evolution and development. We used the R package Monocle 3  to reconstruct the developmental trajectory of the control donor BM single-cell atlas and simulate the evolutionary trajectory of the malignant subclone of MM. Monocle 2 was applied to simulate the developmental trajectory of the immune cells of the BM microenvironment subjected to evolutionary reprogramming of the MM malignant clone. Subsequently, Moran’s I statistic was used to identify genes expressed in complex trajectories under the malignant clonal evolution of MM. These genes may be the molecular driving force for the natural development of MM from a malignant origin and the evolution of drug resistance under drug selection. Additionally, we calculated the RNA rate (time derivative of the gene expression state) to predict the future state and final fate of a single cell, and to analyze its developmental lineage and cell dynamics using velocyto. R  by distinguishing between unspliced and spliced mRNAs in the single-cell atlas.
Gene regulatory network
Single cell regulatory network inference and clustering (SCENIC)  was used to infer gene regulatory networks based on single-cell expression profiles and identifying cell states, providing important biological insight of the mechanism driving cell heterogeneity. To identify the internal transcriptional regulation driving force of control donor BM cell development and MM malignant clone evolution, we used the python module tool pySCENIC to analyze and reconstruct the gene regulatory network with TFs as the core.
The workflow starts by describing the input single-cell expression abundance spectrum matrix and applying a regression per-target approach (GRNBoost2) to infer the co-expression module, from which the indirect targets were pruned based on cis-regulatory motif discovery (cisTarget). Subsequently, AUcell was used to quantify the activity of these regulators by enriching and scoring the target genes of the regulators to obtain the regulon activity score (RAS). The single-cell data were further downscaled using the RAS matrix, and the regulon-specific score (RSS) was calculated based on Jensen–Shannon divergence (JS) to determine the cell cluster-specific regulon. The most specific and significant regulon was mapped to the single-cell cluster atlas and verified using massively parallel sample sequencing. Finally, the connection specificity index (CSI) matrix was calculated and the regulon was hierarchically clustered based on the CSI to define the regulon module to obtain the relationship between the regulon module and regulon and visualized based on the R package ComplexHeatmap .
Three intercellular communication event identification tools were used to identify high-confidence ligand-receptor interactions between cells: CellPhoneDB , iTALK , and our newly developed tool CellCrosstalk. Based on the joint expression of multi-subunit ligand-receptor complexes to infer intercellular communication, CellPhoneDB emphasizes joint positive expression but neglects the importance of interaction in single-arm expression and marker genes. The R package iTALK compensates for this defect by prioritizing the identification of highly expressed or differentially expressed genes (DEGs) in cell clusters, which are matched through the ligand-receptor database to identify important intercellular communication events. However, when genes with high or differential expression are mapped to the ligand-receptor network, the lack of corresponding significance identification may result in a large number of false-positive events. Therefore, we developed an innovative tool for recognizing intercellular communication events named CellCrosstalk (original code: https://github.com/ydlife/CellCrosstalk), which prioritizes the identification of highly expressed genes or DEGs to ensure that independent ligands or receptors are positive in the corresponding cell clusters. Notably, the built-in ligand-receptor interaction network database of CellCrosstalk summarizes the built-in data of CellPhoneDB and iTALK. These genes were mapped to real and 1000 random walk ligand-receptor networks; the ligand-receptor interaction pair mapped between two cell clusters was recorded each time, and the relevant cells were identified based on hypergeometric tested intercellular communication events. Finally, combining the results of the three tools, the ligand-receptor interaction pair identified by any two tools was considered a high-confidence intercellular communication event.
Biological process and pathway enrichment analysis
The R package clusterProfiler  was used to perform functional enrichment analysis on Gene Ontology biological processes and Kyoto Encyclopedia of Genes and Genomes pathways for related genes (P < 0.05).
Genome chromosome structure variation and copy number variation (CNV)
Based on minimap2 software , clean sequencing data of BM aspirates from the control donor and MM patients were compared to the reference genome (Hg19) to obtain the comparison SAM file. Using open-source samtools software (https://sourceforge.net/projects/samtools/), the files were converted into BAM format and sorted. Sniffles software  was used to identify genomic chromosome structural variations of multiple samples, which were merged using SURVIVOR software. The CNV of whole-genome sequencing data in the BAM file was detected using the R package QDNAseq (https://github.com/ccagc/QDNAseq).
Single-nucleotide variant (SNV) analysis
Four different GATK framework processes were used to identify SNVs (https://gatk.broadinstitute.org/hc/en-us). First, the germline short variant discovery (SNPs + Indels) process of the GATK framework was performed on clean Nanopore sequencing data from control donor BM aspirates (https://gatk.broadinstitute.org/hc/en-us/articles/360035535932-Germline-short-variant-discovery-SNPs-Indels-) to obtain a list of SNVs used as controls. Second, the sequencing data of BM aspirates from MM patients were subject to the somatic short variant discovery (SNVs + Indels) process of the GATK framework (https://gatk.broadinstitute.org/hc/en-us/articles/360035894731-Somatic-short-variant-discovery-SNVs-Indels-) to screen MM disease-associated somatic mutation information using control SNVs as a reference. Third, single-cell transcriptional profiles of control donors and MM patient BM tumor microenvironment (TME) cells were subject to the short variant discovery (SNPs + Indels) process (https://gatk.broadinstitute.org/hc/en-us/articles/360035531192-RNAseq-short-variant-discovery-SNPs-Indels-) to obtain a list of SNVs of control single cells and information on germline mutations in MM patients. Fourth, single-cell transcriptional profiles of malignant plasma cells from MM patients were applied to the RNA-seq short variant discovery (SNPs + Indels) process of on somatic mutation calling. The SNVs detected in MM patients and single-cell transcriptional profiles of malignant plasma cells, but not detected in control donors were considered to be MM malignant subclone-associated SNVs.
First, STAR was applied to map the scRNA-seq FASTQ format clean data of malignant plasma cells to the reference. Second, MergeBamAlignment MarkDuplicates was performed for data cleanup, and the data were processed using SplitNCigarReads. Subsequently, base quality recalibration was completed using BaseRecalibrator, Apply Recalibration, and AnalyzeCovariates. Mutect2 was applied to all candidate variants instead of the RNA-seq short variant discovery (SNPs + Indels) process for BM TME cells of control donors and MM patients, where HaplotypeCaller was applied for all variants. Next, GetPileupSummaries and CalculateContamination were applied to calculate contamination and Learn Orientation Bias Artifacts was completed based on the LearnReadOrientation Model. Finally, FilterMutectCalls was used to filter variants, and the implementation of annotated variants was based on Funcotator.
Some SNVs only existed in some malignant subclones, indicating their involvement in the malignant clonal evolution of MM. Although similar strategies for SNV detection have been reported previously [38,39,40], several false-positive and false-negative results were obtained. Therefore, we adopted the latest process recommended by the GATK framework to obtain reliable results. Based on these high-confidence SNVs, we calculated the tumor mutational burden (TMB) of MM patients at the BM aspirate and malignant subclonal levels based on the number of somatic variants per megabase in the genome . The high concordance of the TMB at both levels suggests high reliability and robustness of our SNV-calling strategy.
Estimation of CNV in single cells
Single-cell CNV was estimated using the inferCNV package (inferCNV of the Trinity CTAT Project, http://github.com/broadinstitute/inferCNV), which compares the gene expression of each tumor cell with the average expression or “normal” reference cell gene expression to determine its expression intensity, and displays the relative gene expression on each chromosome in the form of a heatmap. Compared with normal cells, MM malignant subclones always have over- or under-expression of local genome fragments. The resulting inferred malignant subclonal CNV events were highly consistent with the true CNV events detected by Nanopore sequencing of BM aspirates from MM patients, which not only validated the feasibility and reliability of the inferCNV strategy but also localized the CNV in patients to specific malignant subclones.
Analysis of clonal evolution based on malignant cell clusters
A total of 11 MM malignant clonal clusters were identified in this study. Subsequently, based on the proposed temporal analysis, we preliminarily established that these malignant clonal clusters had two potential malignant origins. The R package fishplot was applied to trace the clonal structural evolution of these malignant plasma cell clusters (https://github.com/chrisamiller/fishplot).
To predict the prognostic potential of clonal profiles, we extracted specific marker genes with log fold change (logFC) > 1, which were used as cell characteristics to obtain non-parametric and unsupervised scores in a large-scale clinical cohort of 9574 patients with 24 independent datasets to determine the relative abundance of malignant progenitor cells. The difference in relative abundance was used as a standard score for MM malignant origin dominance to classify patients into origin types. Patients were also staged according to the relative abundance of malignant origin progenitor cells: type I and IX double positive, double negative, type I-specific positive, and type IX-specific positive. The R package survminer (https://github.com/kassambara/survminer) was used to determine the prognostic significance of the relative abundance of progenitor cells of malignant origin, the advantage score of malignant origin, and the classification of malignant origin.
Prognostic model for MM drug resistance-related markers
Drug resistance evolution-related markers were extracted in the evolutionary pattern of MM malignant clones that were significantly highly expressed in the large-scale datasets of patients. These drug resistance markers were used for univariate Cox regression analysis of overall survival (OS) in the training cohort GSE136400, and the characteristic genes significantly correlated to prognosis were obtained (P < 0.01). These prognostic characteristic genes were included in multivariate Cox regression analysis for OS and relapse-free survival (RFS) to establish a prognostic model, which was tested using the R package survminer, and externally validated in independent GSE9782 and MMRF-CoMMpass cohorts (https://themmrf.org/). Time-independent receiver operator characteristic (ROC) curves were then conducted to assess the prediction performance of the prognostic model.
Self-organizing mapping (SOM) analysis
The contribution of malignant subclones to the molecular heterogeneity of BM aspirates from MM was quantified by SOM analysis . First, patient-specific marker genes in malignant subclones were identified using the R package Seurat, and their expression matrices were extracted. Each malignant subclone was divided for differential expression analysis, respectively.
A matrix of Pi values was used as the standard input to the SOM, where the row represents genes and the columns represent malignant subclones, where Pi = −log10 (adjusted P (P.adj)) × logFC. Next, the R package Kohonen (https://github.com/iamciera/SOMexample) was applied for SOM analysis, resulting in marker gene–MM patient-specific association patterns and malignant subclonal contributions of patient heterogeneous molecules. The Kohonen topology-preserving map creates a multidimensional continuous representation of “perceptual space” on a neuronal grid, where input by vector x = (× 1,..., x d). Exemplars of these vectors were repeatedly presented to the two-dimensional neuronal network organization to simulate the various stimuli experienced by the senses, which in turn condensed the labeled genes into the corresponding neural units. Thus, genes with similar patient heterogeneity contributions were organized in the same SOM grid neural unit and similar neural units were clustered in close proximity.
According to the strategy of Ma et al. , the degree of intratumor heterogeneity was measured by calculating diversity scores with gene expression profiles of malignant cells within the tumor. First, the intratumoral heterogeneity of each tumor was measured based on principal components (PCs) rather than the original gene expression profile to obtain the primary information and reduce noise. The top 30 PCs were selected for subsequent calculations based on the alignment of eigenvalues and were significant. This criterion was determined after using different numbers of PCs to determine the robustness of the diversity score calculation. Subsequently, the diversity score for each patient was calculated.
Expression profile data analyzed in this study were obtained from the Gene Expression Omnibus database (https://www.ncbi.nlm.nih.gov/geo/) and EMBL-EBI database (http://www.ebi.ac.uk) with accession numbers GSE136400, GSE9782, and E-TABM-1138. RNA-seq parallel sequencing data were obtained from the Multiple Myeloma Research Foundation CoMMpass Study (https://research.themmrf.org/).
The raw data of scRNA-seq and single-molecule long-read genome sequencing have been deposited in the Genome Sequence Archive in National Genomics Data Center (GSA: HRA001335) and are publicly accessible at https://ngdc.cncb.ac.cn/gsa.
Single-cell landscape of control BM
From scRNA-seq, 64,718 cells were obtained, including 5238 from the control donor BM aspirate and 59,480 from MM patients after quality control (Fig. 1A, Supplementary Table 3). The t-SNE approach captured 16 cell types, including hematopoietic lineages, myeloid lineage cells, and lymphoid cells, all of which belong to common cell groups in the BM (Fig. 1B), consistent with the established phenotypic characteristics of immune cells (Fig. 1C). The pseudotime trajectory of cell development revealed the continuous process of myeloid development and hematopoiesis, that hematopoietic stem cells (HSCs) differentiate into promonocytes to monocytes to finally dendritic cells (DCs), and that HSCs differentiate into erythroblasts to post erythroblasts (Fig. 1D). The pseudo-sequential differentiation trajectory was negatively correlated with the stemness and cycle scores of BM cells, which is in line with the physiological changes of cell differentiation (Fig. 1E).
The single-cell profile of the control donor was constructed to reveal the ecological composition of different hematopoietic cell types in the normal BM, suggesting differentiation trajectory and fate choice, which is consistent with the current hematopoietic concept. Thus, this concept can be applied as a training cohort and reference system for the cellular ecological landscape of MM.
Cellular ecosystem landscape in MM
The malignant plasma cells extracted from NDMM and RRMM patients (CD38+CD56+CD138+CD19−CD20−; Supplementary Fig. 1A, B and Fig. 2A-C) were further clustered, with 11 malignant plasma cell clusters obtained, whereas six malignant plasma cell subsets were obtained for CD20+ RRMM patients (CD38+CD56+CD138+CD19−CD20+; Supplementary Fig. 1C and Fig. 2D, E). Malignant cell clusters shared cancer characteristics and widely significantly overexpressed HLA-A, HLA-B, MCL1, HDAC1, LCK, HSPB1, and IL6R (Fig. 2F), indicating a common tumor origin. Corresponding specific markers were also identified between different malignant cell clusters (Fig. 2G, Supplementary Table 4), providing direct evidence for the formation of different subpopulations.
We identified and annotated the major cell types within the microenvironment in MM patients by investigating the expression patterns of known marker genes, which enabled to accurately define the specific identity of the microenvironment cells in MM patients (Fig. 2H, I, Supplementary Table 5). There was complex and active communication between different BM cells in MM patients, regardless of newly diagnosed and relapsed patients (Supplementary Fig. 1D, E and Fig. 2J). Consequently, the global cell ecological landscape of BM from MM patients can thus be characterized and clarified.
Variation events related to cell clonal evolution
Chromosomal instability is a hallmark of human cancer and tumor heterogeneity. The structural variations catalog (Supplementary Fig. 2A) in the present study showed significant genomic instability in the majority of chromosomes, except for chr4, including deletion (DEL), duplication, inversion, and insertion (INS). Genes related to DEL and INS events were expressed and played a major role in driving malignant cell clustering (Fig. 3A). Furthermore, malignant CNV events (Fig. 3B) and the corresponding transcripts (Fig. 3C), such as TNFSF13B, CD79A, TNFRSF13B, PARP1, IMPDH2, and MYC, were identified, expressed either alone or in combination in different malignant cell subpopulations and contributing to polyclonality. IFITM2 has been proven as an effector gene of the type I interferon response that protects cells against invading viral pathogens ; the shared mutation of IFITM2 (chr11:309127:A > G) in all malignant subpopulations was detected at both the bulk and the single-cell levels, with high expression activity in most malignant cells (Fig. 3D), suggesting its potential as a candidate mutation of malignant origin. Similarly, the shared mutation of ANK1 (chr8:41510767:T > G) was detected in all malignant subpopulations, but only expressed in type I malignant subpopulation. These mutations were concentrated in the patterns of C > T and T > C (Supplementary Fig. 2B), which is consistent with the basic genetic concept that C in CpG dinucleotides tends to mutate to T after methylation [45, 46]. Further analysis showed that other patients had higher TMB levels than that of the CD20+ RRMM patient and had enriched type IX malignant cell clusters (Supplementary Fig. 2C). TMB was consistent at the bulk and the single-cell levels (Supplementary Fig. 2D). In addition to genomic variation, the gene regulatory network with TFs as pivots was organized into six modules (Fig. 3E), such as BCL6, FOXO1, E2F7, and FOXP2, to regulate the specific gene expression (Fig. 3F) and RNA transcription rate (Fig. 3G) of MM malignant cell subpopulations to guide cellular fate choice. This promotes the transformation and differentiation of the core state (Supplementary Fig. 3A), ultimately mediating the formation of a series of clonal phenotypes .
The pseudo-sequential clonal evolution atlas of MM malignant cells was finally constructed (Fig. 3H), which is consistent with the trend of the cell stemness index and cell cycle score (Supplementary Fig. 3B, C), whose corresponding gene expression, signaling pathways, and biological functions are shown in Fig. 3I. According to the cell stemness index score and pseudo-sequential clonal evolution analysis, the malignant origins of the cells were divided into type I and type IX malignant progenitor cell origins with the highest level of cell stemness index score (Fig. 3J). During the natural development of MM cells (Fig. 3K), type I and type IX malignant progenitor cells evolved into type II, III, and IV and type VIII, X, VII, and V malignant cells, respectively, with different drug sensitivity profiles; thus, type IV, VIII, VII, and V subpopulations explosively grew to an occupied advantage. However, with the selection of drugs and the occurrence of chemical carcinogenic variants, new cloned type VI malignant plasma cells were formed and became a superior subpopulation.
Multi-omics abnormal program identification of primitive MM malignant progenitor cells
Subsequently, we focused on primitive MM cell subpopulations, including type I and IX malignant progenitor cells, which promote the occurrence and growth of tumors with high stemness and cell cycle activity (Fig. 4A). Type IX malignant progenitor cells showed higher stem cell characteristics, resulting in a higher degree of malignancy in patients. Their cell markers with logFC > 1 were extracted to verify the expression and relative abundance in a large-scale clinical cohort of 9574 patients with 24 independent datasets (Supplementary Table 6) and the MM pathological classification of malignant origin was performed (Fig. 4B, Supplementary Table 7). A higher relative abundance of malignant progenitor cells was associated with a worse prognosis (P < 0.0001; Fig. 4C) for both types I and IX, especially for the latter (P < 0.0001). We observed similar results at the single-cell level with high robustness. Patients with a double-positive origin of type I and IX had the worst prognosis, followed by those with a double-negative origin (P < 0.0001). Consistent with the relative abundance analysis, the prognosis of patients positive only for type IX origin was worse than that for only type I origin.
We evaluated the malignant forerunner events that were preferentially highly expressed by two malignant origin progenitor cells, including DNA damage repair, cell cycle, proliferation, migration, invasion, and stemness (Supplementary Fig. 3D). In this regard, they showed extensive similarities, with many previously reported hematological tumor-related genes, such as HMGB1, CCND2, CDK1, CDKN2A, and MYC [47,48,49,50,51]. These genes were more significantly dysregulated in type IX malignant cells, which explains the higher stem cell characteristics and malignant degree in the type IX origin. The similarity of malignant origin plasma cells also lies in their significantly dysregulated signaling pathways (Supplementary Fig. 3E). It is noteworthy that the type I malignant origin was significantly enriched in both base excision repair and DNA damage repair to reduce random variation events in the process of malignant proliferation, while type IX only activated the base excision repair, causing more random variation events in the process of clonal proliferation with a richer clonal pattern. Related dysregulated genes were clustered into co-expression modules according to expression similarity, which are mainly regulated by CNV events, involving SNVs and TFs (Fig. 4D), resulting in a shared carcinogenic initiation program and its own specific carcinogenic regulation mechanisms in two malignant origins of MM (Fig. 4E). Homogeneity and specificity are reflected not only in the carcinogenic processes of these origin types but also in the change in the expression pattern of their evolution (Supplementary Fig. 3F) with joint high expression of HBG2, MYC, CD79B, and MCL1 (Supplementary Fig. 3G). By contrast, the type I origin will reduce the expression of CDK6, ITGB1, BCL2L1, and STAT1 during evolution, whereas the type IX origin will promote their expression. Additionally, the type IX origin shows precursory potential for further evolution, in which the local cell population takes the lead in showing the gene expression pattern of subclone type VIII, with higher evolution efficiency. In contrast, the type I origin maintained a similar expression pattern in some cell populations cloned by the offspring, resulting in relatively slower evolution. These different evolutionary transition forms and efficiency indicate that the type IX origin can evolve more cell subclonal clusters than the type I origin within the same timeframe.
We further explored early immune escape mechanisms, as the foundation for immune resistance and evasion, throughout the malignant process of MM. Immune checkpoints PD-1, PD-L1, and CTLA4 were negatively expressed in malignant plasma and microenvironment cells (Fig. 4F). The expression pattern of genes related to tumor immune killing in the microenvironment cells (Fig. 4G) was used to build a high-confidence immune escape intercellular communication network in the early stage of MM malignant transformation (Fig. 4H), providing insight into the potential mechanism (Fig. 4I) from innate to specific immunity. Among them, HLA-DMA and HLA-DRA were generally missing in antigen-presenting cells (B cells, classical dendritic cells (cDCs), and monocytes), which is the first weakness in the immune response cascade against MM (P.adj < 0.0001). CCL5 expression was also inhibited in CTLs and NK cells (P.adj < 0.0001), which significantly reduces the efficiency of recruiting DCs and more CTLs into the nidus . Additionally, defective expression of the perforin PRF1 and the granzyme GZMB in NK cells is the main reason for the inherent antitumor immune inactivation in MM patients (P.adj < 0.0001). Moreover, the expression of CD38 and CD27 in the B cell lineage (B, ProgB-I, ProgB-II, and plasma) was inhibited (P.adj < 0.0001), indicating deficiency of humoral immunity, while low levels of CCL3L1 also suggest low immune activity of monocytes (P.adj < 0.0001) .
Relapse and drug resistance in MM patients at the single-cell level
We evaluated the inhibitory effect of drugs on the expression of the established target genes at the patient level (Fig. 5A) to infer the part of the drug action that failed. The expression levels of target genes for thalidomide, melphalan, lenalidomide, and cyclophosphamide in the RRMM group did not show significant differences from those of the NDMM group, whereas those for dexamethasone, bortezomib, and doxorubicin decreased in sectional cell subpopulations. It is clear that the drug resistance of RRMM1 mainly occurs prior to the action on target genes, whereas the drug resistance of RRMM2 mainly occurs after interaction with the target gene. RRMM3 shares these two drug resistance mechanisms due to the extensive application of drugs. The characteristics of drug sensitivity at the patient level are also reflected at the cell subpopulation level, as demonstrated by systematic comparison between the cell subpopulations of the RRMM and NDMM groups (Fig. 1A, Fig. 5B). In addition, we observed the increased expression of genes such as IGHG3, IGLC2, and IGHG2 and the decreased expression of genes such as IGHG4, ITM2C, and IGHA1, which mediates the formation and development of malignant clones, causing heterogeneity in patients. The signal score of malignant cells was applied to estimate the pathway activity (Fig. 5C), in which apoptosis and the FOXO/p53 signaling pathways were inhibited in RRMM, while the calcium/Rap1/JAK-STAT/VEGF/mTOR signaling pathways related to survival, proliferation, migration, and stem cell characteristics were activated, which are mainly driven by CNV events supplemented by SNVs and TFs (Fig. 5D). This is similar to the pattern of the malignant origin. During clonal evolution of malignant cells to acquire drug resistance, the transformation of gene expression (Fig. 5E) was clustered into two gene expression modules, involving four evolutionary patterns. The acquisition and maintenance of the drug resistance phenotype of malignant cells require high expression of gene expression module 1 (related to drug sensitivity) and relatively low expression of gene expression module 2 (related to the development advantage in the natural state). The genes of both modules are related to the above-mentioned drug resistance signals, which are significantly regulated by CNV events as the main driving force for the clonal evolution of drug resistance.
The global regulatory network of intercellular communication (Fig. 5F) showed that activation of the drug metabolism signal pathway in B cells and cDCs may be reprogrammed by malignant cells, whereas the drug resistance evolution of malignant cells may involve microenvironment cells such as B cells, plasma cells, T cells, and cDCs. Genes in the drug resistance global regulatory network also showed significant prognostic potential in the training cohort (Fig. 5G). Among them, IL5RA, KRAS, and PPP2R5C were independently linked to prognosis (Fig. 5G); thus, a drug resistance-related prognostic model for MM based on multivariate Cox regression analysis was established for OS and RFS, which was verified in an independent cohort (Fig. 5H). Time-independent ROC curves demonstrated that the prognostic model enabled to accurately predict MM patients’ OS and RFS (Fig. 5I).
Intratumoral cellular heterogeneity in MM patients
We applied the tumor cell-specific transcriptional diversity score to measure the intratumoral heterogeneity of MM (Supplementary Fig. 4A), detailed in the single-cell profile map of malignant cells from each patient (Supplementary Fig. 4B), which demonstrates diverse and specific ecological components with vague traces of two malignant origins. Most of the malignant cells in MM patients exhibited the activated antigen processing and presentation (Fig. 6A), although the genes involved were different, indicating that malignant plasma cells retained incomplete antigen presentation ability to a certain extent. In particular, compared with other patients, for patient RRMM3, more distinctive malignant cell ecological components were found, regulated by TFs such as E2F7, E2F8, and SAP30 (Supplementary Fig. 4C) and leading to greater heterogeneity. SOM analysis showed the contribution (Fig. 6B) and concrete manifestations (Supplementary Fig. 4D) of malignant cell subpopulations to patient heterogeneity, of which NDMM and RRMM patients showed distinct SOM neural units, suggesting that different gene clusters respond differently to the drug sensitivity of MM. Compared with malignant cells, patient heterogeneity of microenvironment cells was relatively obscure (Supplementary Fig. 4E), which was also observed at the molecular level (Fig. 6C).
Immune cell subpopulation atlas of MM
The single-cell profiles (Fig. 6D) of B, plasma, CTLs, T, and NK cells were obtained by clustering the immune cell subpopulations, so that the heterogeneity of the immune microenvironment in MM patients was amplified at the subpopulation level. Interestingly, we found that plasma, CTLs, T, and NK cells highly expressed markers of MM malignant clonal evolution as subpopulation-specific markers, reflecting the immune cell heterogeneity of patients. Among them, IGHG3 was mainly expressed in lymphocytes of RRMM1, involving plasma-1, CTL-0, T-2, and NK-0, whereas IGLC2 was specifically expressed only in RRMM2 lymphocytes, related to CTL-1, T-5, and NK-2. Alternatively, IGHG2 was specifically expressed in the lymphocytes of CTL-2 and NK-4 in RRMM3. IGHG4 was specifically expressed in the lymphocytes of NDMM1 and NDMM3, touching upon CTL-3, T-4, and NK-4. The molecular specific markers of these microenvironment lymphocyte subpopulations are consistent with the distinctive expression of patient heterogeneous malignant single cells (Supplementary Fig. 4D). This provides strong evidence that clonal evolution of MM can reprogram lymphocyte-mediated patient heterogeneity, inspiring the novel strategy of individualized treatment in clinics. As this is an unprecedented discovery, we ruled out the possibility of wrong cell type identification for the sake of caution. The CTLs, T, and NK cells showed positive expression of their corresponding specific marker genes, which have significant cell co-localization with the corresponding patient heterogeneous malignant marker immunoglobulins (Supplementary Fig. 5A-C). Importantly, the same phenomenon was observed in cDCs, monocytes (Fig. 6E), and HSCs (Supplementary Fig. 5D). The transformation of these cell phenotypes from normal to depleted was accompanied with transformation of a series of gene expression patterns (Supplementary Fig. 5E, F) with various degrees of cell death signal activation, involving ferroptosis, necroptosis, and apoptosis (Fig. 6F, G). This raises the question of how malignant plasma cell markers appear in microenvironment cells to reprogram the microenvironment. By reviewing the biological signals activated by these depleted cells, we found that they not only activated various cell death signals, but also activated phagocytosis-related pathways, such as phagosome, Fc gamma, R-mediated phagocytosis (Fig. 6F, G). Correspondingly, during the clonal evolution of MM malignant cells, the key gene of exosome synthesis CD63 was expressed and a series of biological pathways related to vesicle synthesis and secretion were activated to varying degrees (Supplementary Fig. 5G, Fig. 3I). It is axiomatic that the MM malignant cells reprogramming microenvironment to augment immune escape is involved in vesicle synthesis and secretion of malignant cells as well as phagocytosis of microenvironment cells during the clonal evolution process, resulting in the emergence of malignant marker mRNA and mediating apoptosis and depletion.
In contrast, the positive expression of GNLY in T-1 emphasizes the potential malignant events that effectively activate the immune response of T cells. T-1 mainly existed in NDMM, concentrated in NDMM2, which may be related to the negative expression of malignant clonal evolution markers in the lymphocytes of NDMM2 patients, further indicating that chemoresistance of MM not only resists drugs but also inhibits immune killing.
Importantly, we found the negative expression of malignant clonal evolution markers in B cells, implying certain defensive means of B cells to resist the reprogramming and transformation during malignant cell clonal evolution, suggesting that a cellular immunotherapy strategy may benefit MM patients. CCL5 and CCL4 were specifically expressed in B-0 (Fig. 6D, F), which activate antigen presentation pathways and PD-L1/PD-1 immune checkpoints in cancer to recover the adverse effect of low CCL5 expression reducing the recruitment efficiency for cells in immune killing cascades (e.g., cDCs and CTLs). Surprisingly, B-0 were widely distributed in MM patients, mostly in the NDMM group (concentrated in NDMM2 patients), but not in the control, suggesting that an immunotherapy strategy of injecting CCL5-positive B cells derived from autologous expansion in vitro may be beneficial. Likewise, T-1 with potential immune killing efficacy highly expressed CCL5 and CCL4, activating antigen processing/presentation and chemokine signaling pathways to augment leukocyte transdermal migration. There are similar patterns of molecular characteristics and gene expression of T-1 in T-3 that were specific to controls, further highlighting the important role of CCL5 and CCL4 in the antitumor immunity of MM.
The unique challenges of complexity of tumor cell ecology have hindered progress in clonal evolution and intratumoral heterogeneity of MM. However, high-resolution single-cell sequencing technology has enabled the study of tumor evolution. In this study, we characterized the tumor cell and microenvironment cell ecosystem for the BM of MM patients, clarifying the potential of malignant clonal evolution patterns and their reprogramming microenvironment cells.
Moreover, we observed that from the beginning of malignant origin, the shared mutation of IFITM2 in all malignant subpopulations enhances its expression. Highly expressed IFITM2 participates in the type I interferon response to protect cells from viral pathogens and activates the expression and secretion of IL-6, which promotes the differentiation of B cells into plasma cells, mediating the growth of myeloma and stimulating bone resorption [44, 54]. Meanwhile, the activation of IFITM2 and IL-6 forms a reciprocally positive feedback loop to make MM malignant progenitor cells act out the expression pattern with virus infection. This builds a suitable host environment and microenvironment for the virus that is conducive to viral infection, such as human herpesvirus type 8 (HHV-8), thus increasing the susceptibility risk by approximately 10-fold [55,56,57]. HHV-8 and other viruses can continuously expand and encode the homolog of IL-6 in the host malignant plasma cells, which further activates the positive feedback loop to promote viral spread and the malignant process of MM . Notably, the viral susceptibility caused by this positive feedback loop is a high-risk factor, leading to the serious consequences and high mortality of MM, similar to the recent outbreak of SARS-CoV-2 [59, 60]. Likewise, ANK1 was also identified to mutate in the early origin stage of MM, resulting in the presence of the mutation in all malignant subpopulations. ANK1 participates in the malignant progression of acute myeloid leukemia in the elderly and in children with Down syndrome, related to erythropoietin, which provides a potential molecular link for MM secondary acute myeloid leukemia and mediating poor prognosis [4, 59, 60]. In particular, according to the tumor cell stemness index score and pseudo-sequential clonal evolution analysis, the malignant origins were divided into type I and IX malignant progenitor cells. Compared with the type I origin, more significant stem cell characteristics were found in type IX malignant progenitor cells, leading to a more severe degree of malignancy and progression. It is becoming increasingly clear that MM is not a single clonal genome, but rather distinct subclones that evolve from one or more origins in the disease period, leading to intratumoral heterogeneity and different performance in environmental adaptability. The increase in clonality also depends on the accumulation of gene mutations in offspring cells to regulate the environmental adaptive competition between malignant cell subpopulations. A portion of malignant cells maintain the ability for self-renewal and long-term clonal maintenance, becoming the dominant population in the current state, while others are in a dormant state [61, 62].
Significantly, the clonal evolution of MM malignant plasma cells is not only under natural selection from the microenvironment, but also capable of reacting to the TME to create more suitable conditions for survival, including the reprogramming of immune cells to inhibit immune killing and promote drug metabolism through vesicle synthesis and secretion of malignant plasma cells as well as the phagocytosis of microenvironment cells, which has been reported in other tumors [63, 64]. Interestingly, the microenvironment reprogramming produced by clonal evolution showed significant patient specificity, mediating the tumor heterogeneity of patients, which is an unprecedented discovery that warrants full consideration of personalized and accurate treatment in clinical immune salvation therapy for MM. Specifically, the BM microenvironment of NDMM2 patients with abundant CCL5-positive B and T cells can be protected from reprogramming with an antagonistic response; thus, relevant drugs and cell therapy schemes may be of great benefit. According to previous studies, CCL5 should also be expressed in malignant cells to promote the secretion of CXCL9 by immune cells through epigenetic regulation, thus jointly recruiting T cells, whose infiltration and activation exert immune killing effects to inhibit tumor progression ; ligand-receptor interaction between CCL5 and CCR1 is notably active in the monocyte-plasma communication in MM patients with stage III . However, we found negative expression of CCL5 in almost all of the malignant plasma cells of MM patients, which proved to be the main cause of insufficient CD8+ T cell infiltration and exhaustion, rather than the role of classical immune checkpoints, such as PD-1, PD-L1, and CTLA4. This emphasizes that the immune escape mediated by the downregulation of CCL5 expression occurs not only in solid tumors but also in MM, offering a new biomarker and candidate target of MM immunotherapy. As a result, a new adjuvant strategy to improve the sensitivity of immune checkpoint inhibitors for MM by jointly reversing CCL5 silencing drugs such as decitabine may be effective . The new adjuvant strategy extends the indications of therapy combined with epigenetic agent and immune checkpoint blockage in MM for further treatment of RRMM. Decitabine combined with PD-1/PD-L1 inhibitors has recently entered clinical trials for patients with hematological tumors such as relapsed/refractory classic Hodgkin lymphoma, which achieved satisfactory benefits with high safety . Alternatively, B cells have not been reprogrammed as other microenvironment immune cells with potential defensive measures, which warrants further research as a novel immunotherapy strategy.
In summary, we explored the clonal evolution of MM occurrence and development, mediating intratumoral cell heterogeneity, and the connectivity from initial treatment to drug resistance. We then identified the interaction and response between malignant plasma cells and the microenvironment during clonal evolution. Furthermore, this study broadens the cognitive boundary of MM; however, due to the limited sample size, a larger prospective cohort is required to obtain more universal and versatile conclusions.
Availability of data and materials
All data generated or analysed during this study are included in this published article [and its supplementary information files].
Single-cell RNA sequencing
Newly Diagnosed MM
Relapsed and/or refractory MM
- CD20+ RRMM:
Unique molecular identifiers
t-distributed stochastic neighbor embedding
DNA methylation-based stemness index
Gene expression-based stemness index
Single cell regulatory network inference and clustering
Regulon activity score
Connection specificity index
Differentially expressed genes
Copy number variation
Tumor mutational burden
Log fold change
Receiver operator characteristic
Hematopoietic stem cells
Classical dendritic cells
Herpesvirus type 8
Teras LR, DeSantis CE, Cerhan JR, Morton LM, Jemal A, Flowers CR. 2016 US lymphoid malignancy statistics by World Health Organization subtypes. CA Cancer J Clin. 2016;66(6):443–59.
Kumar SK, Rajkumar V, Kyle RA, van Duin M, Sonneveld P, Mateos MV, et al. Multiple myeloma. Nat Rev Dis Primers. 2017;3:17046.
Rollig C, Knop S, Bornhauser M. Multiple myeloma. Lancet. 2015;385(9983):2197–208.
Morgan GJ, Walker BA, Davies FE. The genetic architecture of multiple myeloma. Nat Rev Cancer. 2012;12(5):335–48.
Landau HJ, Yellapantula V, Diamond BT, Rustad EH, Maclachlan KH, Gundem G, et al. Accelerated single cell seeding in relapsed multiple myeloma. Nat Commun. 2020;11(1):3617.
Rasche L, Chavan SS, Stephens OW, Patel PH, Tytarenko R, Ashby C, et al. Spatial genomic heterogeneity in multiple myeloma revealed by multi-region sequencing. Nat Commun. 2017;8(1):268.
Johnson DC, Lenive O, Mitchell J, Jackson G, Owen R, Drayson M, et al. Neutral tumor evolution in myeloma is associated with poor prognosis. Blood. 2017;130(14):1639–43.
Maley CC, Aktipis A, Graham TA, Sottoriva A, Boddy AM, Janiszewska M, et al. Classifying the evolutionary and ecological features of neoplasms. Nat Rev Cancer. 2017;17(10):605–19.
Shen YJ, Mishima Y, Shi J, Sklavenitis-Pistofidis R, Redd RA, Moschetta M, et al. Progression signature underlies clonal evolution and dissemination of multiple myeloma. Blood. 2021;137(17):2360–72.
Zeissig MN, Zannettino ACW, Vandyke K. Tumour dissemination in multiple myeloma disease progression and relapse: a potential therapeutic target in high-risk myeloma. Cancers (Basel). 2020;12(12):3643.
Ledergor G, Weiner A, Zada M, Wang SY, Cohen YC, Gatt ME, et al. Single cell dissection of plasma cell heterogeneity in symptomatic and asymptomatic myeloma. Nat Med. 2018;24(12):1867–76.
Cohen YC, Zada M, Wang SY, Bornstein C, David E, Moshe A, et al. Identification of resistance pathways and therapeutic targets in relapsed multiple myeloma patients through single-cell sequencing. Nat Med. 2021;27(3):491–503.
Waldschmidt JM, Kloeber JA, Anand P, Frede J, Kokkalis A, Dimitrova V, et al. Single-cell profiling reveals metabolic reprogramming as a resistance mechanism in BRAF-mutated multiple myeloma. Clin Cancer Res. 2021;27(23):6432–44.
Gerlinger M, Rowan AJ, Horswell S, Math M, Larkin J, Endesfelder D, et al. Intratumor heterogeneity and branched evolution revealed by multiregion sequencing. N Engl J Med. 2012;366(10):883–92.
McGranahan N, Swanton C. Clonal heterogeneity and tumor evolution: past, present, and the future. Cell. 2017;168(4):613–28.
McGranahan N, Swanton C. Biological and therapeutic impact of intratumor heterogeneity in cancer evolution. Cancer Cell. 2015;27(1):15–26.
Bolli N, Avet-Loiseau H, Wedge DC, Van Loo P, Alexandrov LB, Martincorena I, et al. Heterogeneity of genomic evolution and mutational profiles in multiple myeloma. Nat Commun. 2014;5:2997.
Mishima Y, Paiva B, Shi J, Park J, Manier S, Takagi S, et al. The mutational landscape of circulating tumor cells in multiple myeloma. Cell Rep. 2017;19(1):218–24.
Tanay A, Regev A. Scaling single-cell genomics from phenomenology to mechanism. Nature. 2017;541(7637):331–8.
Deamer D, Akeson M, Branton D. Three decades of nanopore sequencing. Nat Biotechnol. 2016;34(5):518–24.
Wick RR, Judd LM, Holt KE. Performance of neural network basecalling tools for Oxford Nanopore sequencing. Genome Biol. 2019;20(1):129.
Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics. 2013;29(1):15–21.
Butler A, Hoffman P, Smibert P, Papalexi E, Satija R. Integrating single-cell transcriptomic data across different conditions, technologies, and species. Nat Biotechnol. 2018;36(5):411–20.
Aran D, Looney AP, Liu L, Wu E, Fong V, Hsu A, et al. Reference-based analysis of lung single-cell sequencing reveals a transitional profibrotic macrophage. Nat Immunol. 2019;20(2):163–72.
van Galen P, Hovestadt V, Wadsworth Ii MH, Hughes TK, Griffin GK, Battaglia S, et al. Single-cell RNA-Seq reveals AML hierarchies relevant to disease progression and immunity. Cell. 2019;176(6):1265–1281 e1224.
Hanzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics. 2013;14:7.
Colaprico A, Silva TC, Olsen C, Garofano L, Cava C, Garolini D, et al. TCGAbiolinks: an R/Bioconductor package for integrative analysis of TCGA data. Nucleic Acids Res. 2016;44(8):e71.
Malta TM, Sokolov A, Gentles AJ, Burzykowski T, Poisson L, Weinstein JN, et al. Machine learning identifies Stemness features associated with oncogenic dedifferentiation. Cell. 2018;173(2):338–354 e315.
Trapnell C, Cacchiarelli D, Grimsby J, Pokharel P, Li S, Morse M, et al. The dynamics and regulators of cell fate decisions are revealed by pseudotemporal ordering of single cells. Nat Biotechnol. 2014;32(4):381–6.
La Manno G, Soldatov R, Zeisel A, Braun E, Hochgerner H, Petukhov V, et al. RNA velocity of single cells. Nature. 2018;560(7719):494–8.
Aibar S, Gonzalez-Blas CB, Moerman T, Huynh-Thu VA, Imrichova H, Hulselmans G, et al. SCENIC: single-cell regulatory network inference and clustering. Nat Methods. 2017;14(11):1083–6.
Gu Z, Eils R, Schlesner M. Complex heatmaps reveal patterns and correlations in multidimensional genomic data. Bioinformatics. 2016;32(18):2847–9.
Efremova M, Vento-Tormo M, Teichmann SA, Vento-Tormo R. CellPhoneDB: inferring cell-cell communication from combined expression of multi-subunit ligand-receptor complexes. Nat Protoc. 2020;15(4):1484–506.
Wang Y, Wang R, Zhang S, Song S, Jiang C, Han G, et al. iTALK: an R Package to Characterize and Illustrate Intercellular Communication. bioRxiv. 2019:507871. https://doi.org/10.1101/507871.
Yu G, Wang LG, Han Y, He QY. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS. 2012;16(5):284–7.
Li H. Minimap2: pairwise alignment for nucleotide sequences. Bioinformatics. 2018;34(18):3094–100.
Sedlazeck FJ, Rescheneder P, Smolka M, Fang H, Nattestad M, von Haeseler A, et al. Accurate detection of complex structural variations using single-molecule sequencing. Nat Methods. 2018;15(6):461–8.
Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, et al. Genome project data processing S: the sequence alignment/map format and SAMtools. Bioinformatics. 2009;25(16):2078–9.
Tirosh I, Venteicher AS, Hebert C, Escalante LE, Patel AP, Yizhak K, et al. Single-cell RNA-seq supports a developmental hierarchy in human oligodendroglioma. Nature. 2016;539(7628):309–13.
Liu F, Zhang Y, Zhang L, Li Z, Fang Q, Gao R, et al. Systematic comparative analysis of single-nucleotide variant detection methods from single-cell RNA sequencing data. Genome Biol. 2019;20(1):242.
Hellmann MD, Ciuleanu TE, Pluzanski A, Lee JS, Otterson GA, Audigier-Valette C, et al. Nivolumab plus Ipilimumab in lung Cancer with a high tumor mutational burden. N Engl J Med. 2018;378(22):2093–104.
Müller B, Reinhardt J, Strickland MT. KOHOMAP: The Kohonen Self-organizing Map. In: Neural Networks: An Introduction. Berlin: Springer Berlin Heidelberg; 1995. p. 291–5.
Ma L, Hernandez MO, Zhao Y, Mehta M, Tran B, Kelly M, et al. Tumor cell biodiversity drives microenvironmental reprogramming in liver Cancer. Cancer Cell. 2019;36(4):418–430 e416.
Schoggins JW, Wilson SJ, Panis M, Murphy MY, Jones CT, Bieniasz P, et al. A diverse range of gene products are effectors of the type I interferon antiviral response. Nature. 2011;472(7344):481–5.
Youk J, An Y, Park S, Lee JK, Ju YS. The genome-wide landscape of C:G > T:A polymorphism at the CpG contexts in the human population. BMC Genomics. 2020;21(1):270.
Zhou Y, He F, Pu W, Gu X, Wang J, Su Z. The impact of DNA methylation dynamics on the mutation rate during human germline development. G3 (Bethesda). 2020;10(9):3337–46.
Yuan S, Liu Z, Xu Z, Liu J, Zhang J. High mobility group box 1 (HMGB1): a pivotal regulator of hematopoietic malignancies. J Hematol Oncol. 2020;13(1):91.
Martin-Garcia D, Navarro A, Valdes-Mas R, Clot G, Gutierrez-Abril J, Prieto M, et al. CCND2 and CCND3 hijack immunoglobulin light-chain enhancers in cyclin D1(−) mantle cell lymphoma. Blood. 2019;133(9):940–51.
Huang X, Di Liberto M, Jayabalan D, Liang J, Ely S, Bretz J, et al. Prolonged early G(1) arrest by selective CDK4/CDK6 inhibition sensitizes myeloma cells to cytotoxic killing through cell cycle-coupled loss of IRF4. Blood. 2012;120(5):1095–106.
Dilworth D, Liu L, Stewart AK, Berenson JR, Lassam N, Hogg D. Germline CDKN2A mutation implicated in predisposition to multiple myeloma. Blood. 2000;95(5):1869–71.
Wen Z, Rajagopalan A, Flietner ED, Yun G, Chesi M, Furumo Q, et al. Expression of NrasQ61R and MYC transgene in germinal center B cells induces a highly malignant multiple myeloma in mice. Blood. 2021;137(1):61–74.
Seo W, Shimizu K, Kojo S, Okeke A, Kohwi-Shigematsu T, Fujii SI, et al. Runx-mediated regulation of CCL5 via antagonizing two enhancers influences immune cell function and anti-tumor immunity. Nat Commun. 2020;11(1):1562.
Shalekoff S, Meddows-Taylor S, Schramm DB, Donninger SL, Gray GE, Sherman GG, et al. Host CCL3L1 gene copy number in relation to HIV-1-specific CD4+ and CD8+ T-cell responses and viral load in south African women. J Acquir Immune Defic Syndr. 2008;48(3):245–54.
Xu L, Zhou R, Yuan L, Wang S, Li X, Ma H, et al. IGF1/IGF1R/STAT3 signaling-inducible IFITM2 promotes gastric cancer growth and metastasis. Cancer Lett. 2017;393:76–85.
Ma HJ, Sjak-Shie NN, Vescio RA, Kaminsky M, Mikail A, Pold M, et al. Human herpesvirus 8 open reading frame 26 and open reading frame 65 sequences from multiple myeloma patients: a shared pattern not found in Kaposi's sarcoma or primary effusion lymphoma. Clin Cancer Res. 2000;6(11):4226–33.
Dhakal B, D'Souza A, Chhabra S, Hari P. Multiple myeloma and COVID-19. Leukemia. 2020;34(7):1961–3.
Blimark C, Holmberg E, Mellqvist UH, Landgren O, Bjorkholm M, Hultcrantz M, et al. Multiple myeloma and infections: a population-based study on 9253 multiple myeloma patients. Haematologica. 2015;100(1):107–13.
Burger R, Neipel F, Fleckenstein B, Savino R, Ciliberto G, Kalden JR, et al. Human herpesvirus type 8 interleukin-6 homologue is functionally active on human myeloma cells. Blood. 1998;91(6):1858–63.
Shaham L, Vendramini E, Ge Y, Goren Y, Birger Y, Tijssen MR, et al. MicroRNA-486-5p is an erythroid oncomiR of the myeloid leukemias of Down syndrome. Blood. 2015;125(8):1292–301.
Whitman SP, Maharry K, Radmacher MD, Becker H, Mrozek K, Margeson D, et al. FLT3 internal tandem duplication associates with adverse outcome and gene- and microRNA-expression signatures in patients 60 years of age or older with primary cytogenetically normal acute myeloid leukemia: a Cancer and leukemia group B study. Blood. 2010;116(18):3622–6.
Kreso A, Dick JE. Evolution of the cancer stem cell model. Cell Stem Cell. 2014;14(3):275–91.
Salehi S, Kabeer F, Ceglia N, Andronescu M, Williams MJ, Campbell KR, et al. Clonal fitness inferred from time-series modelling of single-cell cancer genomes. Nature. 2021;595(7868):585–90.
Ombrato L, Nolan E, Kurelac I, Mavousian A, Bridgeman VL, Heinze I, et al. Metastatic-niche labelling reveals parenchymal cells with stem features. Nature. 2019;572(7771):603–8.
Chen S, Zhu G, Yang Y, Wang F, Xiao YT, Zhang N, et al. Single-cell analysis reveals transcriptomic remodellings in distinct cell types that contribute to human prostate cancer progression. Nat Cell Biol. 2021;23(1):87–98.
Dangaj D, Bruand M, Grimm AJ, Ronet C, Barras D, Duttagupta PA, et al. Cooperation between constitutive and inducible chemokines enables T cell engraftment and immune attack in solid tumors. Cancer Cell. 2019;35(6):885–900 e810.
Zhong L, Yang X, Zhou Y, Xiao J, Li H, Tao J, et al. Exploring the R-ISS stage-specific regular networks in the progression of multiple myeloma at single-cell resolution. Sci China Life Sci. 2022;65(9):1811–23.
Nie J, Wang C, Liu Y, Yang Q, Mei Q, Dong L, et al. Addition of low-dose Decitabine to anti-PD-1 antibody Camrelizumab in relapsed/refractory classical Hodgkin lymphoma. J Clin Oncol. 2019;37(17):1479–89.
This work was funded by National Natural Science Foundation of China (81873450, 82170181); Beijing Municipal Administration of Hospitals’ Youth Program (QMS20200201).
Ethics approval and consent to participate
The study was approved by the Ethics Committee of Beijing Tongren Hospital and Sun Yat-sen University Cancer Center. The research has been carried out in accordance with the World Medical Association Declaration of Helsinki. All methods were carried out in accordance with relevant guidelines and regulations. All patients provided written informed consent.
Consent for publication
All patients provided written informed consent.
The authors declare no conflicts of interest.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
BM intercellular communication events in patients with NDMM and RRMM. (A and B) Expression patterns of CD19 and CD20 in single cells from NDMM (A) and RRMM (B) patients. Top: t-SNE single-cell atlas showing the expression of CD19 and CD20 in BM single cells from MM patients. Bottom: Verification of the expression of CD19 and CD20 in BM cells by flow cytometry. (C) Expression pattern of CD19 in single cells from CD20+ RRMM patients. Top: t-SNE single-cell atlas showing the expression of CD19 in BM single cells from CD20+ RRMM patients. Bottom: Verification of CD19 expression in BM cells from CD20+ RRMM patients by flow cytometry. (D) High-confidence receptor-ligand interaction for BM intercellular communication events in patients with NDMM and RRMM. Bubbles represent the average expression of high-confidence ligand-receptor interaction pairs in the source and target cells. (E) Overview of high-confidence intercellular communication networks in BM patients with NDMM and RRMM. Each bubble represents a cellular identity of the BM in NDMM and RRMM patients, and the coloring is consistent with that of the corresponding single-cell atlas. Bubble size represents active cell communication with other cells. Each arrow represents the interaction between the source cell ligand and the target cell receptor, and its thickness represents the number of ligand-receptor interaction pairs. All communication events were detected using any two tools among CellPhoneDB, iTalk, and CellCrosstalk, which represent the confidence of the reciprocal pair. BM, bone marrow; NDMM, newly diagnosed multiple myeloma; RRMM, relapsed and/or refractory MM; t-SNE, t-distribution and stochastic neighbor embedding.
Nanopore sequencing combined with single-cell transcriptome to identify malignant driving events. (A) Diagram of genomic SVs in patients with NDMM and RRMM. “INS” and “INV” of patients with counts > 5 are displayed, and colors on chromosomes represent the number of patients with DEL (green to yellow) and DUP (green to blue). (B) Patterns of gene mutations in MM patients. (C) Volume and single-cell levels of the TMB in the BM of MM patients. Top: TMB levels of different malignant subclones. Middle: Gene mutation patterns in different malignant subclones. Bottom: TMB horizontal curve of patients with MM. (D) Correlation of TMB in MM patients at the bulk and single-cell levels. The BM TMB of MM patients was significantly correlated at the bulk level and at the single-cell level (P < 0.001), demonstrating the feasibility and reliability of detecting the SNV spectrum and TMB level at the single-cell level. SV, structural variation; INS, insert; INV, inversion; DEL, deletion; DUP, duplication; CNV, copy number variation; SNV, single nucleotide variation; TMB, tumor mutational burden.
Exploring malignant origins from the perspective of clonal evolution. (A) Pseudotime trajectory of myeloid malignant subclones in NDMM and RRMM patients mapped in the t-SNE–based single-cell atlas. (B, C) Stemness score (B) and cell cycle score (C) of myeloid malignant subclones in patients with NDMM and RRMM. Top: Comparison of stemness and cell cycle scores between myeloid malignant subclones from patients with NDMM and RRMM. Bottom: Stemness and cell cycle scores of myeloid malignant subclones in NDMM and RRMM patients mapped using the t-SNE single-cell atlas. (D) Dysregulated expression patterns of malignant leader genes in type I and IX malignant origins compared to normal B cells. (E) Type I and IX malignant origin expression dysregulated genes are significantly involved in malignant precursor biological signals. These biological signals of malignant precursors can be grouped into five modules according to shared gene members. Each module was assigned an independent color. (F) Gene expression patterns of the evolution of malignant origins toward dominant subclones during natural development. (G) Expression patterns of dysregulated genes across Type I, IV, VIII and IX.
Cellular heterogeneity in the myeloma of MM patients. (A) Heterogeneity score of the myeloma in NDMM and RRMM patients. (B) Single-cell atlas based on t-SNE showing heterogeneity of malignant subclonal ecological components in different MM patients. (C) Specific markers of malignant subclones in patients with NDMM and RRMM. Left: Expression pattern of malignant subclonal-specific marker genes in patients. Right: TFs that regulate malignant subclonal-specific markers. (D) Contribution of cell subsets with malignant subclone-specific markers in patients. Left: Expression patterns of malignant subclonal-specific marker genes in patients. The malignant subclonal marker genes illustrated are those that are the most significant and associated with subsequent MM malignant plasma cell clonal evolution reprogramming of immune cells in the microenvironment. Right: Contribution of cell subsets specifically labeled by the top specific marker malignant subclones in patients. (E) The single-cell atlas based on t-SNE showing the heterogeneity of cell ecological components in the myeloma microenvironment in different MM patients. TF, transcription factor.
Expression patterns of immune cell subpopulations in MM. (A-C) Expression patterns of known cell identity-specific markers and their co-location with malignant markers in lymphocytes provide evidence for MM reprogramming of microenvironment cells. The pie chart shows the co-location between known cell identity-specific markers and corresponding patient heterogeneous malignant marker immunoglobulins. (D) HSC subset profiles of control donor BM and the MM patient BM microenvironment. Left: Single-cell subpopulation atlas of HSCs. Middle: Proportion of HSCs in control donors and different MM patients. Right: Expression patterns of cell subset-specific markers mapped on the HSC single-cell atlas. (E, F) MM malignant clonal evolution reprograms pseudotime trajectories and expression pattern changes in lymphoid (E) and myeloid cells (F). Lymphocytes include B cells, plasma, CTLs, T cells, and NK cells, whereas myeloid cells include cDCs and monocytes. For each cell type, the pseudotime trajectory (left), pseudotime value change (middle), and expression pattern change (right) from normal cells to MM during MM malignant clonal evolution reprogramming are shown for lymphocyte and myeloid cell subsets. (G) Expression patterns of the exosome-specific marker gene CD63 in malignant cells. HSCs, hematopoietic stem cells; CTLs, cytotoxic T lymphocytes; NK cells, natural killer cells; cDCs, conventional dendritic cells.
Demographic, clinical, molecular, and diagnostic information of all patients in this study.
The known marker genes used for annotating cell clusters.
Specific marker genes of control donor BM cells. BM, bone marrow.
Specific marker genes of malignant plasma cell subclones in NDMM and RRMM patients. NDMM, newly diagnosed multiple myeloma; RRMM, relapsed and/or refractory MM.
Specific marker genes of BM tumor microenvironment cells in patients with NDMM and RRMM.
Characterization of specific marker genes for type I and type IX malignant progenitors.
Relative abundance estimation and typing of type I and IX malignant progenitors in a large clinical cohort.
About this article
Cite this article
Liang, Y., He, H., Wang, W. et al. Malignant clonal evolution drives multiple myeloma cellular ecological diversity and microenvironment reprogramming. Mol Cancer 21, 182 (2022). https://doi.org/10.1186/s12943-022-01648-z