Abstract
Multiple myeloma, a plasma cell malignancy, is the second most common blood cancer. Despite extensive research, disease heterogeneity is poorly characterized, hampering efforts for early diagnosis and improved treatments. Here, we apply single cell RNA sequencing to study the heterogeneity of 40 individuals along the multiple myeloma progression spectrum, including 11 healthy controls, demonstrating high interindividual variability that can be explained by expression of known multiple myeloma drivers and additional putative factors. We identify extensive subclonal structures for 10 of 29 individuals with multiple myeloma. In asymptomatic individuals with early disease and in those with minimal residual disease post-treatment, we detect rare tumor plasma cells with molecular characteristics similar to those of active myeloma, with possible implications for personalized therapies. Single cell analysis of rare circulating tumor cells allows for accurate liquid biopsy and detection of malignant plasma cells, which reflect bone marrow disease. Our work establishes single cell RNA sequencing for dissecting blood malignancies and devising detailed molecular characterization of tumor cells in symptomatic and asymptomatic patients.
This is a preview of subscription content, access via your institution
Access options
Access Nature and 54 other Nature Portfolio journals
Get Nature+, our best-value online-access subscription
$29.99 / 30 days
cancel any time
Subscribe to this journal
Receive 12 print issues and online access
$209.00 per year
only $17.42 per issue
Buy this article
- Purchase on SpringerLink
- Instant access to full article PDF
Prices may be subject to local taxes which are calculated during checkout
Similar content being viewed by others
Data availability
Raw and processed sc RNA-seq data and ChIP-seq data were deposited to the National Center for Biotechnology Information (NCBI)’s Gene Expression Omnibus with accession number GSE117156. Raw and processed genomic DNA targeted sequencing data were deposited to NCBI’s Sequence Read Archive with accession number SRP165705. Source code used for scRNA-seq analysis can be found at https://bitbucket.org/amitlab/multiple-myeloma-2018/. The participats’ clinical information is available in Supplementary Table 1.
References
Rajkumar, S. V. Multiple myeloma: 2016 update on diagnosis, risk-stratification, and management. Am. J. Hematol. 91, 719–734 (2016).
Kyle, R. A. et al. Monoclonal gammopathy of undetermined significance (MGUS) and smoldering (asymptomatic) multiple myeloma: IMWG consensus perspectives risk factors for progression and guidelines for monitoring and management. Leukemia 24, 1121–1127 (2010).
Morgan, G. J., Walker, B. A. & Davies, F. E. The genetic architecture of multiple myeloma. Nat. Rev. Cancer 12, 335–348 (2012).
Dhodapkar, M. V. MGUS to myeloma: a mysterious gammopathy of underexplored significance. Blood 128, 2599–2606 (2016).
Bolli, N. et al. A DNA target-enrichment approach to detect mutations, copy number changes and immunoglobulin translocations in multiple myeloma. Blood Cancer J. 6, e467 (2016).
Chapman, M. et al. Initial genome sequencing and analysis of multiple myeloma. Nature 471, 467–472 (2011).
Egan, J. et al. Whole-genome sequencing of multiple myeloma from diagnosis to plasma cell leukemia reveals genomic initiating events, evolution, and clonal tides. Blood 120, 1060–1066 (2012).
Lohr, J. G. et al. Widespread genetic heterogeneity in multiple myeloma: implications for targeted therapy. Cancer Cell 25, 91–101 (2014).
Walker, B. et al. Intraclonal heterogeneity and distinct molecular mechanisms characterize the development of t(4; 14) and t(11; 14) myeloma. Blood 120, 1077–1086 (2012).
Laganà, A. et al. Integrative network analysis identifies novel drivers of pathogenesis and progression in newly diagnosed multiple myeloma. Leukemia 32, 120–130 (2018).
Shah, V. et al. Prediction of outcome in newly diagnosed myeloma: a meta-analysis of the molecular profiles of 1905 trial patients. Leukemia 32, 102–110 (2018).
Shaughnessy, J. D. et al. A validated gene expression model of high-risk multiple myeloma is defined by deregulated expression of genes mapping to chromosome 1. Blood 109, 2276–2284 (2007).
Bolli, N. et al. Heterogeneity of genomic evolution and mutational profiles in multiple myeloma. Nat. Commun. 5, 2997 (2014).
Gawad, C., Koh, W. & Quake, S. R. Single-cell genome sequencing: current state of the science. Nat. Rev. Genet. 17, 175–188 (2016).
Patel, A. P. et al. Single-cell RNA-seq highlights intratumoral heterogeneity in primary glioblastoma. Science 344, 1396–1401 (2014).
Giladi, A. & Amit, I. Single-cell genomics: a stepping stone for future immunology discoveries. Cell 172, 14–21 (2018).
Paiva, B. et al. Differentiation stage of myeloma plasma cells: biological and clinical significance. Leukemia 31, 382–392 (2017).
Jaitin, D. A. et al. Massively parallel single-cell RNA-seq for marker-free decomposition of tissues into cell types. Science 343, 776–779 (2014).
Paul, F. et al. Transcriptional heterogeneity and lineage commitment in myeloid progenitors. Cell 163, 1663–1677 (2015).
Juneja, S., Viswanathan, S., Ganguly, M. & Veillette, C. A simplified method for the aspiration of bone marrow from patients undergoing hip and knee joint replacement for isolating mesenchymal stem cells and in vitro chondrogenesis. Bone Marrow Res. 2016, 1–18 (2016).
Halliley, J. et al. Long-lived plasma cells are contained within the CD19−CD38hiCD138+ subset in human bone marrow. Immunity 43, 132–145 (2015).
Chesi, M. et al. Frequent translocation t(4;14)(p16.3; q32.3) in multiple myeloma is associated with increased expression and activating mutations of fibroblast growth factor receptor 3. Nat. Genet. 16, 260–264 (1997).
Pawlyn, C. & Morgan, G. Evolutionary biology of high-risk multiple myeloma. Nat. Rev. Cancer 17, 543–556 (2017).
Combes, A. et al. BAD-LAMP controls TLR9 trafficking and signalling in human plasmacytoid dendritic cells. Nat. Commun. 8, 913 (2017).
Defays, A. et al. BAD-LAMP is a novel biomarker of nonactivated human plasmacytoid dendritic cells. Blood 118, 609–617 (2011).
Fathallah-Shaykh, H., Wolf, S., Wong, E., Posner, J. B. & Furneaux, H. M. Cloning of a leucine-zipper protein recognized by the sera of patients with antibody-associated paraneoplastic cerebellar degeneration. Proc. Natl Acad. Sci. USA 88, 3451–3454 (1991).
Hellström, I. et al. The HE4 (WFDC2) protein is a biomarker for ovarian carcinoma. Cancer Res. 63, 3695–3700 (2003).
Nutt, S. L., Hodgkin, P. D., Tarlinton, D. M. & Corcoran, L. M. The generation of antibody-secreting plasma cells. Nat. Rev. Immunol. 15, 160–171 (2015).
Kumar, S. K. & Rajkumar, S. V. The multiple myelomas—current concepts in cytogenetic classification and therapy. Nat. Rev. Clin. Oncol. 15, 409–421 (2018).
Rajan, A. M. & Rajkumar, S. V. Interpretation of cytogenetic results in multiple myeloma for clinical practice. Blood Cancer J. 5, e365 (2015).
Puig, N. et al. The predominant myeloma clone at diagnosis, CDR3 defined, is constantly detectable across all stages of disease evolution. Leukemia 29, 1435–1437 (2015).
Lefranc, M.-P. et al. IMGT, the international ImMunoGeneTics information system 25 years on. Nucleic Acids Res. 43, D413–D422 (2015).
Azizi, E. et al. Single-cell map of diverse immune phenotypes in the breast tumor microenvironment. Cell 174, 1293–1308 (2018).
Tian, E. et al. The role of the Wnt-signaling antagonist DKK1 in the development of osteolytic lesions in multiple myeloma. N. Engl. J. Med. 349, 2483–2494 (2003).
Zhao, X.-Y. Y. et al. Long noncoding RNA licensing of obesity-linked hepatic lipogenesis and NAFLD pathogenesis. Nat. Commun. 9, 2986 (2018).
Leidi, M., Mariotti, M. & Maier, J. Transcriptional coactivator EDF-1 is required for PPARγ-stimulated adipogenesis. Cell. Mol. Life Sci. 66, 2733–2742 (2009).
Simaite, D. et al. Recessive mutations in PCBD1 cause a new type of early-onset diabetes. Diabetes 63, 3557–3564 (2014).
Chen, X. et al. Prognostic value of diametrically polarized tumor-associated macrophages in multiple myeloma. Oncotarget 8, 112685–112696 (2017).
Dubovsky, J. et al. Lymphocyte cytosolic protein 1 is a chronic lymphocytic leukemia membrane-associated antigen critical to niche homing. Blood 122, 3308–3316 (2013).
Mishima, Y. et al. The mutational landscape of circulating tumor cells in multiple myeloma. Cell Rep. 19, 218–224 (2017).
Lohr, J. et al. Genetic interrogation of circulating multiple myeloma cells at single-cell resolution. Sci. Transl. Med. 8, 363ra147 (2016).
Manier et al. Whole-exome sequencing of cell-free DNA and circulating tumor cells in multiple myeloma. Nat. Commun. 9, 1691 (2018).
Rasche et al. Spatial genomic heterogeneity in multiple myeloma revealed by multi-region sequencing. Nat. Commun. 8, 268 (2017).
Picelli, S. et al. Smart-seq2 for sensitive full-length transcriptome profiling in single cells. Nat. Methods 10, 1096–1098 (2013).
Harris, P. A. et al. Research electronic data capture (REDCap)—a metadata-driven methodology and workflow process for providing translational research informatics support. J. Biomed. Inform. 42, 377–381 (2009).
Keren-Shaul, H. et al. A unique microglia type associated with restricting development of Alzheimer’s disease. Cell 169, 1276–1290 (2017).
Levine, J. H. et al. Data-driven phenotypic dissection of AML reveals progenitor-like cells that correlate with prognosis. Cell 162, 184–197 (2015).
Baran, Y. et al. MetaCell: analysis of single cell RNA-seq data using k-NN graph partitions. Preprint at bioRxiv https://doi.org/10.1101/437665 (2018).
Butler, A., Hoffman, P., Smibert, P., Papalexi, E. & Satija, R. Integrating single-cell transcriptomic data across different conditions, technologies, and species. Nat. Biotechnol. 36, 411–420 (2018).
Blecher-Gonen, R. et al. High-throughput chromatin immunoprecipitation for genome-wide mapping of in vivo protein-DNA interactions and epigenomic states. Nat. Protoc. 8, 539–554 (2013).
Kim, S. et al. Strelka2: fast and accurate calling of germline and somatic variants. Nat. Methods 15, 591–594 (2018).
Flores-Montero, et al. Next generation flow for highly sensitive and standardized detection of minimal residual disease in multiple myeloma. Leukemia 31, 2094 (2017).
Acknowledgements
We thank the participants and their families; clinical coordinators D. Yaish, I. Maman and S. Levy; N. Voskoboinik for iFISH analysis; K. Kogan for help with REDCap installation and A. Giladi for review of the manuscript. We thank Brian Fritz and Tarjei Mikkelsen from 10x Genomics for help and support with the Chromium single cell 5′ and V(D)J kits. This study was partly supported by the Benoziyo Family Fund, Clalit Health Care Services and Mrs. and Mr. Barry Lang. I.A. is supported by the Chan Zuckerberg Initiative (CZI); the Howard Hughes Medical Institute International Scholar award; the European Research Council Consolidator Grant (ERC-COG) 724471-HemTree2.0; Melanoma Research Alliance (MRA) Established Investigator Award (509044); the Israel Science Foundation (703/15); the Ernest and Bonnie Beutler Research Program of Excellence in Genomic Medicine; the Helen and Martin Kimmel award for innovative investigation; a Minerva Stiftung research grant; the Israeli Ministry of Science, Technology, and Space; the David and Fela Shapell Family Foundation; the NeuroMac DFG/Transregional Collaborative Research Center Grant,; and the Abramson Family Center for Young Scientists. A.T. is supported by the European Research Council (ERC) (scAssembly), the CZI, and the Flight Attendant Medical Research Institute. B.P. is supported by an ERC 2015 Starting Grant (MYELOMANEXT).
Author information
Authors and Affiliations
Contributions
G.L. conceived and designed the study, performed and analyzed experiments, and wrote the manuscript. A.W. designed the study, performed bioinformatic analyses and wrote the manuscript. M.Z. performed flow cytometry sorting experiments. S.-Y.W performed analysis of single cell RNA sequencing. Y.C.C., M.E.G., N.S., H.M., M.K.-M., K.H.-T., M.S., D.B.Y., A.N., L.S., I.Avivi. provided input on experimental design and collected clinical data and participant samples. H.K.-S. and I.Y. performed experiments and maintained RNA-seq infrastructure. C.B. performed ChIP-seq experiments. R.R. Analyzed clinical data using REDCap. E.D. and J.J.K. performed bioinformatic analyses. V.Y., E.P., O.L. designed targeted genomic sequencing and analyzed data. A.O.-U., K.B.H. and S.I. performed FISH analysis. S.K., J.S.-M. and B.P. helped with design and FACS analysis. G.I.B. coordinated and supervised clinical data. A.T. supervised the project, designed and analyzed experiments and wrote the manuscript. I.A. supervised the project, designed and analyzed experiments and wrote the manuscript.
Corresponding author
Ethics declarations
Competing interests
The authors declare the following competing interests: a patent application (US Provisional Patent Application No. 62/756,640) has been filed related to this work. O.L. declares receiving funding and/or honoraria from Adaptive, Amgen, Binding Site, BMS, Celgene, Cellectis, Glenmark, Janssen, Karyopharm, Pfizer, Seattle Genetics and Takeda.
Additional information
Publisher’s note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Extended data
Extended Data Fig. 1 Sorting strategy of plasma cells: representative flow cytometry plots showing sorting strategy (CD38+CD138+) for plasma cells after doublet exclusion.
Sorting experiments were performed for 40 individuals with multiple myeloma and controls. For 21 subjects, additional sorting was performed for circulating plasma cells from peripheral blood (PB) (separately). a, Bone marrow plasma cells of subject MM08 (active myeloma). b, Circulating plasma cells of subject MM04 (active myeloma). c, Bone marrow plasma cells of subject MGUS05. d, Bone marrow plasma cells of subject MM13 (active myeloma with minimal residual disease). Plots were generated using FlowJo software (Online Methods). During sorting, surface marker expression for additional markers in the Euroflow panel (CD19, CD81, CD27, CD56, CD117, CD45) was recorded for each single cell (Online Methods).
Extended Data Fig. 2 Quality control metrics of single cell sequencing: bone marrow and peripheral blood single cell quality control metrics.
a, Shown are number of reads, number of UMIs and percentage of cells analyzed per batch of 384 cells (that were pooled for library construction) for all CD38+CD138+ bone marrow single cells from 40 participants (29 newly diagnosed subjects and 11 control donors). Cells were sorted into plates, were sequenced and underwent quality control evaluation and filtering out of non-plasma cell contaminants (Online Methods). b, Shown are number of reads, number of UMIs and percentage of cells analyzed per batch of 384 cells (that were pooled for library construction) for all CD38+CD138+ peripheral blood single cells from 21 participants (19 newly diagnosed subjects and 2 control donors). Cells were sorted into plates, were sequenced and underwent quality control evaluation and filtering out of non-plasma cell contaminants (Online Methods). c, Table summarizing the total CD38+CD138+ cell numbers from bone marrow and blood of newly diagnosed individuals and control donors, that were sorted and sequenced, passed quality control and were analyzed (after filtering for non-plasma cell contaminants; see Online Methods).
Extended Data Fig. 3 Filtering strategy for plasma cells, and normalized single cell expression for selected genes.
a, Plot depicting plasma cell in silico filtering according to immunoglobulin load and UMI count per cell (Online Methods) b, Heat map showing clustering analysis of 3,179 ‘contaminating’ cells that pass quality control but do not express immunoglobulin genes above the cut-off (100 UMIs per cell). Representative genes (mostly unrelated to the plasma cell program) are shown. c, Normalized single cell expression (UMI count, log scale) of 16 representative genes, projected onto a tSNE map of all 20,586 bone marrow plasma cells from 40 participants (29 newly diagnosed individuals and 11 control donors).
Extended Data Fig. 4 Clustering analysis of single bone marrow plasma cells.
a, Heat map showing clustering analysis of 20,586 bone marrow plasma cells sorted from 40 participants (29 newly diagnosed individuals and 11 control donors), featuring normalized single cell expression levels of the 100 most variable genes (Online Methods). b, Dendrogram showing the hierarchical clustering of the average transcription profile for all clusters, C1–29 (related to a). c, Cell-to-cell correlation matrix (Spearman) of 200 randomly selected cells from each subject (with n > 200 cells; left); kNN adjacency matrix (k = 100) showing, for each cell, its 100 nearest neighbors (blue; right).
Extended Data Fig. 5 Genomic FISH examples from selected subjects, and epigenetic analysis of LAMP5.
a, Genomic iFISH images of bone marrow plasma cells magnetically enriched for CD138+ for subjects AL03 (left), MM08 (middle) and SMM01 (right). Shown are fluorescent probes for the immunoglobulin heavy chain (green: IGH) and translocation partners (red: FGFR3, CCND1 and CKS1B, from left to right). Each hybridization was performed twice with similar results. b, FACS plots showing intracellular staining (gated on live cells) for LAMP5 protein (red) compared with isotype control (blue) for KHM1B and RPMI-8226 cell lines (top and bottom, respectively; see Online Methods). Each experiment was performed twice with similar results. c, Genome browser view of normalized H3K4me2 profiles of peaks found in a 1-megabase region in the LAMP5 locus. Data are from two independent biological replicates (Online Methods).
Extended Data Fig. 6 Clonal assessment by inferring immunoglobulin sequence at single cell resolution.
a, Heat map depicting the relative frequency of the immunoglobulin sequences from analyzing 20,586 bone marrow plasma cells of 40 participants (29 newly diagnosed individuals and 11 control donors): top, immunoglobulin light chain variable region (IGKV and IGLV); middle, immunoglobulin light chain constant region (IGKC and IGLC); bottom, immunoglobulin heavy chain constant region (IGHC). b, Chromium 10× single cell BCR clonotype distribution of control donor Hip13, created by Chromium’s Loupe V(D)J browser (Online Methods). Left, heat map showing cell frequency for specific IGHV/IGLV/IGKV variable (V) region sequences (x axis) and IGHJ/IGKJ/IGLJ joining (J) region sequences (y axis). Right, frequency of different BCR clonotypes (inferred from the heat map). c, Plots of single cell BCR data for subject MM02 using the Chromium 10× single cell BCR platform (Online Methods). Left, single IGLV for subject MM02. Right, frequency of different BCR clonotypes for subject MM02, generated by Chromium’s Loupe V(D)J browser.
Extended Data Fig. 7 CNAs deduced from scRNA-seq data.
a, Comparison of CNVKit output from genomic DNA-targeted sequencing panel at 500× coverage (black dots) and sciCNA (blue) for subject MM10. b, Scatter plot of CNVKit versus sciCNA estimation for each DNA segment for subject MM10. c, Cell-to-cell correlation matrix (Pearson correlation on row-normalized log transformed data) for subjects SMM02 (1,664 cells), MM04 (439 cells) and MM09 (882 cells). d, Heat maps of normalized single cell gene expression (UMI count, log scale) for bone marrow plasma cells from subjects SMM01 (left, 650 cells) and AL03 (right, 548 cells), clustered with the same number of normal bone marrow plasma cells from control donors. Representative variable genes are shown.
Extended Data Fig. 8 scRNA-seq of rare malignant plasma cells in individuals with minimal residual disease.
a, Scatter plots of single cell expression (average UMI count, log scale) for of ELKAP2 (left) and LCP1 (right) genes in subjects pre- and post-treatment. b, Chromium 10× 5′ single cell expression estimates (Online Methods) visualized as tSNE plots for subject MM13 (n = 1) bone marrow cells enriched with CD138 magnetic beads (Online Methods). Shown are expression estimates for IGLV3–IGLV25 (left), PDIA2 (middle) and NSD2 (right). c, Heat map of normalized single cell gene expression (UMI, log scale) for subject MM14 (n = 1) bone marrow plasma cells (1,252 cells). Representative variable genes are shown. d, Scatter plots showing average single cell gene expression (log2 scale) of control donors’ bone marrow plasma cells (x axis) compared with either MM14 cells from C1–4 (top scatter plot; related to heat map in b) or MM14 cells from C7 (bottom scatter plot). Selected gene names are shown.
Extended Data Fig. 9 Circulating plasma cells reflect the bone marrow disease.
a, Heat map showing the distribution of circulating plasma cells from 21 individuals with multiple myeloma as well as controls in 12 clusters (cPC1–cPC12). b, Two-dimensional tSNE view of circulating plasma cells from 21 subjects and controls. Each dot represents a single cell. Patients are color coded. c, Bar plots of average UMI count for the CD52 gene (Online Methods) in common cluster cPC4 (left; blue bars) compared with circulating tumor cells (right; red bars). d, Scatter plots of single cell expression estimates (average UMI count, log scale) for subects’ abnormal bone marrow plasma cells (x axis) to circulating tumor cells in the blood (y axis). Top, subject SMM03, (198 bone marrow plasma cells and 16 circulating tumor cells); bottom, subject MGUS05 (188 bone marrow plasma cells and 211 circulating tumor cells). e, Immunoglobulin light chain variable region distribution within each subject’s bone marrow and blood tumor cells. Shown are 15 individuals for whom we were able to reconstruct BCR data from rare circulating tumor cells. Color represents percentage of cells, ranging from white (0) to blue (0.5). In each box, the top panel represents PB and the bottom panel represents bone marrow. Cell numbers are shown on the left.
Extended Data Fig. 10 CD52 segregates polyclonal circulating plasma cells from malignant ones.
a, Flow cytometry plots of circulating plasma cells from individuals with relapsed myeloma (n = 3) post-treatment. Plots were generated using Infinicyt software (Online Methods). Circulating tumor cells (red) are marked by an aberrant surface phenotype compared with normal plasma cells (green). Principal component analysis quantifies the significance (contribution to principal component 1) of each surface marker to separate between ‘circulating tumor cells’ and ‘normal plasma cells’. Each row represents a different subject.
Supplementary information
Supplementary Tables
Supplementary Tables 1–4, 6
Supplementary Table
Supplementary Table 5
Rights and permissions
About this article
Cite this article
Ledergor, G., Weiner, A., Zada, M. et al. Single cell dissection of plasma cell heterogeneity in symptomatic and asymptomatic myeloma. Nat Med 24, 1867–1876 (2018). https://doi.org/10.1038/s41591-018-0269-2
Received:
Accepted:
Published:
Issue Date:
DOI: https://doi.org/10.1038/s41591-018-0269-2
This article is cited by
-
Single-cell sequencing reveals the mechanisms of multiple myeloma progression: clarity or confusion?
Annals of Hematology (2025)
-
Visualizing scRNA-Seq data at population scale with GloScope
Genome Biology (2024)
-
Bone marrow stromal cells induce chromatin remodeling in multiple myeloma cells leading to transcriptional changes
Nature Communications (2024)
-
1q amplification and PHF19 expressing high-risk cells are associated with relapsed/refractory multiple myeloma
Nature Communications (2024)
-
The TGFβ type I receptor kinase inhibitor vactosertib in combination with pomalidomide in relapsed/refractory multiple myeloma: a phase 1b trial
Nature Communications (2024)
