Monitoring of Leukemia Clones in B-cell Acute Lymphoblastic Leukemia at Diagnosis and During Treatment by Single-cell DNA Amplicon Sequencing.

Acute lymphoblastic leukemia (ALL) is characterized by the presence of chromosomal changes, including numerical changes, translocations, and deletions, which are often associated with additional single-nucleotide mutations. In this study, we used single cell-targeted DNA sequencing to evaluate the clonal heterogeneity of B-ALL at diagnosis and during chemotherapy treatment. We designed a custom DNA amplicon library targeting mutational hotspot regions (in 110 genes) present in ALL, and we measured the presence of mutations and small insertions/deletions (indels) in bone marrow or blood samples from 12 B-ALL patients, with a median of 7973 cells per sample. Nine of the 12 cases showed at least 1 subclonal mutation, of which cases with PAX5 alterations or high hyperdiploidy (with intermediate to good prognosis) showed a high number of subclones (1 to 7) at diagnosis, defined by a variety of mutations in the JAK/STAT, RAS, or FLT3 signaling pathways. Cases with RAS pathway mutations had multiple mutations in FLT3, NRAS, KRAS, or BRAF in various clones. For those cases where we detected multiple mutational clones at diagnosis, we also studied blood samples during the first weeks of chemotherapy treatment. The leukemia clones disappeared during treatment with various kinetics, and few cells with mutations were easily detectable, even at low frequency (<0.1%). Our data illustrate that about half of the B-ALL cases show >2 subclones at diagnosis and that even very rare mutant cells can be detected at diagnosis or during treatment by single cell-targeted DNA sequencing.

INTRODUCTION

Acute lymphoblastic leukemia (ALL) is a hematological malignancy arising from malignant transformation of immature B or T lymphoid progenitor cells. It is the most common type of childhood cancer and has a higher incidence in pediatric patients compared with adults. Of all diagnosed pediatric ALL cases, over 80% is of the B-cell lineage (B-ALL).[1]

B-ALL is a genetically heterogeneous disease caused by the accumulation of chromosomal rearrangements (mostly deletions and translocations) and mutations (point mutations and insertions/deletions). Recurrent genetic translocations in B-ALL result in fusion genes, such as ETV6-RUNX1,[2,3] TCF3-PBX1,[3,4] TCF3-HLF,[3,4] or BCR-ABL1.[3] ETV6-RUNX1-like[5] and BCR-ABL1-like[6] cases are B-ALL subtypes that have similar transcriptional patterns as ETV6-RUNX1 or BCR-ABL1 cases but do not carry the translocations in question. Furthermore, KMT2A (MLL),[2-4] DUX4,[3] ZNF384,[2,3] MEF2D,[2,3] and CRLF2[7] rearrangements have been described in B-ALL,[6-9] as well as intrachromosomal amplifications of chromosome 21 (iAMP21).[10-12] Moreover, B-ALL can also be characterized by hypodiploidy (<44 chromosomes)[3,9] or high hyperdiploidy (>51 chromosomes),[3] often also involving trisomy of chromosomes 4 and 10.[13]

Aside from larger chromosomal alterations, B-ALL cases can carry smaller indels or somatic point mutations. Frequently observed mutations are alterations in the cell cycle regulators TP53 and CDKN2A/B, or transcriptional regulators involved in B-cell development, such as PAX5, EBF1, or IKZF1.[4,7] PAX5 is the most common target of genetic alterations in B-ALL, with deletions, amplifications, point mutations, and gene fusions in almost 40% of ALL cases.[4,14] Furthermore, RAS signaling mutations in NRAS, KRAS, FLT3, or NF1, and activating mutations in the JAK/STAT signaling pathway are enriched in certain B-ALL subtypes.[15-17] In particular, RAS and JAK/STAT signaling mutations are present in the majority of PAX5 altered cases.[4] Additionally, mutations can occur in the epigenetic regulators SETD2, KMT2D, or CREBBP.[18]

Sequencing studies have provided a wealth of information on the different genomic lesions present in ALL, but still little is known on the order of acquisition of the various mutations and heterogeneity within the leukemic cell population. In the early era of single-cell genomics, low throughput targeted single-cell DNA sequencing was used to study ALL clonality and was able to provide the first data on heterogeneity and mutational acquisition.[19,20] Although these methods were similar to the one used in this paper, they only allowed a limited number of cells to be analyzed. In a previous study,[21] we used high-throughput single-cell DNA amplicon sequencing to dissect the clonal heterogeneity in T-ALL samples at diagnosis and during treatment, based on a custom amplicon library covering mutational hotspot regions in ALL. These data revealed the presence of minor leukemic subclones and evaluated their therapy response in T-ALL.

Here, we applied the same custom ALL panel, consisting of 305 genomic regions across 110 frequently mutated genes in either T- or B-ALL, to a set of 12 B-ALL samples at diagnosis and during treatment to study heterogeneity and clonality in B-cell leukemia. We performed single-cell DNA sequencing with the Tapestri platform (Mission Bio)[22] to genotype over 7000 single cells per sample and screened for single-nucleotide variants (SNVs) and indels. The Tapestri platform uses droplet-based microfluidics technology to isolate single cells and allows high-throughput targeted sequencing of a set of genomic regions of interest. This approach allowed us to study the population of leukemic cells at the single-cell level and delineate the order of acquisition and co-occurrence of subclonal mutations, providing insight in the heterogeneity of leukemia at diagnosis and the evolution of the clonality during chemotherapy treatment. We combined this approach with fluorescent in situ hybridization (FISH), multiplex ligation-dependent probe amplification (MLPA), karyotyping, optical genome mapping of structural variants with the Saphyr instrument (Bionano Genomics), and bulk-RNA sequencing to provide a comprehensive overview of the genetic and chromosomal variants present in this B-ALL cohort.

METHODS

Patients and samples

Bone marrow (BM) and peripheral blood (PB) samples obtained at diagnosis and during chemotherapy treatment were collected from children and adults diagnosed with B-ALL at the University Hospitals Leuven. This study was approved by the Ethical Committee of the University Hospitals Leuven. All pediatric patients were treated according to European Organisation for Research and Treatment of Cancer (EORTC) study protocol 58081 and the adult patient (XG118) according to HOVON100. Written informed consent was obtained from every patient in accordance with the Declaration of Helsinki. Mononuclear cells were isolated from fresh BM and PB samples using Ficoll-Paque and viably frozen in 95% fetal bovine serum and 5% dimethyl sulfoxide.

Single-cell DNA sequencing (Tapestri platform, Mission Bio)

Cryopreserved BM and PB samples were thawed, washed 3 times with phosphate-buffered saline with 2% fetal bovine serum, filtered with a 40 µm cell strainer, and counted with the LUNA-FL Automated Fluorescence Cell Counter (Logos Biosystems) or TC20 Automated cell counter (Bio-Rad). Cell suspensions with concentrations of approximately 4000 cells/µL were prepared, with cell viabilities of >75%. A total of 50 µL containing on average 187,000 ± 7000 (SEM) cells were loaded on the Tapestri platform (Mission Bio) for cell encapsulation with lysis buffer and proteolytic enzymes. Encapsulated cells were lysed and DNA binding proteins were enzymatically digested by incubating the samples for 60 minutes at 50°C and 10 minutes at 80°C. Next, the droplets containing encapsulated cell lysates were combined with droplets containing barcoding beads. These beads are coupled to oligonucleotides carrying the cell barcode and amplicon-specific primers for targeted amplification. For the multiplex polymerase chain reaction (PCR) reaction, we used a custom ALL panel kit,[21] which targets variant hotspots across 110 genes with a total of 305 amplicons of approximately 192 bp insert length. The panel was custom-designed based on Catalogue Of Somatic Mutations In Cancer and pediatric cancer databases[23,24] and covers >1300 reported mutations in ALL (SNV and small indels).

Before the multiplex PCR, cell emulsions were exposed to UV light to cleave off the oligonucleotides from the barcoding beads. Next, the amplicons were PCR amplified under the following conditions: 10 minutes at 95°C; 20 cycles of 30 seconds at 95°C, 10 seconds at 72°C, 10 minutes 30 seconds at 61°C, and 30 seconds at 72°C; and a final step of 2 minutes at 72°C followed by overnight holding at 4°C. The emulsions containing the barcoded PCR products were then broken using extraction agent, containing perfluoro-1-octanol, releasing the barcoded DNA into the aqueous phase. AMPure XP beads (Beckman Coulter) were used for PCR cleanup, and the purified DNA was quantified with the dsDNA High Sensitivity Qubit Assay (ThermoFisher Scientific). Sample indices and Illumina adapter sequences were subsequently incorporated via a 10-cycle library PCR reaction with following conditions: 3 minutes at 95°C; 10 cycles of 20 seconds at 98°C, 20 seconds at 62°C, and 45 seconds at 72°C; and a final step of 2 minutes at 72°C followed by holding at 4°C. The resulting libraries (2 per sample for V1 chemistry) were purified with a second AMPure XP cleanup, analyzed with a High Sensitivity DNA kit on a Bioanalyzer (Agilent Technologies), and normalized to 5 nM. The samples were sequenced on a NovaSeq 6000 instrument (Illumina), using the S2 or S4 flowcell and Novaseq XP workflow for a 150 bp paired-end run, with a 15% ratio of PhiX DNA.

Data analysis: preprocessing and variant filtering

FastQ files were processed using Tapestri Pipeline (Tapestri 1.8.6) for adapter trimming (Cutadapt), sequence alignment (reference genome hg19), barcode correction, cell finding and variant calling (Genome Analysis Toolkit). Loom files generated by the Tapestri pipeline were loaded into Tapestri Insights v2.2 software for prefiltering and exploratory analysis.

Cells and variants were filtered in Tapestri Insights software according to the following 6 parameters: genotype quality score <30, read depth <10 reads, single-cell variant allele frequency (scVAF) <20%, variants genotyped in <60% of the cells, cells with <50% of the genotypes present, and variants mutated in <0.5% of cells. Variants complying any of these filter parameters were removed, and further filtering was performed with custom R scripts to exclude noisy variants resulting from technical artifacts. Variants annotated at the first or final 5 nucleotides of the amplicon were excluded from the analysis, as well as variants that were detected in 10 or more patients of this cohort as those are most likely artifacts. Furthermore, germline, intronic, and synonymous variants were removed, where germline variants were identified as described in Albertí-Servera et al.[21] The remaining variants were further filtered based on their quality metrics, and annotation errors or redundant annotations were corrected by examining the mapped sequencing reads using Integrative Genomics Viewer.

Clonal architecture and inference of mutational histories

Clonal architectures were initially determined by genotype clustering analysis with the Tapestri Insights software and visualized using custom pipelines in R. Only cells with complete genotypes were considered and the minimum clone size was 10 cells. Allelic dropout rates (ADO) were determined as described previously.[21] Since ADO can generate false-positive homozygous variant calling, we represented the clonal architectures without zygosity information.

To infer the phylogenetic trees of the somatic mutations at diagnosis, we used the infSCITE (single-cell inference of tumor evolution) software,[25] which performs a stochastic search to find the maximum likelihood or maximum a posteriori tree. In this model, we applied the infinite sites assumption and did not allow for reoccurrences or losses of variants. Moreover, we omitted the zygosity information and removed all cells with uncertain variant calls. The infSCITE software was run in 3 repetitions with a chain length of 100 000. Furthermore, we assumed the false-positive rate to be 0.001 and used the following parameters: –fn 0.001 –ad 0.001 0.001 –cc 0.3. Further data analysis to count the number of cells per node in the graph was performed with a custom R script. Clonal evolution during chemotherapy treatment was visualized in fish plots generated with the Tapestri Insights software.

Data sharing statement

The single-cell DNA-seq data (fastQ and vcf files) have been deposited at the European Genome-phenome Archive (EGA), which is hosted by the European Bioinformatics Institute and the Centre for Genomic Regulation, under accession number EGAS00001005029.

RESULTS

Sample characteristics

In this study, we analyzed 23 samples (13 BM and 10 PB mononuclear cells) from 12 B-ALL patients (11 pediatric cases and 1 adult case). Clinical patient characteristics are described in Table . Genetic characterization at diagnosis was performed by karyotyping, FISH, MLPA, and (reverse transcription-)PCR. Furthermore, we used the Saphyr instrument (Bionano Genomics) for optical genome mapping of structural variants and performed bulk-RNA sequencing to detect fusion transcripts. Based on these data, the 12 cases were assigned to the major ALL subtypes (Table ; Figure ).[4] The cases were selected based on the availability of viably frozen cells and not based on specific B-ALL subtypes or risk profiles. Due to the cohort size of 12 patients, not all B-ALL subtypes were represented in this study.

Heatmap of genomic characteristics of all patient samples at diagnosis. Each column represents 1 of the 12 B-ALL patients and each row represents a specific genetic alteration. The upper part of the heatmap (orange gradient; number of variants per patient and gene) shows SNVs or indels identified with single-cell DNA sequencing, with the first column showing the number of variants in each gene and the top row (blue gradient) the number of variants per patient. The lower part of the heatmap (red gradient) displays the major chromosomal alterations that were found by molecular diagnostics, karyotyping, bulk RNA-sequencing, or Bionano Optical Genome Mapping. The top row (green gradient) shows the number of chromosomal alterations per patient.

ETV6-RUNX1 rearrangements were found in 2 cases (XF93 and XF99), a BCR-ABL1 fusion was detected in 1 case (XF104), and a PAX5-JAK2 fusion was detected in 1 case (XG125), identifying the latter as a Philadelphia-like case.[4] Furthermore, we observed high hyperdiploidy in XF106, XG111, and XG115. PAX5 fusion genes or mutations that are typical for the PAX5alt subgroup were detected in 4 cases (PAX5-ZNF318, PAX5-NOL4L, PAX5-ZCCHC7, PAX5G183R). PAX5 was the most frequently altered gene, with 8 out of 12 cases (XF93, XF99, XF101, XF104, XF109, XG110, XG124, and XG125) showing PAX5 aberrations, either as deletions, point mutations or fusions with other genes. The adult patient XG118 did not harbor a clear aberration to assign it to any of the ALL subtypes, but single-cell sequencing later revealed the presence of the N159Y mutation in IKZF1. This indicates that this case could be classified as IKZF1N159Y, a recently defined rare subtype of ALL with intermediate prognosis.[8] In addition to these subtype-defining alterations, we also detected additional lesions, as shown in Table .

A custom amplicon panel to detect genetic variants at single-cell level in B-ALL patients

The microfluidics-based MissionBio Tapestri platform allows high-throughput single-cell DNA sequencing to detect SNVs and small indels across a set of predefined amplicons.[22] We used a custom ALL panel targeting 305 amplicons of 192 bp median length across 110 genes, covering mutational hotspot regions in ALL (T-ALL and B-ALL). We optimized this panel previously and validated it in a first cohort of 8 T-ALL cases.[21] In the current study, we used this amplicon panel to characterize the clonal architecture of 12 B-ALL cases at the single-cell level.

Eight cases were analyzed at diagnosis (T0) only (BM for XF99, XF101, XF104, XF106, XG110, and XG125; PB for XG118; BM and PB for XF93). The other 4 cases were analyzed both at diagnosis and during chemotherapy treatment (XF109, XG111, XG115, and XG124), with longitudinal samples collected during pretreatment with glucocorticoids (T1, first week), after the first dose of chemotherapy (T2, second week) and at clinical assessment of MRD (T3, 5 weeks after diagnosis). BM aspirations were performed on T0 and T3, for the other timepoints there were only PB samples available. We have previously shown that the mutational heterogeneity of BM and PB samples is comparable.[21] An overview of the samples sequenced per patient is provided in Suppl. Table S1.

A total of 176 280 cells were sequenced across all samples, with a median of 7973 cells per sample (interquartile range [IQR]: 4431–10 026) and a median coverage of 132X (IQR: 104X–153X) (Figure ). We noticed an uneven distribution of sequencing reads across the different amplicons, with 15 amplicons having insufficient coverage (average read depth across 23 samples under 10 reads) and being discarded from further analysis (Figure ; Suppl. Table S2). Similar findings concerning amplicon coverage were reported previously.[21] The estimated median ADO rate was 7.8% (IQR: 6.1–8.5) (Figure ; Suppl. Table S3).

Quality metrics and amplicon coverage of B-ALL patient samples based on single-cell DNA amplicon sequencing. (A) Total number of sequenced cells (left panel), number of reads per amplicon per cell (middle panel), and mean allele dropout rates (right panel). Details on sample and sequencing metrics are in Suppl. Table S1 and ADO calculations in Suppl. Table S3. Each dot is one of the 23 sequenced samples of the 12 B-ALL patients and the horizontal line represents the median value. (B) Bar plot of average amplicon coverage per cell per sample. Each bar represents 1 of the 305 amplicons, and the height of the bar represents the average number of reads per cell. The 15 lowest performing amplicons with an average coverage below 10 reads (represented by horizontal dashed red line) are listed on the right. Detailed data of amplicon coverage per sample can be found in Suppl. Table S2.

Single-cell DNA sequencing at diagnosis reveals subclonal mutations in a subset of cases

To investigate clonal architecture in the different subtypes of B-ALL, we analyzed the single-cell DNA amplicon sequencing data of the 12 diagnostic samples. Low-quality and low-frequency variants were excluded from the analysis by applying the 6 filtering parameters in Tapestri Insights software, resulting in a total of 239 different candidate gene variants, with a median of 73 (IQR: 61–77) variants per patient. All variants passing this initial filtering are listed in Suppl. Table S4. The prefiltering step was followed by exclusion of variants identified in 10 or more patients of this cohort and of falsely called variants located at the first or last 5 nucleotides of the sequence reads since these often contain artifacts.[21] A total of 137 variants were retained, with a median of 29 (IQR: 24–32) per patient. From the remaining variants, germline, synonymous, and intronic variants were further excluded from downstream analysis, where germline variants were identified as described in Alberti-Servera et al.[21] Finally, variants were removed that were genotyped in lower cell numbers or with large discrepancy between VAF by read count and VAF by cell count. This resulted in a total of 28 good-quality variants in 13 genes across the 12 B-ALL cases (Median per patient: 1.5, IQR: 0–4), as listed in Suppl. Table S4. We detected varying numbers of mutations, ranging from 0 variants in XF93, XF104, and XG125 up to 7 variants in XG111 and XG115 (Figure ).

We observed that the patients with an ETV6-RUNX1 rearrangement (XF93 and XF99) had a low mutational burden. In XF99, we found 1 splicing-related intronic SETD2 mutation in only 1% of the cells, while we did not find any significant variants in XF93. Similarly, the cases with the BCR-ABL1 fusion (XF104), PAX5-JAK2 fusion (XG125), or the IKZF1N159Y mutation (XG118) did not have any other relevant variants. From a mutational point of view, these 5 B-ALL cases were homogeneous, as we did not detect any relevant subclones with SNVs or indels in the genes that we studied.

In contrast to these 5 cases lacking clonal heterogeneity at mutation level, the B-ALL cases classified as high hyperdiploid (n = 3) showed a high mutational heterogeneity while in the cases with a PAX5alt (n = 4), we identified an intermediate number of relevant mutations.

High hyperdiploid cases (XF106, XG111, and XG115)

In 2 of the 3 cases with high hyperdiploidy, we detected several subclones with signaling mutations. In XG111, 2 point mutations were present in the NRAS gene (G12D and G13D). Furthermore, we found a BRAF mutation, as well as 2 PTPN11 variants (A72V and E76Q), a FLT3N676K and KMT2DQ3735* mutation. XG115 showed a high number of different subclones each with their own alterations, including variants in KRAS (G12D and G13D), NRAS (G12D and G12S), and FLT3 (I836del, A680V, and Y589-V592del). Each of these clones represented 0.6% to 17.9% of the cells, with the KRASG13D clone being the largest subclone. Two CREBBP (S1680del in 78% and L1681del in 1% of cells) and 2 OVGP1 variants (P496T in 2% and T495N in 2% of cells) were found in XF106.

PAX5alt cases (XF101, XF109, XG110, and XG124)

In the diagnostic sample of XF101, we found a large fraction of cells with a SETD2K2511Nfs*10 insertion (88% of the cells), while a small portion of cells had a EZH2R690L mutation (1% of the cells). In patient XF109, characterized by the presence of a PAX5G183R mutation in 84% of cells, we found 2 additional mutations in the JAK1 gene (L910P and Y652H), each in a different subclone. XG110 showed a minor CREBBP mutant clone accounting for 4% of the cells. In XG124, we detected 2 NRAS mutations (G12S and G12D), with the G12D mutation present in only a small subset of the cells with NRASG12S.

Clonal architecture and mutational history of B-ALL cases with multiple clones at diagnosis

To get more information on the clonal architecture and order at which mutations were acquired, we studied the 4 B-ALL cases in which we found mutations in the JAK/STAT, FLT3, and RAS signaling pathways. Of these 4 patients, 2 were high hyperdiploid (XG111 and XG115), 2 had PAX5 alterations (XF109 and XG124), and 3 had a deletion of CDKN2A/B (XF109, XG115, and XG124), often in combination with other rearrangements or deletions.

We observed multiple clones at diagnosis and used the InfSCITE algorithm to study the mutational history by inferring the most likely order of acquisition of SNVs and indels (Figure ). In XG111 and XG115, a branched pattern of clonality was predicted, with many minor clones emerging in parallel, whereas XF109 and XG124 showed a more sequential pattern with stepwise acquisition of a smaller number of somatic events (Figure ).

Clonal composition and mutational history based on single-cell sequencing data using the InfSCITE algoritm. (A) Heatmap of clonal architecture of diagnostic samples of patients XF109, XG111, XG115, and XG124. Because of potential allele dropout, zygosity information is not displayed (wild type [WT] = light gray, mutant [Mut] = dark gray). The heatmap shows the genotype consequence of each relevant somatic mutation and the colored horizontal bar at the top shows the clonal composition of the sample. (B) Phylogenetic trees showing the most likely order of variant acquisition during B-ALL development. Mutational history was inferred for diagnostic samples of patients XF109, XG111, XG115, and XG124 using the InfSCITE algorithm. Circles represent the clones and the outline color of the circles corresponds with the clone names in (A). The size of each circle represents the relative clone size in the sample, which is also indicated by the frequency next to the circles.

In XF109, the majority of the cells harbored a PAX5G183R mutation (84%), with most of these cells acquiring an additional mutation in JAK1. Remarkably, the cells either acquired a JAK1L910P or JAK1Y652H mutation, representing 46% and 13% of the total number of cells, respectively. Most of the leukemic cells in XG124 had a NRASG12S mutation (87%), with some cells (3%) acquiring an additional NRASG12D mutation.

In XG115, multiple mutually exclusive RAS signaling mutations in NRAS (G12D; G12S), KRAS (G12D; G13D), and FLT3 (I836del; A680V; Y589-V592del) were acquired in distinct small subclones, but the majority of cells did not harbor additional mutations in the genes that we investigated (labeled as C0 in Figures ). Thus, for FLT3, we detected 3 different mutations in 2.6% to 15.3% of the cells (Figure ), with one of these mutations difficult to detect in the Tapestri software, as this was a 12 nucleotide deletion in the juxtamembrane domain. To improve the detection of indels in the juxtamembrane region, we mapped all reads of this exon and inspected these for indels. In this way, we found a range of other in-frame deletions present in minor subclones, typically representing <0.5% of the cells (Figure ). The fact that so many cells could be detected with various in-frame deletions in the juxtamembrane region of FLT3 strongly suggests that there is high genomic instability in that region.

In-frame alterations in Overview of the alterations in the juxtamembrane domain of FLT3 in XG115. Alterations were identified based on inspection of sequencing reads in Integrative Genomics Viewer: yellow denotes deletions, purple denotes insertions, and green denotes mutations. All alterations identified were in-frame.

Clonal evolution during chemotherapy treatment. Fish plots showing clonal evolution during chemotherapy treatment in patients XF109, XG111, XG115, and XG124. At each timepoint, the total number of cells passing the filtering parameters is indicated. The somatic mutations (which are covered by this amplicon panel) per clone and the number of cells and frequencies are shown in Suppl. Table S5.

In patient XG111, variant PTPN11A72V defined a major clone and was present in 70% of cells, with 2 distinct subclonal mutations FLT3N767K and NRASG12D, alongside 4 minor clones with mutations in PTPN11, KMT2D, NRAS, or BRAF. Inspection of this case for the presence of additional FLT3 mutations in the juxtamembrane region did not identify additional mutations or indels.

Clonal evolution during chemotherapy treatment

Treatment of the pediatric B-ALL patients in this study consisted of intensive combination chemotherapy regimens according to the EORTC 58081 protocol. First, patients were pretreated with glucocorticoids (prednisone) for 1 week, followed by intensive intravenous and intrathecal chemotherapy during the 4-week induction phase. The induction phase consisted of glucocorticoids, vincristine, asparaginase, anthracycline (daunorubicine), and intrathecal chemotherapy, all intended to achieve complete remission defined by the morphological assessment of the BM (<5% blasts). This is followed by intensive chemotherapy during several months in the consolidation and reinduction/reintensification phase, and eventually followed by less intensive maintenance therapy.[1,26] To study the clonal evolution during induction treatment, we analyzed longitudinal samples of PB during glucocorticoid treatment (T1), after the first dose of chemotherapy (T2), and BM at the end of the induction phase (T3) (hematological assessment of minimal residual disease [MRD]). MRD was assessed by PCR for Ig rearrangement, which is an important prognostic marker and is used for risk stratification of the patients. All patients achieved complete remission, with 4 patients remaining MRD positive at day 35 by PCR who became MRD negative at day 90 (Table ).

We studied patients XF109, XG111, XG115, and XG124 at multiple timepoints during chemotherapy treatment and used single-cell sequencing to assess the presence and evolution of the leukemic clones. We analyzed 2 follow-up samples for XF109 and XG111 (T1 and T2) and 3 additional samples were sequenced for XG115 and XG124 (T1, T2, and T3). In all patients, we noticed a rapid decrease in the number of cells with mutations during chemotherapy treatment, in line with clinical observations.

The number of cells of clone 1 (NRASG12S) of XG124 did not decrease much after 3 days of treatment, but declined by T2 (only 2% mutated cells), and there were no residual mutated cells left in the 4670 cells of the BM sample that we analyzed at T3. In XF109, we noticed a slower decline of the leukemic clones with 21% of blood cells still containing the diagnostic mutations after 6 days of treatment, which decreased to 1% at T2 (8 cells with PAX5G183R, 40 cells with JAK1L910P/PAX5G183R, and 37 cells with JAK1Y652H/PAX5G183R remaining). Cells with the single PAX5G183R mutation disappeared faster than the double mutant cells with PAX5 and JAK1 mutation. In XG111, we observed a strong decrease in the number of mutated cells, with 3.5% of cells belonging to the major PTPN11A72V clone at T1, and only 11 residual mutated cells at T2.

Also in patient XG115, the leukemic burden was already greatly decreased by T1, with <2% cells left with detectable mutations. However, we detected that there were still mutant cells of 5 clones present at time of MRD measurement (T3): we detected 13 residual leukemia cells of the total 8886 cells analyzed based on the presence of the SNVs (Figure ): 6 cells with KRASG13D, 1 cell with FLT3A680V, 2 cells with KRASG12D, 2 cells with NRASG12D, and 2 cells with NRASG12S, indicating that many clones were still detectable, with some clones showing higher frequency compared with the diagnostic sample. These data show the sensitivity of single-cell DNA amplicon sequencing and its ability to detect the presence of small numbers of mutated cells during follow-up of the patient.

DISCUSSION

In recent years, single-cell sequencing technologies have emerged and are rapidly evolving. Droplet-based technologies offer the advantage that they can simultaneously analyze thousands of cells in a short amount of time, but this comes at the price of a lower sequencing depth. While this technology has mostly been applied to single-cell RNA sequencing, investigating the heterogeneous transcriptional profiles of complex mixtures, we used the droplet-based single-cell targeted DNA sequencing technique to characterize the mutational heterogeneity in B-ALL, both at diagnosis and during treatment.

We previously designed a custom amplicon panel to target those genes that are frequently mutated in ALL.[21] By using single-cell DNA amplicon sequencing, we aimed to study the clonal heterogeneity of ALL in more detail than can be achieved by whole genome sequencing. The technology used here can detect small subclones that represent <1% of the total cells, which would have been undetectable in routine whole genome sequencing analysis.

A first question we addressed is if clonal heterogeneity is observed in B-ALL at diagnosis and whether this is present in all cases or more restricted to specific subtypes of ALL. Compared with T-ALL, where the majority of cases accumulate additional mutations in subclones,[21] not all B-ALL cases show such pattern. From the 12 B-ALL cases that we analyzed in this study, 4 cases (1 ETV6-RUNX1, 1 BCR-ABL1 fusion, 1 BCR-ABL1-like, and 1 IKZF1N159Y subtype) did not harbor additional subclonal mutations and 1 additional ETV6-RUNX1 positive case had only 1 minor subclone (1% of the cells). These findings are consistent with other bulk sequencing studies that also detected a limited number of mutations in cases from these subtypes.[4,27] Genomic data suggest that patients in these subgroups often harbor deletions of typical B-ALL genes, such as PAX5, TCF3, or ETV6, but more rarely point mutations in signaling proteins or transcriptional regulators.[27] Single-cell sequencing of larger cohorts of B-ALL patients belonging to these subtypes is needed to validate our findings and might allow further fine tuning of these conclusions.

In contrast, cases from the PAX5alt or high hyperdiploid subtypes had a higher mutational burden with multiple subclones that could be identified by single-cell sequencing. We typically found a larger number of alterations in high hyperdiploid cases (4–7 mutations) compared to PAX5alt cases (1–3 mutations). In these 2 subgroups, we observed mutations in the FLT3, RAS, or JAK/STAT signaling pathways, and some of these only in a low number of cells. Indeed, these data further extend previous studies that had identified signaling mutations in these subgroups. Our data show that B-ALL cases with signaling mutations typically have multiple subclones, often with various mutations affecting the same pathway. These findings suggest that there is a higher mutation rate and a higher need for additional mutations in PAX5alt and high hyperdiploid cases compared with other B-ALL subtypes.

These observations also have clinical implications. Our data indicate that B-ALL cases with signaling mutations are likely to harbor additional mutations in often very small subclones. Therefore, it seems better not to use such signaling mutations as markers for MRD studies, since minor subclones are likely to harbor other mutations and those could also contribute to positive MRD or relapse. A second clinical implication is on the future use of targeted therapies for the treatment of B-ALL cases with FLT3, NRAS/KRAS, or JAK kinase point mutations. Our data point out that in B-ALL cases with signaling mutations, there is typically a diverse set of subclones present with various mutations, and that targeted therapy is unlikely to be targeting each of the subclones. Similar observations were made in AML cases with FLT3 mutations, where rapid enrichment of NRAS or KRAS mutant subclones was observed upon treatment with FLT3 inhibitors.[28] Activating mutations in FLT3 are common alterations in AML, and can also occur in ALL, most frequently in the high hyperdiploid subtype.[29] Also in the current study, the 2 B-ALL cases with FLT3 mutations were of the high hyperdiploid subtype and FLT3 mutations were detected together with other RAS signaling mutations in other subclones. In one B-ALL case, we observed multiple (>10) different in-frame alterations in the FLT3 juxtamembrane domain in small subclones, indicative of a high genome instability and a high selection of cells with such signaling mutations during leukemia development in such cases. This is similar to the high mutation rate of NOTCH1 in T-ALL cases.[21]

We observed multiple alterations in the PAX5alt and high hyperdiploid subtypes of ALL, with the highest number of subclones in high hyperdiploid cases. Despite the fact that these subtypes, which are among the most common subtypes of ALL, have a good to intermediate prognosis and the risk for relapse is low, the absolute number of patients who relapse for these 2 subgroups combined is significant. Targeted therapies should be considered for these relapsed patients and single-cell sequencing to map the variant clones and mutations present at relapse could support the choice of the best agents. Sequencing relapse samples from the patients in this cohort would be especially interesting, as this would allow us to study the clonal composition at relapse and identify treatment-resistant clones. However, only one of the patients from this study has relapsed recently.

A second question we wanted to address is on the pattern of mutations and the origin of the subclones. Our data indicate that there are various modes of mutation accumulation in B-ALL. Some cases show a consecutive acquisition of the mutations, while other cases seem to acquire multiple mutations in parallel. While in T-ALL patients a clear hierarchical mutational tree could be identified with multiple subclones emerging from a single parent clone,[21] we see a more parallel hierarchy in B-ALL, where clones are appearing in parallel with only little sequential events. This branched pattern of clonality with absence of a major initial lesion suggests that small indels or somatic SNVs are not the founding lesions for leukemic transformation and progression in B-ALL. Surprisingly, some cases acquired various mutations in the same signaling pathway in parallel. This indicates that similar events happened in different cells independent from each other, suggesting a pressure on the developing leukemia cells to acquire mutations in that specific pathway to allow further progression. We assume that major chromosomal aberrations are at the basis of leukemic transformation in B-ALL and are then enforced by additional subclonal somatic mutations in signaling pathways, which are acquired later during the progression of the leukemia.

A third question is whether the different subclones are equally sensitive to chemotherapy. We therefore performed single-cell DNA sequencing on blood samples taken during the first days/weeks of treatment in the 4 cases with the highest clonal heterogeneity. While we observed that all leukemia clones rapidly decreased during chemotherapy treatment, we also noticed that some clones declined with a slower rate. Although this does not provide direct evidence of resistance or lower sensitivity to chemotherapy, our observations warrant further investigation and suggest that some clones with a specific mutational pattern within a single patient may show lower sensitivity towards therapy. This could result either in a slower response or also to resistance and relapse. It will be important to consolidate these findings with studies in larger ALL cohorts covering all subtypes.

ACKNOWLEDGMENTS

We thank Mission Bio for early access to Tapestri and for help with data analysis.

AUTHOR CONTRIBUTIONS

SM and LAS designed the study, performed experiments, analyzed experimental and bioinformatics data, and wrote the article. OG, ME, and TS performed experiments and analyzed data. JDB, LM, BD, NB, AU, JM, and HS provided samples and analyzed clinical and genetic data. KDK analyzed data. JC, HS, and SD designed and supervised the study, analyzed data, and wrote the article. All authors provided input for the article.

DISCLOSURES

JC is the editor of HemaSphere and member of the speakers bureau of MissionBio. The other authors have no conflicts of interest to disclose.

SOURCES OF FUNDING

This project was funded by grants from Kom op tegen Kanker (Stand up to Cancer) the Flemish cancer society (JC, HS), the VIB Technology Watch Fund (JC), the KU Leuven (C14/18/104, JC, KDK, JM, NB, HS), FWO and EU MSCA fellowships to LAS, a Foundation against cancer fellowship to SD, a KU Leuven and FWO fellowship to SM. This study was supported by the VIB Nucleomics Core and the VIB Technology Innovation Lab.

Table 1.

Clinical Characteristics of the 12 B-ALL Cases at Diagnosis

Patient	Gender	Age (y)	WBC (109/L)	% BM Blasts	% PB Blasts	BM Morphology (Day 35)	MRD PCR (Day 35)	MRD PCR (Day 90)
XF93	Female	1–10	92.0	82	77	CR	<1 × 10–4	Negative
XF99	Male	1–10	17.5	59	25	CR	Negative	Negative
XF101	Male	11–16	42.6	61	53	CR	Negative	Negative
XF104	Female	11–16	78.7	87	86	CR	<1 × 10–4	Negative
XF106	Male	1–10	11.4	87	79	CR	Negative	Negative
XF109	Female	1–10	31.6	80	56	CR	Negative	Negative
XG110	Male	1–10	28.3	90	37	CR	<1 × 10–4	Negative
XG111	Female	1–10	14.3	95	74	CR	Negative	Negative
XG115	Female	1–10	71.2	98	59	CR	8 × 10–4	Negative
XG118	Female	>16	8.8	n.a.	8	CR	Negative	Negative
XG124	Female	1–10	216.6	90	95	CR	Negative	Negative
XG125	Male	1–10	92.4	97	74	CR	Negative	Negative

BM = bone marrow; CR = complete remission; MRD = minimal residual disease; PB = peripheral blood; WBC = white blood cell count.

Table 2.

Genomic Data

Patient	ALL Subtype[4]	Karyotype	Structural Variants	Mutations Detected by Single-cell Sequencing (% Mutated Cells)
XF93	ETV6-RUNX1	46,XX[22]	ETV6-RUNX1 fusionPAX5 exon 2–6 deletionTCF3 deletionETV6 deletion	/
XF99	ETV6-RUNX1	46,XY[15]	ETV6-RUNX1 fusionPAX5 exon 1–6 deletionETV6 exon 1–5 deletion	SETD2 intronic, splicing (1%)
XF104	BCR-ABL1	46,XX,t(9;22)(q34;q11)[8]/46,sl,add(16)(p12)[2]	BCR-ABL1 fusionCDKN2A/B deletionIKZF1 exon 4–7 deletionPAX5 deletionRB1 exon 19–26 deletionBTG1 exon 2 deletion	/
XG125	BCR-ABL1-like	46,XY[23]	PAX5-JAK2 fusion	/
XG118	IKZF1 N159Y	47,XX,+21[9]/47,sl,add(3)(q13),del(4)(q13q24),add(6)(p11),–12,add(14)(p11),+mar[3]/46,XX[8]	/	IKZF1N159Y (43%)
XF101	PAX5alt	91<4n>,XXYY,–17,idic(17)(p11)[8]/92<4n>,XXYY[1]/46,XY[5]	PAX5-ZNF318 fusionTP53 deletion	SETD2K2511Nfs*10 (88%)EZH2R690L (1%)
XF109	PAX5alt	46,XX,idic(9)(p1?3)[7]/46,XX[3]	9p deletion (PAX5, CDKN2A/B)	PAX5G183R (84%)JAK1L910P (47%)JAK1Y652H (14%)
XG110	PAX5alt	46,XY,–2,?t(9;17)(p13;q21),del(11)(q21q23),?der(14)t(2;14)(?p11;q32),+mar[8]/46,XY[1]	KMT2A deletionPAX5-ZCCHC7 fusion	CREBBPH1487Q (4%)
XG124	PAX5alt	46,XX,+8,–20[1]/46,XX,+8,der(9)t(9;20)(p13;p11),–20[5]/46,XX[1]	PAX5-NOL4L fusionPAX5 exon 10 deletion CDKN2A/B deletion	NRASG12S (90%)NRASG12D (3%)
XF106	High hyperdiploid	55,XY,+X,+4,+6,+8,+11,+14,add(17)(p13),+18,+21,+21,inc[7]/46,XY[7]	CDKN2A/B deletionIL3RA duplicationP2RY8 duplication	CREBBPL1681del (1%)CREBBPS1680del (78%)OVGP1P496T (2%)OVGP1T495N (2%)
XG111	High hyperdiploid	55,XX,+4,+6,+8,+10,+14,+17,+18,+21,+21[4]/46,XX[6]	/	PTPN11A72V (78%)FLT3N676K (6%)PTPN11E76Q (3%)NRASG12D (2%)BRAFK601N (1%)NRASG13D (1%)KMT2DQ3735* (1%)
XG115	High hyperdiploid	56,XX,+X,+4,+6,+8,+10,+17,+18,+21,+21,+mar[7]/55,sl,t(15;15)(q21;q24),–17[2]/46,XX[2]	CDKN2A/B deletionETV6 exon 1 deletionPAR1 duplicationIGH rearrangement	KRASG13D (17%)FLT3I836del (15%)KRASG12D (3%)FLT3A680V (3%)FLT3Y589-V592del (2%)NRASG12D (2%)NRASG12S (1%)

Genomic lesions found in the 12 B-ALL cases, based on karyotyping, molecular diagnostics, bulk RNA-sequencing, and Bionano Optical Genome Mapping. Variants present at diagnosis were determined using Mission Bio targeted single-cell DNA sequencing, with indication of the percentages of cells in which the mutations were found in brackets.