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

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.¹  This dysregulation results from primary immunoglobulin heavy chain (IgH) translocations with partners MMSET/FGFR3, CCND1, CCND3, cMAF, and MAFB^2-7  or an alternative pathway of chromosomal hyperdiploidy.⁸  Beyond these initiating events, common mutations that accompany progression include: multiple chromosomal losses^9,10  including deletion of RB1^11-13  and TP53⁷ , amplification of chromosome 1q,⁷  mutations in the NF-κB signaling pathway^14,15  or in K- and N-RAS.¹ PTEN deletions have been reported in 20%-33% of sPCL cases,^16,17  which is significantly higher than the incidence in myeloma (6%)¹⁷  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.¹⁸ 

The Multiple Myeloma Research Consortium (MMRC) recently reported the whole genome and the exome sequences of 38 myeloma patients.¹⁹  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.¹⁹  In addition, a significant clustering of mutations was observed in histone-modifying enzymes.¹⁹  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.¹⁹ 

Although these recent findings from the MMRC¹⁹  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 disease^20,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 cohort¹⁹  and/or the COSMIC database²²  but do not share the same point mutations. PolyPhen-2²³ /SIFT²⁴  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.

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,²⁵  the SNV were validated by capillary sequencing in all samples. The MMRC cohort¹⁹  of 38 myeloma genomes and the COSMIC database²²  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.

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 (http://www.broadinstitute.org/mmgp/home), the GlaxoSmithKline Cancer Cell Line Genomic Profiling Dataset (https://array.nci.nih.gov/caarray/project/woost-00 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.²⁶ 

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.²⁷  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.

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 × 10⁹ single-end reads/sample were sequenced with an average of 1.76 × 10⁹ reads/sample (71.8%) aligning to the hg18 reference (Table 1). The sPCL sample had 1.04 × 10⁹ paired-end reads/sample with 0.97 × 10⁹ 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 genomes¹⁹  and/or the COSMIC database²²  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.

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 cohort¹⁹  and/or the COSMIC database,²²  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.¹⁹  AFF1 is a fusion partner with MLL in both pediatric and adult acute leukemias²⁸  that has been implicated in the alteration of histone methylation signatures in mice.²⁹  In this patient, the AFF1 gene contains a missense mutation (R163H) that Mutation Taster³⁰  predicts may result in the loss of phosphoserine sites downstream of a putative altered splice site.

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 cohort¹⁹  and/or the COSMIC database²²  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 cohort¹⁹  and the COSMIC database.²²  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 Taster³⁰  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.

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 database²²  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 cohort¹⁹  and/or the COSMIC database,²²  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.²⁷  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.

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 leukemia³¹  and overexpression has been associated with poor outcome in pediatric acute myeloid leukemia and adult myelodysplastic syndrome.³²  A translocation was identified between NOTCH2NL and LCLAT1 in the MMRC cohort,¹⁹  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.

The comprehensive study of myeloma genomes in the MMRC cohort¹⁹  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)³³  that regulates transcription.^34,35  The AEP complex colocalizes with MLL at the HOXA9 promoter.³³  In myeloma cell lines, HOXA9 expression was inversely related to H3K24me3 levels, and cells lines depleted of HOXA9 demonstrated a growth disadvantage.¹⁹  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.³⁶  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 G₁/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 myelomagenesis^11-13  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 G₁/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.

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.

National Institutes of Health

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@mayo.edu.