Key Points

Abstract

Childhood acute lymphoblastic leukemia (cALL) is the most frequent pediatric cancer. Refractory/relapsed cALL presents a survival rate of ∼45% and is still one of the leading causes of death by disease among children. Mechanisms, such as clonal competition and evolutionary adaptation, govern treatment resistance. However, the underlying clonal dynamics leading to multiple relapses and differentiating early (<36 months postdiagnosis) from late relapse events remain elusive. Here, we use an integrative genome-based analysis combined with serial sampling of relapsed tumors (from primary tumor to ≤4 relapse events) from 19 pre–B-cell cALL patients (8 early and 11 late relapses) to assess the fitness of somatic mutations and infer their ancestral relationships. By quantifying both general clonal dynamics and newly acquired subclonal diversity, we show that 2 distinct evolutionary patterns govern early and late relapse: on one hand, a highly dynamic pattern, sustained by a putative defect of DNA repair processes, illustrating the quick emergence of fitter clones, and on the other hand, a quasi-inert evolution pattern, suggesting the escape from dormancy of leukemia stem cells likely spared from initial cytoreductive therapy. These results offer new insights into cALL relapse mechanisms and highlight the pressing need for adapted treatment strategies to circumvent resistance mechanisms.

Introduction

Childhood acute lymphoblastic leukemia (cALL) is the most frequent pediatric cancer, accounting for ∼25% of all cases.1  Despite improved treatment strategies, ∼20% of patients are either refractory or ultimately relapse, making cALL the second cause of cancer-related mortality among children and adolescents.2-5 

Tumor plasticity, based on the random accumulation of somatic mutations combined with external pressures, such as chemotherapy, is a breeding ground for clonal competition and evolutionary adaptation leading to treatment resistance.6,7  Over the past few years, next-generation sequencing (NGS) has provided meaningful information toward the characterization of these evolutionary processes and the identification of somatic events providing selective advantages. In the context of hematological malignancies, mutations in the RAS pathway genes have been associated with early relapse and chemoresistance in cALL.8  Hyperdiploid leukemia cells mutated in the histone acetyltransferase (HAT) domain of CREB-binding protein (CREBBP) were identified in glucocorticoid-resistant patients,9,10  whereas the mutated cytosolic 5′-nucleotidase II gene (NT5C2) was shown to be involved in resistance to treatment using nucleoside analog therapies.11,12  Such studies indicated that treatment could either lead to the eradication of nonresistant clones and the emergence of fitter and more aggressive subclones or, in rarer cases, to the persistence of an already existing resistant and dominant clone.7,13-17 

Although important, these studies generally offered only snapshots of tumors and therefore failed to capture the complexity of their temporal evolution.6  The latter is essential to decipher the full spectrum of relapse determinants. Furthermore, early cALL relapses, occurring during treatment, and late relapses, occurring at least 36 months after the initial diagnosis, are likely to present different evolutionary mutational processes. However, large scale genomic studies have thus far focused on early events only13  and have not addressed the yet elusive origin of long-term cALL relapses.

To decipher the full spectrum of somatic events that lead to treatment failure in cALL and shed light on the spatiotemporal clonal evolution of cALL relapse, we performed both whole-genome sequencing (WGS) and whole-exome sequencing (WES) on a series of primary tumor (PT) samples and multiple relapses (≤4 relapses [R1-R4]) from 19 pre–B-cell cALL patients. Based on variant allele frequency (VAF) analysis, we inferred the global clonal population frequencies as well as evolutionary or ancestral relationships. By dividing our cohort into early (<36 months postdiagnosis, n = 8) or late bone-marrow relapse events (n = 11), we captured significant different clonal dynamics and showed a variable clonal origin of late relapse events. This study is the first to use serial sampling to show the difference of clonal evolution in late vs early pre–B-cell cALL relapses.

Methods

Study subjects

All study subjects were French Canadians of European descent. Incident cases were diagnosed in the Division of Hematology-Oncology at the Sainte-Justine Hospital (Montreal, Canada) as part of the Quebec cALL cohort (Table 1).18  Bone marrow samples were collected at diagnosis and/or relapse. Paired normal DNA specimens were collected from peripheral blood samples without blast cells. The CHU Sainte-Justine Research Ethics Board approved the protocol. Informed consent was obtained from the parents of the patients to participate in this study and for publication of this report and any accompanying images.

Table 1.

Clinical characteristics of the pre–B-cell cALL relapse patients

Patient IDSexAge at diagnosis, moYear of diagnosisRisk groupWBC, × 109 /LPlatelet, × 109 /LTreatment protocolRadiationBMTPre-BMT therapyRelapse event numberTumor purity, %TypeEvent-free survival, dDeathSequenced events
34 129 1995 High 2.8 323.0 DFCI 91-01 Y (after R1) — NA|94 Relapse 2147 PT|R1 
62 29 1996 Standard 22.3 85.0 DFCI 91-01 Y (after R1) — 95|97 Relapse 1025 PT|R1 
64 67 1991 Standard 3.8 30.0 DFCI 91-01 — NA|95 Relapse 1600 PT|R1 
217 23 1997 High 99.2 16.0 DFCI 95-01 Y (after R1) — 94|87 Relapse 259 PT|R1 
325* 160 1999 High 7.7 12.0 DFCI 95-01 — 98|94|95|81|99 Relapse 306|106|309|41 PT|R1|R2|R3|R4 
382 60 2000 NA 4.7 NA DFCI 95-01 — 96|70 Relapse 534|105 R1|R2 
391 35 2000 Standard 10.0 47.0 DFCI 95-01 — 98|92 Relapse 1143 PT|R1 
394 92 2000 Standard 10.5 11.0 DFCI 95-01 — 99|95 Relapse 1189 PT|R1 
445 45 1997 NA NA NA DFCI 2000-01 Y (after R1) — 97|93 Refractory (at 2nd induction) 2227|106 R1|R2 
447 129 2002 High 61.1 94.0 DFCI 2000-01 — 98|95 Relapse 982 PT|R1 
579* 75 2002 High 183.1 29.0 DFCI 2000-01 Y (after R1) Antithymocyte globulin, busulfan, fludarabine 96|90|91 Relapse 1240|496 PT|R1|R2 
659 35 2004 Standard 29.91 24.0 DFCI 2000-01 — 100|97 Relapse 1420 PT|R1 
670 33 2005 Standard 5.48 20.0 DFCI 2005-01 — 93|94 Relapse 1199 PT|R1 
684* 74 2005 NA 8.49 317.0 DFCI 2005-01 Y (after R1) Y (after R1) Full-body radiation, cyclophosphamide, etoposide 84|42|48 Relapse 1153|588 PT|R1|R2 
717 73 2006 Standard 5.55 24.0 DFCI 2005-01 — 100|65 Relapse 1380 PT|R1 
764 55 2008 High 59.2 21.0 DFCI 2005-01 — 96|73 Relapse 1107 PT|R1 
772* 39 2008 Standard 29.72 4.0 DFCI 2005-01 Y (after R1) Y (after R1) Full-body radiation, antithymocyte globulin, cyclophosphamide, etoposide 93|45 Refractory (at 2nd induction) 736|228 R1|R2 
786 206 2009 High 120.9 57.0 DFCI 2005-01 — 97|98|95 Refractory 35|22 PT|R1|R2 
808* 15 2010 High 265.7 47.0 DFCI 2005-01 Y (after PT) Y (after PT) Full-body radiation, antithymocyte globulin, cyclophosphamide, etoposide 93|38 Refractory 199 PT|R1 
Patient IDSexAge at diagnosis, moYear of diagnosisRisk groupWBC, × 109 /LPlatelet, × 109 /LTreatment protocolRadiationBMTPre-BMT therapyRelapse event numberTumor purity, %TypeEvent-free survival, dDeathSequenced events
34 129 1995 High 2.8 323.0 DFCI 91-01 Y (after R1) — NA|94 Relapse 2147 PT|R1 
62 29 1996 Standard 22.3 85.0 DFCI 91-01 Y (after R1) — 95|97 Relapse 1025 PT|R1 
64 67 1991 Standard 3.8 30.0 DFCI 91-01 — NA|95 Relapse 1600 PT|R1 
217 23 1997 High 99.2 16.0 DFCI 95-01 Y (after R1) — 94|87 Relapse 259 PT|R1 
325* 160 1999 High 7.7 12.0 DFCI 95-01 — 98|94|95|81|99 Relapse 306|106|309|41 PT|R1|R2|R3|R4 
382 60 2000 NA 4.7 NA DFCI 95-01 — 96|70 Relapse 534|105 R1|R2 
391 35 2000 Standard 10.0 47.0 DFCI 95-01 — 98|92 Relapse 1143 PT|R1 
394 92 2000 Standard 10.5 11.0 DFCI 95-01 — 99|95 Relapse 1189 PT|R1 
445 45 1997 NA NA NA DFCI 2000-01 Y (after R1) — 97|93 Refractory (at 2nd induction) 2227|106 R1|R2 
447 129 2002 High 61.1 94.0 DFCI 2000-01 — 98|95 Relapse 982 PT|R1 
579* 75 2002 High 183.1 29.0 DFCI 2000-01 Y (after R1) Antithymocyte globulin, busulfan, fludarabine 96|90|91 Relapse 1240|496 PT|R1|R2 
659 35 2004 Standard 29.91 24.0 DFCI 2000-01 — 100|97 Relapse 1420 PT|R1 
670 33 2005 Standard 5.48 20.0 DFCI 2005-01 — 93|94 Relapse 1199 PT|R1 
684* 74 2005 NA 8.49 317.0 DFCI 2005-01 Y (after R1) Y (after R1) Full-body radiation, cyclophosphamide, etoposide 84|42|48 Relapse 1153|588 PT|R1|R2 
717 73 2006 Standard 5.55 24.0 DFCI 2005-01 — 100|65 Relapse 1380 PT|R1 
764 55 2008 High 59.2 21.0 DFCI 2005-01 — 96|73 Relapse 1107 PT|R1 
772* 39 2008 Standard 29.72 4.0 DFCI 2005-01 Y (after R1) Y (after R1) Full-body radiation, antithymocyte globulin, cyclophosphamide, etoposide 93|45 Refractory (at 2nd induction) 736|228 R1|R2 
786 206 2009 High 120.9 57.0 DFCI 2005-01 — 97|98|95 Refractory 35|22 PT|R1|R2 
808* 15 2010 High 265.7 47.0 DFCI 2005-01 Y (after PT) Y (after PT) Full-body radiation, antithymocyte globulin, cyclophosphamide, etoposide 93|38 Refractory 199 PT|R1 

Information concerning the different relapse events are separated by a bar in the cells.

–, not applicable; BMT, bone marrow transplant; DFCI, Dana-Farber Cancer Institute; F, female; M, male; N, no; NA, not available; Y, yes.

*

Hypermutable case.

Genome sequencing and variant identification

Exomes (Nextera Rapid Capture Exome Enrichment) and genomes were sequenced on Illumina HiSequation 2000/2500 systems to reach an average coverage of 200× and 80×, respectively. Resulting reads were aligned to the hg19 (Human Genome version 19) reference genome using Bowtie2 (version 2.2.3)19  and Casava software, respectively. Single nucleotide variants (SNVs) and small indels were called using Strelka20  (supplemental Figure 1). To uncover any missed subclonal mutations, the variations identified in each tumor event were screened directly in pileup files of previous and/or next events using an adapted version of SNooPer.21  Copy number variants (CNVs) were called using the Varscan2 copyCaller option22  (supplemental Figure 1). Additional structural variations (SVs) were called from WGS using the Casava software (supplemental Figure 1). Variants were prioritized as previously described.23  See supplemental Data (“Variant annotation and prioritization of cancer driver gene mutations”) for details.

Results

The genomic landscape of pre- and posttreatment cALL

To investigate the temporal evolution of the tumors, we performed serial sampling of 19 pre–B-cell cALL cases who suffered 1 or multiple relapses (R1-R4) and obtained a total of 63 samples (Table 1). The cytogenetic information regarding genome ploidy and large chromosome abnormalities is summarized in Table 2. Using deep WES and/or WGS associated with a strict data reduction strategy (supplemental Methods; supplemental Figure 1), we generated a comprehensive repertoire of somatic mutations, including SNVs and insertions/deletions (indels) (Figure 1; supplemental Table 1), as well as structural variations (SVs), such as CNVs and cryptic inversions or translocations (Figure 1; supplemental Table 2). VAFs of both SNVs and indels were calculated from the distribution of reads either supporting the mutated or the reference allele. The comparison of VAFs obtained for samples for which both WES and WGS were performed indicated that sequence capture had limited influence on read distributions (supplemental Figure 2).

Table 2.

Karyotypes of the pre–B-cell cALL relapse patients

Patient IDKaryotype (PT)Karyotype (R events)
34 46,XY 46,XY 
62 54,XY,+4,+6,+8,+10,+14,+17,+18,+21 55,XXY,+4,+6,+8,+10,+14,+der(17)del(17)(p12),+18,+21[12]/46,XY[7] 
64 46,XX 46,XX,t(7;11)(q32;q14),+der(11)t(1;7;11)(q21;q32;q14),-21[17]/46,XX,idem,t(9;22)(p24;q22.1)[8] 
217 46,XX,inv(4)(p14q27),del(9)(p22) 46,XX,inv(4)(p14q28),del(9)(p22),t(12;14)(q12;q11)[8]/46,XX,inv(4)(p14q28),del(9)(p22)[3] 
325* 46,XY NA|NA|NA|NA 
382 47,XY,t(9;14)(p13;q32),+mar NA | 47,XY,del(2)(p22?),add(3)(p25),t(9;14)(p13;q32),t(16;Y)(q24;q11.2),+mar 
391 52,XY,+3,+6,+11,+17,+21,+22 52,XY,+X,+6,+14,+17,+21,+21 
394 46,XX,-12,-21,+mar1,+mar2 46,XX,-12,-21,+mar1,+mar2[1]/46,XX,der(12),-21,+mar3[5] 
445 NA 47∼48,-X,der(1),del(1)x2,-6,der(8),der(17),+21,+21,+mar | NA 
447 47,XX,add(6)(q27),der(15)t(?13;15)(p11;p11),add(20)(?q11.2),+21 47,XX,add(1)(q32),add(6)(q27),der(15)(p11),+21 
579* 46,XY 46,XY | 46,XY 
659 55,XXY,+2,+?3,+4,+6,+10,+14,+18,+21 54,XY,+4,+6,del(9)(p21),+10,del(12)(p13),+17,+21 
670 51∼56,XY,+X,+2,+4,+6,+9,+10,+14,+17,+21 53,XY,+4,+6,+8,+10,+14,+17,+21 
684* 46,XX 46,XX | 46,XX,dup(1)(q21q23) 
717 46,XX,hsr(21)(q22.1) 46,XX,hsr(21)(q22.1) 
764 56,XY,+X,+4,+5,+6,del(9p),+10,+14,+17,+18,+21,+21 NA 
772* 23∼39,X,Y,+3,+5,+14,+16,+17,+20,[9]/52,XY,+X,+Y,+14,+20,+21,+21[13] 26∼27,XY,+14,+18,+21 
786 NA NA | NA 
808* 46,XY,t(11;19)(q23;p13.1) 46,XY,t(11;19)(q23;p13.1) 
Patient IDKaryotype (PT)Karyotype (R events)
34 46,XY 46,XY 
62 54,XY,+4,+6,+8,+10,+14,+17,+18,+21 55,XXY,+4,+6,+8,+10,+14,+der(17)del(17)(p12),+18,+21[12]/46,XY[7] 
64 46,XX 46,XX,t(7;11)(q32;q14),+der(11)t(1;7;11)(q21;q32;q14),-21[17]/46,XX,idem,t(9;22)(p24;q22.1)[8] 
217 46,XX,inv(4)(p14q27),del(9)(p22) 46,XX,inv(4)(p14q28),del(9)(p22),t(12;14)(q12;q11)[8]/46,XX,inv(4)(p14q28),del(9)(p22)[3] 
325* 46,XY NA|NA|NA|NA 
382 47,XY,t(9;14)(p13;q32),+mar NA | 47,XY,del(2)(p22?),add(3)(p25),t(9;14)(p13;q32),t(16;Y)(q24;q11.2),+mar 
391 52,XY,+3,+6,+11,+17,+21,+22 52,XY,+X,+6,+14,+17,+21,+21 
394 46,XX,-12,-21,+mar1,+mar2 46,XX,-12,-21,+mar1,+mar2[1]/46,XX,der(12),-21,+mar3[5] 
445 NA 47∼48,-X,der(1),del(1)x2,-6,der(8),der(17),+21,+21,+mar | NA 
447 47,XX,add(6)(q27),der(15)t(?13;15)(p11;p11),add(20)(?q11.2),+21 47,XX,add(1)(q32),add(6)(q27),der(15)(p11),+21 
579* 46,XY 46,XY | 46,XY 
659 55,XXY,+2,+?3,+4,+6,+10,+14,+18,+21 54,XY,+4,+6,del(9)(p21),+10,del(12)(p13),+17,+21 
670 51∼56,XY,+X,+2,+4,+6,+9,+10,+14,+17,+21 53,XY,+4,+6,+8,+10,+14,+17,+21 
684* 46,XX 46,XX | 46,XX,dup(1)(q21q23) 
717 46,XX,hsr(21)(q22.1) 46,XX,hsr(21)(q22.1) 
764 56,XY,+X,+4,+5,+6,del(9p),+10,+14,+17,+18,+21,+21 NA 
772* 23∼39,X,Y,+3,+5,+14,+16,+17,+20,[9]/52,XY,+X,+Y,+14,+20,+21,+21[13] 26∼27,XY,+14,+18,+21 
786 NA NA | NA 
808* 46,XY,t(11;19)(q23;p13.1) 46,XY,t(11;19)(q23;p13.1) 

Karyotypes of each relapse event per patient are separated by a vertical bar.

NA, not available.

*

Hypermutable cases.

Figure 1.

Putative driver mutations (somatic SNVs, small indels, and CNVs) identified among 19 pre–B-cell cALL relapse patients and grouped into relevant signaling pathways. Genes (top) are ordered according to the number of events identified in the cohort (descending order from left to right). Pathways are ordered according to the number of genes identified as mutated (descending order from left to right). The mutation color scale indicates the VAFs of SNVs and indels (from light green [VAF ≤ 0.1] to dark blue [VAF = 1]). CNVs are indicated in gray. The mutation rate color scale indicates to which dynamic or subclonal quartile the event belongs to (from light yellow [Q1, low values] to red [Q4, high values]). Numbers in cells indicate the number of multiclonal mutations co-occurring in the same gene. μ, mutation rate.

Figure 1.

Putative driver mutations (somatic SNVs, small indels, and CNVs) identified among 19 pre–B-cell cALL relapse patients and grouped into relevant signaling pathways. Genes (top) are ordered according to the number of events identified in the cohort (descending order from left to right). Pathways are ordered according to the number of genes identified as mutated (descending order from left to right). The mutation color scale indicates the VAFs of SNVs and indels (from light green [VAF ≤ 0.1] to dark blue [VAF = 1]). CNVs are indicated in gray. The mutation rate color scale indicates to which dynamic or subclonal quartile the event belongs to (from light yellow [Q1, low values] to red [Q4, high values]). Numbers in cells indicate the number of multiclonal mutations co-occurring in the same gene. μ, mutation rate.

A mean of 145 (range: 5-1852, median = 21) and 891 (range 8-9,985, median = 31) coding somatic mutations were identified per PT and first relapse event, respectively (supplemental Table 3). Sixty-eight percent (in PT) and 52% (in relapse) were considered subclonal mutations (VAFs adjusted for tumor purity <0.30; supplemental Table 3). Of note, local or complete deletions of the short arm of chromosome 9 was one of the most recurrent SVs observed with 16% and 21% of cases presenting a variation at diagnosis and relapse, respectively (Table 2; supplemental Table 2). Nonsynonymous SNVs, including missense and nonsense variants, were the most prevalent class of mutation in all events and for all cases (supplemental Figure 3). On average, 26% and 46% were persistent between PT/R1 and R1/R2, respectively (supplemental Table 3). All patients presented shared mutations between the PT and subsequent events, suggesting a common ancestral origin of lymphoblasts, as previously reported.13,24  Notably, all deletions of chomosome 9p persisted to relapse (Table 2; supplemental Table 2). Five relapse events (case 325 R3, case 579 R2, case 684 R2, case 772 R2, and case 808 R1) carried >10 times the number of coding somatic mutations identified in the precedent event (supplemental Table 3). These patients were considered hypermutable and were either refractory or suffered multiple relapse events. Importantly, 4 of these (579 R2, 684 R2, 772 R2, and 808 R1) were consecutive to a BMT (Table 1), suggesting an association between their aggressive chemotherapy regimens and the phenotype (P = .004, 2-tailed Fisher's exact test). Of note, we identified no association between this phenotype and irradiation (P = .262). The mutation spectra of the 5 hypermutable cases revealed a substantial increase of transition mutations [A(T)>G(C) and C(G)>T(A)] from the PT (mean percentage of total spectrum = 29%) to the hypermutated event (72%) (supplemental Figure 4), previously identified as a chemotherapy-related mutational signature in other cancer types and reinforced by defective DNA repair mechanisms.25-27  Except for these cases, no significant shift in mutation number was observed (supplemental Table 3). The hypermutated cases 684, 772, and 808 were the only patients from this cohort to harbor a p.G39E germ line polymorphism (rs1042821) in the DNA repair gene MSH6, a known predisposing gene in Lynch syndrome28-30  that is also associated with an increased risk of several cancer types.31  As for patient 325, who did not undergo transplant, the acquisition at R3 of a somatic homozygous loss of function mutation in the mismatch repair endonuclease PMS2 (p.R134*, VAF = 0.77), which presented a modest positive drift at R4 (VAF = 0.85; Figure 1; supplemental Table 1), could have contributed to the hypermutator phenotype. Case 808 also carried the p.R20Q polymorphism in PMS2 that has been associated with reduced apoptotic function.32 

Putative driver mutations were grouped into relevant signaling (“Methods”; Figure 1; supplemental Tables 1 and 2) and therapy-related pathways (“Methods”; supplemental Figure 5; supplemental Tables 1 and 2). Six signaling pathways were identified as recurrent targets (Figure 1; supplemental Table 1). Of these, epigenetic regulation and MAP kinase–RAS signaling showed the highest mutation rates with, overall, 84.2% and 78.9% of patients presenting ≥1 altered gene, respectively. Interestingly, these pathways showed numerous putative driver events, yet counter-selected at relapse (13 of 47 mutations vs 3 of 36 mutations in other pathways; Figure 1; supplemental Table 1). Although each pathway showed a significant proportion of relapse-specific mutations (range: 22.2%-75.0%), most alterations identified in MAP kinase–RAS signaling were already present at diagnosis (8.7% relapse-specific mutations; Figure 1; supplemental Table 1). Multisubclonal hits in 5 genes (KDM6A, SETD2, NRAS, NR3C1, and TP53) were identified in 4 early relapse cases (cases 325, 382, 447 and 772) and 1 late relapse patient (case 659) (Figure 1; supplemental Table 1). Of note, mutations in NR3C1 (p.D724E and p.Y641H) and KDM6A (p.R1351*) were selected at the expense of co-occurring mutations (p.P626S in NR3C1 and p.R519* in KDM6A).

Overall, a slightly higher proportion of early relapse cases compared with late ones were identified as mutated for 4 of the 6 pathways: 100% vs 72.7% of mutated cases for the epigenetic regulation pathway, 75.0% vs 54.5% for the hemopoiesis–immune system development pathway, 87.5% vs 72.7% for the MAP kinase–RAS signaling pathway, and 62.5% vs 36.4% for the cell cycle pathway. A strong yet not significant effect was observed for the DNA repair signaling pathway, with 62.5% of early relapse cases mutated compared with only 18.2% of patients with late relapses (P = .07, 2-tailed Fisher's exact test). Of note, 4 of the 5 early relapse cases harbored these DNA repair pathway mutations within the dominant clonal population. Finally, 80% of patients who presented a hypermutator phenotype harbored ≥1 mutation in the DNA repair pathway vs 21% for other cases (P = .04). No association was identified between relapse types and cytogenetic subgroups, such as hyperdiploid (P = 1.0) and normal karyotypes (P = .3).

Early and late relapse cALL cases show distinct clonal dynamics

To evaluate the influence of tumor plasticity on outcome, we divided our patients into groups according to their progression-free survival (early vs late event), their response to treatment (nonrefractory vs refractory), and their survival status. A “dynamic mutation rate” (number of dynamic mutations per Mb per day), considering only mutations showing significant VAF shifts, was calculated for each relapse (supplemental Data [“Dynamic and subclonal mutation rates”]; supplemental Figure 6). Each value was representative of the clonal dynamics of a tumor between 2 events (either PT vs R1 or R[n] vs R[n + 1]) (Figures 1 and 2).

Figure 2.

Early and late relapses show distinct clonal dynamics. The “dynamic mutation rate” (number of dynamic mutations per Mb per day) was calculated from WES data by considering mutations showing a significant VAF shift only (see “Methods”). Each value represents the clonal dynamics of a tumor between 2 events (either PT vs R1 or R[n] vs R[n + 1]). Patients are divided into 6 panels depending on their survival or refractory status and on the delay before relapse events. R1 stands for the analysis of PT vs R1, R2 for R1 vs R2, R3 for R2 vs R3, and R4 for R3 vs R4. The dynamic mutation rate was significantly reduced in late relapses compared with early events (P = .0094, Mann-Whitney U test) and in nonrefractory compared with refractory events (P = .0473, Mann-Whitney U test).

Figure 2.

Early and late relapses show distinct clonal dynamics. The “dynamic mutation rate” (number of dynamic mutations per Mb per day) was calculated from WES data by considering mutations showing a significant VAF shift only (see “Methods”). Each value represents the clonal dynamics of a tumor between 2 events (either PT vs R1 or R[n] vs R[n + 1]). Patients are divided into 6 panels depending on their survival or refractory status and on the delay before relapse events. R1 stands for the analysis of PT vs R1, R2 for R1 vs R2, R3 for R2 vs R3, and R4 for R3 vs R4. The dynamic mutation rate was significantly reduced in late relapses compared with early events (P = .0094, Mann-Whitney U test) and in nonrefractory compared with refractory events (P = .0473, Mann-Whitney U test).

This way, we observed marked reduction in clonal dynamics in late relapses (>36 months postdiagnosis) compared with early relapse events (WES data, P = .0094, Mann-Whitney U test) and for nonrefractory compared with refractory events (P = .0473) (Figures 2 and 3). The same pattern was observed using WGS data for which first relapses (R1) only were available (P = .0079; supplemental Figure 7). Given that R1 samples were obtained at the moment these relapses were diagnosed, before any relapse-associated treatment, this therapy cannot explain the observed difference in clonal dynamics between early and late subgroups. Of note, the absolute number of somatic mutations was not significantly different between early and late relapse groups (P = .9, Mann-Whitney U test), and no difference in the distribution of mutation classes (eg, synonymous, missense, nonsense) was observed (supplemental Figure 3). These results illustrate the plasticity of the earliest events, with rapid clonal switches over short periods of time (eg, case 808 with a mean dynamic mutation rate 1.4e-01 mutations/Mb per day and case 325 with 2.3e-02 mutations/Mb per day) compared with those associated with long-term and quasi-inert relapses (eg, case 34 with 1.3e-04 mutatons/Mb per day). Interestingly, a similar trend of reduced clonal dynamics was observed in WGS data for relapse events of patients who survived compared with those who did not (P = .0159, supplemental Figure 7). We also observed a progressive increase of dynamics at each new relapse event for the same patient (eg, case 325; Figure 1 and 2).

Figure 3.

Higher clonal dynamics correlate with shorter event-free survival. Pre–B-cell cALL cases were divided into groups of low and high clonal dynamics based on the median dynamic mutation rate value. A Kaplan-Meier curve estimating the relapse probability over time was constructed for each group. Solid and dashed lines stand for clonal dynamics over and under the median value, respectively. Cases with higher clonal dynamics had significantly shorter event-free survival (P = .0067, χ2 test).

Figure 3.

Higher clonal dynamics correlate with shorter event-free survival. Pre–B-cell cALL cases were divided into groups of low and high clonal dynamics based on the median dynamic mutation rate value. A Kaplan-Meier curve estimating the relapse probability over time was constructed for each group. Solid and dashed lines stand for clonal dynamics over and under the median value, respectively. Cases with higher clonal dynamics had significantly shorter event-free survival (P = .0067, χ2 test).

To assess the predictive value of the clonal dynamics in regards to relapse-free survival (first relapse event) while taking variables such as therapy differences into account, we conducted a multivariable analysis (Cox proportional hazard model) adjusted for sex, age at diagnosis, white blood cell (WBC) count, cytogenetic risk, as well as treatment protocol (Dana-Farber Cancer Institute 1991, 1995 2000, and 2005). We showed that a high dynamics (greater than the median dynamic mutation rate value) was independently associated with a significantly lower relapse-free survival (adjusted hazard ratio = 9.92 [95% confidence interval, 1.48-66.52], P = .01). Treatment protocols and clinical characteristics showed no significant influence (P = .4, .4, .7, .2, and .6 for sex, age at diagnosis, clinical risk group, WBC, and treatment protocol, respectively).

Assuming that subclonal diversity plays a central role in therapeutic failure and relapse,3,6,7,33  we quantified the number of newly acquired subclonal mutations representative of the general expansion of subclones at each event. The estimated “subclonal mutation rate” revealed the exact same pattern as reported above with a significantly reduced subclonal expansion in late events compared with early events (P = .0053 and P = .0317 for WES and WGS data, respectively; Figure 1; supplemental Figures 8 and 9) and in nonrefractory compared with refractory events (P = .0022, WES data; supplemental Figure 8). Again, cases with higher subclonal mutation rates were also the ones who did not survive (P = .0159, WGS data; supplemental Figure 9).

Clonal evolution of early and late relapses

Based on the analysis of somatic VAFs in loss of heterozygosity–free regions of the genome, we inferred the number of subclones and tracked their dynamics over time (supplemental Data [“Spatial and temporal analysis of clonal populations”]).34  We also built clonal trees to determine the ancestral relationships between somatic mutations (supplemental Data [“Spatial and temporal analysis of clonal populations”]).35  Contrasting clonal evolution phenotypes are illustrated by 3 cases, 34, 325, and 62 (Table1; Figure 4). Case 34 presented the longest period before relapse (close to 6 years), but also showed the slowest evolving tumor of our cohort (Figure 3). As for case 325, one of the patients with a hypermutator phenotype, he presented a first and early relapse event less than a year after PT (day 306) followed by multiple relapse events (R1-R4) (Table 1). With a relapse occurring 1025 days after PT, case 62 was the last patient of the early event group (Table 1). He showed the slowest evolving tumor of this group (supplemental Figure 7) and presented an overall number of mutations comparable to patient 34 (supplemental Table 3). Case 34 showed a near perfect clonal equilibrium with the same tumor architecture reoccurring 6 years postdiagnosis (Figure 4A-B). Six clusters of clones were detected. The evolution of the mutation p.W2006* in MLL2/KMT2D (VAF-PT = 0.31, VAF-R1 = 0.55) was the only difference observed between the 2 events. The re-emergence of subclonal mutations at relapse with similar frequencies (eg, KRAS p.K117N, VAF-PT = 0.03, VAF-R1 = 0.05) further highlighted the inertia of this tumor. Clonal evolution in this patient was reminiscent of an escape from long-term dormancy of premalignant neoplastic stem cells that also initiated the primary leukemia. Although this case was the most extreme example of long-term clonal equilibrium and the only one that presented the exact same clonal architecture at diagnosis and relapse, all late relapses presented common mutations with PTs (supplemental Table 3). On the other hand, with >25 SNVs/indels identified in tumor suppressors, oncogenes, and therapy-related genes, case 325 presented a complex clonal architecture composed of 9 clusters of clones spatially evolving over the 5 time points (Figures 1 and 4C), including the hypermutable event R3. Although case 325 did not present the highest evolutionary dynamics within this cohort (Figure 2), these multiple relapses in a short period of time illustrated the step-by-step replacement of the dominant clonal population to the benefit of a fitter clonal progeny under selective pressures, such as chemotherapy. Finally, case 62 showed an intermediate phenotype with ∼50% of the mutations identified at diagnosis that persisted at relapse (supplemental Table 3). Four clusters were identified (Figure 4E-F) and an accumulation of newly acquired mutations was observed at relapse. Although modest in comparison with other patients of the early relapse group, the clonal evolution pattern of patient 62 was still clearly distinguishable from the late relapse cases.

Figure 4.

Clonal architectures and ancestral relationships of early and late relapses. Clusters were obtained using Sciclone (variational Bayesian mixture model)34  for patients 34 (A), 325 (B), and 62 (C) based on VAFs at PT (x-axis) and R1 (y-axis) or R(n) and R(n + 1). Each dot stands for a mutation (SNVs or indels, WES coverage ≥40×). Different clusters are depicted by different colors and dot shapes (6, 9, and 4 clusters identified for case 34, 325, and 62, respectively). Mutations of interest are annotated. Case 34: The ancestral clone (clone 1), predominant at PT, contained the driver mutation NRAS p.G12D. Clone 2, descendant of clone 1, harbored a nonsense mutation (p.W2006*) in the tumor suppressor gene MLL2/KMT2D and became predominant at relapse (VAF_PT = 0.31, VAF_R1 = 0.55). A low-frequency subclonal mutation (VAF ≤ 0.05) in the oncogene KRAS was detected at diagnosis and rose again at late relapse with a similar frequency (VAF_PT = 0.03, VAF_R1 = 0.05). Case 325: The PT architecture showed a very close structure to that observed at R1 with a persistent ancestral subclone (clone 4) and a similar distribution of subclones. Clone 4 harbored MAPK pathway–activating mutation p.A72D in the SH2 domain of the oncogene PTPN11 subclonal at diagnosis (VAF = 0.23), selected at R1 (VAF = 0.49), and dominant at the 3 subsequent time points (VAF = 0.52, 0.51, and 0.47 in R2, R3, and R4, respectively). Subclonal mutations in driver genes were counter-selected after treatment: KRAS p.K117N (VAF-PT = 0.16 and VAF-R1 = 0.003) and SETD2 p.R456* (VAF-PT = 0.24 and VAF-R1 = 0). A missense mutation in the glucocorticoid receptor NR3C1 (p.P626S) emerged at R1 (VAF = 0.22) and expanded at R2 to become predominant (VAF = 0.42). NR3C1 p.P626S was counter-selected at R3 (VAF = 0) during the most important changes in this tumor history and was replaced by 2 co-occurring mutations in NR3C1 (p.Y641H, VAF = 0.50 and p.D724E, VAF = 0.53) limited to the newly dominant clone 5. This latter clone, descendant of ancestral clone 4, acquired a homozygous loss-of-function mutation p.R134* in the DNA repair gene PMS2 and carried several mutated driver genes: CREBBP (p.Y1503H), IKZF1 (p.F173L), MLL3/KMT2C (p.R2609*), and SRSF2 (p.G12S). This led to the expansion of this clone at R3, which constituted the predominant population at both R3 and R4. Multiple minor subclones, descendant of clone 5 and harboring mutations in TP53 (p.R273H, p.R273C and p.Y236C), MSH6 (p.F573L), and WHSC1 (p.A457T), also emerged at R3 and were maintained as minor clones at R4. Case 62: The dominant population at relapse, probably pulled up by the driver mutation KRAS p.G12D (VAF-PT = 0.39, VAF-R1 = 0.55), harbored a series of newly acquired mutations partly belonging to cluster 3. A subclonal mutation (VAF-PT = 0.10), predicted to alter splicing in the HAT domain of the histone acetyl transferase CREBBP, was lost at relapse (VAF-R1 = 0). Clonal trees depicting the evolutionary history of tumors based on VAFs of somatic mutations predicted as putative drivers in PTs and/or the subsequent relapse(s) of cases 34 (D), 325 (E), and 62 (F). Ancestral relationships between clones were inferred using AncesTree35  under the infinite sites (perfect phylogeny) assumption and are represented by black solid lines. Posterior probabilities of the ancestral relationships are indicated on the side of the black lines. AncesTree uses a probabilistic model for the observed read counts: if Xpj and Xpk are random variables describing the VAFs of mutation j and k in the sample p (PT or the subsequent relapse[s]), Pr[XpjXpk] denote the posterior probability that XpjXpk, and the sample with the smallest probability minpPr[XpjXpk] represents the weakest evidence that mutation j preceded mutation k (if close to 1, j is likely to be ancestral to k). Dashed lines show ancestral clones that existed at the time of sequencing. Each sample is indicated in a colored box at the bottom of the trees (A, B, C, D, and E stand for PT, R1, R2, R3, and R4, respectively). Colored lines indicate the inferred composition of clones and their fraction in each sample (only edges with usage of at least 0.05 are shown).

Figure 4.

Clonal architectures and ancestral relationships of early and late relapses. Clusters were obtained using Sciclone (variational Bayesian mixture model)34  for patients 34 (A), 325 (B), and 62 (C) based on VAFs at PT (x-axis) and R1 (y-axis) or R(n) and R(n + 1). Each dot stands for a mutation (SNVs or indels, WES coverage ≥40×). Different clusters are depicted by different colors and dot shapes (6, 9, and 4 clusters identified for case 34, 325, and 62, respectively). Mutations of interest are annotated. Case 34: The ancestral clone (clone 1), predominant at PT, contained the driver mutation NRAS p.G12D. Clone 2, descendant of clone 1, harbored a nonsense mutation (p.W2006*) in the tumor suppressor gene MLL2/KMT2D and became predominant at relapse (VAF_PT = 0.31, VAF_R1 = 0.55). A low-frequency subclonal mutation (VAF ≤ 0.05) in the oncogene KRAS was detected at diagnosis and rose again at late relapse with a similar frequency (VAF_PT = 0.03, VAF_R1 = 0.05). Case 325: The PT architecture showed a very close structure to that observed at R1 with a persistent ancestral subclone (clone 4) and a similar distribution of subclones. Clone 4 harbored MAPK pathway–activating mutation p.A72D in the SH2 domain of the oncogene PTPN11 subclonal at diagnosis (VAF = 0.23), selected at R1 (VAF = 0.49), and dominant at the 3 subsequent time points (VAF = 0.52, 0.51, and 0.47 in R2, R3, and R4, respectively). Subclonal mutations in driver genes were counter-selected after treatment: KRAS p.K117N (VAF-PT = 0.16 and VAF-R1 = 0.003) and SETD2 p.R456* (VAF-PT = 0.24 and VAF-R1 = 0). A missense mutation in the glucocorticoid receptor NR3C1 (p.P626S) emerged at R1 (VAF = 0.22) and expanded at R2 to become predominant (VAF = 0.42). NR3C1 p.P626S was counter-selected at R3 (VAF = 0) during the most important changes in this tumor history and was replaced by 2 co-occurring mutations in NR3C1 (p.Y641H, VAF = 0.50 and p.D724E, VAF = 0.53) limited to the newly dominant clone 5. This latter clone, descendant of ancestral clone 4, acquired a homozygous loss-of-function mutation p.R134* in the DNA repair gene PMS2 and carried several mutated driver genes: CREBBP (p.Y1503H), IKZF1 (p.F173L), MLL3/KMT2C (p.R2609*), and SRSF2 (p.G12S). This led to the expansion of this clone at R3, which constituted the predominant population at both R3 and R4. Multiple minor subclones, descendant of clone 5 and harboring mutations in TP53 (p.R273H, p.R273C and p.Y236C), MSH6 (p.F573L), and WHSC1 (p.A457T), also emerged at R3 and were maintained as minor clones at R4. Case 62: The dominant population at relapse, probably pulled up by the driver mutation KRAS p.G12D (VAF-PT = 0.39, VAF-R1 = 0.55), harbored a series of newly acquired mutations partly belonging to cluster 3. A subclonal mutation (VAF-PT = 0.10), predicted to alter splicing in the HAT domain of the histone acetyl transferase CREBBP, was lost at relapse (VAF-R1 = 0). Clonal trees depicting the evolutionary history of tumors based on VAFs of somatic mutations predicted as putative drivers in PTs and/or the subsequent relapse(s) of cases 34 (D), 325 (E), and 62 (F). Ancestral relationships between clones were inferred using AncesTree35  under the infinite sites (perfect phylogeny) assumption and are represented by black solid lines. Posterior probabilities of the ancestral relationships are indicated on the side of the black lines. AncesTree uses a probabilistic model for the observed read counts: if Xpj and Xpk are random variables describing the VAFs of mutation j and k in the sample p (PT or the subsequent relapse[s]), Pr[XpjXpk] denote the posterior probability that XpjXpk, and the sample with the smallest probability minpPr[XpjXpk] represents the weakest evidence that mutation j preceded mutation k (if close to 1, j is likely to be ancestral to k). Dashed lines show ancestral clones that existed at the time of sequencing. Each sample is indicated in a colored box at the bottom of the trees (A, B, C, D, and E stand for PT, R1, R2, R3, and R4, respectively). Colored lines indicate the inferred composition of clones and their fraction in each sample (only edges with usage of at least 0.05 are shown).

Discussion

Although the observed effects require validations in larger cohorts to generalize our findings, to our knowledge, this is the first study capturing 2 clonal dynamics governing late and early relapses through the study of somatic allele frequencies. However, some limitations have to be noted: first, although the mean coverage on the targeted region of WES reached 200×, bulk cell sequencing is limiting and excludes an exhaustive characterization of the subclonal architecture. Secondly, RNA sequencing was not performed, which did not allow for the full characterization of cases and a comprehensive analysis of early and late relapse-specific expression profiles.

Five of the 6 pathways identified as frequently targeted in our cALL cohort (epigenetic regulation, hemopoiesis–immune system development, MAP kinase–RAS signaling, DNA repair, and cell cycle) confirmed their relapse-driving potential.13  The JAK-STAT/hemopoiesis SH2B3 gene was found mutated in 2 cases (Figure 1; supplemental Table 1). However, we did not observe the recently reported enrichment of mutations in the JAK-STAT pathway in relapsed cALL patients.13  Alterations of this pathway are common genomic features of Ph-like acute lymphoblastic leukemia (ALL),36  and although we missed expression data allowing a classification of cases in this particular subgroup, we suspect the observed divergence to rely on a lack of representatives in our cohort. With the rise and fall of numerous mutated clones, epigenetic regulation and MAP kinase–RAS signaling were the 2 most represented (84.2% and 78.9% of mutated patients, respectively) and dynamic signaling pathways. They perfectly illustrated the clonal competition that occurs under external selective pressures, such as chemotherapy. Of particular interest, we also identified a mutational turnover of the glucocorticoid receptor NR3C1 with a succession and a co-occurrence of missense mutations (p.P626S, p.Y641H, and p.D724E) in the dominant clones of case 325. Despite its limitation to a single case in this cohort, the identified reduction of transactivation by glucocorticoids due to a mutation of the residue P626 in CV-1 cells37  suggests a possible involvement in resistance to therapy and thus could explain the selection of mutations. Further suggesting a role of NR3C1 in glucocorticoids resistance and relapse, a very recent study characterized a novel transcriptional program in B-lymphoid cells controlling malignant transformation through the repression of glucose and energy supply and identified the product of NR3C1 as a central effector of this restriction.38 

We characterized 2 general evolution patterns: a highly dynamic/adaptive clonal pattern specific to early relapses, illustrating the quick emergence of fitter clones and the eradication of other subclones, and a slow-to-quasi-inert evolution pattern associated with late events. By quantifying both the general clonal dynamics and the newly acquired subclonal diversity, we captured the underlying differences of resistance mechanisms between early and late pre–B-cell ALL relapse. In the most plastic tumors, the accumulation of somatic mutations allows rapid subclonal diversification and increases the chance of acquisition of resistant mutations. It leads to a postchemotherapy turnover of predominant populations and to the selection of the fitter clone. Interestingly, we found an overrepresentation of cases from this early relapse group with mutation in DNA repair genes (62.5% vs 18.2% in the late relapse group). Furthermore, most early relapses presented dominant mutations in this pathway. Although this observation requires validations in larger cohorts, a defect in the DNA repair process might explain the fast evolutionary adaptability of early relapse cases. On the other hand, the limited occurrence of hypermutator phenotypes (26% of the cohort) suggests that a strong accumulation of somatic mutations is not a prerequisite mechanism for the acquisition of resistance, but obviously increases the odds of acquisition of advantageous mutations and leads to a higher tumor plasticity, likely giving rise to resistant clones. In regards to these hypermutable patients, 4 of them (2 early and 2 late relapsers) underwent a BMT and the associated aggressive chemotherapy regimens, including high doses of alkylating agents, such as busulfan or cyclophosphamide. In other cancer types, somatic mutations in the DNA repair pathway, such as the ones observed here for hypermutable patients, were suspected to confer resistance to alkylator therapy and to generate hypermutations because of strong exposure to mutagenic agents. 25-27  Similar to our observations, these cases also presented an excess of transition mutations. Importantly, given that as many early as late relapsers were affected, this particularity is unlikely to have contributed to the captured difference of clonal evolution. As for case 325, who presented a hypermutator phenotype at his R3, which was not preceded by a transplant, the loss of function of the DNA repair gene PMS2 combined with repetitive treatments were likely responsible for this phenotype. This confirmed the recently identified hypermutator potential of the somatic loss of function of the mismatch repair endonuclease PMS213, that, until recently, was mostly associated with Lynch syndrome.39  Although 3 of the 5 hypermutables cases were the only patients of the cohort identified as carriers of germ line mutations in MSH6 and PMS2, these polymorphisms alone can’t be considered responsible for the phenotype. This is further illustrated by the absence of a significant difference in the number of mutations harbored at the PT by these cases compared with the nonhypermutable ones. However, as pointed out previously, an association between the pretransplant intensive treatment regimen and a defect in DNA repair mechanisms, caused either by somatic mutations or, less likely, by germ line variants, could have led to this phenotype. Although further investigation is required to substantiate the role of these variants, in this context, a predisposition for a destabilization of the DNA repair functions should not be disregarded as a contributing factor.

As for the significantly reduced mutational dynamics of late events, it suggests a different mechanism of relapse than early ones, less likely relying on mutation-acquired chemoresistance. Until recently, data that attributed late recurrences to the emergence of clones derived from the original population at the time of diagnosis were generally limited to the use of rearrangements of immunoglobulin genes (IGH/IGK) as specific markers of clonal populations.40-42  Although little is known regarding cell-of-origin and the process of long-term cALL relapses, recent studies failed to identify shared mutations between diagnosis and late relapse, except for the most ancestral events (eg, fusion BCR-ABL1).43  This led to the suggestion that late relapse cells likely derive from stem cells of a minor clone at diagnosis.43,44  Here, we showed that all late events harbored somatic alterations that are shared with the PT. Some of these cases presented no variation of allele frequency between primary and relapse leukemias for most of the identified mutations, illustrating the conservation of the same clonal architectures. This provides evidence of a possible re-emergence of the predominant population at the time of diagnosis. The absence of clonal evolution implies that the leukemic cells have not undergone any biological changes during their prolonged dormancy. If so, it is likely that the reappearance of leukemia has been caused by a relaxation of the host's immune surveillance or an environmental change allowing the reactivation/re-expansion of dormant cells spared from the initial cytoreductive therapy41  rather than a genetic evolution enabling an escape from immune surveillance, as previously suggested.43  This is a clinical asset because the leukemia retains its chemosensitivity. The treatment quickly leads to a new remission, as observed in our cohort. Altogether, these characteristics suggested that late events could originate from a subpopulation of relapse-inducing cells exhibiting properties of long-term dormancy and stemness, as recently isolated from pediatric and adult patients with minimal residual disease (MRD).45  These cells also displayed in vivo drug resistance due to limited proliferation when located in their niche and were sensitive when dissociated from this protective environment.45 

Overall, if confirmed in other cohorts, our results have direct clinical implications considering that rapid subclonal expansion and high tumor plasticity are key determinants of rapid cALL evolution and predictive factors of early therapeutic escapes. This is in line with previous findings associating shorter remissions with the early presence of subclonal driver mutations in leukemia samples7,46  or with MRD levels on day 19 of remission induction.47  In this context, appropriate clinical use of novel genomic technologies could be an asset and, as previously suggested,13  would allow a systematic postdiagnosis driver mutation–oriented survey of MRD to detect early re-emergence of cancer cells. As for the prevention of long-term relapses that are less likely relying on mutation-based chemoresistance, patients would benefit from therapeutic strategies that release MRD cells from the niche.

The full-text version of this article contains a data supplement.

Acknowledgments

The authors thank the patients and their parents for participating in this study. Patient tissue samples were provided by the CHU Sainte-Justine Pediatric Cancer biobank. Whole-exome sequencing was performed at the Integrated Clinical Genomics Centre In Pediatrics, CHU Sainte-Justine. Whole-genome sequencing was performed by Illumina, Inc. (San Diego, CA). Whole-genome genotyping was performed at the McGill University and Genome Quebec Innovation Centre. Computations were made on the supercomputer Briarée from the University of Montreal, managed by Calcul Québec and Compute Canada.

This work was supported by research funds provided by the Canadian Institutes for Health Research. D.S. holds the François-Karl-Viau Research Chair in Pediatric Oncogenomics. The operation of the Calcul Québec supercomputer is funded by the Canada Foundation for Innovation, NanoQuébec, Réseau de médecine génétique appliquée, and the Fonds de Recherche du Québec–Nature et technologies.

Authorship

Contribution: D.S. is the principal investigator and takes primary responsibility for the paper; J.-F.S. and D.S. contributed to the conception and design of the study; P.C., C.R., and M.O. were involved in sample and library preparation for the whole-exome sequencing; J.-F.S. performed bioinformatics analysis, data integration, and analysis; J.-F.S. interpreted the data and wrote the paper; J.H. and D.S. were involved in critical revision of the manuscript; and all authors approved the final version.

Conflict-of-interest disclosure: The authors declare no competing financial interests.

Correspondence: Daniel Sinnett, Division of Hematology-Oncology, CHU Sainte-Justine Research Center, 3175 Côte Sainte-Catherine, Montreal, QC H3T 1C5, Canada; e-mail: daniel.sinnett@umontreal.ca.

References

References
1.
Institute of Medicine and National Research Council
.
Hewitt
M
,
Weiner
SL
,
Simone
JV
, eds.
Childhood Cancer Survivorship: Improving Care and Quality of Life
.
Washington, DC
:
National Academies Press
;
2003
.
2.
Siegel
DA
,
King
J
,
Tai
E
,
Buchanan
N
,
Ajani
UA
,
Li
J
.
Cancer incidence rates and trends among children and adolescents in the United States, 2001-2009
.
Pediatrics
.
2014
;
134
(
4
):
e945
-
e955
.
3.
Hunger
SP
,
Mullighan
CG
.
Acute lymphoblastic leukemia in Children
.
N Engl J Med
.
2015
;
373
(
16
):
1541
-
1552
.
4.
Smith
MA
,
Seibel
NL
,
Altekruse
SF
, et al
.
Outcomes for children and adolescents with cancer: challenges for the twenty-first century
.
J Clin Oncol
.
2010
;
28
(
15
):
2625
-
2634
.
5.
Linabery
AM
,
Ross
JA
.
Trends in childhood cancer incidence in the U.S. (1992-2004)
.
Cancer
.
2008
;
112
(
2
):
416
-
432
.
6.
Greaves
M
,
Maley
CC
.
Clonal evolution in cancer
.
Nature
.
2012
;
481
(
7381
):
306
-
313
.
7.
Landau
DA
,
Carter
SL
,
Stojanov
P
, et al
.
Evolution and impact of subclonal mutations in chronic lymphocytic leukemia
.
Cell
.
2013
;
152
(
4
):
714
-
726
.
8.
Irving
J
,
Matheson
E
,
Minto
L
, et al
.
Ras pathway mutations are prevalent in relapsed childhood acute lymphoblastic leukemia and confer sensitivity to MEK inhibition
.
Blood
.
2014
;
124
(
23
):
3420
-
3430
.
9.
Inthal
A
,
Zeitlhofer
P
,
Zeginigg
M
, et al
.
CREBBP HAT domain mutations prevail in relapse cases of high hyperdiploid childhood acute lymphoblastic leukemia
.
Leukemia
.
2012
;
26
(
8
):
1797
-
1803
.
10.
Mullighan
CG
,
Zhang
J
,
Kasper
LH
, et al
.
CREBBP mutations in relapsed acute lymphoblastic leukaemia
.
Nature
.
2011
;
471
(
7337
):
235
-
239
.
11.
Meyer
JA
,
Wang
J
,
Hogan
LE
, et al
.
Relapse-specific mutations in NT5C2 in childhood acute lymphoblastic leukemia
.
Nat Genet
.
2013
;
45
(
3
):
290
-
294
.
12.
Tzoneva
G
,
Perez-Garcia
A
,
Carpenter
Z
, et al
.
Activating mutations in the NT5C2 nucleotidase gene drive chemotherapy resistance in relapsed ALL
.
Nat Med
.
2013
;
19
(
3
):
368
-
371
.
13.
Ma
X
,
Edmonson
M
,
Yergeau
D
, et al
.
Rise and fall of subclones from diagnosis to relapse in pediatric B-acute lymphoblastic leukaemia
.
Nat Commun
.
2015
;
6
:
6604
.
14.
Landau
DA
,
Carter
SL
,
Getz
G
,
Wu
CJ
.
Clonal evolution in hematological malignancies and therapeutic implications
.
Leukemia
.
2014
;
28
(
1
):
34
-
43
.
15.
Green
MR
,
Gentles
AJ
,
Nair
RV
, et al
.
Hierarchy in somatic mutations arising during genomic evolution and progression of follicular lymphoma
.
Blood
.
2013
;
121
(
9
):
1604
-
1611
.
16.
Welch
JS
,
Ley
TJ
,
Link
DC
, et al
.
The origin and evolution of mutations in acute myeloid leukemia
.
Cell
.
2012
;
150
(
2
):
264
-
278
.
17.
Walter
MJ
,
Shen
D
,
Ding
L
, et al
.
Clonal architecture of secondary acute myeloid leukemia
.
N Engl J Med
.
2012
;
366
(
12
):
1090
-
1098
.
18.
Healy
J
,
Bélanger
H
,
Beaulieu
P
,
Larivière
M
,
Labuda
D
,
Sinnett
D
.
Promoter SNPs in G1/S checkpoint regulators and their impact on the susceptibility to childhood leukemia
.
Blood
.
2007
;
109
(
2
):
683
-
692
.
19.
Langmead
B
,
Salzberg
SL
.
Fast gapped-read alignment with Bowtie 2
.
Nat Methods
.
2012
;
9
(
4
):
357
-
359
.
20.
Saunders
CT
,
Wong
WS
,
Swamy
S
,
Becq
J
,
Murray
LJ
,
Cheetham
RK
.
Strelka: accurate somatic small-variant calling from sequenced tumor-normal sample pairs
.
Bioinformatics
.
2012
;
28
(
14
):
1811
-
1817
.
21.
Spinella
JF
,
Mehanna
P
,
Vidal
R
, et al
.
SNooPer: a machine learning-based method for somatic variant identification from low-pass next-generation sequencing
.
BMC Genomics
.
2016
;
17
(
1
):
912
.
22.
Koboldt
DC
,
Zhang
Q
,
Larson
DE
, et al
.
VarScan 2: somatic mutation and copy number alteration discovery in cancer by exome sequencing
.
Genome Res
.
2012
;
22
(
3
):
568
-
576
.
23.
Spinella
JF
,
Cassart
P
,
Richer
C
, et al
.
Genomic characterization of pediatric T-cell acute lymphoblastic leukemia reveals novel recurrent driver mutations
.
Oncotarget
.
2016
;
7
(
40
):
65485
-
65503
.
24.
Mullighan
CG
,
Phillips
LA
,
Su
X
, et al
.
Genomic analysis of the clonal origins of relapsed acute lymphoblastic leukemia
.
Science
.
2008
;
322
(
5906
):
1377
-
1380
.
25.
Alexandrov
LB
,
Nik-Zainal
S
,
Wedge
DC
, et al
;
ICGC PedBrain
.
Signatures of mutational processes in human cancer [published correction appears in Nature. 2013. 502(7470):258]
.
Nature
.
2013
;
500
(
7463
):
415
-
421
.
26.
Hunter
C
,
Smith
R
,
Cahill
DP
, et al
.
A hypermutation phenotype and somatic MSH6 mutations in recurrent human malignant gliomas after alkylator chemotherapy
.
Cancer Res
.
2006
;
66
(
8
):
3987
-
3991
.
27.
Tomita-Mitchell
A
,
Kat
AG
,
Marcelino
LA
,
Li-Sucholeiki
XC
,
Goodluck-Griffith
J
,
Thilly
WG
.
Mismatch repair deficient human cells: spontaneous and MNNG-induced mutational spectra in the HPRT gene
.
Mutat Res
.
2000
;
450
(
1-2
):
125
-
138
.
28.
Miyaki
M
,
Konishi
M
,
Tanaka
K
, et al
.
Germline mutation of MSH6 as the cause of hereditary nonpolyposis colorectal cancer
.
Nat Genet
.
1997
;
17
(
3
):
271
-
272
.
29.
Kolodner
RD
,
Tytell
JD
,
Schmeits
JL
, et al
.
Germ-line msh6 mutations in colorectal cancer families
.
Cancer Res
.
1999
;
59
(
20
):
5068
-
5074
.
30.
Wu
Y
,
Berends
MJ
,
Mensink
RG
, et al
.
Association of hereditary nonpolyposis colorectal cancer-related tumors displaying low microsatellite instability with MSH6 germline mutations
.
Am J Hum Genet
.
1999
;
65
(
5
):
1291
-
1298
.
31.
Li
Z
,
Kong
L
,
Yu
L
, et al
.
Association between MSH6 G39E polymorphism and cancer susceptibility: a meta-analysis of 7,046 cases and 34,554 controls
.
Tumour Biol
.
2014
;
35
(
6
):
6029
-
6037
.
32.
Marinovic-Terzic
I
,
Yoshioka-Yamashita
A
,
Shimodaira
H
, et al
.
Apoptotic function of human PMS2 compromised by the nonsynonymous single-nucleotide polymorphic variant R20Q
.
Proc Natl Acad Sci USA
.
2008
;
105
(
37
):
13993
-
13998
.
33.
McGranahan
N
,
Swanton
C
.
Clonal heterogeneity and tumor evolution: past, present, and the future
.
Cell
.
2017
;
168
(
4
):
613
-
628
.
34.
Miller
CA
,
White
BS
,
Dees
ND
, et al
.
SciClone: inferring clonal architecture and tracking the spatial and temporal patterns of tumor evolution
.
PLOS Comput Biol
.
2014
;
10
(
8
):
e1003665
.
35.
El-Kebir
M
,
Oesper
L
,
Acheson-Field
H
,
Raphael
BJ
.
Reconstruction of clonal trees and tumor composition from multi-sample sequencing data
.
Bioinformatics
.
2015
;
31
(
12
):
i62
-
i70
.
36.
Iacobucci
I
,
Mullighan
CG
.
Genetic basis of acute lymphoblastic leukemia
.
J Clin Oncol
.
2017
;
35
(
9
):
975
-
983
.
37.
Bledsoe
RK
,
Montana
VG
,
Stanley
TB
, et al
.
Crystal structure of the glucocorticoid receptor ligand binding domain reveals a novel mode of receptor dimerization and coactivator recognition
.
Cell
.
2002
;
110
(
1
):
93
-
105
.
38.
Chan
LN
,
Chen
Z
,
Braas
D
, et al
.
Metabolic gatekeeper function of B-lymphoid transcription factors
.
Nature
.
2017
;
542
(
7642
):
479
-
483
.
39.
Senter
L
,
Clendenning
M
,
Sotamaa
K
, et al
.
The clinical phenotype of Lynch syndrome due to germ-line PMS2 mutations
.
Gastroenterology
.
2008
;
135
(
2
):
419
-
428
.
40.
Vora
A
,
Frost
L
,
Goodeve
A
, et al
.
Late relapsing childhood lymphoblastic leukemia
.
Blood
.
1998
;
92
(
7
):
2334
-
2337
.
41.
Frost
L
,
Goodeve
A
,
Wilson
G
,
Peake
I
,
Barker
H
,
Vora
A
.
Clonal stability in late-relapsing childhood lymphoblastic leukaemia
.
Br J Haematol
.
1997
;
98
(
4
):
992
-
994
.
42.
Levasseur
M
,
Maung
ZT
,
Jackson
GH
,
Kernahan
J
,
Proctor
SJ
,
Middleton
PG
.
Relapse of acute lymphoblastic leukaemia 14 years after presentation: use of molecular techniques to confirm true re-emergence
.
Br J Haematol
.
1994
;
87
(
2
):
437
-
438
.
43.
Ford
AM
,
Mansur
MB
,
Furness
CL
, et al
.
Protracted dormancy of pre-leukemic stem cells
.
Leukemia
.
2015
;
29
(
11
):
2202
-
2207
.
44.
Ford
AM
,
Fasching
K
,
Panzer-Grümayer
ER
,
Koenig
M
,
Haas
OA
,
Greaves
MF
.
Origins of “late” relapse in childhood acute lymphoblastic leukemia with TEL-AML1 fusion genes
.
Blood
.
2001
;
98
(
3
):
558
-
564
.
45.
Ebinger
S
,
Özdemir
EZ
,
Ziegenhain
C
, et al
.
Characterization of rare, dormant, and therapy-resistant cells in acute lymphoblastic leukemia
.
Cancer Cell
.
2016
;
30
(
6
):
849
-
862
.
46.
Barrio
S
,
Shanafelt
TD
,
Ojha
J
, et al
.
Genomic characterization of high-count MBL cases indicates that early detection of driver mutations and subclonal expansion are predictors of adverse clinical outcome
.
Leukemia
.
2017
;
31
(
1
):
170
-
176
.
47.
Pui
CH
,
Pei
D
,
Coustan-Smith
E
, et al
.
Clinical utility of sequential minimal residual disease measurements in the context of risk-based therapy in childhood acute lymphoblastic leukaemia: a prospective study
.
Lancet Oncol
.
2015
;
16
(
4
):
465
-
474
.

Author notes

The datasets supporting the conclusions of this article have been deposited in the NCBI BioProject data repository (accession number PRJNA393314).

Supplemental data