Genomic landscape and chronological reconstruction of driver events in multiple myeloma

The multiple myeloma (MM) genome is heterogeneous and evolves through preclinical and post-diagnosis phases. Here we report a catalog and hierarchy of driver lesions using sequences from 67 MM genomes serially collected from 30 patients together with public exome datasets. Bayesian clustering defines at least 7 genomic subgroups with distinct sets of co-operating events. Focusing on whole genome sequencing data, complex structural events emerge as major drivers, including chromothripsis and a novel replication-based mechanism of templated insertions, which typically occur early. Hyperdiploidy also occurs early, with individual trisomies often acquired in different chronological windows during evolution, and with a preferred order of acquisition. Conversely, positively selected point mutations, whole genome duplication and chromoplexy events occur in later disease phases. Thus, initiating driver events, drawn from a limited repertoire of structural and numerical chromosomal changes, shape preferred trajectories of evolution that are biologically relevant but heterogeneous across patients.

The genome of multiple myeloma (MM) is complex and heterogeneous, with a high frequency of structural variants (SVs) and copy-number abnormalities (CNAs) [1][2][3] . Translocations between the immunoglobulin heavy chain (IGH) locus and recurrent oncogenes are found in ~40% of patients. Cases without IGH translocations often have a distinctive pattern of hyperdiploidy affecting oddnumbered chromosomes, where the underlying target genes remain mysterious.
These SVs and recurrent CNAs are considered early drivers, being detectable also in pre-malignant stages of the disease 1,2,4 . Cancer genes are also frequently altered by driver point mutations, with MAPK and NF-κB signaling as major targets 5-8 .
Many blood cancers develop along preferred evolutionary trajectories. Early driver events, drawn from a restricted set of possible events, differ in which subsequent cancer genes confer clonal advantage, leading to considerable substructures of co-operativity and mutual exclusivity among cancer genes. These subtypes vary in chemosensitivity and survival, suggesting that although patients share a common histological and clinical phenotype, the underlying biology is distinctly heterogeneous. Preliminary studies have suggested that these patterns exist in MM as well 5,7,9-12 , but have not yet been systematically defined in large cohorts with broad sequencing coverage.
We performed whole genome sequencing (WGS) of 67 tumor samples collected at different time points from 30 MM patients, together with matched germline controls (Supplementary Fig. 1, Supplementary Table 1, Methods). We also included in our analyses published whole exome data from 804 patients 13,14 . To discover significant cancer genes, we analyzed the ratio of non-synonymous to synonymous mutations, correcting for mutational spectrum and covariates of mutation density across the genome using a published algorithm 15,16 . Overall, 55 genes were significantly mutated with a false discovery rate of 1% (Fig. 1a,   Supplementary Table 2). A significant fraction of these driver mutations was detected at subclonal level, suggesting a major role in late phases of cancer development (Supplementary Fig. 2a). Beyond well-known myeloma genes such as KRAS,NRAS, Figure 2f). MAX, a DNA-binding partner of MYC, showed an interesting pattern of start-lost mutations, nonsense and splice site mutations, together with hotspot missense mutations at residues Arg35, Arg36 and Arg60, known to abrogate DNA binding 7 (Supplementary Fig. 2g). Genes with rather more mysterious function were also significant: the zinc finger ZNF292, recently described as mutated in chronic lymphocytic leukemia and diffuse large B-cell lymphoma 19,20 , showed an excess of protein-truncating variants (Supplementary Fig. 2h); the uncharacterized TBC1D29 gene showed two hotspots of missense mutations 21 ( Supplementary Fig. 2i).
In pairwise comparisons, these cancer genes showed distinct patterns of comutation and mutual exclusivity (Supplementary Fig. 2j). To define the logic rules underpinning the conditional dependencies of driver events, we employed Bayesian networks and the hierarchical Dirichlet process (Fig. 1b-c, Supplementary Fig. 3).
This confirmed that the strongest determinants of genomic substructure in myeloma are IGH translocations and recurrent CNAs. Some co-operating genetic lesions were non-randomly distributed across these main groups: a significant fraction of patients without IGH translocations were generally enriched for 1q gain and deletions on 1p13, 1p32, 13q, TRAF3 and CYLD (Cluster 1). RAS signaling mutations, especially NRAS and KRAS, were associated with hyperdiploidy and MYC translocations (Cluster 2). A significant fraction (33%) of patients with t(11;14) were characterized by low genomic complexity and high prevalence of IRF4 and CCDN1 mutations (Cluster 3). Patients harboring TP53 bi-allelic inactivation were clustered in an independent sub group (Cluster 4). A significant fraction of MMSET (51%) and CCND1 (19%) translocated patients were characterized by multiple cytogenetic aberrations, with low prevalence of CYLD/FAM46C deletions and TRAF3 deletion respectively (Cluster 5). A second fraction of patients with MMSET translocation (46%) were grouped with deletion of 13q14, gain of 1q21, DIS3 and FGFR3 mutations (Cluster 6); finally, patients with either MAF/MAFB translocations or no IGH translocations was characterized by a high driver mutation rate (Cluster 7).
Thus, there is evidence for at least 7 distinct genomic subtypes of myeloma, each with distinct combinations of driver mutations and recurrent SVs (Figure 1c).
In 6/30 (20%) patients, we found a novel complex pattern characterized by cycles of templated insertions. Here, several low-amplitude copy number gains on different chromosomes were linked together through SVs demarcating the region of duplication (Fig. 2b,d; Supplementary Fig. 5d-h). The most plausible explanation for this pattern is that the templates are strung together into a single chain, hosted within one of the chromosomes. We have observed such events in some solid cancers 22 , where the copy number gains suggest a mutational process that is replication-based, rather than the break-and-ligate processes generating chromothripsis and chromoplexy 23 .
This complex landscape of SVs in myeloma often involved known driver genes: MYC (14/30 cases; 46%), CCND1 (7/30; 23%) and MMSET (3/30; 10%) were common targets (Supplementary Fig. 4b) of these non-canonical events. The juxtaposition of CCND1 to the IGH locus was caused by either unbalanced translocations or insertional events in 5/8 patients (Supplementary Fig. 6). Similarly, MYC translocations showed unanticipated complexity, with four cases of templated insertions involving MYC or its regulatory regions (Supplementary Fig. 7). Such events are the structural basis of oncogene amplification observed by FISH in many cases of t(11;14) and t(8;14) 28,29 . Interestingly, many of the MYC SVs involved the immunoglobulin light chain loci, IGK or IGL, rather than the heavy chain IGH locus, and were seen in patients with hyperdiploidy (Supplementary Fig. 4b). Although sometimes occurring late, these events were under strong selective pressure: we identified a striking case of convergent evolution where a subclone bearing an IGL:MYC translocation was lost and one bearing an IGH:MYC was acquired at relapse (Supplementary Fig. 8). SVs also led to loss of tumor suppressor genes such as BIRC2/3, CDKN2A/B, CDKN2C, TRAF3 30 and FAM46C, either within focal deletions or more complex events. These data confirm that SVs, accessing both simple and complex mechanisms of genome rearrangement, is a major force shaping the myeloma genome.
In our cohort of serial WGS samples we could evaluate the temporal evolution of the disease, integrating SVs, CNAs and point mutations. This included not only defining clonal and subclonal mutations and CNAs, but also timing the relative order of acquisition of clonal chromosomal gains, which we estimated with "molecular time" analysis based on the fraction of duplicated to single-copy point mutations (Methods) 31 .
Hyperdiploidy was found in 18/30 patients. Applying our "molecular time" analysis, we found that chromosomal gains were not necessarily acquired simultaneously (Fig. 3a). In 11/18 cases, we found evidence of several independent gains over time, although the other 7 patients acquired all chromosome gains in a much narrower time window (Fig. 3b). As might be expected if gains occur as independent events, subclonal evolution within the myeloma cells between diagnosis and relapse could lead to considerable diversification of chromosome complements.
In fact, we observed significant karyotype changes within the same patient over time, including loss of some trisomies in hyperdiploid patients (Supplementary Fig. 9ab) 32 . At the extreme end of this cytogenetic dynamism, we found 4 patients acquiring a whole genome duplication at relapse, suggesting its potential role in relapsed/refractory stages (Supplementary Fig. 9c-g).
By pairwise comparison of the relative timings of copy number alterations we reconstructed the preferred chronological order of CNAs acquisition (Methods).
Gains of chromosomes 19, 11, 9 and 1q were amongst the earliest in our series ( Fig. 3c) and recurrent chromosome losses were generally acquired later than trisomies, consistent with the proposition that hyperdiploidy is an early driver event in MM.
Translocations involving CCND1 and MMSET were always fully clonal, similarly confirming their early driver role in MM pathogenesis. Most chromothripsis events were clonal and conserved during evolution (17/22; 77%), suggesting they occurred early in MM pathogenesis. However, a small fraction of patients showed some evidence of subclonal or late chromothripsis (5/22; 23%), implying a potential involvement in drug resistance and late cancer progression (Supplementary Fig.   9h).
We integrated all extracted chronological data on SVs, hyperdiploidy and point mutations to generate phylogenetic trees for each sample (Methods and Supplementary Fig. 10) 5,33 . The methodology is worked through for one illustrative patient carrying i) several chromosome gains, ii) 3 separate chromothripsis events and iii) a whole genome duplication ( Fig. 3d-f). One chromothripsis involved chromosomes 8 and 15, duplicating the long arm of chromosome 15. Because few mutations were present on chromosome 15 at the time it duplicated, this must have occurred early in molecular time (Fig. 3e). Gain of chromosome 3 and copy-neutral loss-of-heterozygosity of small arm of chromosome 1 (chr1p) occurred not long after, and were followed by a second chromosomal crisis involving chromosomes 3, 5, 13 and 22. This chromothripsis must have occurred in one of the two duplicated alleles of chromosome 3 (and therefore after the acquisition of a chromosome 3 trisomy) because losses within the chromothripsis region had copy number 2 and SNPs were heterozygous. Within the same time window, a separate chromothripsis event occurred on chr1p after copy-neutral loss-of-heterozygosity. Finally, this patient underwent whole genome duplication after two therapy lines (Supplementary Fig.   9f).
Applying this approach to all patients, we observed that the trunks of the phylogenetic trees of 29/30 (97%) patients were characterized by few genomic events generally acquired during different time windows of the MM life history before emergence of the most recent common ancestor (Fig. 4). Overall, chromothripsis, cycles of templated insertions, chromosomal gains and other SVs accounted for most of the earliest events, emerging as key early drivers of the disease and paving the way for subsequent driver mutations that would confer further selective advantage to the clone.
Taken together, these data suggest that MM development follows preferred evolutionary trajectories, with stuttering accumulation of driver events in keeping with its insidiously progressive but unpredictable clinical course. Critical early events include immunoglobulin translocation with MMSET and CCND1; hyperdiploidy and focal complex structural variation processes hitting key myeloma genes. These early driver mutations shape the subsequent evolution of myeloma, each with preferred sets of co-operating cancer genes.