We studied 51 childhood B-ALL patients diagnosed in Shenzhen Children’s Hospital. BM specimens and (or) PB were obtained at diagnosis and during treatment, and in total 169 specimen were collected (Table S1 in Supplementary Material). The study was carried out in accordance with the recommendations of Declaration of Helsinki and was approved by BGI-IRB. Written informed consent was obtained from the parent(s) or guardian(s) of each child. All the BM and PB specimens were collected in heparin and stored at −80°C until analysis was conducted. High-quality gDNA were extracted from the frozen BM and PB samples using DNA Blood mini kit (QIAGEN, Cat. no.51106).

B-ALL patients were stratified into three risk groups according to the following criteria:

Standard risk (SR): (1) age at diagnosis between 1 and 6 years; (2)WBC < 20 × 10⁹/L; (3) good prednisone respond (GPR) at 7 days treatment, peripheral blasts < 1.0 × 10⁹/L at day 8; (4) BM aspiration results M1 (blasts < 5%) or M2 (blasts 5% ~25%) at day 15 post induction; and (5) BM aspiration results M1 at day 33 post induction.

Intermediate risk (IR): (1) age at diagnosis <1 year or ≥6 years; (2) WBC ≥ 20 × 10⁹/L; (3) GPR; (4) BM aspiration results M1 or M2 at day 15 post induction; (5) BM aspiration results M1 at day 33 post induction; (6) reach to SR but BM aspiration results M3 (blasts > 25%) on day 15 of induction therapy, and BM aspiration results M1 on day 33 of induction therapy.

High risk (HR): (1) IR but BM aspiration results M3 on day 15 of induction therapy; (2) poor prednisone respond (PPR), peripheral blasts > 1.0 × 10⁹/L at day 8; (3) BM aspiration results M2 or M3 at day 33 of induction; (4) t(9:22) (BCR/ABL) positive or t(4:11) (MLL/AF4); (5) testicle leukemia at diagnosis and did not clinically resolve by day 33 of chemotherapy; (6) large mediastinum mass at diagnosis, and did not completely resolve by day 33; (7) central nervous system leukemia (CNSL) at diagnosis.

Immunophenotyping and MRD detection by FCM was performed at KINGMED CENTER FOR CLINICAL LABORATORY. The red blood cells of BM samples were lysed to get nucleated cells. For MRD detection, BM nucleated cells were stained with the following monoclonal antibody combinations: (1) CD58 FITC, CD34 PE, 7AAD PerCP-CY5.5, CD10 PE-Cy7, CD19 APC, CD38 V450, CD45 V500; and (2) CD66c FITC, CD13+CD33 PE, CD34 PerCP-CY5.5, CD10 PE-Cy7, CD19 APC, CD15 V450, CD45 V500. The labeled cells were analyzed using a BD CANTO II flow cytometer (Becton Dickinson, San Jose, CA, USA) with at least 100,000 events were acquired. Data were analyzed using BD FACSDiva (Becton Dickinson) and FCS Express (De Novo Software, Los Angeles, CA, USA) software. The compensation matrix was set up using BD CompBeads (Becton Dickinson) for fluorochrome-conjugated antibodies. Quality control was performed using BD Cytometer Setup and Tracking Beads (Becton Dickinson).

The complementarity determining region 3 (CDR3) of the variable regions of IGH was amplified by multiplex PCR. Concretely, the complete VDJ rearrangements of IGH were amplified from 1200 ng gDNA with 12 degenerate forward primers annealed to the 55 functional variable genes (V) and 4 reverse primers annealed to the 6 functional joining genes (J) listed in the IMGT database. The primers have been carefully evaluated to minimize PCR bias, and were listed in Table S4 in Supplementary Material. PCRs (50 μL) were set up at 25 μL of 2× QIAGEN Multiplex PCR master mix, 5 μL of QIAGEN Q solution, 1 μL of 10 μM forward primer pool, 1 μL of 10 μM reverse primer pool, and 18 μL of 67 ng/μL gDNA. The reaction cycling conditions were: 95°C for 15 min, 30 cycles at 94°C for 30 s, 60°C for 90 s, and 72°C for 30 s, followed by a final extension at 72°C for 5 min. The target amplified products (120–200 bp) was purified by electrophoresis on 2% agarose gel, and the Illumina Hiseq sequence adaptors were ligated. Then the sequencing libraries were sequenced with standard 2 × 150 paired end reads on Illumina Hiseq2500 platform.

Sequencing data were analyzed by the TCR and BCR repertoire analyzing pipeline IMonitor (25). About 2–5 million mapped reads (include correct V and J genes) remained after the above pipeline for each sample, and we used a random subset of 2 million mapped reads for continued analysis. In order to remove the false IGH clones derived from sequencing error or contamination from other samples in the same Illumina sequencing lane, we filtered the clonotypes with supported reads lower than the inferred average sequencing depth for every B cells in each sample. The average sequencing depth was inferred as following: first, the total input cell number for amplification was 2 × 10⁵ cells, equating to 1200 ng input gDNA. Second, we calculated the B cell number by multiplying the total input cell number by the B cell percentage (detected by FCM according the CD19 marker). Third, the average sequencing depth was inferred by dividing the total used reads (2 million) by the B cell number used in the sample. The CDR3s were classified as non-functional if frame-shift mutations or stop codons existing in CDR3 region.

Similar to previous studies, we defined the CDR3s in diagnostic samples with frequency higher than 10% as leukemic clones (17, 18). In the post-treatment samples, we determined the MRD by calculating the frequencies of identical CDR3 sequences with the leukemic clones in the level of total nucleated cells. The MRD is calculated by multiplying the frequencies of leukemic clones by B cell percentage in the total nucleated cells, which is determined by FCM.

The method used by Gawad et al. (24) was used to identify the evolved IGH clones from the leukemic clones with minor modification. Concretely, the evolved IGH clones were defined as: (1) identical J sequence with the leukemic clone; (2) share at least eight bases of identical NDN sequence with the leukemic clone; (3) different V gene from the leukemic clone; (4) more than three mismatches in CDR3 sequence if the evolved IGH clone had the same length as the leukemic clone.

We randomly, without any prior knowledge of the clonal IGH rearrangement existence, collected 51 pediatric B-ALL patients in total. According to the disease risk stratification criteria, 19 of them were classified as SR, 18 as IR, and 14 as HR (Table S1 in Supplementary Material, details in the section “Materials and Methods”). The characteristics and clinical information of the patients, including diagnostic age, gender, and cytogenetics were collected and listed in Table S1 in Supplementary Material. Ninety-six percent of patients (49/51) were aged between 1 and 10, and 2 patients were more than 10 years old. Among the patients, 35.3% (18/51) were females, and 23.5% (12/51) were diagnosed with TEL-AML1 gene fusion. All the 51 patients showed good treatment result after 2–3 months of chemotherapy.

According to the definition of leukemic clones (see Materials and Methods), we identified leukemic IGH CDR3s with complete V–D–J rearrangement in a frequency of above 10% in the diagnostic BM samples in 47 of 51 patients (92.2%). Of the 47 patients with leukemic IGH clones, 24 had two or more leukemic clones (Table 1; Table S2 in Supplementary Material). Therefore, we identified 77 leukemic clones from the 47 patients in total (Table S2 in Supplementary Material). The four patients without identification of leukemic clones did not show any significantly distinct immunophenotypes. They may harbor incomplete VDJ rearrangement, or use the pseudo genes to recombine VDJ which cannot be amplified by our multiplex PCR primer set. We then investigated if the number of leukemic clones related with the disease risk. The result showed that more patients in SR group have just one disease clone; however, more patients in HR group have two or more disease clones (Figure 1).

Because B-ALL is a kind of cancer originating from the immature stage of B cells, we investigated if the leukemic IGH rearrangements underwent positive selection. Of the 77 leukemic clones, 28 (about one-third) were functional with correct open reading frame (ORF) in the CDR3, and the other 49 were classified as non-functional due to the occurrence of frame-shift or stop-codon mutations in the CDR3 (Table S2 in Supplementary Material). The ratio of functional to non-functional clones was proportional with that expected by chance. The data above demonstrated that these B cells were immature precursor B cells without experiencing positive selection. We next compared the mRNA transcription of the functional and non-functional leukemic clones by sequencing the IGH repertoire of RNA extracted from the same BM specimen of nine patients. All the leukemic clones, including the non-functional clones, were detected in RNA with high frequency. Surprisingly, the transcription of non-functional leukemic clones was not eliminated or even decreased compared with the functional clones (Figure S1 in Supplementary Material), implying the resistance to the nonsense mRNA decay pathway (26, 27). However, compared with the frequencies in DNA samples, the frequencies of most leukemic clones in RNA decreased largely, both for functional and non-functional ones (Table 2), which could be due to decreased transcription of the pre-B cells (28), including B-ALL leukemic cells, and it is consistent with other report (29).

With the precise IGH CDR3 sequence information of the leukemic clones, we determined the frequencies of the leukemic clones in the follow-up BM samples. MRD based on the total nucleated cells was calculated by multiplying the frequency of the leukemic clone by the percentage of the B cells in total nucleated cells, which was determined by FCM for those samples. MRD levels and leukemic clone frequencies decreased as treatment progressed (Table S2 in Supplementary Material and Figure S2 in Supplementary Material). With the purpose of evaluating the reproducibility of MRD detection by HTS, we performed two replicated MRD detections from equal BM specimen for three samples of two patients (P049-15, P049-33, and P051-15) and the result showed good reproducibility (Table S3 in Supplementary Material).

In order to assess the accuracy of the MRD detection by IGH repertoire HTS, we compared the detection results of HTS with the results detected by FCM, which were performed as the golden standard for MRD detection. All 34 MRD positive samples by FCM were detected as positive by HTS with similar MRD levels. Twenty-eight samples were detected as MRD positive by HTS, but MRD negative by FCM. The other 19 samples were MRD negative by both HTS and FCM. No MRD positive samples by FCM were detected negative by HTS, which demonstrated the higher sensitivity of HTS (Figure 2).

Because BM aspiration can be painful, using PB instead of BM was evaluated for MRD detection. For 15 patients, we collected both diagnostic BM and PB samples (Table S1 in Supplementary Material). We investigated if the leukemic clones identified in BM also existed in PB with high frequencies. The results revealed that all leukemic clones could be detected in PB with frequency of above 5%, except one clone, which was 18.405% in BM, but was just 0.810% in PB (Figure 3A). The Pearson correlation coefficient between the frequencies of leukemic clones in diagnostic PB and BM was 0.7860. Therefore, PB can be used instead of BM in disease diagnosis and in leukemic clone identification for most patients.

To test MRD detection in PB specimens and to compare it with BM, post-treatment BM and PB samples were collected simultaneously in 15 pairs from 11 patients (Table S1 in Supplementary Material). We determined the frequencies of leukemic clones using the same method. In total, 22 leukemic clones were identified in pre-treatment BM samples for the 11 patients, and 18 clones were positive in post-treatment BM. When detecting the 18 clones in the corresponding PB samples, 12 clones were detected and other 6 were PB negative. When evaluating their frequencies in BM, 11/12 PB positive clones had relatively higher frequency than the 6 PB negative clones. Three leukemic clones were MRD negative in both PB and BM; one was positive in PB with clonal frequency of less than 0.01%, but was negative in BM. We then assessed if the leukemic clone frequency in PB correlated positively with that in BM, and a weak correlation was found with a Pearson correlation coefficient of 0.322 (Figure 3B). Overall, although PB could be used to detect MRD to some extent, the sensitivity of MRD detection using PB was still lower than using BM and it resulted in false negative for some samples.

On-going change at the IGH gene is one of the most important ways leading to the diversity of cancer cells in B-ALL. After finishing IGH variable (V) gene, diversity (D) gene, and joining (J) gene rearrangement during B cell development (30), the V segments upstream of the used V gene are still reserved, which make the V gene replacement possible by rearranging an upstream V gene with the complete V–D–J rearrangement. Because the loss of 5′ recombination signal sequences (RSS) in D segment of the rearranged V–D–J exon, cryptic RSS within the used V gene is involved in V gene replacement, which effectively increases the length of the V–D–J junction of the new rearrangement (31). The identified evolved IGH rearrangements of leukemic clones in this paper could be generated by three different mechanisms. The first and second mechanisms are based on V replacement in IGH gene. The first mechanism is that the leukemic clone is the ancestral clonotype and the evolved IGH clones derived from continuous V replacement. In this situation, the V segment of the evolved IGH clones will be located at the upstream of the leukemic clone, and the CDR3 length of evolved IGH clones would be longer than that of leukemic clone (31). The second mechanism is that the leukemic clone is produced by V replacement from a more ancestral pre-leukemic clonotype. In this situation, the V segment of leukemic clone will located at the downstream of the ancestral clone, and the CDR3 length of the leukemic clone would be longer than the ancestral clone. The third mechanism is that the ancestral (pre-leukemic) clone finished an incomplete D–J rearrangement of IGH locus, and the leukemic clone and other evolved IGH clones underwent independent IGH V rearrangement. Therefore, the CDR3 length of the leukemic clone and the evolved IGH clones would be comparable.

To explore the roles of these mechanisms, we assessed and compared the V segments used and the CDR3 length between the leukemic clones and the evolved IGH clones. In total, 51.62% of evolved IGH clones were mainly replaced by upstream V genes compared to 19.67% by downstream V segments (Figures 4A,B). For the evolved IGH clones using upstream V segments, most of the CDR3s were longer than the corresponding leukemic clones; however, the CDR3 lengths of the evolved IGH clones using downstream V segments were comparable with the leukemic clones (Figure S3 in Supplementary Material). The above result demonstrated that the evolved IGH clones we identified were produced mainly by the first mechanism; in addition, the third mechanism also contributed to the evolved IGH clones. Because the pre-leukemic clones were also important in leukemic development and could lead to relapse of disease, we included all those evolved IGH clones in subsequent analysis.

In the diagnostic BM samples, we identified evolved IGH clones in most of the patients, and only two patients did not have evolved IGH clones. We also found that four leukemic clones were produced due to V gene replacement of the other four leukemic clones in three patients (Table S2 in Supplementary Material). The number of identified evolved IGH clones varied from 0 to 4558 for each leukemic clone (Figure 4A), with 14 (18.2%) of the leukemic clones did not experience IGH evolution. When investigating the ORF of the evolved IGH clones, only 11.5% of those clones were functional with correct ORF, which is much lower than the proportion of the leukemic clones. The total frequency of these clones varied from 0 to 38.55% for each leukemic clone (Figure 4C). There was a positive correlation between the total frequency and the number of evolved IGH clones for each leukemic clone (Pearson’s r = 0.6661; Spearman’s r = 0.9851). Whether the number and frequencies of the evolved IGH clones were determined by the total number of the leukemic cells is an interesting question. We found that both the number (Pearson’s r = −0.0490) and total frequencies (Pearson’s r = −0.0985) of evolved IGH clones was not associated with the frequency of the cancer clone (Figure 4C). The association of evolved IGH clones with clinical features is another important question; however, we did not identify their correlation with disease characteristics in our dataset (Figures S4A,B in Supplementary Material).

We next investigated the dynamic change of the evolutionary IGH clones following chemotherapy treatment in BM. Forty-six percentage of the MRD positive samples contained evolved IGH clones, and we did not detect evolved IGH clones in MRD negative samples (Figure 4D). We also investigated the ORF of these evolved IGH clones, and the results showed only 12.7% of those clones were functional with correct ORF, which was similar to the proportion of evolved IGH clones in diagnostic samples. The total frequency of these evolved IGH clones showed very weak correlation with the frequency of leukemic clones (Pearson’s r = 0.2203, Figure 4D). Tracking the fluctuation of frequency ratios of the evolved IGH clones to the leukemic clones revealed the imbalanced response of leukemic clones and evolved IGH clones to chemotherapy. The ratios of total frequencies of evolved IGH clones to the frequencies of leukemic clones decreased during treatment for most of clones; however, for several clones, the day 15 samples showed higher ratio than the diagnostic samples (Figure S5 in Supplementary Material), which may imply positive selection of the evolved IGH clones upon the pressure of chemotherapy in the begin of treatment. Interestingly, we found a significant number of evolved IGH clones newly emerged in the post-treatment samples (Figure S6 in Supplementary Material), even though these clones eventually disappeared with the negative detection of MRD after chemotherapy.

We compared the number and total frequency of evolved IGH clones in PB and BM of the same patients. The result showed that the number and total frequency of evolved IGH clones in PB were correlated positively with that in BM (Pearson’s r = 0.9941 for number; Pearson’s r = 0.7154 for total frequency). However, for most samples, the total frequency of evolved IGH clones in PB was lower than that in BM (Figure 5). This result revealed that PB has the potential to serve as a non-invasive specimen to monitor clonal IGH evolution during therapy, but it may underestimate the frequency of evolved IGH clones.