Darwinian evolution of tumor cells remains underexplored in childhood cancer. We here reconstruct the evolutionary histories of 56 pediatric primary tumors, including 24 neuroblastomas, 24 Wilms tumors, and 8 rhabdomyosarcomas. Whole-genome copy-number and whole-exome mutational profiling of multiple regions per tumor were performed, followed by clonal deconvolution to reconstruct a phylogenetic tree for each tumor. Overall, 88% of the tumors exhibited genetic variation among primary tumor regions. This variability typically emerged through collateral phylogenetic branching, leading to spatial variability in the distribution of more than 50% (96/173) of detected diagnostically informative genetic aberrations. Single-cell sequencing of 547 individual cancer cells from eight solid pediatric tumors confirmed branching evolution to be a fundamental underlying principle of genetic variation in all cases. Strikingly, cell-to-cell genetic diversity was almost twice as high in aggressive compared with clinically favorable tumors (median Simpson index of diversity 0.45 vs. 0.88; P = 0.029). Similarly, a comparison of multiregional sampling data from a total of 274 tumor regions showed that new phylogenetic branches emerge at a higher frequency per sample and carry a higher mutational load in high-risk than in low-risk tumors. Timelines based on spatial genetic variation showed that the mutations most influencing relapse risk occur at initiation of clonal expansion in neuroblastoma and rhabdomyosarcoma, whereas in Wilms tumor, they are late events. Thus, from an evolutionary standpoint, some high-risk childhood cancers are born bad, whereas others grow worse over time.

Significance:

Different pediatric cancers with a high risk of relapse share a common generic pattern of extensively branching evolution of somatic mutations.

Cancer is one leading cause of death among children in developed countries. Based on the histology and the genetic profile of the primary tumor, most pediatric cancers can be split into subtypes predicting their risk of relapse and progression. For example, neuroblastomas can be broadly subdivided into one low-risk group of tumors affecting young children (<18 months of age) with only whole-chromosome (numerical) aberrations, one high-risk group in older children signified by MYCN oncogene amplification, and another high-risk group occurring in older children with tumors carrying structural chromosome changes, ATRX deletions, and/or mechanisms leading to telomere dysregulation (1, 2). Wilms tumors, on the other hand, are subdivided according to the current European protocol into the favorable intermediate-risk histologic subtype, the blastemal histology subtype with a moderate risk of relapse, and the high-risk diffuse anaplastic histology (3), which correlates closely to somatic inactivation of the TP53 tumor suppressor (4, 5). Finally, rhabdomyosarcomas are largely divided into the highly aggressive alveolar subtype having rearrangements of the FOXO1 gene along with complex structural rearrangements and the more favorable embryonal subtype, which has a genotype dominated by numerical aberrations (6).

Despite this abundance of genetic prognostic markers, the causal link between somatic mutations on the one hand and the risk for treatment resistance and relapse remains unclear. The course leading up to fatal disease typically goes via a tumor that is at least partially sensitive to chemotherapy, followed by one or more treatment-resistant relapses. This suggests that chemotherapy provides a selection pressure that favors enrichment of treatment-resistant clones over time through a Darwinian process. A key requirement of selection is heritable variation, in cancers manifested by intratumoral genetic heterogeneity. A small number of recent publications have revealed that intratumoral genetic diversity is indeed a common phenomenon in pediatric tumors (5, 7–12). However, we still know little about how the capacity of cancer cells to evolve relates to their clinical features and the presence of specific prognostic markers. We recently demonstrated how different cancer cell clones in the same tumor can grow together or in distinct domains by following one of four distinct evolutionary trajectories (5). Here, we employ multiregional sampling (MRS) followed by whole-genome copy-number analysis supplemented by whole-exome sequencing (WES) and whole-genome sequencing of single cells to trace ancestral relationships through classical phylogenetic methods in prognostically distinct childhood tumor subtypes. We also use MRS to answer the question of when in each tumor's evolutionary history the genetic hallmarks for each prognostic subtype emerge.

Study design

An extension of the MRS-based dataset published by Karlsson and colleagues (5) was employed to construct tumor cell phylogenies (see Supplementary Materials and Methods). Building robust cancer phylogenies requires an approach that detects a sufficient amount of genetic variation among tumor samples while still allowing the inclusion of as many samples as possible per patient. Because childhood cancers are dominated by large-scale chromosomal alterations rather than functional point mutations, we focused our investigations on copy-number aberrations and copy-neutral genomic imbalances. Such mutations can be detected robustly using high-resolution SNP arrays (SNP-A) on both frozen and formalin-fixed paraffin-embedded tumor material, hence unlocking information from diagnostic samples stored in pathology archives (Supplementary Fig. S1A–S1E). Inclusion criteria for the present study was a histopathologically confirmed diagnosis of Wilms tumor, neuroblastoma or rhabdomyosarcoma with prognostic subtyping performed according to present guidelines, availability of at least two primary tumor samples with >50% tumor cells according to histopathologic assessment, and a documented distance between samples of at least 10 mm. In total, we analyzed whole-genome copy-number and allelic imbalance profiles from 274 tumor samples (240 from primary, 34 from metastatic sites) from 24 neuroblastomas, 24 Wilms tumors, and 8 rhabdomyosarcomas (Supplementary Fig. S1F; Supplementary Table S1). We obtained a median of 4 samples per patient without any significant difference among tumor types (3.5 for neuroblastoma, 4.0 for Wilms tumor, and 4.5 for rhabdomyosarcoma). We also performed WES of 35 biopsies from 19 tumors, followed by targeted deep sequencing (TDS) of 76 samples, allowing complementary scoring of single-nucleotide variants (SNV) and indels. Risk assessment was performed according to standard pediatric oncology protocols (see Supplementary Materials and Methods). Genetic analysis of human tumors was approved by the Regional Ethics Board of Southern Sweden (permit nos. LU289-11 and LU119-03) with written informed consent obtained from parents or legal guardians of the patients. Sequencing and array data are deposited in the European Genome-phenome Archive (www.ebi.ac.uk/ega/home) under accession number EGAS00001002662.

Clonal deconvolution

For statistical robustness, clonal deconvolution requires many data points for each biopsy. This makes sequencing-based methods hard to adapt to the paucity of mutations typical of many pediatric cancers. However, by providing multiple data points (probe read outs) for each detected allelic imbalance and copy-number aberration, SNP-A data can be readily used for clone size calculations and clonal deconvolution using long-established and comprehensively validated methods (7, 13, 14). Accordingly, clone size estimations were here performed by calculating the prevalence of chromosomal imbalances as a function of copy-number and log2 ratio or allele frequency deviation relative to the largest clone (5). For sequence alterations, clone prevalence in each sample was calculated from variant allele frequencies combined with allele copy numbers (see Supplementary Materials and Methods). Genome profiles for individual clones (clonal deconvolution) were then inferred by clustering chromosomal alterations and sequence alterations with similar prevalence levels (Supplementary Table S2). Unbalanced segments <0.1 Mb or covered by <100 markers and clones <10% were excluded from estimations to ensure precision.

Phylogenetic analyses

Following deconvolution, phylogenies modeling the ancestral relationship among detected clones were constructed using a custom-made R-script especially adapted for SNP-A data. As input data, we used an event matrix detailing for each clone/subclone (MRS data) or cell [single-cell sequencing (SCS); see below] the absence or presence of each detected aberration. To validate the robustness of phylogenetic analysis, trees based on maximum likelihood was compared with trees based on maximum parsimony for eight SCS datasets and 12 datasets based on MRS. For 18 of 20 datasets, the two methods produced identical phylogenies. The two remaining datasets (SCS data on RMS6 and MRS data on RMS1) failed to produce stable optima for maximum likelihood trees. Based on this, we concluded maximum parsimony to be a more stable method for the present study. Branches leading to clones that encompassed all (≥90%) cancer cells in a sample were regarded as clonal, whereas those encompassing <90% were regarded as subclonal.

Single-cell sequencing

To determine single-cell DNA copy numbers, shallow genome sequencing was performed on single nuclei (15). To establish cutoff values for scoring copy-number aberrations in each of these cells, we used the copy-number profiles of accompanying noncancer cells present in the samples as a reference. These cells were defined by not showing any of the imbalances detected by SNP array in the corresponding tumor bulk sample. Comparing chromosomes that appeared normal disomic in both tumor cells and noncancer cells, there was no difference in the average alteration rate per 1 Mb bin of sequence reads (0.13% for tumor cells; 0.56% for noncancer cells; P = 0.16; two-sided Student t test). At a cutoff of five or more consecutively altered 1 Mb bins, only one copy-number change was detected over >46,000 bins in these disomic chromosomes, resulting in false positive rate of <2:100,000. We then used a series of interstitial deletions and duplications with known breakpoints from bulk sample array analysis to analyze breakpoint precision by SCS, again resulting in a cutoff of five or more 1 Mb bins for scoring significant differences in breakpoints. Based on these criteria, tumor cells that did not show significant differences in copy-number profiles or breakpoints were grouped together into clones, whereas cells with unique genotypes were kept as single cells at phylogenetic analyses.

Temporal analysis

Phylogenetic branching and the accumulation of mutations over time are not necessarily linked in a straightforward fashion. This is of particular concern in pediatric cancer, having many chromosomal rearrangements that may emerge through sudden bursts (e.g., chromothripsis) or single catastrophic mitotic events (e.g., aneuploidization). Also, pediatric tumors rarely show whole-genome duplication events that can be used in conjunction with sequencing data to derive timelines of mutation. This makes it difficult to employ methodology developed for adult cancers to analyze the temporal sequence of mutations (16). To circumvent this problem, we employed our MRS data on spatial genetic variation to make an estimation of whether individual chromosomal aberrations and mutations occurred preferentially early or late in tumor development. In summary, we dichotomized mutations into early and late events, corresponding to those present in the phylogenetic stem (>90% of tumor cells, all regions) versus those present in less than all regions, respectively (see Supplementary Materials and Methods).

Branching evolution is a common feature of pediatric cancer

Following analysis of allelic imbalances and copy-number aberrations (all cases), supplemented by data on SNVs and indels (19 cases), intratumoral genetic heterogeneity was found in 88% (49/56) of the investigated tumors (Supplementary Table S3). Only in one of these cases (WT5) was genetic variation ascertained only by WES and TDS in the absences of variation by SNP-A, underscoring that the latter method is highly informative while at the same time being employable on archived pathology samples. The genetic variation across geographic regions was then subjected to clonal deconvolution and phylogenetic trees created to link the detected clones by ancestry. In principle, tumor phylogenies can have three different configurations based on the relationship of the branches leading up to genetically distinct clones and subclones (Fig. 1A). In linear evolution, branches are arranged in a straight lineage from the stem toward a series of clones that are direct descendants of each other, reflecting a cumulative accumulation of mutations (Fig. 1B). Alternatively, branches radiate in a collateral fashion from the founder genome, reflecting the presence of private mutations in each detected clone. Finally, there can be a combination of these two scenarios. A total of 212 branches were detected in samples from primary tumors in our dataset (relapses/metastases excluded). Phylogenies with collateral branching (mixed or collateral only) was the typical pattern, being present in 59% to 67% of cases from each of the three tumor types (Fig. 1C). Linear evolution was observed as the only scenario in approximately one third of primary tumors, with an equal distribution (33%–41%) in Wilms tumors, neuroblastomas, and rhabdomyosarcomas. However, tumors where no branching or only linear evolution was detected had fewer anatomic regions available for analysis (Fig. 1D). This indicated that a failure to detect collateral branching in a specific case could be due to a relative undersampling because of tumor necrosis or other reactive changes, leading to loss of viable tumor cells that could represent collateral lineages.

Figure 1.

Sample numbers, tumor size, and phylogenetic patterns. A, Three fundamental types of phylogenetic branching illustrated by fish plots showing genetic variation over evolutionary time (top row) and the resulting phylogenetic trees (bottom row). Linear branching infers a straight line of descent from an ancestral tumor population (founder genome detected in region P1) via mutations x to an intermediate descendant (detected in region P2) and further via the additional mutations y to an ultimate descendant population (detected in region P3). Collateral branching implies evolution of parallel branches from a common ancestor via private mutations x, y, and z, distinguishing populations detected in P1, P2, and P3, respectively. Mixed branching refers to combinations of linear and collateral branches in any given order. B, Examples of linear and mixed branching from two Wilms tumors (WT). Sampled regions are marked on macroscopic images adapted from Karlsson and colleagues (5). In WT16, the three sampled areas (P1–P3) surround a cystic necrotic center and share a common deletion in the long arm of chromosome 11 (11q-), the only aberration detected in all cells in P1. P1 also contains a subclone (P1s) sharing a series of additional structural changes (3q+, 4q-) and trisomies (+7, +9, etc.) with all cells in P2 and P3. These two regions in addition contain a clonal gain of the long arm of chromosome 1 (1q). There is also a subclone in P3 (P3s), having gain of 6p (6p+) and homozygous loss of the tumor suppressor JADE1. Contrasting this linear accumulation of mutation in WT07, the five tumor regions (P1–P5) available for sampling exhibit collateral branching into three subclones (P4s, P5s, P1s1) with unique features including two different deletions of 11q distinguished by different breakpoints (a and b; details in Supplementary Table 2) and loss of the AMER1 tumor suppressor. In P1, this is followed by linear accumulation of extra whole chromosomes in the subclone P1s2. Branch lengths correspond to the number of detected copy-number aberrations (CNA). C, Relative distribution of the three types of branching among primary tumors from Wilms tumor (WT), neuroblastoma (NB), and rhabdomyosarcoma (RMS). Percentages within brackets correspond to the fraction of tumors with branching evolution either as the only mode or together with linear evolution. D, Tumors where collateral branching (solely or mixed) was detected had more viable tumor regions available for analysis than tumors where only linear branching or no branching was detected. P values by two-sided Mann–Whitney U test. Medians are demarcated by black bars and colored numerals. E, Confidence intervals of correlation coefficients from linear regression (r values as black horizontal bars surrounded by boxes of 95% confidence interval) showing the correlation between sample numbers and maximum tumor diameters on the one hand, and mean branch length, the number of detected subclonal branches, and the number of detected clonal aberrations on the other hand. F, The number of detected branches in each tumor leading up to subclonal copy number aberrations shown as a function of the number of regions available for analysis; medians are denoted by black bars and colored numerals.

Figure 1.

Sample numbers, tumor size, and phylogenetic patterns. A, Three fundamental types of phylogenetic branching illustrated by fish plots showing genetic variation over evolutionary time (top row) and the resulting phylogenetic trees (bottom row). Linear branching infers a straight line of descent from an ancestral tumor population (founder genome detected in region P1) via mutations x to an intermediate descendant (detected in region P2) and further via the additional mutations y to an ultimate descendant population (detected in region P3). Collateral branching implies evolution of parallel branches from a common ancestor via private mutations x, y, and z, distinguishing populations detected in P1, P2, and P3, respectively. Mixed branching refers to combinations of linear and collateral branches in any given order. B, Examples of linear and mixed branching from two Wilms tumors (WT). Sampled regions are marked on macroscopic images adapted from Karlsson and colleagues (5). In WT16, the three sampled areas (P1–P3) surround a cystic necrotic center and share a common deletion in the long arm of chromosome 11 (11q-), the only aberration detected in all cells in P1. P1 also contains a subclone (P1s) sharing a series of additional structural changes (3q+, 4q-) and trisomies (+7, +9, etc.) with all cells in P2 and P3. These two regions in addition contain a clonal gain of the long arm of chromosome 1 (1q). There is also a subclone in P3 (P3s), having gain of 6p (6p+) and homozygous loss of the tumor suppressor JADE1. Contrasting this linear accumulation of mutation in WT07, the five tumor regions (P1–P5) available for sampling exhibit collateral branching into three subclones (P4s, P5s, P1s1) with unique features including two different deletions of 11q distinguished by different breakpoints (a and b; details in Supplementary Table 2) and loss of the AMER1 tumor suppressor. In P1, this is followed by linear accumulation of extra whole chromosomes in the subclone P1s2. Branch lengths correspond to the number of detected copy-number aberrations (CNA). C, Relative distribution of the three types of branching among primary tumors from Wilms tumor (WT), neuroblastoma (NB), and rhabdomyosarcoma (RMS). Percentages within brackets correspond to the fraction of tumors with branching evolution either as the only mode or together with linear evolution. D, Tumors where collateral branching (solely or mixed) was detected had more viable tumor regions available for analysis than tumors where only linear branching or no branching was detected. P values by two-sided Mann–Whitney U test. Medians are demarcated by black bars and colored numerals. E, Confidence intervals of correlation coefficients from linear regression (r values as black horizontal bars surrounded by boxes of 95% confidence interval) showing the correlation between sample numbers and maximum tumor diameters on the one hand, and mean branch length, the number of detected subclonal branches, and the number of detected clonal aberrations on the other hand. F, The number of detected branches in each tumor leading up to subclonal copy number aberrations shown as a function of the number of regions available for analysis; medians are denoted by black bars and colored numerals.

Close modal

To probe this potential confounder from sample numbers further, we then subdivided phylogenetic branches into those leading up to clonal and subclonal aberrations, respectively. Of the detected branches, 65% (137/212 also including linear branches) led to subclones only, i.e., tumor cell populations encompassing less than 90% of the entire cancer cell population in all samples where it was found. The remaining 35% (75/212) of branches led to populations that were clonal in at least one sample, i.e., present in 90% or more of tumor cells in that sample. The number of subclonal branches found in each tumor showed a significant but modest (r2 = 0.16) positive correlation to the number of regions available for sampling (Fig. 1E and F), whereas neither the number of clonal branches nor the mean branch length (aberration burden) was influenced by sample numbers. Tumor diameter did not correlate to phylogenetic parameters. We conclude that collateral branching, often mixed with linear features, appears to be the typical evolutionary mode of the three tumor types analyzed. Its absence may reflect undersampling, and comparison of branch numbers between tumor phylogenies must therefore be strictly normalized against sample numbers to avoid methodological bias.

Low- and high-risk tumors have distinct phylogenetic features

There was a wide variation in phylogenetic features observed among the 49 tumors where intratumoral genetic variation was detected (Fig. 2AL). To assess differences in phylogeny between tumor entities free from biases due to sampling and analysis platforms, we restricted the data to allelic imbalances and copy-number aberrations detected in primary tumor regions of the three main diagnostic tumor types. This analysis showed that Wilms tumors overall have a shorter stem length than neuroblastomas and rhabdomyosarcomas (Supplementary Table S1; P = 0.009; two-sided Mann–Whitney U test), well in accordance with a predominance of epigenetic rather than genetic factors underlying their inception (17). There were no other significant differences in phylogenetic features among the three diagnostic entities. Instead there was a marked variation in phylogenetic features within each group.

Figure 2.

Phylogenetic trees and iteration of canonical events. A–D, Representative phylogenetic trees based on maximum parsimony from Wilms tumors (WT) of intermediate-risk (IRH), blastemal type (BTH), and diffuse anaplastic (DAH) histology. P, samples from primary tumors; M samples from metastases; R, samples from local relapses. Main clones are annotated as red circles and subclones as black circles within blue circles as described in Supplementary Fig. S1. Scale bars indicate the stem/branch length of a single mutation/allelic imbalance/copy-number alteration (M). Plus (+) and minus (–) signs denote gains and losses of chromosomes or chromosomal segments. ICAs are denoted by the same gene/chromosome segment denoted by blue and red text in stem and branches, respectively. E–H, Trees representing two low-risk favorable (FAV; E–F), one high-risk MYCN-amplified (MNA; G), and one high-risk neuroblastoma (NB) dominated by structural (STR; H) rearrangements. H, NB19 exhibits chromothripsis-like (CHL) changes in chromosome 5, with one set confined to the stem (:1) and two other subsets (:2 and :3) confined to each of two phylogenetic branches. One branch also has CHL involving chromosome 8, and the other branch has a homozygous deletion (–/–) of CDKN2A. I, A high-risk MYCN-amplified neuroblastoma with metastases (M1–M3) at presentation, followed by relapses (R1, R2). There is an activating ALK mutation in the stem, followed by amplification of the mutated oncogene in R2. J–L, Trees representing embryonal (EMB) and alveolar (ALV) rhabdomyosarcoma (RMS). RMS1 (J) and RMS6 (K) each exhibit a subclone (green line continued from stem) whose descendants span all sampled regions. In RMS6, these subclonal descendants show parallel evolution of CDKN2A deletions (colored circles). L, RMS7 shows multiple driver mutations in the stem, followed by tetraploidization (2n→4n) and clonal sweeps leading up to each sampled region.

Figure 2.

Phylogenetic trees and iteration of canonical events. A–D, Representative phylogenetic trees based on maximum parsimony from Wilms tumors (WT) of intermediate-risk (IRH), blastemal type (BTH), and diffuse anaplastic (DAH) histology. P, samples from primary tumors; M samples from metastases; R, samples from local relapses. Main clones are annotated as red circles and subclones as black circles within blue circles as described in Supplementary Fig. S1. Scale bars indicate the stem/branch length of a single mutation/allelic imbalance/copy-number alteration (M). Plus (+) and minus (–) signs denote gains and losses of chromosomes or chromosomal segments. ICAs are denoted by the same gene/chromosome segment denoted by blue and red text in stem and branches, respectively. E–H, Trees representing two low-risk favorable (FAV; E–F), one high-risk MYCN-amplified (MNA; G), and one high-risk neuroblastoma (NB) dominated by structural (STR; H) rearrangements. H, NB19 exhibits chromothripsis-like (CHL) changes in chromosome 5, with one set confined to the stem (:1) and two other subsets (:2 and :3) confined to each of two phylogenetic branches. One branch also has CHL involving chromosome 8, and the other branch has a homozygous deletion (–/–) of CDKN2A. I, A high-risk MYCN-amplified neuroblastoma with metastases (M1–M3) at presentation, followed by relapses (R1, R2). There is an activating ALK mutation in the stem, followed by amplification of the mutated oncogene in R2. J–L, Trees representing embryonal (EMB) and alveolar (ALV) rhabdomyosarcoma (RMS). RMS1 (J) and RMS6 (K) each exhibit a subclone (green line continued from stem) whose descendants span all sampled regions. In RMS6, these subclonal descendants show parallel evolution of CDKN2A deletions (colored circles). L, RMS7 shows multiple driver mutations in the stem, followed by tetraploidization (2n→4n) and clonal sweeps leading up to each sampled region.

Close modal

When comparing prognostic subtypes with high and low risks of relapse, respectively, two differences consistently emerged across all three tumor types (Fig. 3A; see Supplementary Materials and Methods for risk classification). First, branches were found to be shorter in the more prognostically favorable subtypes of Wilms tumor, neuroblastoma, and rhabdomyosarcoma, than in the corresponding high-risk types (Fig. 3B). Second, low-risk subtypes of Wilms tumor, neuroblastoma, and rhabdomyosarcoma exhibited a significantly lower number of clonal branching events per tumor sample than did the corresponding high-risk tumors (Fig. 3C and D). These differences were also present at comparison between prognostically favorable (intermediate-risk histology) and moderate-risk (blastemal type) Wilms tumors (Supplementary Fig. S2A–S2C). Even though the analysis was normalized to remove bias due to sampling variation among tumors, subclonal branch frequency still did not correlate consistently to prognostic subgroups.

Figure 3.

Phylogenetic features and prognostic subtypes. A, Generic phylogenetic trees for each prognostic tumor subtype, based on mean stem length, mean branch length, and branch numbers normalized by sample numbers, denoted by branch color to reflect the proportion of subclonal (green) to clonal (red) branches. Mean lengths of phylogenetic stems (s, gray lines) and branches (b, green lines for subclonal, red lines for clonal) are denoted by numbers as well as their relative scale. The median number of branches (n) is standardized (x4/sample number) to reflect that the overall median sample number was 4 per primary tumor. For Wilms tumor (WT), generic trees are presented for the low-risk “intermediate risk histology” (IRH), moderate-risk “blastemal type histology” (BTH), and the high-risk “diffuse anaplastic histology” (DAH); for neuroblastoma (NB), the prognostically low-risk favorable tumors (FAV), high-risk MYCN-amplified tumors (MNA), and high-risk tumors dominated by other structural chromosome rearrangements (STR); for rhabdomyosarcoma (RMS), the low-risk embryonal subtype (EMB) and the high-risk alveolar subtype (ALV). BD, Dot plots covering the aberration burden (absence of branching = 0) in each detected branch (B), as well as the number of subclonal (C) and clonal (D) phylogenetic branches normalized by the number of samples analyzed in each case in low- and high-risk tumors, respectively. Red bars reflect mean values, and P values were obtained by two-tailed Mann–Whitney U test. The significance threshold was Bonferroni adjusted to 0.017 (0.05/3), with significant differences denoted by red type P values. Only copy-number aberrations and allelic imbalances (SNP-A data) from primary tumors were included in scatter plots and significance testing to avoid bias due to skewed availability of material for WES. E, Relative distribution of canonical aberrations in the phylogenies of all tumors. P values for differences in the proportion of canonical aberrations present only in phylogenetic branches were obtained by two-sided Fisher exact test. F, SNP array log2 plots showing ICAs in NB19 and NB23, comprising a stepwise loss of CDKN2A through loss/deletion and a stepwise ALK gain of function through mutation/amplification, respectively. Plots for chromosomes 10 (NB19) and 1 (NB23) are included as references. G, Panorama of ICAs in Wilms tumor, neuroblastoma, and rhabdomyosarcoma, with the number of tumors exhibiting such events given below each pie chart. m, mutation; +, gain; −, loss; cnni, copy-number neutral event. H, Results of enrichment analysis of all copy-number aberrations involving stem gain followed by branch gain (gain–gain) or stem loss followed by branch loss (loss–loss). Of tumors with this series of events, around half (19/39) showed overlap between stem and branch aberrations, but only five cases showed a significant enrichment indicative of selection. I, Enrichment analyses results from six tumors, with open bars showing the expected frequency distribution of stem-branch overlap sizes in base pairs (bp) versus the actual size of overlap (red vertical lines). Blue/red text denotes segments/genes affected in tumors with a higher degree of stem-branch overlap than expected from a random distribution. P values were calculated directly from frequency distribution and denote the probability that an overlap of at least the detected size would occur by random distribution of gains or losses. Losses in NB20 exemplify an actual degree of overlap well in accordance with a random distribution, whereas the remaining tumors illustrate significantly enriched overlap of losses or gains.

Figure 3.

Phylogenetic features and prognostic subtypes. A, Generic phylogenetic trees for each prognostic tumor subtype, based on mean stem length, mean branch length, and branch numbers normalized by sample numbers, denoted by branch color to reflect the proportion of subclonal (green) to clonal (red) branches. Mean lengths of phylogenetic stems (s, gray lines) and branches (b, green lines for subclonal, red lines for clonal) are denoted by numbers as well as their relative scale. The median number of branches (n) is standardized (x4/sample number) to reflect that the overall median sample number was 4 per primary tumor. For Wilms tumor (WT), generic trees are presented for the low-risk “intermediate risk histology” (IRH), moderate-risk “blastemal type histology” (BTH), and the high-risk “diffuse anaplastic histology” (DAH); for neuroblastoma (NB), the prognostically low-risk favorable tumors (FAV), high-risk MYCN-amplified tumors (MNA), and high-risk tumors dominated by other structural chromosome rearrangements (STR); for rhabdomyosarcoma (RMS), the low-risk embryonal subtype (EMB) and the high-risk alveolar subtype (ALV). BD, Dot plots covering the aberration burden (absence of branching = 0) in each detected branch (B), as well as the number of subclonal (C) and clonal (D) phylogenetic branches normalized by the number of samples analyzed in each case in low- and high-risk tumors, respectively. Red bars reflect mean values, and P values were obtained by two-tailed Mann–Whitney U test. The significance threshold was Bonferroni adjusted to 0.017 (0.05/3), with significant differences denoted by red type P values. Only copy-number aberrations and allelic imbalances (SNP-A data) from primary tumors were included in scatter plots and significance testing to avoid bias due to skewed availability of material for WES. E, Relative distribution of canonical aberrations in the phylogenies of all tumors. P values for differences in the proportion of canonical aberrations present only in phylogenetic branches were obtained by two-sided Fisher exact test. F, SNP array log2 plots showing ICAs in NB19 and NB23, comprising a stepwise loss of CDKN2A through loss/deletion and a stepwise ALK gain of function through mutation/amplification, respectively. Plots for chromosomes 10 (NB19) and 1 (NB23) are included as references. G, Panorama of ICAs in Wilms tumor, neuroblastoma, and rhabdomyosarcoma, with the number of tumors exhibiting such events given below each pie chart. m, mutation; +, gain; −, loss; cnni, copy-number neutral event. H, Results of enrichment analysis of all copy-number aberrations involving stem gain followed by branch gain (gain–gain) or stem loss followed by branch loss (loss–loss). Of tumors with this series of events, around half (19/39) showed overlap between stem and branch aberrations, but only five cases showed a significant enrichment indicative of selection. I, Enrichment analyses results from six tumors, with open bars showing the expected frequency distribution of stem-branch overlap sizes in base pairs (bp) versus the actual size of overlap (red vertical lines). Blue/red text denotes segments/genes affected in tumors with a higher degree of stem-branch overlap than expected from a random distribution. P values were calculated directly from frequency distribution and denote the probability that an overlap of at least the detected size would occur by random distribution of gains or losses. Losses in NB20 exemplify an actual degree of overlap well in accordance with a random distribution, whereas the remaining tumors illustrate significantly enriched overlap of losses or gains.

Close modal

As expected from a higher frequency of clonal but not subclonal branching in high-risk tumors, the interregional but not the intraregional genetic variation was higher in high-risk than in low-risk tumor subtypes (Supplementary Fig. S2D–S2G). Taking the whole dataset of 56 tumors into account, tumors that relapsed at least once had a higher frequency of clonal branching per sample and a higher mean branch length, whereas there was no difference in stem length or in subclonal branch frequency (Supplementary Fig. S2H). Univariate analyses accordingly revealed lower survival rates for tumors with more clonal branches and longer branches (Supplementary Fig. S2I–S2J). Multivariable survival analysis was then performed, including risk group (blastemal type Wilms tumors treated as high-risk), diagnostic tumor type, age, tumor diameter, stem length, branch frequencies, and mean branch length per tumor. Of these parameters, risk group and mean branch length were independent predictors (P ≤ 0.01) of relapse-free survival (Supplementary Fig. S2K). In summary, our data showed that tumors in which new clones were more frequent and had more de novo mutations were also associated with a higher risk of relapse, independent of diagnostic tumor type.

Phylogenetic branching contributes to molecular pathogenesis

We then evaluated whether phylogenetic branching could contribute to the acquisition of genetic aberrations of importance to pathogenesis or was merely a side effect of an elevated mutation rate. Because this analysis did not include case-to-case comparisons, we used all available genetic profiles from the 56 tumors, including those from metastases and relapses (Supplementary Table S2) and compared it with a set of genomic imbalances and point mutations that have been reported as highly recurrent in each of Wilms tumor, neuroblastoma, and rhabdomyosarcoma (termed canonical aberrations; see Supplementary Materials and Methods for a precise list). Out of these canonical aberrations, a total of 13, 11, and 5, respectively, were identified in at least one tumor in the present dataset (Supplementary Table S4). Of 173 identified instances of canonical aberrations, 77 were confined to the phylogenetic stem and 64 were found in branches only (Fig. 3E), exemplified by TP53 mutation/loss in Fig. 2D and CDKN2A deletions in Fig. 2K. Well in accordance with the finding of relatively short stems in Wilms tumor, the proportion of canonical aberrations confined to branches was higher in this tumor type than in the other two. Notably, 32 of 173 canonical aberrations were present twice in the same lineage of stem and branches. Such iterated canonical aberrations (ICA) were found in 50% to 38% of tumors of each diagnostic subtype (Fig. 3F and G).

In total, 22% (7/32) of ICAs involved additional gain or loss of well-established driver genes, including TP53, AMER1, MYCN, ALK, and CDKN2A (examples in Fig. 2H and Fig. 3F NB19; Fig. 2I and Fig. 3F NB23). However, 78% (25/32) of detected ICAs could not be linked directly to specific canonical driver genes (Fig. 3G), but consisted of successive gains or losses of whole or parts of chromosomes, e.g., a successive accumulation of copies of 7q in Wilms tumor (Fig. 2C; ref. 18), and 17q in neuroblastoma (Fig. 2G; ref. 19). To assess whether these ICAs were coincidental results of chromosomal instability during branching evolution (genetic drift) or truly statistically overrepresented events in the specific tumor regions, we employed a permutation-based analysis. For each tumor, chromosomal segments of the same length as those gained or lost in branches were stochastically reshuffled 10,000 times with respect to genomic location, after which, the relative frequency of different degrees of overlap was calculated (see Supplementary Materials and Methods). This approach successfully identified canonical and noncanonical bona fide driver events as statistically enriched by overlap in 5 of 19 tumors with ICAs produced by repeated aberration of the same segment in the stem and subsequent branches (Fig. 3H and I). However, none of the canonical large-scale chromosomal rearrangements were more frequently iterated than expected from a stochastic distribution of rearranged segments in branches compared with the stem. We conclude most repeated rearrangement of the same large-scale chromosomal segment within a tumor phylogeny can be explained by genetic drift, but that phylogenetic branching may still significantly contribute to pathogenesis by either providing the first emergence of a canonical aberration (64/173 instances of canonical aberrations, 36%) or by further modifying the copy number of driver genes (5/19 stem-branch ICAs, 26%).

SCS corroborates collateral branching as the main mode of evolution

Estimating phylogenetic features at the single-cell level circumvents potential biases from sample numbers when creating clonal phylogenies. In addition, it does not suffer from the methodological weaknesses present in clonal deconvolution of data from bulk samples (20), e.g., ambiguity in whether small subclones are nested within each other or not (Supplementary Fig. S1C). We therefore subjected eight representative tumor samples from eight tumors (three neuroblastomas, three Wilms tumors, and two rhabdomyosarcomas) to single-cell flow-sorting followed by single-cell low-pass (0.01–0.02x) whole-genome sequencing. In total, 547 libraries were of sufficient quality for conclusive SCS analysis (50–93 cells per case). This revealed a broad range of intercellular genetic diversity among the tumors (Supplementary Figs. S3, S4A, and S4B). Notably, all analyzed tumors revealed collateral branching at construction of phylogenies based on single cells, including copy-number variation in canonical aberrations and oncogenes (Fig. 4A and B; ref. 21).

Figure 4.

Single-cell phylogenies. A, Maximum parsimony trees from high-risk tumors including three neuroblastomas (NB18, NB19, and NB20) and one rhabdomyosarcoma (RMS7). The number of copy-number aberrations (CNA) in the stem is specified. Lower case letters at the end of branches correspond to the location of detected clones, i.e., two or more cells with identical CNA profiles, whereas branches not ending in letters correspond to genotypes detected in single cells. Stem aberrations characteristic of each tumor type are given in brown text, tetraploidization by 4n, and ICAs in branches are specified in blue type. Copy numbers (n =) for specific chromosomes or segments are exemplified by plots of sequence read numbers. In the TP53-mutated RMS7, where no cells show identical profiles, the complex scenario of variation is exemplified by variable copy numbers of the FNDC3B oncogene (21). See Supplementary Fig. S1 for examples of full single-cell copy-number aberration profiles. B, Maximum parsimony trees from low-risk tumors including one rhabdomyosarcoma (RMS6) and three Wilms tumors (WT6, WT9, and WT14), with annotations as in A. C, Branch lengths measured by the total number of copy-number aberrations included in each branch (circles), with medians denoted by red lines and P values by two-tailed Mann–Whitney U test.

Figure 4.

Single-cell phylogenies. A, Maximum parsimony trees from high-risk tumors including three neuroblastomas (NB18, NB19, and NB20) and one rhabdomyosarcoma (RMS7). The number of copy-number aberrations (CNA) in the stem is specified. Lower case letters at the end of branches correspond to the location of detected clones, i.e., two or more cells with identical CNA profiles, whereas branches not ending in letters correspond to genotypes detected in single cells. Stem aberrations characteristic of each tumor type are given in brown text, tetraploidization by 4n, and ICAs in branches are specified in blue type. Copy numbers (n =) for specific chromosomes or segments are exemplified by plots of sequence read numbers. In the TP53-mutated RMS7, where no cells show identical profiles, the complex scenario of variation is exemplified by variable copy numbers of the FNDC3B oncogene (21). See Supplementary Fig. S1 for examples of full single-cell copy-number aberration profiles. B, Maximum parsimony trees from low-risk tumors including one rhabdomyosarcoma (RMS6) and three Wilms tumors (WT6, WT9, and WT14), with annotations as in A. C, Branch lengths measured by the total number of copy-number aberrations included in each branch (circles), with medians denoted by red lines and P values by two-tailed Mann–Whitney U test.

Close modal

ICAs were a common feature (6/8 cases) also in the SCS phylogenies, including accumulation of 17q copies in the branches of 3 of 3 neuroblastomas and further copy-number enrichment in 3 of 4 Wilms tumors and embryonal rhabdomyosarcomas of 7q+, +8, and +12. Similar to the MRS data, the four high-risk tumors had longer phylogenetic branches than did the low-risk tumors (Fig. 4C). SCS allowed a precise estimate of the cell numbers present in each clone in the sequenced cell population (Fig. 5A and B), making it possible to assess the Simpson index of diversity, i.e., the likelihood that two randomly sampled cells would have different copy-number profiles in each of the tumors (Fig. 5C), as well as the fraction of cells having a genotype only detected once (Fig. 5D). Both these estimates showed that cell-to-cell genetic diversity was significantly larger in the high- than in the low-risk tumors. The SCS data thus replicated the impression from MRS analysis that collateral branching rather than just the linear evolution is the main evolutionary modus operandi in the analyzed pediatric tumors and that high-risk tumors exhibit more frequent as well as more extensive phylogenetic branching than do low-risk tumors.

Figure 5.

Population structure based on SCS. A and B, Genotype distribution among cells in high- and low-risk tumors as derived by phylogenetic analysis. Diameters of blue circles correspond to relative clone sizes, whereas red lines correspond to single-cell genotypes. The total number of analyzed cells in each tumor is specified (N =), and the number of cells detected in each clone (a, b, c, etc.) is given in parentheses. Black lines mark the phylogenetic relationship between each clone, the stem, and the other clones as detailed in Fig. 4. C, The Simpson index of diversity of each tumor (circles) where medians are denoted by red lines and the P value was calculated by two-tailed Mann–Whitney U test. D, The proportion of cells in each tumor whose genotype was only detected once (singlet cells) as compared with the proportion of cells whose genotype was detected in two or more cells (clonal cells). P value was calculated by two-tailed Fisher exact test from absolute cell numbers.

Figure 5.

Population structure based on SCS. A and B, Genotype distribution among cells in high- and low-risk tumors as derived by phylogenetic analysis. Diameters of blue circles correspond to relative clone sizes, whereas red lines correspond to single-cell genotypes. The total number of analyzed cells in each tumor is specified (N =), and the number of cells detected in each clone (a, b, c, etc.) is given in parentheses. Black lines mark the phylogenetic relationship between each clone, the stem, and the other clones as detailed in Fig. 4. C, The Simpson index of diversity of each tumor (circles) where medians are denoted by red lines and the P value was calculated by two-tailed Mann–Whitney U test. D, The proportion of cells in each tumor whose genotype was only detected once (singlet cells) as compared with the proportion of cells whose genotype was detected in two or more cells (clonal cells). P value was calculated by two-tailed Fisher exact test from absolute cell numbers.

Close modal

Some childhood tumors are born bad, others grow bad

Mutations that occur either before or at the time of the last complete clonal sweep are likely to be present in all tumor cells in all tumor regions at the time of presentation/sampling, thus placing these in the stem of each tumor's phylogeny (Fig. 6A). Conversely, mutations that occur later during tumor growth are likely to be confined to a few samples or regions. Based on this, we made a simple dichotomy of mutations into early and late events, corresponding to those present in the phylogenetic stem versus those present in less than all regions, respectively. To prevent any bias due to differences in sample numbers among subgroups, we normalized the number of detected late aberrations by sample numbers.

Figure 6.

Reconstructing evolutionary history. A, Fish plot illustrating the subdivision into early and late aberrations, with the former being those present clonally in all regions (stem mutations) and the latter being aberrations present clonally or subclonally in less than all regions. Global branch mutations, i.e., those present in all regions but not always clonally, were assigned to an intermediate stage and filtered out from further analysis. B–D, The number of copy-number aberrations (CNA) in Wilms tumor, neuroblastoma, and rhabdomyosarcoma assigned to early and late phases, subdivided into charts covering whole-chromosome anomalies (WCA) and structural rearrangements (STR), respectively. Each chart is in turn is subdivided into a field for each prognostic subgroup including for Wilms tumors intermediate-risk histology (IRH), the moderate-risk blastemal type histology (BTH), and the high-risk diffuse anaplastic histology (DAH); for neuroblastoma, low-risk favorable tumors (FAV), high-risk MYCN-amplified (MNA), and high-risk tumors with structural rearrangements without MYCN amplification (STR); for rhabdomyosarcoma, the low-risk embryonal (EMB) and high-risk alveolar (ALV) subtypes. Box plots, 25 and 75 percentiles; whiskers, 95% confidence intervals; single points, outliers. P values were derived by two-sided Mann–Whitney U test between low/moderate and high-risk groups. The significance threshold was Bonferroni adjusted to 0.004 (0.05/12), with significant differences denoted by red type P values. Late anomalies were normalized to a standard of four samples per tumor. SNVs/indels are not included due to their sparsity in the data. Pie charts flanking each box plot display the frequency of cases having (blue) or not having (gray) one or more canonical prognostic aberration, with the location of the pie indicating whether the aberration in question is a WCA or STR, which prognostic subgroups it signifies, and its distribution between early and late anomalies. For Wilms tumors, the only canonical prognostic aberration was deletion/mutation (STR) of TP53, which was found only in DAH and only among late copy-number aberrations. For neuroblastomas and rhabdomyosarcomas, the canonical prognostic anomalies are detailed in Supplementary Materials and Methods. E–G, Distribution among early and late phases of cases having canonical anomalies, where green/red bars denote markers of prognostic subtypes and brown/beige bars denote characteristic diagnostic but prognostically insignificant aberrations. Aberrations for which there is a significant overrepresentation among early or late anomalies are denoted by the absolute number of tumors (black type) where the aberration was present early and late, and the P value (red type) from two-tailed Fisher exact test. Chromosomal aberrations are denoted as copy-number neutral imbalances (cnni), gains (+), losses (–), whereas SNV/indels affecting genes are denoted by the name of the gene. t., translocations causing gene fusions; MNA, MYCN amplification. Only aberrations occurring in two or more tumors were included.

Figure 6.

Reconstructing evolutionary history. A, Fish plot illustrating the subdivision into early and late aberrations, with the former being those present clonally in all regions (stem mutations) and the latter being aberrations present clonally or subclonally in less than all regions. Global branch mutations, i.e., those present in all regions but not always clonally, were assigned to an intermediate stage and filtered out from further analysis. B–D, The number of copy-number aberrations (CNA) in Wilms tumor, neuroblastoma, and rhabdomyosarcoma assigned to early and late phases, subdivided into charts covering whole-chromosome anomalies (WCA) and structural rearrangements (STR), respectively. Each chart is in turn is subdivided into a field for each prognostic subgroup including for Wilms tumors intermediate-risk histology (IRH), the moderate-risk blastemal type histology (BTH), and the high-risk diffuse anaplastic histology (DAH); for neuroblastoma, low-risk favorable tumors (FAV), high-risk MYCN-amplified (MNA), and high-risk tumors with structural rearrangements without MYCN amplification (STR); for rhabdomyosarcoma, the low-risk embryonal (EMB) and high-risk alveolar (ALV) subtypes. Box plots, 25 and 75 percentiles; whiskers, 95% confidence intervals; single points, outliers. P values were derived by two-sided Mann–Whitney U test between low/moderate and high-risk groups. The significance threshold was Bonferroni adjusted to 0.004 (0.05/12), with significant differences denoted by red type P values. Late anomalies were normalized to a standard of four samples per tumor. SNVs/indels are not included due to their sparsity in the data. Pie charts flanking each box plot display the frequency of cases having (blue) or not having (gray) one or more canonical prognostic aberration, with the location of the pie indicating whether the aberration in question is a WCA or STR, which prognostic subgroups it signifies, and its distribution between early and late anomalies. For Wilms tumors, the only canonical prognostic aberration was deletion/mutation (STR) of TP53, which was found only in DAH and only among late copy-number aberrations. For neuroblastomas and rhabdomyosarcomas, the canonical prognostic anomalies are detailed in Supplementary Materials and Methods. E–G, Distribution among early and late phases of cases having canonical anomalies, where green/red bars denote markers of prognostic subtypes and brown/beige bars denote characteristic diagnostic but prognostically insignificant aberrations. Aberrations for which there is a significant overrepresentation among early or late anomalies are denoted by the absolute number of tumors (black type) where the aberration was present early and late, and the P value (red type) from two-tailed Fisher exact test. Chromosomal aberrations are denoted as copy-number neutral imbalances (cnni), gains (+), losses (–), whereas SNV/indels affecting genes are denoted by the name of the gene. t., translocations causing gene fusions; MNA, MYCN amplification. Only aberrations occurring in two or more tumors were included.

Close modal

We then compared the total numbers of early and late aberrations among the three tumor types and found no consistent differences. In contrast, there were again clear differences among prognostic subgroups within each tumor type (Fig. 6BG). In Wilms tumor, it is well known that high-risk diffuse anaplastic Wilms tumors compared with moderate- and low-risk tumors contain overall a higher burden of copy-number alterations (22). Our temporal analysis showed that this difference in aberration burden, including both structural anomalies and whole chromosome aberrations, emerged late and was not present in the founder genome (Fig. 6B). Analysis of canonical aberrations was remarkably consistent with this scenario. Mutations in TP53 and deletions of the corresponding part of 17p—phenomena closely associated with high-risk tumors and anaplasia—were late events in all five diffuse anaplastic tumors investigated (Fig. 6E). In contrast, copy-number neutral imbalance of 11p containing the WT1 and H19/IGF2 loci, a common aberration in all prognostic subtypes of Wilms tumors, typically emerged early.

In contrast to the scenario observed in Wilms tumors, low-risk neuroblastomas differed from the two high-risk neuroblastoma subtypes (MYCN amplified and structurally rearranged nonamplified) by having their characteristic profile of whole-chromosome anomalies occur early (Fig. 6C). Similarly, the structural aberrations signifying high-risk neuroblastomas emerged early in high-risk cases, whereas such aberrations only appeared as late anomalies in low-risk tumors. Well in accordance with this scenario, whole chromosome 17 gain and partial gain of 17q through structural rearrangement, signifying low-risk and high-risk neuroblastomas respectively (23), were overrepresented among early aberrations (Fig 6F). So was MYCN amplification, a long-established, strong prognostic predictor (24).

Rhabdomyosarcomas showed a scenario similar to neuroblastoma. Here, the embryonal subtype is denoted by having multiple alterations in whole chromosome copy numbers, whereas the alveolar subtype is denoted by multiple structural aberrations, including translocations involving the FOXO1 gene (25). This contrasting prevalence of numerical versus structural aberrations was reflected among early aberrations, with only 1 of 6 embryonal rhabdomyosarcomas containing any early structural aberrations but all 6 containing early numerical aberrations (Fig 6D and G). Similarly, the two cases of alveolar rhabdomyosarcoma contained FOXO1 rearrangements as well as a large number of structural changes among their early anomalies. In summary, our analysis showed that the genomic signatures that differentiate prognostic subgroups in neuroblastomas and rhabdomyosarcomas emerge earlier than they do in Wilms tumors, where the complex genomic features of diffuse anaplastic tumors occur closer to the point of surgical sampling.

The present study explores the phylogenetic features of the three most common solid malignancies in children outside the central nervous system. It is extensively based on a dataset from which we previously reported that pediatric cancers display a repertoire of four evolutionary trajectories, i.e., subclonal variation, clonal coexistence, regional clonal sweeps, and regional evolutionary explosions (5). These trajectories describe how different cancer cell clones in the same patient can compete or share space with each other, but they do not inform on the ancestral relationship between clones. In contrast, the present study analyzes clonal ancestry and the evolutionary history of mutations important from a clinical diagnostic and/or prognostic perspective. In tumor phylogenetics, it is critical to guard against bias due to sample numbers. Analysis of too few samples may lead to an underestimation of the number of phylogenetic branches and, as a consequence, produce a bias toward longer stems. Indeed, the present study highlighted that sample numbers especially influence the number of subclonal branches detected and the probability of missing collateral branching. As countermeasures against sample bias, we here normalized our data against sample numbers and also included an equal average number of samples from each tumor type. For further validation, we also performed single-cell karyotyping based on which phylogenetic analysis can be performed without the additional step of clonal deconvolution (26). Although only eight cases could be analyzed by SCS, this revealed results consistent with the overall conclusion from MRS data, with more frequent and more divergent branching in high-risk than in low-risk tumors. An alternative approach to our sample number standardization and single-cell validation would be to include strictly the same number of samples for each tumor subtype. However, this would create other biases, such as a failure to include small tumors that provide few samples or, conversely, missing biologically essential phylogenetic branches in large tumors due to relative undersampling.

Understanding what factors drive a high-risk phenotype in pediatric solid tumors is essential for finding new ways to cure these neoplasms for which the greatest threat to patient survival is relapse after chemotherapy. Relapse tumors have, in the vast majority of cases, a genomic profile distinct from that of their corresponding primary tumors (5, 7–12). This indicates that the founder lineage for any such relapse has been filtered through a Darwinian selection process that enriches for phenotypes that can withstand cytotoxic drugs. From this point of view, it is reasonable to assume that a factor associated to the risk of relapse would be evolvability per se, i.e., a cancer's capacity to adapt to changing circumstances through natural selection of genetic variants. The present study highlights that very simple phylogenetic parameters, such as branch length and branch frequency, could act as proxy markers for such overall cancer cell evolvability. Notably, evolvability is not equal to genetic instability. Genetically unstable tumor cells are defined by a high mutation rate only, not taking into account the fitness of generated genetic variants. In contrast, our phylogenetic approach is set to detect novel genetic variants sufficiently fit to form clones (bulk analysis based on MRS) or at least survive long enough to be sequenced (SCS analysis).

One should remember that our approach is based on an inference of events over time from samples obtained at just one or a few time points. It employs standard assumptions from population genetics such as low rates of back mutations (revertants) and scarcity of parallel (convergent) evolution. The trees presented here are thus models of past events, not objective point-by-point measurements. Nevertheless, our data enabled a systematic phylogenetic inventory of genetic aberrations characteristic of the tumor types in question, many of which are used for diagnostic or prognostic purposes. We found that a significant number of such canonical aberrations were confined to branches and thus were often detected only in some tumor regions. This was most pronounced in Wilms tumors where only around one third (27/74) of the detected canonical aberrations were present in all cells in all samples (stem anomalies). This high intratumoral variability genetic markers could explain why relatively few clinically useful molecular prognosticators have been detected in Wilms tumor despite decades of biomarker research over large datasets (27). The present study suggests that quantifying genetic variation per se could be a more successful strategy in Wilms tumor than focusing on specific genetic markers. This approach based on analyzing a generic pattern of cancer cell evolution, rather than a specific genetic alteration, may be suitable to accommodate the high diversity of somatic mutations present also in other pediatric cancers.

Even though our data indicated that the investigated high-risk tumors share a common phenotype of high genetic variability, there was a difference in the time point at which the mutational patterns predictive of relapse emerged. For neuroblastomas and rhabdomyosarcomas, the critical period in this respect was the time point leading up to the last clonal sweep, as reflected by the early appearance of characteristic high-risk aberrations such as MYCN amplification and structural rearrangement of 17q for the former, and the FOXO1 rearrangements for the latter. These mutations were concomitant to a large number of other structural chromosomal rearrangements, typical of high-risk tumors of both subtypes. In contrast, the anaplastic Wilms tumors started out largely as other Wilms tumors, often with the characteristic but nonprognostic copy-number neutral imbalance of the IGF2/H19 and WT1 loci in 11p. In high-risk tumors, this was subsequently followed by a regional emergence of complex chromosome changes that signify diffuse anaplasia. This concurrence in time between bursts of chromosomal instability in high-risk tumors of all three tumor types indicates that a set of specific driver mutations could be directly responsible for the shared feature of frequent and highly divergent phylogenetic branching. Indeed, the canonical aberrations associated to relapse have connections to DNA damage response and repair. This is most obvious for TP53 mutations in Wilms tumors that have been strongly linked to the emergency of high-risk histology (4), but also true for MYCN amplification and structural gain if 17q found in high-risk neuroblastomas, as well as the PAX3/7-FOXO1 gene fusion of alveolar rhabdomyosarcoma (28–34). However, one must also consider that children with high-risk disease of all three tumor types are generally older than those with low-risk disease. Hence, more extensive branching could to some extent also be a consequence of high-risk tumors having passed through more mitotic divisions since the onset of monoclonal expansion.

In all, our data indicate that evolvability of the cancer cell genome must be taken into account at tumor sampling for the purpose of personalized medicine and must be anticipated when designing novel therapeutic strategies for high-risk childhood cancers.

No potential conflicts of interest were disclosed.

Conception and design: N. Andersson, A. Valind, D. Gisselsson

Development of methodology: N. Andersson, B. Bakker, A. Valind, D.C.J. Spierings, D. Gisselsson

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): B. Bakker, D.C.J. Spierings, F. Foijer, D. Gisselsson

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): N. Andersson, B. Bakker, J. Karlsson, A. Valind, L. Holmquist Mengelbier, D.C.J. Spierings, D. Gisselsson

Writing, review, and/or revision of the manuscript: J. Karlsson, A. Valind, L. Holmquist Mengelbier, D.C.J. Spierings, F. Foijer, D. Gisselsson

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): B. Bakker, J. Karlsson, L. Holmquist Mengelbier, D. Gisselsson

Study supervision: F. Foijer, D. Gisselsson

The authors acknowledge technical support from the Science for Life Laboratory, the Knut and Alice Wallenberg Foundation, the National Genomics Infrastructure founded by the Swedish Research Council, and Uppsala Multidisciplinary Center for Advanced Computational Science for assistance with massively parallel sequencing and access to the UPPMAX computational infrastructure. We would also like to thank the Swegene Centre for Integrative Biology at Lund University (SCIBLU) for assistance. This study was supported by grants to D. Gisselsson from the Swedish Research Foundation (2016-01022), the Swedish Cancer Society (CAN2015/284), the Swedish Childhood Cancer Foundation (PR2016-024, NCP2015-0035), the Crafoord Foundation, the Royal Physiographic Society, and the Gunnar Nilsson Cancer Foundation.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.

1.
Valentijn
LJ
,
Koster
J
,
Zwijnenburg
DA
,
Hasselt
NE
,
van Sluis
P
,
Volckmann
R
, et al
TERT rearrangements are frequent in neuroblastoma and identify aggressive tumors
.
Nat Genet
2015
;
47
:
1411
4
.
2.
Ackermann
S
,
Cartolano
M
,
Hero
B
,
Welte
A
,
Kahlert
Y
,
Roderwieser
A
, et al
A mechanistic classification of clinical phenotypes in neuroblastoma
.
Science
2018
;
362
:
1165
70
.
3.
Vujanic
GM
,
Gessler
M
,
Ooms
A
,
Collini
P
,
Coulomb-l'Hermine
A
,
D'Hooghe
E
, et al
The UMBRELLA SIOP-RTSG 2016 Wilms tumour pathology and molecular biology protocol
.
Nat Rev Urol
2018
;
15
:
693
701
.
4.
Wegert
J
,
Vokuhl
C
,
Ziegler
B
,
Ernestus
K
,
Leuschner
I
,
Furtwangler
R
, et al
TP53 alterations in Wilms tumour represent progression events with strong intratumour heterogeneity that are closely linked but not limited to anaplasia
.
J Pathol Clin Res
2017
;
3
:
234
48
.
5.
Karlsson
J
,
Valind
A
,
Holmquist Mengelbier
L
,
Bredin
S
,
Cornmark
L
,
Jansson
C
, et al
Four evolutionary trajectories underlie genetic intratumoral variation in childhood cancer
.
Nat Genet
2018
;
50
:
944
50
.
6.
Skapek
SX
,
Anderson
J
,
Barr
FG
,
Bridge
JA
,
Gastier-Foster
JM
,
Parham
DM
, et al
PAX-FOXO1 fusion status drives unfavorable outcome for children with rhabdomyosarcoma: a children's oncology group report
.
Pediatr Blood Cancer
2013
;
60
:
1411
7
.
7.
Mengelbier
LH
,
Karlsson
J
,
Lindgren
D
,
Valind
A
,
Lilljebjorn
H
,
Jansson
C
, et al
Intratumoral genome diversity parallels progression and predicts outcome in pediatric cancer
.
Nat Commun
2015
;
6
:
6125
.
8.
Eleveld
TF
,
Oldridge
DA
,
Bernard
V
,
Koster
J
,
Daage
LC
,
Diskin
SJ
, et al
Relapsed neuroblastomas show frequent RAS-MAPK pathway mutations
.
Nat Genet
2015
;
47
:
864
71
.
9.
Schramm
A
,
Koster
J
,
Assenov
Y
,
Althoff
K
,
Peifer
M
,
Mahlow
E
, et al
Mutational dynamics between primary and relapse neuroblastomas
.
Nat Genet
2015
;
47
:
872
7
.
10.
Padovan-Merhar
OM
,
Raman
P
,
Ostrovnaya
I
,
Kalletla
K
,
Rubnitz
KR
,
Sanford
EM
, et al
Enrichment of targetable mutations in the relapsed neuroblastoma genome
.
PLoS Genet
2016
;
12
:
e1006501
.
11.
Schleiermacher
G
,
Javanmardi
N
,
Bernard
V
,
Leroy
Q
,
Cappo
J
,
Rio Frio
T
, et al
Emergence of new ALK mutations at relapse of neuroblastoma
.
J Clin Oncol
2014
;
32
:
2727
34
.
12.
Brady
SW
,
Ma
X
,
Bahrami
A
,
Satas
G
,
Wu
G
,
Newman
S
, et al
The clonal evolution of metastatic osteosarcoma as shaped by cisplatin treatment
.
Mol Cancer Res
2019
;
17
:
895
906
.
13.
Staaf
J
,
Lindgren
D
,
Vallon-Christersson
J
,
Isaksson
A
,
Goransson
H
,
Juliusson
G
, et al
Segmentation-based detection of allelic imbalance and loss-of-heterozygosity in cancer cells using whole genome SNP arrays
.
Genome Biol
2008
;
9
:
R136
.
14.
Rasmussen
M
,
Sundstrom
M
,
Goransson Kultima
H
,
Botling
J
,
Micke
P
,
Birgisson
H
, et al
Allele-specific copy number analysis of tumor samples with aneuploidy and tumor heterogeneity
.
Genome Biol
2011
;
12
:
R108
.
15.
van den Bos
H
,
Bakker
B
,
Taudt
A
,
Guryev
V
,
Colome-Tatche
M
,
Lansdorp
PM
, et al
Quantification of aneuploidy in mammalian systems
.
Methods Mol Biol
2019
;
1896
:
159
90
.
16.
Jolly
C
,
Van Loo
P
. 
Timing somatic events in the evolution of cancer
.
Genome Biol
2018
;
19
:
95
.
17.
Gadd
S
,
Huff
V
,
Walz
AL
,
Ooms
A
,
Armstrong
AE
,
Gerhard
DS
, et al
A children's oncology group and TARGET initiative exploring the genetic landscape of Wilms tumor
.
Nat Genet
2017
;
49
:
1487
94
.
18.
Hawthorn
L
,
Cowell
JK
. 
Analysis of wilms tumors using SNP mapping array-based comparative genomic hybridization
.
PLoS One
2011
;
6
:
e18941
.
19.
Decaesteker
B
,
Denecker
G
,
Van Neste
C
,
Dolman
EM
,
Van Loocke
W
,
Gartlgruber
M
, et al
TBX2 is a neuroblastoma core regulatory circuitry component enhancing MYCN/FOXM1 reactivation of DREAM targets
.
Nat Commun
2018
;
9
:
4866
.
20.
Alves
JM
,
Prieto
T
,
Posada
D
. 
Multiregional tumor trees are not phylogenies
.
Trends Cancer
2017
;
3
:
546
50
.
21.
Cai
C
,
Rajaram
M
,
Zhou
X
,
Liu
Q
,
Marchica
J
,
Li
J
, et al
Activation of multiple cancer pathways and tumor maintenance function of the 3q amplified oncogene FNDC3B
.
Cell Cycle
2012
;
11
:
1773
81
.
22.
Williams
RD
,
Al-Saadi
R
,
Natrajan
R
,
Mackay
A
,
Chagtai
T
,
Little
S
, et al
Molecular profiling reveals frequent gain of MYCN and anaplasia-specific loss of 4q and 14q in Wilms tumor
.
Genes Chromosomes Cancer
2011
;
50
:
982
95
.
23.
Bown
N
,
Cotterill
S
,
Lastowska
M
,
O'Neill
S
,
Pearson
AD
,
Plantaz
D
, et al
Gain of chromosome arm 17q and adverse outcome in patients with neuroblastoma
.
N Engl J Med
1999
;
340
:
1954
61
.
24.
Brodeur
GM
,
Seeger
RC
,
Schwab
M
,
Varmus
HE
,
Bishop
JM
. 
Amplification of N-myc in untreated human neuroblastomas correlates with advanced disease stage
.
Science
1984
;
224
:
1121
4
.
25.
Arnold
MA
,
Anderson
JR
,
Gastier-Foster
JM
,
Barr
FG
,
Skapek
SX
,
Hawkins
DS
, et al
Histology, fusion status, and outcome in alveolar rhabdomyosarcoma with low-risk clinical features: a report from the children's oncology group
.
Pediatr Blood Cancer
2016
;
63
:
634
9
.
26.
Davis
A
,
Navin
NE
. 
Computing tumor trees from single cells
.
Genome Biol
2016
;
17
:
113
.
27.
Dome
JS
,
Graf
N
,
Geller
JI
,
Fernandez
CV
,
Mullen
EA
,
Spreafico
F
, et al
Advances in Wilms tumor treatment and biology: progress through international collaboration
.
J Clin Oncol
2015
;
33
:
2999
3007
.
28.
Bell
E
,
Premkumar
R
,
Carr
J
,
Lu
X
,
Lovat
PE
,
Kees
UR
, et al
The role of MYCN in the failure of MYCN amplified neuroblastoma cell lines to G1 arrest after DNA damage
.
Cell Cycle
2006
;
5
:
2639
47
.
29.
Agarwal
S
,
Milazzo
G
,
Rajapakshe
K
,
Bernardi
R
,
Chen
Z
,
Barberi
E
, et al
MYCN acts as a direct co-regulator of p53 in MYCN amplified neuroblastoma
.
Oncotarget
2018
;
9
:
20323
38
.
30.
Newman
EA
,
Chukkapalli
S
,
Bashllari
D
,
Thomas
TT
,
Van Noord
RA
,
Lawlor
ER
, et al
Alternative NHEJ pathway proteins as components of MYCN oncogenic activity in human neural crest stem cell differentiation: implications for neuroblastoma initiation
.
Cell Death Dis
2017
;
8
:
3208
.
31.
Newman
EA
,
Lu
F
,
Bashllari
D
,
Wang
L
,
Opipari
AW
,
Castle
VP
. 
Alternative NHEJ pathway components are therapeutic targets in high-risk neuroblastoma
.
Mol Cancer Res
2015
;
13
:
470
82
.
32.
Richter
M
,
Dayaram
T
,
Gilmartin
AG
,
Ganji
G
,
Pemmasani
SK
,
Van Der Key
H
, et al
WIP1 phosphatase as a potential therapeutic target in neuroblastoma
.
PLoS One
2015
;
10
:
e0115635
.
33.
Mercado
GE
,
Xia
SJ
,
Zhang
C
,
Ahn
EH
,
Gustafson
DM
,
Lae
M
, et al
Identification of PAX3-FKHR-regulated genes differentially expressed between alveolar and embryonal rhabdomyosarcoma: focus on MYCN as a biologically relevant target
.
Genes Chromosomes Cancer
2008
;
47
:
510
20
.
34.
Gryder
BE
,
Yohe
ME
,
Chou
HC
,
Zhang
X
,
Marques
J
,
Wachtel
M
, et al
PAX3-FOXO1 establishes myogenic super enhancers and confers BET bromodomain vulnerability
.
Cancer Discov
2017
;
7
:
884
99
.