Whole-genome analysis identifies novel drivers and high-risk double-hit events in relapsed/refractory myeloma

Despite tremendous advance in therapeutics, multiple myeloma (MM) remains largely incurable, characterized by periods of remission and relapse, with most patients eventually dying from the disease. Currently, combinations of immunomodulatory agents (IMiDs; eg, lenalidomide and pomalidomide), proteasome inhibitors, steroids, anti-CD38 antibodies, and more recently B-cell maturation antigen (BCMA)-targeted therapies are standards of care in newly diagnosed MM (ndMM) and/or relapsed/refractory MM (rrMM) settings. The complex array of combinatorial therapeutic options creates selection pressure that dictates the evolution of genomic resistance in patients through the acquisition of novel and/or expansion of existing genetic drivers. Owing to its higher resolution in accurately detecting genomic events,¹ whole-genome sequencing (WGS) analysis of a large cohort of rrMM genomes may uncover evidence regarding the relationship between therapeutic exposure and the selective pressure driving the evolution of therapy-exposed and resistant subclones.

Previously, our ndMM analysis identified an array of genomic aberrations as well as oncogenic dependencies indicating that the genomic landscape in ndMM is characterized by primary events followed by acquisition of additional genetic events.² The WGS landscape of rrMM has only previously been described in small data sets with limited investigation into the effects of specific therapeutic pressures on tumor genetics.^3-5 We recently demonstrated a highly significant enrichment of genomic CRBN aberrations, and loss of the 2q37 region, which houses COP9 signalosome members following acquisition of lenalidomide or pomalidomide resistance (designated as LENR and POMR, respectively).^6,7 Gene expression signatures have also been shown to transition to those that indicate higher risk as disease progresses.⁸ 

Here, we present a genome-wide unsupervised analysis of somatic mutations, copy number aberrations (CNAs), and structural variants (SVs) along with mutational signatures in a set of 418 rrMM WGS samples from 386 patients. Drug-specific resistance-associated tumor genomic variants have previously been described (eg, cereblon aberration associated with IMiD resistance,⁷ proteasome subunit PSMB5 aberrations associated with protease inhibitor resistance,⁹ and BCMA locus loss after BCMA-directed therapy¹⁰), as well as immune cell abnormalities associated with resistance to immunotherapies. However, deciphering how drug-specific genomic evolution enhance the acquisition or enrichment of “high-risk” genetic features during sequential therapeutic exposures has not been explored systematically. By utilizing the latest stage-available samples from rrMM patients, we compared them with ndMM WGS samples from 198 unrelated patients¹¹ to enable an in-depth population-level case-control analysis of how the MM genomes evolve in response to therapeutic exposures of lenalidomide- and pomalidomide-based regimens (Figure 1). We identified mutually exclusive rare driver genes (EZH2, DUOX2, and PIGO), enrichment of existing high-risk features such as 1q Gain and TP53/17p events (suggesting a transition to high-risk biology in rrMM) and enrichment of a novel subclonal DNA mismatch repair signature in high mutation burden tumors.

This study is based on an analysis of a set of unrelated ndMM and rrMM cases with WGS data (supplemental Methods and Table 1 for details; available on the Blood website). For our analysis, patients with rrMM were classified based on refractory status LENR or POMR. Because most patients were from pomalidomide trials, which required either exposure or refractoriness to lenalidomide, the LENR cohort (n = 269/386 patients, 70%) is larger than the POMR cohort (n = 55/386 patients, 14%) (Table 1, supplemental Table 1). The ndMM data set comprised WGS data from patients in IFM2009 (n = 198).¹¹ 

CD138⁺ plasma cells and matched germline controls (peripheral blood mononuclear cells) were stored in RLT buffer (Qiagen) and DNA was extracted using the Qiagen AllPrep DNA Mini Kit. WGS was performed on tumor/normal pairs at an average of 60×/30× depth. All WGS raw data processing and analyses are described in detail in the supplemental Methods. Briefly, FASTQ files were aligned to the human genome assembly hg38 using BWA-mem. Duplicate reads were removed using Picard MarkDuplicates and base recalibration of alignments was performed using Base Quality Score Recalibration according to Genome Analysis Toolkit Best Practices. Somatic variants were called using MuTect2 and annotated by ANNOVAR.¹² Genome-wide copy number aberration analysis was undertaken in the entire data set by using Battenberg.¹³ In addition to calling clonal and subclonal allele-specific CNAs, it was also used to estimate purity and average ploidy of each tumor. To estimate clonality, we calculated the cancer cell fraction (CCF) of each variant by adjusting the variant allele frequency (VAF) for copy number status at the given locus, multiplicity of the variant, and tumor purity.¹⁴ To assign somatic mutations as clonal or subclonal, the observed VAF of each mutation was modeled using a binomial distribution to obtain its 95% confidence interval (CI) and CCF computed for the 95% boundaries. Variants with an upper CCF boundary above and below 1 were considered to be clonal and subclonal, respectively. To identify candidate drivers of the relapse/refractory stage of MM, we used cDriver.¹⁵ Assignment of identified drivers into tumor suppressor gene (TSG) and oncogene (ONC) categories was undertaken using the 20/20 principle.¹⁶ We analyzed all drivers for gene essentiality as reported in DepMap.¹⁷ De novo extraction of single-base substitution (SBS) signatures was independently implemented using SigProfilerExtractor¹⁸ and signeR,¹⁹ and assigned to COSMIC signatures¹⁸ by computing cosine similarities. To fit the mutational spectrum of each tumor with the identified COSMIC signatures, we used mmsig.²⁰ This was implemented to not only minimize mutational signature bleeding across samples²¹ but to also accurately estimate the relative contribution of each signature. SVs from whole-genome data were detected using MANTA.²² All statistical analyses were implemented in R.²³ 

The demographics and patient characteristic of total, LENR, and POMR cohorts are described in Table 1. Most LENR/POMR patients were also refractory to bortezomib (83% in LENR; 70% in POMR) and nearly all patients had received dexamethasone; therefore, whereas they are classified by their refractory status to LEN or POM in this analysis, these patients were also exposed or refractory to proteasome inhibitors, and/or a corticosteroid. In addition, 66% of the rrMM population received a stem cell transplant, which exposed them to high-dose melphalan. Only a small number of patients exposed (n = 38) or refractory (n = 21) to daratumumab were present. This data set was generated from trials of IMiDs, thus data from patients that were non-IMiD treated are not represented.

We classified the rrMM and ndMM genomes into MM subgroups based on canonical immunoglobulin heavy-chain translocations and hyperdiploidy (HRD, ≥3 copy number gain of odd-number chromosomes). To address the aim of this study and given that no significant difference in the distribution of MM subtypes was observed between ndMM and rrMM, we focused on subgrouping patients based on IMiD resistance (LENR and POMR, Figure 1). Because of sample size limitations, for a number of differential analyses we first compared the rrMM and ndMM patient data sets, and then sought differential enrichment in LENR and POMR subgroups.

We initially sought ndMM-specific drivers (Walker et al²; N = 1273) in the rrMM data set and detected 34 (out of 63) of these drivers (recurrence ≥2%). Mutational driver analysis on the rrMM data set also detected an additional 10 previously unreported promising candidate driver genes (recurrence ≥2%; false discovery rate [FDR] < 0.05), including a COSMIC cancer gene (https://cancer.sanger.ac.uk/census), EZH2. Analysis of clonality of drivers based on their CCF^14,24 distribution showed that novel driver genes were generally subclonal (Figure 2A). Analysis of the distribution of nonsilent mutations showed distinct variant types in TDG and MAML3 where the predicted functional consequence of most mutations was splicing defect and frameshift, respectively (Figure 2B, supplemental Table 2). All 44 drivers were examined to see whether they could be classified as TSG or ONC based on their mutation spectrum (Figure 2C). Among the novel genes, MAML3, TDG, and PIGO showed the strongest TSG signals, whereas NBPF15 showed a strong ONC signal similar to BRAF and NRAS/KRAS. We also analyzed all drivers for gene essentiality using DepMap.¹⁷ Novel drivers categorized as TSG such as DUOX2, PIGO, MAML3, and TDG showed overlapping distributions to well-known TSGs such as RB1, SETD2, and TRAF3. In addition, we saw a strong CRISPR effect score signal (median = −1.46) for LILRA6, a driver categorized as an ONC, similar to well-known ONCs such as IRF4 and SF3B1 (Figure 2D). Overall, when comparing ONC and TSG genes, we saw a highly significant elevated CRISPR effect score in ONC genes (mean of −0.63 vs mean of −0.18; P = 3.3 × 10⁻¹⁷).

Interestingly, a total of 46 LENR/POMR tumors (14.2%; Figure 2B) were negative for all mutational drivers. No significant difference was observed between LENR and POMR for the prevalence of these tumors (P = .2). These tumors were further assessed for translocations and CNA. In total, 40 tumors (87%) had at least 1 such major event in their genome (supplemental Figure 1), however, the rest were devoid of any known genomic drivers.

To identify which drivers are unique or enriched in rrMM and thus relevant to therapy-driven evolution, we compared their recurrence and clonality between the rrMM and ndMM data sets. Interestingly, although novel candidate drivers (except PIGO) were not rrMM-specific, when compared with ndMM, the most extreme changes in CCF and recurrence were observed in these drivers (Figure 3A). Specifically, DUOX2, EZH2, and PIGO had the highest enrichment (fold-change > 4.5). Interestingly, these 3 drivers were virtually mutually exclusive, suggesting alternate trajectories of therapeutic resistance (Figure 3B). To ensure robustness of these findings, we included the MMRF CoMMpass²⁵ whole-exome sequencing data set (N = 1001), making a combined ndMM data set (N = 1001 + 198 = 1199). We showed that the frequency of almost all drivers is overestimated in the WGS ndMM data set (supplemental Figure 2A) and, in turn, most enrichment fold-changes are higher when the combined ndMM is used. Interestingly, in addition to the top-enriched novel drivers remaining highly enriched (log₂FC ≥ 2; PIGO, DUOX2, and EZH2), MAML3 also became more enriched in rrMM (log₂FC = 2.9, supplemental Figure 2B).

We then sought to identify driver genes enriched across IMiD therapy stages (ie, ndMM to LENR to POMR). A positive trend in prevalence was observed for 16 driver genes (Figure 3C), among which TSGs were enriched (5.5-fold vs oncogenes). TP53 showed a statistically significant increase (proportion trend test, P = .016). Although the significance of DUOX2 and IRF4 was borderline owing to lower recurrence (P < .1), both were enriched at the POMR stage when compared with ndMM (7.3- and 2.9-fold increase, respectively). Interestingly, although there was no difference between ndMM and LENR in MAML3 frequency, POMR showed a 2.75-fold higher recurrence of this driver gene (Figure 3C).

To gain more insight into the biology of rrMM and to increase statistical power, we also analyzed the prevalence of nonsilent mutations in 10 canonical cancer pathways²⁶ across IMiD therapy stages. Of 3 pathways related to cell cycle, MYC and TP53 showed an overall positive trend (Figure 3D), however, only TP53 showed a statistically significant increase (proportion trend test, P = .049). Overall, rrMM genomes showed acquisition of novel TSG drivers such as EZH2, PIGO, and DUOX2 or enrichment of high-risk recurrent drivers such as TP53²⁷ and DIS3.²⁸ Furthermore, we undertook network and pathway enrichment analysis based on the 44 drivers to investigate possible mechanisms of disease progression potentially driven by the novel candidate drivers. First, we generated a functional interaction network among all driver genes using STRING version 11.5.²⁹ Of the 10 promising candidate drivers, TDG, MAML3, and EZH2 had 2, 3, and 5 interactions, respectively, with known MM drivers (supplemental Figure 3A), suggesting potential functional association. All 3 also clustered within the same subnetwork with known driver genes such as RB1, CCND1, and CREBBP (supplemental Figure 3B). Next, we used Enrichr³⁰ to identify enriched pathways unique to promising candidate drivers. Among the enriched pathways (supplemental Table 3), base excision repair (TDG), notch signaling pathway (MAML3), and B-cell receptor signaling pathway (LILRA6) were notable.

At the individual tumor level, we sought to assess the extent of biallelic events given their association with relapse in MM.³¹ We observed biallelic inactivation of 15 driver genes (Figure 3E), of which one (CDKN2A) was solely inactivated by homozygous deletion (HD) events (supplemental Figure 4). We showed that biallelic CDKN2A and CREBBP are unique to rrMM. This is consistent with that reported previously for ndMM where both were observed at only 0.1% and 0.2%, respectively (N = 1074).^2,TP53 was not only the most frequent biallelic event in rrMM (9.3%) but it was also enriched compared with ndMM (1.8 fold). Among the novel drivers, notably, TDG showed enriched biallelic inactivation (2.6 fold) at 1.3% prevalence in rrMM. Given its high TSG signal (Figure 2B) and recent evidence supporting its tumor suppressive role,³² this observation further corroborates that TDG acts as a TSG in MM.

Leveraging the WGS breadth of our rrMM and ndMM data sets, we sought to identify novel rrMM drivers by detecting differentially enriched somatic noncoding mutational hotspots flanking protein-coding genes (refer to Ansari-Pour et al³³). For this analysis, we used a construct similar to that of a genome-wide association study and looked for regions with a higher somatic mutation recurrence in the rrMM data set. A total of 10 genome-wide signals (FDR < 0.05) were identified (Figure 3F). Of note, TP53BP1, which is a haploinsufficient TSG,³⁴ has been previously shown to be involved in double-strand DNA break repair in MM³⁵ and its loss has also been reported to induce chemoresistance in colorectal cancer.³⁶ Given that these mutations are noncoding, we hypothesize that, in the context of DNA damaging drugs, acquiring such variants may lead to the dysregulation of genes such as TP53BP1, which may, in turn, confer resistance to tumors that harbor them.

SV analysis was undertaken to specifically detect translocations with a focus on those that involved the IGH, IGL, and MYC loci in the rrMM data set. IGH events (canonical or otherwise) were found in 42% of tumors (Table 1), similar to the ndMM data set (44%) and to that previously reported (41%),³⁷ with the most common IGH events being t(11;14), t(4;14), and t(8;14), whereas t(6;14) and t(14;16) were present in a few cases (Figure 4A). As expected, the first 2 were found to be early events in MM (ie, clonal), however t(8;14) was found to be at subclonal levels across the data set (Figure 4A). IGL- and MYC-centric translocations were present in 11% and 29.5% of tumors, respectively, and were mainly subclonal (Figure 4B-C), consistent with these events being secondary genomic aberrations.^38,39 The prevalence of these secondary translocations was then examined across clinical settings. MYC translocations showed an increase in prevalence from LENR to POMR (27% to 44%, odds ration [OR], 2.1; P = .023) but no difference was observed in VAF distribution after adjusting for purity (P = .157). A significant increase in prevalence was also observed for IGL events from LENR to POMR (10% to 22%, OR, 2.5; P = .021) but with no change in clonality (P = .68). This increase in IGL prevalence is consistent with a recent study showing that no significant survival benefit was provided by IMiD-containing regimens for patients with IGL translocations.³⁷ 

We compared the rrMM landscape of CNA events against those in ndMM (Figure 4D) and identified the differentially enriched CNAs. We observed significant overrepresentation of 1q Gain, 6q LOH, and 17p LOH (FDR < 0.01; Figure 4E). To trace CNA enrichment across therapy stages, significantly enriched copy number regions (supplemental Methods) were identified in the rrMM data set across the genome and their recurrence frequency was compared (Figure 4F). Notably, we observed a significant positive cline of 1q Gain and 17p LOH recurrence from ndMM to LENR to POMR (proportion trend test, P < .05). We also observed a positive cline of recurrence of LOH events at 1p and 14q but they did not reach statistical significance. The higher prevalence of these events with later therapy stage is consistent with these events being associated with disease progression and poor prognosis.^2,40,41 Although the rrMM patient population is extremely complex with concurrent layers of genomic, immune, clinical, and treatment heterogeneity, we nonetheless examined the prognostic effect of 1q Gain and 17p LOH events per trial. We only observed a significant prognostic effect in the MM-010 trial for 1q Gain based on progression-free survival (P = .031; supplemental Figure 5).

The associations among all mutational drivers, enriched CNAs, WGD, and recurrent translocations (TLs) were analyzed together to identify novel patterns of somatic interaction in the form of significant cooccurrence and mutual exclusivity in the rrMM data set. Genomic events with at least 1 significant pairwise association (FDR < 0.1) were identified (Figure 5A). The most significant interactions were the cooccurrence of TP53-LOH17p and mutual exclusivity of HRD-WGD, of which the latter has not been previously reported. When compared with ndMM (supplemental Figure 6), we observed a significant association between IRF4 and t(11:14) only in the rrMM data set (OR,10.4; P = 8 × 10⁻⁶).

We also sought to identify interactions in the form of double-hit events (ie, 2 independent events occurring at 2 different loci in the same genome at the patient level. Double-hit events were identified in 171 (44.3%) patients with rrMM (supplemental Figure 7). The rate of double-hit events was markedly lower in the ndMM cohort (30.8%), thus showing a significant enrichment in rrMM (OR, 1.8; P = 1.7 × 10⁻³). We calculated the frequency of each double-hit event in rrMM (Figure 5B) and ndMM (supplemental Figure 8). We observed 2 double-hit events enriched significantly in rrMM (FDR < 0.05; Figure 5C), namely Amp1q-ISS3 (12.4 fold) and 1qGain-17pLOH (3.7 fold), of which the latter was the most frequent double-hit event in rrMM (11%; Figure 5B). Interestingly, Amp1q-ISS3, which shows the highest enrichment, has been previously shown to be a high-risk subgroup in ndMM.⁴² In addition, 1qGain-17pLOH is the cooccurrence of the 2 most frequent enriched CNA in rrMM (Figure 4E). Across clinical settings, the overall rate of double-hit events was not significantly different between POMR and LENR (P = .17), and the same was true for each individual event (supplemental Figure 9).

Finally, we analyzed the mutational signatures to detect an etiological basis for exogenous and endogenous factors associated with relapse/refractory status. De novo extraction of mutational signatures detected 8 SBS COSMIC signatures¹⁸ (supplemental Figure 10A). All signatures, except SBS12, had been previously observed in ndMM^11,43 and rrMM⁴³ tumors. When fitting the identified COSMIC signatures to each tumor, we also included SBS-MM1^5,24,43 to see whether its inclusion better reconstructs the observed mutational spectrum. We observed an increase in both prevalence and activity of SBS2, SBS13, SBS12, and SBS-MM1 (Figure 6A), all of which, except SBS-MM1, are either associated with APOBEC (SBS2/SBS13) or defective mismatch repair (DMR; SBS12). SBS-MM1 is a relatively flat signature that has been associated with melphalan treatment and generally with rrMM.⁴³ Our analysis suggests that SBS-MM1 is likely to be associated with rrMM owing to higher prevalence and activity. However, we should note that the highly significant co-occurrence of SBS-MM1 with SBS8 (Figure 6B) is elusive and suggests that fine-tuned functional experiments are required to establish the underlying etiology of SBS-MM1 and how it may relate or differ from SBS8. Analysis of signature activity against mutation burden showed strong positive correlations with APOBEC and DMR signatures (Figure 6C). The activity of mutational signatures was also analyzed at clonal and subclonal stages of evolution. Of note, SBS12 (DMR) showed the strongest increase in activity among all signatures at the subclonal level (Figure 6D, supplemental Figure 10B). We therefore examined whether mutations in DMR-related genes led to the increase in DMR signature activity. We curated a list of 6 well-known mismatch repair genes (MSH2, MLH1, PMS2, MSH3, MSH6, and MLH3) and identified tumors with nonsilent mutations in them. We observed a slight enrichment (OR, 1.7) of DMR-related genes in tumors with increased DMR activity but owing to the small sample size of mutant tumors (N = 13), this was not statistically significant (P = .42). We further explored the DMR signature to see whether more genomic aberrations are present in tumors with increased DMR activity. We computed 2 genome-wide metrics, namely (a) proportion of genome altered (PGA) and (b) mutation burden. We observed no significant increase in PGA (mean of 35.1% in increased DMR vs mean of 32.7% in other tumors, P = .48). However, interestingly, a highly significant increase in mutation burden was observed in increased DMR tumors (mean of 35 286 vs mean of 17 998, 2-fold, P = 7.8 × 10⁻⁸), which was driven by a 2.2-fold increase in single nucleotide variants (mean of 31 884 vs mean of 14 409, P = 4.8 × 10^-9). This suggests that DMR may, in part, be an endogenous factor associated with expanding subclones under therapeutic selection pressure.

Recent studies characterizing the genomic features of rrMM have either analyzed small numbers of samples^3,5,44 or lacked the breadth of WGS.⁴⁵ In combination with deep WGS ndMM data,¹¹ we performed an unsupervised analysis of the largest deep WGS rrMM data set to date (total N = 616), thus providing, to the best of our knowledge, the most comprehensive analysis of the rrMM genomic landscape to date. We identified novel rrMM driver genes through coding and noncoding variation in the genome. Interestingly, all 10 coding-variant drivers showed either a significant clonal expansion or strong enrichment in recurrence in rrMM compared with ndMM. This observation suggests that these driver genes are associated with therapeutic resistance, either already present before therapeutic exposure but providing clonal advantage upon exposure to therapy or arising de novo and enriching at the rrMM stage at the time of acquisition of resistance. Functional analysis of these promising candidate drivers is nonetheless required to corroborate these findings and provide insight into the exact role of these novel drivers in therapeutic resistance. There is a possibility that we may have not detected these driver genes at the ndMM stage owing to their very low CCF (ie, tiny subclones). However, this is unlikely, given that the average sequencing coverage of the ndMM cohort was higher than that in the rrMM cohort (94× vs 83×; Kolmogorov-Smirnov test, P = 7 × 10⁻¹⁴). In addition to the enrichment of TP53 (and TP53-17pLOH biallelic events, along with its pathway in general) and IRF4, we also observed a positive cline in the prevalence of a novel highly enriched driver (DUOX2) from ndMM to LENR to POMR, which is associated with multiple cancers and involved in the regulation of the AKT pathway.⁴⁶ Furthermore, by detecting enrichment of somatic variants in noncoding hotspots,³³ a total of 10 previously unreported drivers in MM such as TP53BP1 and BLM were identified, which need to be examined at the transcriptomic level to assess whether they dysregulate these bona fide drivers.

When comparing CNA enrichment patterns, 1q Gain and 17p LOH were not only enriched in the rrMM data set, their occurrence as double-hits was also enriched. In addition, they had a significant increasing trend from LENR to POMR. We also observed an increase, albeit statistically nonsignificant, in the prevalence of LOH events at 1p and 14q with rrMM therapy stages. The nonsignificance of the increase may be because of the relatively small sample size of the POMR cohort (N = 55) in this study. Although we found 3 genomic features that distinguished POMR from LENR tumors (MAML3 driver gene and IGL/tMYC translocations), a larger POMR cohort would be beneficial to examining more robustly whether LEN and POM select subclones with distinct genomic drivers or POM generally selects the same subclones already expanding in LENR tumors. For instance, we identified 2 driver genes (NF1 and ADGRL3), which increased in prevalence from ndMM (0.015 and 0.02) to LENR (0.048 and 0.037) but were completely absent from POMR. It is not possible to determine whether their absence is because of favorable response to POM or simply not observed owing to the small size of this subgroup. To evaluate the latter possibility further, we calculated the expected upper limit frequency for an unobserved event⁴⁷ in POMR (q = 0.053). Because this frequency is above that observed in LENR, the absence of those 2 drivers may be owing to low statistical power.

Mutational signature analysis identified a novel enriched subclonal signature at rrMM (SBS12) with its etiology likely to be DMR. DMR tumors accumulate high levels of mutations (consistent with our findings; Figure 6C), a feature linked to acquisition of drug resistance,⁴⁸ including in myeloma.⁴⁹ Alternative therapies such as those based on synthetic lethal approaches^50,51 that could exploit DMR in mismatch repair–deficient rrMM tumors may be fruitful.

Our analysis provides insight into the nature of the evolution of myeloma genome under therapeutic pressure. Although in rrMM, therapeutic clonal selection might be expected to reduce the intratumoral heterogeneity, we instead find that ongoing clonal evolution provides genomes with highly diverse enrichment of LOH events, mutational signatures, genetic aberrations, and subclonal complexity. It is therefore likely that ongoing acquisition of novel features is occurring alongside enrichment of those already present. The increase in high-risk changes associated with shortened survival (such as the double-hit events described) alongside direct selection of therapeutic resistance markers are thus unsurprising. However, this reconstruction of the rrMM genome, which blends acquisition of the functional high-risk myeloma state with drug-specific resistance features (such as CRBN loss) pose significant challenges in designing new therapies. Specifically, there must be an even greater imperative to consider clonal evolution in pinpointing which events are connected to disease relapse (ie, present in the most recent common ancestor of the expanding subclone[s]). A limitation of our study was the unrelatedness of the ndMM and rrMM cohorts. We⁶ and others^5,52,53 have shown that sequential sampling is the ideal approach for analyzing clonal dynamics. Detection of subclones that predict relapse might enable investigation of change of therapy and thus their appearance or expansion must be carefully considered as we plan the trials and therapeutic approaches of the future.

The authors thank Marissa Hirst, Hillary Mosso, and Dan Rozelle from Rancho Biosciences for data curation, cleaning, and summary of patient demographics for this analysis. The authors also thank Garth McGrath from Bristol Myers Squibb for significant assistance identifying samples that made this analysis possible and for assistance with data curation. B.W. is partly funded by a grant from the Leukemia & Lymphoma Society.

Contribution: E.F. and A.T. conceived and designed the project, and acquired funding; E.F., F.T., and M.O.E., performed project administration; E.F., M.O.E., N.S., S.G., K.M., N.A.-P., F.T. H.A.-L., M.S., N.M., G.M., B.W., and A.T performed oversight and management of resources (data generation, collection, transfer, infrastructure, data processing); N.A.-P. performed computational and statistical analyses; N.A.-P., S.G., M.O.E., N.S., F.T., E.F., B.W., H.A.-L., G.M., N.M., M.S., and A.T. performed analyses and interpretation; N.A.-P., M.S., M.O.E., F.T., N.S., S.G., E.F., and A.T. performed and structured data visualization; A.T. provided supervision and scientific direction; N.A.-P., S.G., M.S., M.O.E., E.F., and A.T. wrote the manuscript; and all authors critically reviewed the manuscript and figures.

Conflicts-of-interest disclosure: Bristol Myers Squibb (BMS) Corporation has employment and equity ownership; funding for data generation, processing, and storage was provided by BMS Corporation; F.T. and A.T. were BMS employees at the time this analysis was conducted; B.W. receives research funding from BMS; and N.A.-P. is a BMS consultant. The remaining authors declare no competing financial interests.

Correspondence: Anjan Thakurta, Nuffield Department of Orthopaedics, Rheumatology and Musculoskeletal Sciences, University of Oxford, Windmill Road, Oxford, OX3 7LD United Kingdom; e-mail: anjan.thakurta@ndorms.ox.ac.uk.