Whole-genome sequencing of multiple myeloma from diagnosis to plasma cell leukemia reveals genomic initiating events, evolution, and clonal tides

Jan B. Egan, Chang-Xin Shi, Waibhav Tembe, Alexis Christoforides, Ahmet Kurdoglu, Shripad Sinari, Sumit Middha, Yan Asmann, Jessica Schmidt, Esteban Braggio, Jonathan J. Keats, Rafael Fonseca, P. Leif Bergsagel, David W. Craig, John D. Carpten and A. Keith Stewart


The longitudinal evolution of a myeloma genome from diagnosis to plasma cell leukemia has not previously been reported. We used whole-genome sequencing (WGS) on 4 purified tumor samples and patient germline DNA drawn over a 5-year period in a t(4;14) multiple myeloma patient. Tumor samples were acquired at diagnosis, first relapse, second relapse, and end-stage secondary plasma cell leukemia (sPCL). In addition to the t(4;14), all tumor time points also shared 10 common single-nucleotide variants (SNVs) on WGS comprising shared initiating events. Interestingly, we observed genomic sequence variants that waxed and waned with time in progressive tumors, suggesting the presence of multiple independent, yet related, clones at diagnosis that rose and fell in dominance. Five newly acquired SNVs, including truncating mutations of RB1 and ZKSCAN3, were observed only in the final sPCL sample suggesting leukemic transformation events. This longitudinal WGS characterization of the natural history of a high-risk myeloma patient demonstrated tumor heterogeneity at diagnosis with shifting dominance of tumor clones over time and has also identified potential mutations contributing to myelomagenesis as well as transformation from myeloma to overt extramedullary disease such as sPCL.


Multiple myeloma is a treatable mature B-cell malignancy, and progress in improving survival has been rapid in recent years. Despite this, therapy is usually not curative and further insights into the genomic basis for disease and its progression are critical. An early event in myelomagenesis is the dysregulation of Cyclin D.1 This dysregulation results from primary immunoglobulin heavy chain (IgH) translocations with partners MMSET/FGFR3, CCND1, CCND3, cMAF, and MAFB27 or an alternative pathway of chromosomal hyperdiploidy.8 Beyond these initiating events, common mutations that accompany progression include: multiple chromosomal losses9,10 including deletion of RB11113 and TP537, amplification of chromosome 1q,7 mutations in the NF-κB signaling pathway14,15 or in K- and N-RAS.1 PTEN deletions have been reported in 20%-33% of sPCL cases,16,17 which is significantly higher than the incidence in myeloma (6%)17 suggesting PTEN deletion may contribute to leukemic transformation. The presence of t(4;14), t(14;16), chromosome 13 deletion, deletion of 17p, amplification of chromosome 1, and hypodiploidy are associated with high risk/poor prognosis.18

The Multiple Myeloma Research Consortium (MMRC) recently reported the whole genome and the exome sequences of 38 myeloma patients.19 Along with common mutations in known genes such as NRAS, KRAS, and TP53, novel mutations were detected relatively frequently in several genes including, for example, FAM46C and DIS3.19 In addition, a significant clustering of mutations was observed in histone-modifying enzymes.19 Of clinical translational relevance for that study was the first report of BRAF kinase domain mutations in myeloma, with obvious implications for the treatment of disease.19

Although these recent findings from the MMRC19 have provided a detailed characterization of the myeloma genome at a single time point, the genomic basis of the progression of myeloma over time from diagnosis to end-stage disease20,21 remains unclear. Therefore, to further understand the natural history of disease genome progression from diagnosis of myeloma to terminal secondary plasma cell leukemia (sPCL), we conducted a longitudinal study of a single myeloma patient for whom 4 tumor samples were obtained throughout the disease course. Whole-genome sequencing (WGS) conducted on these samples has identified 36 validated novel single-nucleotide variants (SNVs), 27 of which are located in genes that are mutated in the MMRC cohort19 and/or the COSMIC database22 but do not share the same point mutations. PolyPhen-223/SIFT24 predict 18 of these 27 SNVs are damaging. In addition, 10 of the 27 SNVs are shared at all tumor time points while 5 SNVs are unique to the sPCL sample and may be associated with leukemic transformation. Furthermore, mate-pair sequencing (MPS) of the sPCL has identified 79 structural variants.


Detailed methods are available in supplemental Methods (available on the Blood Web site; see the Supplemental Materials link at the top of the online article). Methods in brief are presented here.


Patient samples were acquired with patient consent in accordance with the Declaration of Helsinki with approval from the Mayo Clinic Institutional Review Board. Iliac crest bone marrow (BM) samples were acquired at diagnosis, and first and second relapse. Peripheral blood was obtained at the time of second relapse for use as germline control and again at the time of sPCL diagnosis for circulating tumor cells. BM and peripheral blood samples were treated with ACK lysis buffer to remove red cells and CD138+ cell populations were isolated using anti-CD138 Abs on a Robocept (StemCell Technologies). Tumor cell purity was estimated using a slide-based κ/λ assay. Purified tumor cells were preserved either in RNALater (QIAGEN) or as dry pellets at −80°C. DNA was isolated from preserved samples with the Puregene kit (QIAGEN) following the manufacturer's protocol.

Whole-genome sequencing

Single-end WGS was completed for the constitutional, diagnostic, first relapse, and second relapse samples on the Life Technologies SOLiD platform. Paired-end WGS and mate-pair sequencing (MPS) was conducted for the sPCL sample on the Illumina HiSeq 2000 platform. All sequencing reads were aligned to NCBI36/hg18; SNVs were identified and pairwise comparison made between the tumor samples and the germline sample for the identification of somatic variants. After visual confirmation in the Integrative Genomics Viewer,25 the SNV were validated by capillary sequencing in all samples. The MMRC cohort19 of 38 myeloma genomes and the COSMIC database22 of multiple tumor types were then queried to determine, within larger tumor populations, the frequency of all somatic mutations in the genes containing validated somatic SNVs in this case study. Comparative pathway analysis was conducted to identify potentially shared pathways among the SNV unique to specific time points. Mate-pair analysis was conducted for the identification of structural rearrangements.

Microarrays and FISH

DNA samples from all tumor time points were fragmented, labeled, and hybridized to the 244A Human Genome CGH Microarray (Agilent) according to the manufacturer's recommendation. Aberration calling was performed in Genomic WorkBench Version 6.5 (Agilent Technologies) using the ADM-2 algorithm with a 5.5 threshold and 0.2 log2 and 3 probe filters.

RNA samples from diagnosis and sPCL were submitted to the Mayo Clinic genomics core for gene expression studies. Samples were labeled and hybridized to the HG-U133Plus 2.0 genechip (Affymetrix). Expression estimates were extracted using the MAS5.0 algorithm implemented in the Expression Console using default parameters.

The expression level of identified candidate genes was analyzed across public databases including: the MMRC database (, the GlaxoSmithKline Cancer Cell Line Genomic Profiling Dataset ( 041), and our Mayo Clinic myeloma research laboratory database (J.J.K., E.B., unpublished data, November 14, 2011). FISH assays were conducted at the time of diagnosis in a clinical laboratory in compliance with federal regulatory standards. The microarray datasets associated with this manuscript have been deposited in GEO as a super series under accession GSE36825; the individual platforms are available under GSE36822, GSE36823, and GSE36824.26


Case history

A 67-year-old woman was diagnosed with multiple myeloma and presented with BM containing 25% plasma cells harboring a t(4;14)(p16;q32) and chromosome 13 deletion detected by FISH. The patient was initially treated with lenalidomide and low-dose dexamethasone with very good response, but progression of disease was first observed at 22 months. BM was collected for the first relapse sample at 23 months. The patient then received treatment with carfilzomib on a clinical trial for the next 4 months with good response. Treatment was then discontinued to pursue stem cell mobilization that was ultimately unsuccessful. At 33 months, disease progression was again observed and the patient was treated with bortezomib in combination with an experimental agent, SGN-40, for 2 months; however, therapy was discontinued because of toxicity. Collection of the BM for the second relapse sample, as well as peripheral blood for constitutional genome analysis, was conducted at 39 months. Treatment with melphalan, prednisone, and bortezomib was then begun and continued until month 43 and was again discontinued because of toxicity. For the next 4 months, the patient received various combination chemotherapies with limited success. At month 53, the myeloma progressed to sPCL and treatment was discontinued; the patient was ultimately referred to hospice. Peripheral blood was collected for the sPCL sample at month 54.

Array comparative genomic hybridization (aCGH) of the diagnostic and the first relapse samples in this patient contained numerous copy number abnormalities which in total encompassed 3968 genes, of which 1235 were expressed at diagnosis.27 We wished to interrogate the specific genes that may be involved in myelomagenesis, disease progression, and drug sensitivity or resistance. We therefore decided to conduct WGS to further investigate the individual genetic events associated with the natural history of progressive myeloma.

Whole-genome sequencing

WGS was conducted for germline DNA and at 4 time points during the tumor evolution including: diagnosis, first relapse (23 months), second relapse (39 months), and sPCL (54 months). For the first 3 tumor time points and germline samples, an average of 2.45 × 109 single-end reads/sample were sequenced with an average of 1.76 × 109 reads/sample (71.8%) aligning to the hg18 reference (Table 1). The sPCL sample had 1.04 × 109 paired-end reads/sample with 0.97 × 109 reads/sample (93.1%) aligning to the hg18 reference. The average coverage of the aligned sequence was 30 times before duplicate removal and 20 times after removal of duplicates. More than 2 million constitutional SNVs were identified in the germline sample, and tumor samples had > 14 000 additional somatic SNVs with the sPCL sample containing a much larger number of variants than the other samples. The majority of the SNVs were located within introns and intergenic space. From the somatic SNV located in exons, 124 candidate, nonsynonymous SNVs were identified that were tumor specific, 36 of which validated with capillary sequencing (Table 2, supplemental Table 1). Of the 36 SNVs that validated, 27 were located in genes containing somatic mutations in the MMRC 38 myeloma genomes19 and/or the COSMIC database22 of multiple tumor types. While the genes containing the 27 SNVs were also somatically mutated in these other tumor populations, no somatic SNVs were shared between these populations and this work. Eighteen (68%) of these 27 SNVs were predicted as potentially damaging and 3 (11%) were predicted as truncating amino acid changes. Of those 27 genes, 9 (33%) were located in genes identified as mutated in both the MMRC cohort and COSMIC database including: AFF1, C12orf42, CSMD3, LRRC4C, PCDH7, PTPRD, PPFIBP1, RB1, and ZKSCAN3 (Table 2, Figure 1). As DNA was limiting, mate-pair analysis could only be performed on the sPCL sample. More than 2000 small insertions and deletions (indels) were detected, as well as 79 structural variants (supplemental Table 2), in this plasma cell leukemia illustrating a profoundly disturbed genome.

Table 1

Summary of whole-genome sequencing metrics

Table 2

Summary of somatic SNV validated by capillary sequencing

Figure 1

Frequency of somatic mutations in genes of interest within other tumor cohorts. Genes containing SNV identified in this case study patient were compared against the MMRC cohort of 38 myeloma genomes and the COSMIC database of other tumor types to determine the frequency of all somatic mutations in these genes within larger populations. This graph illustrates the frequency (percentage) of all somatic mutations in the MMRC cohort and in other tumors present in the COSMIC database. Frequencies with < 0.5% were counted as absent and were not included in the graph.

Identification of SNV-driving myelomagenesis

SNVs shared at all tumor time points are likely to be tumor-initiating events. Fifteen of the validated SNVs were shared at all tumor time points, with 10 of these 15 presenting in somatically mutated genes in the MMRC cohort19 and/or the COSMIC database,22 including: AFF1, ATXN1, CNGA3, COL2A1, CSMD3, KRT9, LRRC4C, MAGI1, MYPN, and RNF145 (Table 2). Of particular interest is AFF1, which is widely expressed in myeloma, and was also identified as mutated in the MMRC cohort.19 AFF1 is a fusion partner with MLL in both pediatric and adult acute leukemias28 that has been implicated in the alteration of histone methylation signatures in mice.29 In this patient, the AFF1 gene contains a missense mutation (R163H) that Mutation Taster30 predicts may result in the loss of phosphoserine sites downstream of a putative altered splice site.

Identification of SNV-driving leukemic transformation

SNVs unique to the sPCL sample may contribute to leukemic transformation from myeloma to sPCL. Seven SNVs were unique to sPCL, of which 5 were present in genes somatically mutated in the MMRC cohort19 and/or the COSMIC database22 and include: RB1, TNN, TUBB8, ZKSCAN3, and ZNF521 (Table 2). Two genes of particular interest are RB1 and ZKSCAN3, which were somatically mutated in both the MMRC cohort19 and the COSMIC database.22 A deletion of RB1 was detected in this patient at diagnosis, thus the SNV in the sPCL sample is hemizygous as it occurred in the only remaining copy of RB1. The mutations in RB1 (E287*) and ZKSCAN3 (E28*) are predicted to result in truncated proteins.23,30 Although RB1 is expressed at both diagnosis and sPCL, the RB1 truncation at the sPCL time point is predicted by Mutation Taster30 to result in a protein lacking critical binding pockets and interaction domains. The well-conserved KRAB and SCAN domains within ZKSCAN3 are predicted by Mutation Taster to be eliminated with this truncation. Furthermore, no expression of ZKSCAN3 was detected in this patient.

Clonal evolution over time

Although 10 common SNVs were shared at all tumor time points, we also observed variants which were only detectable at alternating time points—that is, diagnosis and second relapse, or first relapse and sPCL—suggesting the waxing and waning of different clones with time and treatment (Figure 2, supplemental Table 1). At diagnosis, one SNV, PDE4DIP, was unique, not found at later time points. Biologically, this is difficult to explain and may be artifact although the mutation was identified both in libraries prepared from genomic DNA and from capillary sequencing which used whole-genome amplified DNA. This SNV is different from the somatic mutations present in other tumor types within the COSMIC database22 including a single chronic lymphocytic leukemia case. Six SNV were shared by only the diagnostic and the second relapse samples, of which 5 SNV were located in genes with somatic mutations in the MMRC cohort19 and/or the COSMIC database,22 including: C12orf42, DOK5, PARD3B, PPFIBP1, and ZNF557. The first relapse and sPCL samples did not share any unique SNV beyond the 10 observed at all tumor time points. Array CGH and FISH data, described in detail in a companion manuscript by Keats et al in this issue, indicate, however, that these time points share a common progenitor that is not detectable with these validated, coding SNV.27 The first relapse sample had 7 unique SNVs not found at any other time point, of which 6 were present in genes somatically mutated in the MMRC cohort and/or the COSMIC database, including: ATXN1, CACNA1S, DSC1, PCDH7, PTPRD, and TLR9. The sPCL had 7 SNVs unique only to the sPCL sample, of which 5 were present in genes containing somatic mutations in the MMRC cohort and/or the COSMIC database, including: RB1, TNN, TUBB8, ZKSCAN3, and ZNF521. Comparative pathway analysis revealed no shared pathways between the genes containing SNV unique to specific time points. While the evolutionary divergence of the diagnostic and second relapse samples is relatively small with only one unique SNV, the evolutionary divergence of the first relapse and sPCL samples is much larger with the samples acquiring 6 and 5 mutations, respectively.

Figure 2

Clonal divergence of validated variant alleles. Graph illustrating the presence of shared and unique SNV at each tumor time point. There are 15 variants common to all time points and shared by a common ancestor. Six variants are common to only the diagnostic and second relapse, while no variants are common to the first relapse and sPCL. The greatest divergence is observed between the first relapse and sPCL, with 7 unique variants detected in each sample.

Identification of structural variants

The sPCL mate-pair analysis revealed the presence of 79 structural variants including the diagnostic t(4;14), 3 rearrangements involving chromosomes 10 and 12, as well as complex structural rearrangement in regions on chromosomes 11 and 12 (Figure 3, supplemental Table 2). All of the translocation breakpoints in chromosome 12 occur in regions with multiple breakpoints indicating complex structural rearrangement. This suggests intrachromosomal rearrangement likely occurred within chromosome 12 before a translocation event with chromosome 10. Interpretation of these complex events is beyond the scope of this study and is, we suggest, not pertinent to the main message. Thirteen inversions were detected, 10 of which occurred in gene-rich chromosomal regions including MLLT11 and NOTCH2NL which are expressed in myeloma. MLLT11 is a translocation partner in acute leukemia31 and overexpression has been associated with poor outcome in pediatric acute myeloid leukemia and adult myelodysplastic syndrome.32 A translocation was identified between NOTCH2NL and LCLAT1 in the MMRC cohort,19 suggesting this NOTCH signaling pathway protein may be prone to structural rearrangement. One-half of the 10 biallelic deletions observed occurred in noncoding regions while the remaining deletions included deletion of one or more genes including: BIRC2, BIRC3, and TP53.

Figure 3

Circos plot depicting the summary of structural variation in the sPCL. The center, line plots indicate the presence of discordant read pairs. The middle ring contains the array CGH plot, and the outermost ring indicates mutated or deleted genes.


The comprehensive study of myeloma genomes in the MMRC cohort19 has confirmed the heterogeneity of myeloma tumors, and identified potential pathways and mutations of interest for further study. However, the study was unable to address genomic changes over time and specifically changes that contribute to the leukemic transformation leading to sPCL. This work reports several potential drivers of myeloma as well as initiators of leukemic transformation from myeloma to sPCL.

Potential initiating mutations would be expected to present at all time points. One of these genes, AFF1, expressed in myeloma and also mutated in the MMRC cohort, was identified in our study to share the same SNV at all tumor time points. AFF1 (also referred to as AF4) is involved in RNA polymerase II–mediated transcription. Specifically, it associates with positive transcription elongation factor b (P-TEFb), and ENL as part of the AEP complex (AF4/ENL/P-TEFb)33 that regulates transcription.34,35 The AEP complex colocalizes with MLL at the HOXA9 promoter.33 In myeloma cell lines, HOXA9 expression was inversely related to H3K24me3 levels, and cells lines depleted of HOXA9 demonstrated a growth disadvantage.19 Therefore, a damaging mutation in AFF1, which is predicted to disrupt phosphoserine sites downstream of the mutation, could result in dysregulation of HOXA9 leading to dysregulated cell growth and contributing to myelomagenesis.

Variants that do not present until the final sPCL time point are potential candidates for leukemic transformation. One of these, ZKSCAN3, has been demonstrated to modulate CCND2 in myeloma cells with increased levels of CCND2 observed with ZKSCAN3 overexpression in U266 myeloma cells.36 In contrast, the opposite is expected in the sPCL as the truncating SNV in ZKSCAN3 observed in the sPCL is predicted to result in a loss-of-function protein. Thus, in the absence of functional ZKSCAN3 protein, we would expect decreased levels of CCND2 with subsequent increased levels of other cyclin D proteins potentially leading to dysregulation of the G1/S checkpoint.

Another gene with a truncating variant that does not present until sPCL is RB1, a critical regulator of the cell cycle. Monoallelic deletion of RB1 has been implicated in myelomagenesis1113 and was present in the chromosome 13 deletion detected at diagnosis in this patient. The SNV in the remaining RB1 allele of this patient is predicted to result in a loss-of-function protein, as all critical binding and interaction domains are lost. Consequently, the functional loss of RB1 in the remaining allele is significant as the resulting protein is expected to be unable to regulate the critical G1/S cell-cycle checkpoint. Furthermore, a biallelic deletion in TP53, a key tumor suppressor that regulates cell-cycle arrest, was also identified only in the sPCL sample resulting in the absence of 2 key regulators of the cell cycle, RB1 and TP53. Taken together, truncating mutations in both ZKSCAN3 and RB1 together with deletion of TP53 could result in catastrophic dysregulation of cell-cycle checkpoints. This combined with a mutation in AFF1 could confer the growth advantage necessary for leukemic transformation from myeloma to sPCL.

The presence of variants detectable only at alternating time points suggests that all clonal progenitors were present at diagnosis, but selection pressures from treatment and clonal evolution caused dominance of these clones to rise and fall over time. Our work suggests the diagnostic and the second relapse clones demonstrate little divergence over time as only one SNV differentiates the 2 clones. In contrast, the first relapse and the sPCL clones have 6 and 5 unique SNV, respectively, a total of 11 SNVs distinguishing the 2 clones, suggesting greater evolutionary divergence with time and aggressiveness of disease. Taken together, this suggests common initiating events represented by the 10 shared SNVs, but divergence in clonal evolution with at least 2 dominant clones at diagnosis that then evolve in clonal tides detected only at alternating time points. This suggests a selective pressure of therapeutic drugs. Of interest, correlation of genomic change with therapy suggests that the dominant diagnostic clone was sensitive to lenalidomide but as that clone declined, a new resistant clone emerged resulting in the first relapse. Furthermore, the first relapse clone was responsive to carfilzomib but not the original diagnostic clone that re-emerged. In addition, the presence of unique SNVs at the final sPCL time point that were undetectable at earlier time points suggests tumor evolution that involved the acquisition of novel mutations that enabled leukemic transformation to occur. In addition, myeloma can present as multiple, discrete foci that are not adequately represented in serial iliac crest BM samplings; thus, the clonal heterogeneity of myeloma may in fact be even greater than that reported here.

Although the findings reported here are limited to a single case study, they are the first to study genomic evolution of a myeloma over time, and demonstrate tumor heterogeneity and clonal dynamics of myeloma. Furthermore, these findings suggest possible driver mutations that may contribute to myelomagenesis and to leukemic transformation. To further understand the natural history and progression of myeloma to sPCL, additional studies in larger cohorts of serially collected patient samples are warranted.


Contribution: J.B.E. conducted sequencing, analysis, and manuscript preparation; C.-X.S., J.S., E.B., and J.J.K. conducted analysis and validation; W.T., A.C., A.K., S.S., S.M., and Y.A. conducted sequencing alignment and postalignment analysis; R.F. and P.L.B. contributed to manuscript preparation; D.W.C. contributed to analysis and manuscript preparation; and J.D.C. and A.K.S. designed the study, and contributed to analysis and manuscript preparation.

Conflict-of-interest disclosure: R.F. has received a patent for the prognostication of MM based on genetic categorization of the disease; has received consulting fees from Medtronic, Otsuka, Celgene, Genzyme, BMS, Lilly, Millennium, and Amgen; and has also sponsored research from Cylene and Onyx. The remaining authors declare no competing financial interests.

Correspondence: A. Keith Stewart, Mayo Clinic Collaborative Research Bldg, 13400 E Shea Blvd, Scottsdale, AZ 85259-5494; e-mail: stewart.keith{at}


The authors thank Greg Ahmann for his assistance with processing samples for this project. They thank GlaxoSmithKline for the use of the Cancer Cell Line Genomic Profiling Dataset. They are indebted to the Translational Genomics Research Institute High-performance Bio-computing Division for making available the supercomputing resources funded by National Institutes of Health (NIH) grants 1S10RR25056-01 and 1S10RR023390-01.

This work was supported by NIH grant CA133115-01 and the Mayo Clinic Cancer Center.


  • There is an Inside Blood commentary on this article in this issue.

  • The online version of this article contains a data supplement.

  • The publication costs of this article were defrayed in part by page charge payment. Therefore, and solely to indicate this fact, this article is hereby marked “advertisement” in accordance with 18 USC section 1734.

  • Submitted January 21, 2012.
  • Accepted April 17, 2012.


View Abstract