A total of eight patients with myeloid neoplasia and clonal evolution were included in this study. The patients were analyzed between 2011 and 2021 at our department. The ethical review boards of Hannover Medical School approved this study and all patients gave their written consent. Our cohort was analyzed at one or more time points. All results of performed analyses as well as clinical data are shown in Table 1 and in more detail in Table S1 .

Cytogenetic and molecular cytogenetic analyses were performed on bone marrow aspirates or peripheral blood cultures. Chromosome preparation and fluorescence R-banding were performed as described previously (21, 22). Altogether, whenever possible, 15 to 25 metaphases were examined per patient. The karyotype was described according to guidelines of the International System for Human Cytogenetic Nomenclature (23). Depending on the cytogenetic aberrations in each patient, FISH analyses on interphase nuclei were performed using: (1) break apart probes for the loci 3q26 (EVI1 (MECOM) Break Apart; Cytocell, Cambridge, UK), 11q23 and 12p13 (Vysis MLL and ETV6 FISH Probe Kit; Abbott, IL, USA); (2) dual color probes for the loci 5p15.2, 5q31, CEP7, 7q31, 9q34.1, 17p13.1, CEP17, 22q11.2 (Vysis EGR1/D5S23, D5S721, D7S486/CEP7, BCR/ABL and TP53/CEP 17 FISH Probe Kit; Abbott); and (3) single color probe for the locus 20q12 (Vysis D20S108 FISH Probe Kit; Abbott). At least 200 interphase nuclei were analyzed for each probe.

According to the manufacturer’s protocol (e-Array design 84704, Agilent Technologies, Waldbronn, Germany), three patients (patients c1, 6 and 7) were screened for CNVs by high-resolution array-CGH. Microarray slides were scanned immediately using an Agilent microarray scanner at a resolution of 2 μm. Fluorescence ratios were calculated using Feature Extraction Software and copy number states analyzed using the CGH data analysis software Genomic Workbench (Agilent Technologies).

DNA sequencing was performed for all patients using an IDT custom panel (Integrated DNA Technologies Inc, Iowa, USA) including 148 leukemia-associated genes. Sequencing was performed on an Illumina NextSeq sequencer. All samples were aligned against the reference genome GRCh37 using BWA mem (24). The pipeline uses eight open-source tools for independent variant calling. Subsequently, raw calls are automatically combined and characterized with respect to data quality. Default settings were used (minimum number of reads 50, minimum number of reads with the alternate allele 20, minimum variant allele frequency VAF 0.01). Excluding low-quality calls, the remaining variants are further characterized (e.g. presence in common databases and in silico effect prediction). Combining all information, an artifact- and a polymorphism score is calculated for every SNV and small indel. The scores allow for automatic classification of the calls as true variants, polymorphisms and artifacts. Additionally, all detected variants were manually checked using the Integrative Genomics Viewer (IGV) to confirm automatic classification (25, 26). All detected variants were classified according to the standards and guidelines of the American College of Medical Genetics and Genomics (ACMG) (27). Here, we report only variants of unknown significance (VUS), likely pathogenic and pathogenic variants.

SNVs and short indels were detected using targeted NGS and variant calling pipeline appreci8 (28). SVs and CNVs were detected performing karyotyping, FISH and array-CGH analysis. Clonal evolution was reconstructed manually, integrating all mutational information on SNVs, indels, SVs and CNVs as described in Reutter et al. (29). The results were visualized by fishplots (30). In short, percentages of cancer cell fractions (CCFs) for SNVs and indels were estimated based on VAFs, assuming heterozygous variants (2*VAF=CCF). Percentages of CCFs for SVs and CNVs were estimated based on cell counts reported for karyotyping and FISH analyses. For SVs and CNVs that were only detected in array-CGH percentages of CCF were estimated based on log₂Ratio of fluorescence intensities for reference versus tumor probes. In case of CNVs overlapping the position of an SNV or indel, the calculation of CCF is less straightforward. Altogether, we differentiate between three possible scenarios:

We defined w as the ratio of cells featuring a CNV and an SNV/indel. Analogously, x, y and z were defined in the same manner. The sum over all cell ratio was always w+x+y+z=1.0. Furthermore, the known CCF for the CNV was defined as CCF_CNV=w+x. Additionally, two formulas were derived from the model (see Figure 1A ):

We defined cnv_value as a parameter to model the number of alleles with the SNV/indel depending on the CNV. For scenarios 1 and 3 the cnv_value=1 and for scenario 2 the cnv_value can be 0, 1 or 2. Dependent on the three scenarios we considered, the two equations may be further simplified: For scenario 1, the CNV occurring prior to the SNV/indel, y=0 (see Figure 1B ). Thus, CCF_CNV=w+x and CCF_SNV/indel=w. As VAF_SNV/indel is also known, w and CCF_SNV/indel can be easily determined. For scenario 2, a distinction of cases has to be made. Both, the allele with or without SNV/indel may be affected by CNV. Regarding deletions, cnv_value=0 if the allele with SNV/indel is deleted, and cnv_value=1 if the allele without SNV/indel is deleted. Regarding duplications, cnv_value=2 if the allele with SNV/indel is duplicated, and cnv_value=1 if the allele without SNV/indel is duplicated. Additionally, x=0, so CCF_CNV=w. The unknown value of y and thus CCF_SNV/indel=w+y can be determined using VAF_SNV/indel , following the distinction of cases. For the third scenario of SNV/indel and CNV occurring in parallel, independent of each other, w=0. Thus, CCF_CNV=x and CCF_SNV/indel=y. As VAF_SNV/indel is known, the unknown value of CCF_SNV/indel can be determined. Although CCF_SNV/indel can always be determined using the above-mentioned formulas, the result may not always be distinct. Information on whether the CNV occurred prior or after the SNV/indel is in general not available. Therefore, for application of our approach, all scenarios have to be considered and four possible values for CCF_SNV/indel have to be calculated. Scenarios may be excluded if the basic assumption w+x+y+z=1.0 is violated.

In order to show that clonal evolution can be reconstructed and visualized even when only one time point is available, we examined three patients (patient 1, 2 and 3) with the aforementioned bioinformatic approach. The total number of aberrations in a patient ranged from 3 (patient 3) to 11 (patient 1). Due to this approach it was possible to determine clones by clustering mutations. These clones were sorted according to CCF, in order of their appearance. Through reconstruction and visualization of clonal evolution it was possible to identify the initial clone (stemline) as well as the passenger aberrations for all patients. For example, patient 1 (AML) showed a deletion in 5q (del(5)(q14q33)) in the initial clone, patient 2 (MDS) showed a trisomy 8 as well as a pathogenic variant in DNMT3A in the initial clone, and patient 3 (therapy related myeloid neoplasm) showed a translocation t(12;22) in combination with a pathogenic variant in NF1 in the initial clone (see Figure 2 and Table 1 / S1 ). The total number of clones in a patient ranged from 2 (patient 3) to 7 (patient 1). In summary, our approach was successful in identifying linear clonal evolution in all patients.

Next, we analyzed five patients where data was available for more than one time point (2 to 5 time points of analysis) to show that the bioinformatic approach can reconstruct different models of clonal evolution. At each time point, all patients were analyzed by a combination of methods, including karyotyping, FISH and a custom NGS-Panel. The estimated therapy effect was included in the reconstruction of clonal evolution. The approach reconstructed a distinct model of clonal evolution for patients 5, 7 and 8 (see Figures 2C–E ). Patient 5 was analyzed at two time points – initial and progression. The bioinfomatic approach reconstructed a linear clonal evolution pattern. The initial clone showed a translocation t(4;14;11) and a structural aberration in chromosome arm 10p (add(10)(p14)). At this time point, pathogenic variants in KRAS, ASXL1 and NRAS were identified with a very low percentage of CCF (close to 1%). However, the pathogenic variants were only identified by deep sequencing. At the second time point (progression), additional chromosome aberrations (+der(4)(4;14;11), +8, +9, +19 and +21) were identified and this subclone became the mainline (see Figure 2D ). Patient 7 was analyzed at four time points – initial, progression, remission and relapse. Stem cell transplantation was performed between the third and fourth time point of analysis. The initial clone showed an insertion of a segment of chromosome 12 into chromosome 9 (ins(9;12)) and a trisomy 12. At time point two, the percentage of cells only harboring the stemline’s mutations were greatly reduced (~5%) and a subclone developed. This clone showed a deletion in 12p (del(12)(p13)) and a trisomy 11, 19 and X. For the subclone, the approach calculated linear clonal evolution. At time point three, the subclone was no longer detectable, but a new subclone was identified at time point four. This subclone presented a derivative chromosome 7 with a translocation t(7;9) (der(7)t(7;9)(q35;q21)). At this time point, the clonal evolution pattern of patient 7 changed from linear to branched (see Figure 2E ). Patient 8 was analyzed at five time points – initial, remission, progression, relapse and remission. Stem cell transplantation took place between the first and second time point. The initial clone showed a deletion in 9q (del(9)(q21q31)) and a subclone with a likely pathogenic variant in WT1 as well as a pathogenic variant in NRAS. At this time point, a linear evolution pattern was reconstructed. At time point two, after stem cell transplantation, no variants were detected – neither by karyotyping nor deep sequencing. However, the aforementioned variants were present at time point three (progression). At this time point, a linear clonal evolution pattern was reconstructed in one subclone (light red). Furthermore, four additional subclones were observed. All subclones developed independently of the initial clone with a deletion in 9q (del(9)(q21q31)). At the fourth time point of analysis, three additional subclones were detected. The clonal evolution pattern of patient 8 changed from linear to neutral. At the last time point of analysis (remission), no variants were detected (see Figure 2F ). In summary, a linear evolution pattern was reconstructed for patient 5, a branched evolution pattern was reconstructed for patient 7 and a neutral evolution pattern was reconstructed for patient 8.

Reconstruction of clonal evolution in patients 4 and 6 showed more than one possible clonal evolution pattern. All possible versions of reconstruction are shown in Figure 3 . Patient 4 was analyzed at two time points (initial and progression). The initial clone showed a deletion in 5q (del(5)(q14q34)). Additionally, four subclones developed. At this time point, patient 4 showed a linear clonal evolution pattern. At time point two, four additional subclones were detected and the evolution pattern changed from a linear to a branched evolution pattern. Although the model of clonal evolution might be clear, it was not possible to reconstruct the precise subclonal relationship. The smallest subclone with (likely) pathogenic variants in KRAS, NRAS and PTPN11 showed a low percentage of CCF (5%). Based on the available data, the parent clone of the smallest subclone cannot be definitely determined and, therefore, patient 4 showed three possible versions of branched evolution patterns (see Figure 3A ). Patient 6 was analyzed at three time points – initial and two progression time points. The initial clone showed a pathogenic variant in JAK2. At the second time point of analysis, we identified two events: (1) a likely pathogenic variant in TP53 and (2) a derivative chromosome 17 resulting in loss of TP53 (der(17)t(13;17)(q21;p12)) with further aberrations. Based on the available data, it was not possible to reconstruct which of the two events took place first. Furthermore, it was not possible to calculate if both events affected the same allele or different alleles. As described in Materials and Methods, four different cases have to be considered (see Figure 1A ). Two out of four cases were in line with the basic assumption w+x+y+z=1.0. Both cases were confirmed by the observed percentages of CCF and VAFs at the time points one and two of analysis. According to the first case, the CNV event took place first, followed by the likely pathogenic variant in TP53. Both aberrations affected different alleles. For this case, a branched clonal evolution pattern was reconstructed (see Figure 3B ). Mathematically, it is also possible that the point mutation in TP53 took place first and subsequently the mutated allele was affected by the deletion. For this case, a linear clonal evolution pattern was reconstructed. However, a biallelic inactivation of TP53 frequently occurs during disease progression of hematologic neoplasms (e.g. MDS) (31). In summary, three possible branched evolution patterns were reconstructed for patient 4, and two models of clonal evolution (linear and branching) were reconstructed for patient 6.