University of Groningen Evolution of karyotype landscapes in cancer Bakker, Bjorn

(1)

University of Groningen

Evolution of karyotype landscapes in cancer

Bakker, Bjorn

DOI:

10.33612/diss.166886747

IMPORTANT NOTE: You are advised to consult the publisher's version (publisher's PDF) if you wish to cite from it. Please check the document version below.

Document Version

Publisher's PDF, also known as Version of record

Publication date: 2021

Link to publication in University of Groningen/UMCG research database

Citation for published version (APA):

Bakker, B. (2021). Evolution of karyotype landscapes in cancer. University of Groningen. https://doi.org/10.33612/diss.166886747

Copyright

Other than for strictly personal use, it is not permitted to download or to forward/distribute the text or part of it without the consent of the author(s) and/or copyright holder(s), unless the work is under an open content license (like Creative Commons).

Take-down policy

If you believe that this document breaches copyright please contact us providing details, and we will remove access to the work immediately and investigate your claim.

Downloaded from the University of Groningen/UMCG research database (Pure): http://www.rug.nl/research/portal. For technical reasons the number of authors shown on this cover page is limited to 10 maximum.

(2)

553979-L-bw-Bakker 553979-L-bw-Bakker 553979-L-bw-Bakker 553979-L-bw-Bakker Processed on: 25-2-2021 Processed on: 25-2-2021 Processed on: 25-2-2021

(3)

Processed on: 25-2-2021 PDF page: 188PDF page: 188PDF page: 188PDF page: 188

188 | Chapter 8

ABSTRACT

Darwinian evolution of tumour cells remains underexplored in childhood cancer. By whole-genome allelic copy number profiling of 274 tumour regions from 24 neuroblastomas, 24 Wilms tumours and eight rhabdomyosarcomas we find that 88% of these neoplasms exhibit genetic variation among primary tumour territories. This variability can affect diagnostically informative mutations and typically emerges through collateral phylogenetic branching that in subset of cases leads to amplification of oncogenes and loss of tumour suppressors. New branches emerge more frequently and carry a higher mutational load in high-risk than in low-risk tumours, a scenario corroborated by single cell sequencing of eight cases. By creating timelines from spatial variation, we show that the mutations most influencing relapse risk occur at the initiation of clonal expansion in neuroblastoma and rhabdomyosarcoma, while in Wilms tumour they are late events. Thus, from an evolutionary standpoint, some childhood cancers are born bad while others grow worse over time.

(4)

Extensive clonal branching shapes the evolutionary history of high-risk pediatric cancers | 189

INTRODUCTION

Cancer is one of the leading causes of death in children in developed countries. For most types of childhood cancer, the course leading up to fatal disseminated disease goes via a tumour that is first sensitive to chemotherapy, followed by relapsing disease with an insufficient response to treatment. The common denominator in these situations is an alteration of the tumour cell phenotype that enables neoplastic cells to survive cytotoxic treatment. 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 intratumoural genetic heterogeneity. A small number of recent publications have revealed that intratumoural genetic diversity is indeed a common phenomenon in pediatric tumours1-7_.

However, the principles of clonal dynamics resulting from such diversity in pediatric cancer is largely unknown. Most prominently, we still know very little about how the capacity of cancer cells to evolve relates to their clinical features.

Based on the histology and the genetic profile of the primary tumour, 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 tumours affecting young children (<18 months of age) with 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 tumours carrying structural chromosome changes often including deletions in 1p and 11q, gain of 17q, ATRX deletions, TERT rearrangements, or other mechanisms leading to telomere dysregulation8,9_{. Wilms tumours,}

on the other hand, are subdivided according to the current European protocol into the favorable intermediate-risk histological subtype, the blastemal histology subtype with a moderate risk of relapse, and the high-risk diffuse anaplastic histology10_{. While the two}

former subtypes have no specific genetic markers, diffuse anaplasia correlates closely to somatic inactivation of the TP53 tumour suppressor and the emergence of complex chromosomal alterations7,11_{. Finally, rhabdomyosarcomas are largely divided into the highly}

aggressive alveolar subtype having fusions of the FOXO1 gene to either PAX3 or PAX7 along with complex structural rearrangements and the more favorable embryonal subtype, which is fusion-negative and has a genotype dominated by altered copy numbers of whole chromosomes12_{. Despite this abundance of genetic markers for relapse risk, little is known}

on how they relate mechanistically to the process of clonal evolution under treatment. In fact, low risk tumours of neuroblastoma, Wilms tumour and rhabdomyosarcomas all

(5)

190 | Chapter 8

typically exhibit multiple whole chromosome aberrations (aneuploidy), which in adult tumours has been strongly linked to progression, metastasis and drug resistance13_{. Their}

generally excellent response to chemotherapy therefore remains enigmatic.

In the present study, we employ multiregional sampling (MRS) followed by whole-genome allelic copy number analysis supplemented by whole exome sequencing to investigate tumour cell phylogenies in prognostic childhood tumour subtypes. In particular we focus on each subtype’s capacity to generate tumour cells with heritable variation as measured by the frequency of phylogenetic branching into divergent clones and the extent to which such clones differ genetically from each other. We also use MRS to answer the question of when in each tumour’s evolutionary history the genetic hallmarks for each prognostic subtype emerge. Because bulk samples have clear limitations in correctly deconvoluting small clones, data from MRS were validated against phylogenies based on single cell whole genome sequencing.

RESULTS

From geographic variation to phylogeny

We recently employed MRS to resolve how four evolutionary trajectories determine intratumour genetic variation in childhood cancer7_{. Here, we use an extension of the same}

dataset (see Methods) to compare tumour cell phylogenies among prognostic tumour subtypes. Building robust cancer phylogenies requires an approach that detects sufficient amount of genetic variation among tumour 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 single nucleotide polymorphism (SNP) arrays on both frozen and formalin fixed paraffin embedded tumour material, hence unlocking information from diagnostic samples stored in pathology archives. In total, we analyzed whole-genome copy number and allelic imbalance profiles from 274 tumour samples (240 from primary, 34 from metastatic sites) from 24 neuroblastomas, 24 Wilms tumours and eight rhabdomyosarcomas, all with a tumour cell content >50% (Fig. 1; Supplementary Table 1). We obtained a median of 4 samples per patient without any significant difference among tumour types (3.5 for neuroblastoma, 4.0 for Wilms tumour, 4.5 for rhabdomyosarcoma). As a complement to copy number profiling, we performed

(6)

Extensive clonal branching shapes the evolutionary history of high-risk pediatric cancers | 191 whole-exome sequencing of 35 biopsies from 19 tumours, followed by targeted deep sequencing of 76 samples, allowing complementary scoring of single nucleotide variants (SNVs) and indels in a total of 19 tumours. Clone size estimations were performed by calculating the prevalence of chromosomal imbalances as a function of copy number and log2 ratio deviation relative to the largest clone7_{. For sequence alterations, clone prevalence}

in each sample was calculated from variant allele frequencies combined with allele copy numbers (see Methods). Genome profiles for individual clones (clonal deconvolution) were then inferred by clustering chromosomal alterations and sequence alterations with similar prevalence levels (Supplementary Table 2). This revealed intratumoural genetic heterogeneity in 88% (49/56) of the investigated tumours (Supplementary Table 3). Based on the distribution of genetic aberrations among clones and subclones, we then constructed tumour phylogenies, modelling the evolutionary history of each tumour.

Branching evolution is a common feature of pediatric cancer

In principle, tumour phylogenies can have three different configurations based on the relationship of the branches leading up to genetically distinct clones and subclones (Fig. 2a-b). In linear evolution (linear branching), branches are arranged in a straight lineage from the stem towards a series of clones that are direct descendants of each other, reflecting a cumulative accumulation of mutations. Alternatively, branches radiate in a collateral fashion from the founder genome, reflecting the presence of private mutations in each detected clone (collateral branching). Finally, there can be a combination of these two scenarios (mixed branching). A total of 212 branches were detected in samples from primary tumours in our dataset. Linear branching was observed as the only scenario in approximately one third of tumours, with an equal distribution (33%-41%) in Wilms tumours, neuroblastomas and rhabdomyosarcomas (Fig. 2c). However, phylogenies containing some aspect of collateral branching was the most common scenario.

Tumours where no branching or only linear branching was detected had fewer anatomic regions available for analysis (Fig. 2d). This indicated that an absence of collateral branching was due to a relative undersampling because of tumour necrosis or other reactive changes leading to loss of viable tumour cells that could represent collateral lineages. Of the detected branches, 65% (137/212) led to subclones only, i.e. tumour 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 tumour cells in that sample. The number of subclonal branches found in each tumour showed a significant but modest (r2_=0.16)

(7)

192 | Chapter 8

positive correlation to the number of regions available for sampling (Fig. 2e and f), while neither the number of clonal branches nor the mean branch length (aberration burden) was influenced by sample numbers (Fig. 2e). Tumour diameter did not correlate to phylogenetic parameters (Fig. 2e). These results led us to focus further analyses on aberration burden and clonal branching as these parameters appeared to be little influenced by the number of tumour regions available for analysis as compared to tree configurations and subclonal branching.

Low and high-risk tumours have distinct phylogenetic features

There was a wide variation in phylogenetic features observed among the 49 tumours where intratumoural genetic variation was detected (Fig. 3). When comparing phylogenies among the three main diagnostic tumour types, Wilms tumours were found to have shorter mean stem length than had neuroblastomas and rhabdomyosarcomas (Supplementary Table 1; p=0.009; two-sided Mann-Whitney U test), well in accordance with a predominance of epigenetic rather than genetic factors underlying their inception14_{. 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. When comparing prognostic subtypes with high and low risks of relapse, respectively, two differences consistently emerged over all three tumour types (Fig. 4a). First, branches were found to be shorter in the more prognostically favourable subtypes of Wilms tumour, neuroblastoma and rhabdomyosarcoma, than in the corresponding high-risk types (Fig. 4b; Supplementary Fig. 1a). Second, low-risk subtypes of Wilms tumour, neuroblastoma and rhabdomyosarcoma exhibited a significantly lower number of clonal branching events per sample than did the corresponding high-risk tumours (Fig. 4c-d; Supplementary Fig. 1b-c). Even though the analysis was normalized to remove bias due to sampling variation among tumours, subclonal branch frequency still did not correlate consistently to prognostic subgroups. In summary, we found that the emergence of new clones was more frequent and was associated with a higher load of de novo mutations in high-risk than in low-risk tumours.

(8)

Fig. 1 Methodological overview. a From each patient, multiple samples from the primary tumour (P1, P2..) as

well as samples from metastatic sites (M1, M2..) were obtained. b A median of four samples per tumour from 24 neuroblastomas (NB), 24 Wilms tumours (WT) and eight rhabdomyosarcomas (RMS) were subjected to whole genome genotyping by SNP-array (SNP-A), resulting in copy number aberration (CNA) profiles for all samples, complemented by whole exome sequencing (WES) to detect single nucleotide variants (SNVs) and small insertions/deletions (indels). c CNA data were combined with the variant allele frequencies of SNVs/ indels to perform clonal deconvolution (see Methods), where smaller clones (e.g. yellow and brown territories) were always assumed, to be nested within larger clones so as not to overcall phylogenetic branching7_{. d After}

deconvolution, phylogenetic analysis was performed, allocating mutations found in all tumour cells in all samples to the stem, and the remaining mutations to branches. Branches leading to clones that encompassed all (≥90%) cells in a sample were regarded as clonal, denoted by red large circles, whereas those encompassing <90% were regarded as subclonal, denoted by black small circles within larger blue circles. For example, the mutations forming the branch marked by a green point are clonal in the region P1 and subclonal in P2.

e The geographical distribution of mutations across samples was also used to infer a time line of mutational

acquisition. f Number of neuroblastomas (NB), Wilms tumours (WT) and rhabdomyosarcomas (RMS) from which multiple samples from each primary tumour (T) as well as from metastases (M) were analyzed by SNP array to obtain copy number aberration (CNA) profiles. A subset of these samples was also analyzed by whole exome sequencing followed by targeted deep sequencing to detect non-neutral single nucleotide variants (SNV) and small insertions/deletions (indels).

(9)

194 | Chapter 8

Both stems and branches contribute to putative driver events

Histology is often insufficient for sub-classification of pediatric tumours. Therefore, frequent and typical genomic aberrations have been systematically inventoried over many decades to aid in diagnostic work. However, the extent to which stems and branches, respectively, contribute to such diagnostic aberrations has been little investigated in pediatric tumours. To analyse this, we first used database searches supplemented with recent broad sequencing studies to identify a set of genomic imbalances and point mutations (here termed canonical aberrations) that have been repeatedly linked to each of Wilms tumour, neuroblastoma, and rhabdomyosarcoma (see Methods). Out of these canonical aberrations, a total of 13, 11 and 5, respectively, were identified in at least one tumour in the present dataset (Supplementary Table 4). This subset included both large-scale chromosomal alterations such as loss of heterozygosity of 11p in Wilms tumour, gain of 17q/trisomy 17 in neuroblastoma and trisomy 13 in rhabdomyosarcoma, as well as single-gene rearrangements such as AMER1 mutation/deletion in Wilms tumour, ALK mutations in neuroblastoma and FOXO1 rearrangements in rhabdomyosarcoma.

a P1 founder genome tumor regions 0 +0.6 -0.6 r (95% confidence interval)

samples - branch length samples - subcl. branches samples - clonal branches t. diameter - branch length t. diameter - subcl. branches t. diameter - clonal branches

Sampling, tumor size and branching

16

2

no. of sampled regions

Sampling and branching type

collat/mixed branching linear branching no branching p = 0.009 p = 0.04 d 1 cm 1 cm P2 P3 P1

tumor regionsP2 P3 P1tumor regionsP2 P3

Linear

branching Collateral branching branchingMixed

P1 P2 P3 P1 P2 P3 P1 P2 P3 b WT16 WT07 c mixed branching collateral branching linear branching Sampled regions WT NB RMS Branching types 33% 57% 55% 41% 38% 63% 1 CNA P1 P1s P2 P3 P3s 11q-Linear branching 3q+ +12 4q- +13 +7 +18 +9 1q+ JADE1 -/-6p+ Sampled

regions branchingMixed

P1-P5 P4s P5s P1 s1 P1s2 1CNA LOH 19p 3q-11q-A 11q-B 9q- AMER1-+6 +12x2 +13 e 4 ₃ 2 x x+y x y z x x+y x+z Subclonal branches f

no. of sampled regions2 3 4 5 6

detected subclonal branches 1

2 2

(10)

Extensive clonal branching shapes the evolutionary history of high-risk pediatric cancers | 195 ◀ Fig. 2 Sample numbers, tumour 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 ensuing phylogenetic trees (bottom row). Linear branching infers a straight line of descent from an ancestral tumour 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 tumours (WT). Sampled regions are marked on macroscopic images adapted from Karlsson et al.7_{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 tumour suppressor JADE1. Contrasting this linear accumulation of mutation, in WT07 the five tumour 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 tumour 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 Wilms tumour (WT), neuroblastoma (NB) and rhabdomyosarcoma (RMS). d Tumours where collateral branching (solely or mixed) was detected had more viable tumour regions available for analysis than tumours 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 tumour 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 tumour leading up to subclonal copy number aberrations showed as a function of the number of regions available for analysis; medians are denoted by black bars and colored numerals.

Out of a total of 173 identified instances of canonical aberrations, 77 were confined to phylogenetic stems and 64 were found in branches only (Figure 5a). Well in accordance with the finding of relatively short stems in Wilms tumour, the proportion of canonical aberrations confined to branches were higher in this tumour type than in the other two. Notably, 32/173 canonical aberrations were present twice, i.e. iterated, in the same lineage of stem and branches. Such iterated canonical aberration (ICAs) were found in 50-38% of tumours of each diagnostic subtype (Fig. 5b-c). A total 22% (7/32) of ICAs involved additional gain or loss of well-established driver genes, including TP53, AMER1, MYCN, ALK and CDKN2A. Examples of this type of ICA include the sequential elimination of both CDKN2A copies in NB19: first through loss of chromosome 9 in the stem and then by interstitial deletion in a branch (Fig. 3h and Fig. 5b NB19). Another example is the mutation ALK p.F1174L in the stem of NB23 followed by amplification in a branch (Fig. 3i and Fig. 5b NB23). This escalation of ALK copy number occurred under conventional chemotherapy treatment, not under ALK inhibitor treatment. There were in fact no trend

(11)

196 | Chapter 8

for differences in ICAs when the comparing samples obtained before or after treatment in any patient (ICAs ascertained before, after and over the treatment period in 9, 13, and 2 patients respectively). WT6 a 1q+ -1 _{1 M} 1 CDK12 I1209T -21 7q- 2q-B P2 P1 P3 1 M b WT9 LOH 1 1p +7,+8 B1 B2 d 1 M P2 P1 M P3 P4 TP53 p.R273H -17 WT21 P2 P4 P1 P3 1 M c WT14 MYCN p.P44L 2p24+ 7q+ 7q++ [DAH] [IRH] [IRH] [BTH] e 1 M NB3 [FAV] B1 B2 +1 +2 +17x2 +10 OR4K2 p.A163S f NB4 [FAV] 1 M 2n-4n -4 -9 1p-+17 17q+ B1 B2 B2 1 M g NB18 [MNA] MYCN amp 17q+ B P4 P3 P1 P2 M1M2 h NB19 [STR] 1 M M1 M2 B8 B7 B1-B6 CHL(5):1 -9 CDKN2A -/-CHL(5):2 CHL(8) CHL(5):3 main clone sub-clone j 2n-4n +2 1 M B2 B1 +15 +17 +7 +8 +13 RMS1 [EMB] k RMS6 [EMB] 1 M LOH 11p l RMS7 [ALV] 1 M P3 P4 B1A B1B B3 B2 P1 P5 P6 P2 M-LN P M-LU B CDKN2A deletions WNT8B p.C324R CDKN2A del ho TP53 p.R116W 17p- 2n->4n +2 +8 +10 +12 17p-17q+ +20 stem retained clonal sweep Subclone present in all samples i NB23 [MNA] 1 M MYCN amp P2 P1 R1M1-M3 R2 ALK p.F1 174L ALK amp Subclone present in all samples 17q++ Subclone present in all samples 17q++

Fig. 3 Phylogenetic trees and iteration of canonical events. a-d Representative phylogenetic trees based

on maximum parsimony from Wilms tumours (WTs) of intermediate risk (IRH), blastemal type (BTH), and diffuse anaplastic (DAH) histology. Main clones (red circles) and subclones (black circles within blue circles) are annotated as described in Fig 1. Scale bars indicate the stem/branch length of a single mutation (M). Plus (+) and minus (-) signs denote gains and losses of chromosomes or chromosomal segments. e-h Trees representing two low-risk favorable (e-f, FAV), one high-risk MYCN-amplified (g, MNA), and one high-risk neuroblastoma (NB) dominated by structural (h, STR) rearrangements. NB19 (h) 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 NB 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). RMS7 (l) shows multiple driver mutations in the stem, followed by tetraploidization (2n->4n) and clonal sweeps leading up to each sampled region.

(12)

a s=4.6 b=1.3 n=0.5 x4 WT Low-risk (IRH) s=4.2 b=2.1 n=1.0 x4 Moderate (BTH) High-risk(DAH) s=1.3 n=2.0 x4 b=3.8 NB s=13.4 b=1.1 n=0.5 x4 Low-risk

(FAV) High-risk(MNA) High-risk(STR)

s=14.0 b=3.3 n=1.5 x4 s=14.3 b=3.7 n=1.7 x4 RMS s=4.2 b=5.7 n=0.85 x4 Low-risk

(EMB) High-risk(ALV)

s=25.5 b=8.8 n=2.5 x4 stem aberrations subclonal branch clonal branch b WT NB RMS

copy number aberrations 0 35 0 15 0 42

IRH DAH FAV MNA/STR EMB ALV

0 1.5 0 4.0 0 2.5 0 1.5 0 0 2.5 p= 0.002 1.5 c d Branch lengths WT NB RMS

Subclonal branch frequency

Clonal branch frequency

WT NB RMS

branches per sample

p= 0.0001 p= 0.001

p= 0.009 p= 0.25 p= 0.71

p= 0.0001 p= 0.0001 p= 0.016

Fig. 4 Phylogenetic features and prognostic subtypes. a Generic phylogenetic trees for each prognostic

tumour 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, grey 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 tumour. For Wilms tumour (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 tumours (FAV), high-risk MYCN amplified tumours (MNA), and high-risk tumours dominated by other structural chromosome rearrangements (STR); for rhabdomyosarcoma (RMS) the low-risk embryonal subtype (EMB) and the high-risk alveolar subtype (ALV). b-d Scatter 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 tumours, 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 were included in scatter plots and significance testing to avoid bias due to skewed availability of material for whole exome sequencing.

In total, 78% (25/32) of detected ICAs could not be linked directly to specific canonical driver genes (Fig. 5c), but consisted of successive gains or losses of whole or parts of chromosomes, e.g. a successive accumulation of copies of 7q in Wilms tumour (Fig. 3c)15_{, and 17q in neuroblastoma (Fig. 3g)}16_{. To assess whether these ICAs were}

coincidental results of chromosomal instability during branching evolution (genetic drift) or actively selected for in the specific tumour regions, we employed a simple

(13)

198 | Chapter 8

permutation-based analysis. For each tumour, chromosomal segments of the same bp 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 Methods and URLs for software access). This approach successfully identified canonical and non-canonical bona fide driver events as selectively enriched by overlap in 5/19 tumours (Fig. 5d-e). 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 to the stem. Neither were ICAs overall overrepresented in clonal compared to subclonal branching events (p=0.55; two-sided Fisher’s exact test). We conclude that while ICAs involving specific driver genes are most probably non-random events favoured by selection, ICAs involving larger chromosomal segments can be explained by genetic drift.

Some childhood tumours are born bad, others grow bad

That canonical aberrations showed a highly diverse distribution among stems versus branches could indicate that some of them generally occurred at an earlier stage of tumour development than others. This prompted us to investigate the temporal sequence of aberrations further, with a particular focus on those canonical aberrations that define good or bad prognostic subgroups for the three studied tumour entities (see Methods for specification). 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 few sequence mutations and at the same time many chromosomal rearrangements that may emerge through sudden bursts (e.g. chromothripsis) or show back-mutation (e.g. reversion from trisomy to disomy). Also, the tumour types in question 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 mutations17_{. We therefore instead employed our data on spatial}

genetic variation to make a simple estimation of whether individual chromosomal aberrations and mutations occurred preferentially early or late in tumour development. Mutations that occur either before or at the time of the last complete clonal sweep are likely to be present in all tumour cells in all tumour regions at the time of presentation/sampling, thus placing these in the stem of each tumour’s phylogeny (Fig. 6a). Conversely, mutations that occur later during tumour growth are likely

(14)

Extensive clonal branching shapes the evolutionary history of high-risk pediatric cancers | 199 to be confined to a few samples or regions. Finally, we reasoned that the timing of aberrations detected in all samples but not clonally would have an ambiguous status. 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. We then compared the total numbers of early and late aberrations among the three tumour types and found no consistent differences. In contrast, there were clear differences among prognostic subgroups within each tumour type (Fig. 6b-g). In Wilms tumour, it is well known that high-risk diffuse anaplastic Wilms tumours compared to moderate and low-risk tumours contain overall a higher burden of copy number alterations18_{. 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 specific canonical aberrations were remarkably consistent with this scenario. Mutations in TP53 and deletions of the corresponding part of 17p – phenomena closely associated with high-risk tumours and anaplasia – were late events in all five diffuse anaplastic tumours 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 tumours, typically emerged early.

In contrast to the scenario observed in Wilms tumours, low-risk neuroblastomas differed from the two high-risk neuroblastoma subtypes (MYCN amplified and structurally rearranged non-amplified) by having their characteristic profile of whole-chromosome anomalies confined to the early group of mutations (Fig. 6c). Similarly, the structural aberrations signifying high-risk neuroblastomas emerged early in high-risk cases, while such aberrations only appeared as late anomalies in low-risk tumours. 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, respectively19_{, were overrepresented among early aberrations (Fig. 6f). So was MYCN}

(15)

200 | Chapter 8

Fig. 5 Distribution of canonical genetic aberrations. For each of the tumour types Wilms tumour

(WT), neuroblastoma (NB) and rhabdomyosarcoma (RMS) a specific set of characteristic and highly recurrent canonical genetic aberrations of presumed pathogenic importance were defined (see Methods and Supplementary Table 4). a Relative distribution of canonical aberrations in the phylogenies of all tumours; p-values for differences in the proportion of canonical aberrations present only in phylogenetic branches were obtained by two-sided Fisher’s exact test. b SNP array log2 plots showing iterated canonical aberrations (ICAs) in NB19 and NB23, comprising a step-wise loss of CDKN2A through loss/deletion and a step-wise ALK gain-of-function through mutation/amplification, respectively. Plots for chromosomes 10 (NB19) and 1 (NB23) are included as references. c Panorama of ICAs in WT, NB and RMS with the number of tumours exhibiting such events given below each pie chart; m, mutation; +, gain; -,loss; cnni, copy number neutral event. d Results of enrichment analysis of all CNAs involving stem gain followed by branch gain (gain-gain) or stem loss followed by branch loss (loss-loss). Of tumours with this series of events, around half (19/39) showed overlap between stem and branch aberrations, but only five cases a significant enrichment indicative of selection. e Enrichment analyses results from six tumours, 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 tumours with a higher degree of stem-branch overlap than expected from a random distribution; p-values were calculated directly from frequency distribution and denotes the probability that an overlap of at least the detected size would occur by random distribution of gains or losses. Losses in NB20 exemplifies an actual degree of overlap well in accordance with a random distribution, while the remaining tumours illustrate significantly enriched overlap of losses or gains.

(16)

Extensive clonal branching shapes the evolutionary history of high-risk pediatric cancers | 201 Rhabdomyosarcomas showed a scenario similar to neuroblastoma. Here, the embryonal subtype is denoted by having multiple alterations in whole chromosome copy numbers, while the alveolar subtype is denoted by multiple structural aberrations, including translocations involving the FOXO1 gene21_{. This contrasting prevalence of numerical}

versus structural aberrations was reflected among early aberrations, with only 1/6 embryonal rhabdomyosarcomas containing any early structural aberrations but all 6 containing early numerical aberrations typical of the subtype including trisomy and copy number neutral imbalance of chromosome 11 (Fig 6d and 6g). 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 tumours, where the complex genomic features of diffuse anaplastic tumours occur closer to the point of tumour presentation and surgical sampling.

Single cell sequencing corroborates extensive clonal branching in high-risk tumours

Clonal deconvolution based on data from bulk samples suffers from several limitations in resolving clonal hierarchies22_{. For example, it is often impossible to deduce whether}

small subclones are nested within each other or not, i.e. if one is a descendant of the other (Fig. 1c). Also, subclones with chromosomal imbalances leading to opposite copy number shifts (e.g. gain and loss of the same chromosome) can cancel each other out and be missed at analysis. To improve the resolution of phylogenetic analysis, we therefore subjected eight tumour samples to single cell flow-sorting followed by single cell low-pass (0.01-0.02x) whole genome sequencing23_{. The samples were chosen to represent an equal}

number of high-risk and low-risk tumours (Fig. 7a and 7b; Supplementary Fig. 2). In total, 547 libraries were of sufficient quality for conclusive single cell sequence (SCS) analysis (50-93 cells per case) . To establish cut-off values for scoring copy number aberrations in each of these cells, we used the copy number profiles of accompanying non-cancer 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 tumour bulk sample. Comparing chromosomes that appeared normal disomic in both tumour cells and non-cancer cells, there was no difference in the average alteration rate per 1 Mb bin of sequence reads (0.13% for tumour cells; 0.56% for non-cancer cells; p=0.16; two-sided Student’s t-test). At a cut-off 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

(17)

202 | Chapter 8

breakpoints from bulk sample array analysis to analyze breakpoint precision by SCS, again resulting in a cut-off of five or more 1 Mb bins for scoring significant differences in breakpoints. Based on these criteria, tumour cells that did not show significant differences in copy number profiles or breakpoints were grouped together into clones, while cells with unique genotypes were kept as single cells at phylogenetic analyses.

Cell-to-cell variation was readily assessable in all eight cases (Supplementary Fig. 3). This revealed a broad range of intercellular genetic diversity among the tumours. Similar to the scenario found by MRS, ICAs were a common feature (6/8 cases) also in the SCS phylogenies, including accumulation of 17q copies in the branches of 3/3 neuroblastomas and further copy number enrichment in 3/4 Wilms tumours of chromosomes/chromosome arms typically gained in this neoplasm (7q+, +8, and +12)15_{. There was again a systematic}

difference between high- and low-risk cases. The former had longer phylogenetic branches than did the latter, leading up to either clones or single cells (Fig. 7c). Because the single cell data were derived from only one sample per tumour, it could not be used to calculate the number of phylogenetic branches detected per sample. On the other hand, SCS allowed a precise estimate of the cell numbers present in each clone in the sequenced cell population (Fig. 8a-b). This made 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 tumours (Fig. 8c), as well as the fraction of cells having a genotype only detected once (Fig. 8d). Both these estimates showed that cell-to-cell genetic diversity was significantly larger in the high- than in the low-risk tumours. Overall, the SCS data thus replicated the impression from MRS analysis that high-risk tumours exhibit more frequent as well as more extensive phylogenetic branching than do low-risk tumours.

Single sample single cell sequencing detects transregional overlap of clones

Even though SCS is advantageous for the deconvolution of clones in a specific sample, little is known about the extent to which it reflects the clonal landscapes of samples from different areas of the primary tumour. If clonal territories overlap inside sampled regions it should be possible to identify, by SCS analysis of a certain region, not only the clones that dominate the region in question but also other clones, which have transgressed into the sampled area from neighbouring areas where they in turn predominate (Fig. 9a-b). According to this scenario, one or more cells constituting a phylogenetic branch in the sample subjected to SCS would then show a genetic profile identical to that of a clone ascertained as the major clone by bulk sample MRS analysis from another area of the same tumour (Fig. 9c-d). Among the eight tumours subjected to SCS, only NB19 showed

(18)

Extensive clonal branching shapes the evolutionary history of high-risk pediatric cancers | 203 evidence of such transregional overlap. This neuroblastoma exhibited highly complex chromothripsis-like rearrangements of chromosome 5 with variation among phylogenetic branches at both MRS (Fig. 3h) and SCS (Fig. 7a) analysis. Two anatomic regions not analysed by SCS but covered by MRS (samples B7 and B8) shared with two samples from bone marrow metastases (M1 and M2) a whole genome profile with a specific chromothripsis-like complex rearrangement in chromosome 5 and multiple deletions in chromosome 8 (Fig. 9e). An identical profile of aberrations was detected at SCS analysis in a specific branch leading to the small subclone A (4/64 cells; Figure 7a NB19). In the remaining cases with regional copy number variation by MRS (NB18, NB20, RMS6, RMS7 and WT6), SCS did not detect specific clones found by MRS in other regions. Thus, SCS analyses of a single region may be used to detect small cell populations, which come to dominate other regions of the same primary tumour and/or metastatic sites in the same patient. However, the phenomenon appears to be rare and the wide intercellular genetic variation observed by SCS will make it difficult to predict from one sampled region which, if any, subclones will come to dominate other foci. We conclude that MRS of the type performed in the present study is necessary to comprehensively cover genetic variation in solid pediatric tumours.

(19)

204 | Chapter 8 Early a A founder genome tumor regions normal genome B C stem mutations global branch mutations Late regional clonal mutations regional subclonal mutations c FAV(9) MNA(7) early late Neuroblastoma WCA 0 25 15 late CNAs d 0 40 EMB(6) ALV(2) 160 earlyRhabdomyosarcomalate

EMB(6) ALV(2) 0 10 10 STR(8) STR 0 30 30 early CNAs b early late Wilms tumor WCA STR IRH(14) BTH(5) DAH(5) IRH(14) BTH(5) DAH(5) 0 10 40 0 10 _early _late 15 p = 0.0001 p = 0.001 p = 0.17 p = 0.97 p = 0.001 p = 0.03 p = 0.0001 p = 0.096 late CNAs early CNAs early late

prognostic aberration prevalence

FAV(9) MNA(7) STR(8) late CNAs early CNAs early late WCA STR

cases with canonical prognostic aberration cases without canonical prognostic aberration

e

Other canonical aberrations Canonical prognostic marker

early late Wilms tumor 0 50% 30% 10_p=0.0044 1 0 p=0.0495 early late early late early late Neuroblastoma 9 2 p=0.036 p=0.023 8 1 7 0 p=0.0094 0 50% 30% early late Rhabdomyosarcoma 0 80% 20% 6 0 2 0 2 0 f g aberration prevalence aberration prevalence aberration prevalence prognostic aberration prevalence

(20)

Extensive clonal branching shapes the evolutionary history of high-risk pediatric cancers | 205 ◀ Fig. 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 (CNAs) in Wilms tumour, 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 subdivided into a field for each prognostic subgroup including for Wilms tumours intermediate risk histology (IRH), the moderate-risk blastemal type histology (BTH) and the high-risk diffuse anaplastic histology (DAH); for neuroblastoma low-risk favorable tumours (FAV), high-risk MYCN amplified (MNA) and high-risk tumours with structural rearrangements without MYCN amplification (STR); for rhabdomyosarcoma the low-risk embryonal (EMB) and high-risk alveolar (ALV) subtypes. Box plots denote the 25 and 75 percentiles, whiskers 95% confidence intervals and 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 4 samples per tumour. 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 color) or not having (gray color) 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 tumours, the only canonical prognostic aberration was deletion/mutation (STR) of TP53, which was found only in DAH and only among late CNAs. For neuroblastomas and rhabdomyosarcomas, the canonical prognostic anomalies are detailed in Methods. e-f 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 tumours (black type) where the aberration was present early and late, and the p-value (red type) from two-tailed Fisher’s exact test. Chromosomal aberrations are denoted as copy number neutral imbalances (cnni), gains (+), losses (-), while SNV/indels affecting genes are denoted by the name of the gene. Translocations causing gene fusions are denoted by (t.) and MYCN amplification by MNA. Only aberrations occurring in two or more tumours were included.

(21)

206 | Chapter 8 A 1 CNA NB18 B C F E MYCN amp 17q+ 8 CNA 17q++ 3 CNA NB19 10 CNA 1p-17q+ D E C A B CHL(5)C CHL(5)D 17q++ 4n NB20 3 CNA 25 CNA 4n -- 1p--11 17q+ E C A F 17q++ 5 CNA 25 CNA PAX7/FOXO1 FNDC3B amp TP53-/-4n RMS7 RMS6 1 CNA +2 +5 +10 +12x2 +20 9 CNA A B C 1 CNA 9 CNA 1p-1q+ -11 A C B WT6 WT9 1 CNA 4 CNA +7 +8x2 F D A B C E 1 CNA WT14 A 2 CNA 7p-7q+ 7q++ B 7q++ a b c HIGH-RISK TUMORS

LOW-RISK TUMORS branch length (CNA)

NB18 NB19 NB20 RMS7 RMS6 WT6 WT9 WT14 HIGH RISK LOW RISK p<0.0001 p<0.0001 p<0.0001 p=0.002 n=4 n=2 n=5 n=2 p q 17 p q 17 CHL(5)E p 17q p 17q n=5 n=4 n=6 n=4 6 n=3 n=2 +8x3 7 8 n=5 7 8 n=3 n=4 +12x312 13 12 13 n=2 n=4 n=5 q chr. 3 p q chr. 3 p FNDC3B amp n=13 n=8 n=4 6 7 n=2 n=1 n=5 6 7 n=2 n=1 n=3 CHL(5)A CHL(5)B n=4 n=2 p q 17 n=5 n=2 p q 17 p q 17 n=3 +17

(22)

Extensive clonal branching shapes the evolutionary history of high-risk pediatric cancers | 207 ◀ Fig. 7. Single cell phylogenies. a Maximum parsimony trees from high-risk tumours including three neuroblastomas (NB18, NB19, NB20) and one rhabdomyosarcoma (RMS7). The number of copy number aberrations in the stem (CNA) is specified. Capital letters at the end of branches correspond to the location of detected clones, i.e. two or more cells with identical CNA profiles, while branches not ending in letters correspond to genotypes detected in single cells. Stem aberrations characteristic of each tumour type are given in brown text, tetraploidization by 4n, and iterated canonical aberrations (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 oncogene50_{. See Supplementary Fig. 1 for}

examples of full single cell CNA profiles. b Maximum parsimony trees from low risk tumours including one rhabdomyosarcoma (RMS6) and three Wilms tumours (WT6, WT9, WT14) with annotations as in panel a. c Branch lengths measured by the total number of CNAs included in each branch (circles) with medians denoted by red lines and P values by two-tailed Mann-Whitney U-test.

A(36) B(2) C(2) E(2) F(2) NB18 N = 66 A(4) B(7) NB19 N = 64 C(13) E(2) D(13) NB20 N = 45 A(15) E(4) C(2) F(2) RMS7 N = 84 8 36 40 A(41) RMS6 N = 50 B(2) C(2) A(74) B(10) C(2) WT6 N = 93 A(50) B(3) WT9 N = 75 C(4) D(2) E(3)F(2) A(41) B(17) D(2) WT14 N = 70 HIGH RISK LOWRISK 1 0 Simpson index Genetic diversity p=0.029 c HIGH-RISK LOW -RISK a

b _{Proportion of cells with unique genotype}

NB18 NB19 NB20 RMS7

RMS6 WT6 WT9 WT14

d

p<0.00001

clonal cells singlet cells 50% clone size singlet descendants HIGH-RISK LOW-RISK

Fig. 8 Population structure based on single cell sequencing. a,b Genotype distribution among cells in high-

and low-risk tumours as derived by phylogenetic analysis. Diameters of blue circles correspond to relative clone sizes while red lines correspond to single cell genotypes. The total number of analyzed cells in each tumour is specified (N=) and the number of cells detected in each clone (A, B, C etc.) given in parentheses. Black lines mark the phylogenetic relationship between each clone, the stem, and the other clones as detailed in Fig. 7.c The Simpson index of diversity of each tumour (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 tumour whose genotype was only detected once (singlet cells) as compared to the proportion of cells whose genotype was detected in two or more cells (clonal cells). P-value was calculated by two-tailed Fisher’s exact test from absolute cell numbers.

(23)

208 | Chapter 8

Multiregional sampling (MRS)

target region for SCS (T) other sampled regions (O)

primary tumor

Single cell phylogeny of T

stem cell T1 cell T2 cell T3 cell T4 cell T5 cell T7 cell T6 region T clone B clone A Transregional overlap a Multiregional phylogeny T O1 O4 O5 branch O2 region O1 b c O3 cell T8

minor clone invading T from O1 identical CNAs

stem

branch

5 5 8 8 5 5 8 8

Bone marrow met (M1) Primary tumor (B4=T) log 2 ratio -2 2 log 2 ratio -2 2 1 5 8 X 1 5 8 X MRS read numbers 0 450

Cell from subclone A

read numbers

1 5 8 X

Cell from subclone C

1 5 8 X 0 500 SCS e d NB19 NB19 overlap

Fig. 9. Phylogenies based on multiregional versus single cell analysis. a Multiregional sampling (MRS) was

arranged to cover the target region (T) of single cells sequencing (SCS) as well as other regions of the primary tumour (O). b SCS analysis of a sample from region T can yield information on the major clone for this region (A) as well as smaller clones spilling over into T from neighboring regions, here exemplified by clone (B) dominating region O1. c-d If such transregional overlap of clones occurs, groups of cells identified by SCS

in branches of the phylogenetic tree from T (c; cell T5) can be identical (double-headed arrow) to a branch

(d; O1) other than T in the MRS tree. e This scenario was verified in NB19. Copy number aberration profiles

determined by MRS identified only one clone in the target region for SCS (B4). A second clone was detected by MRS in a bone marrow metastasis (M1). The two clones differed by complex structural rearrangements of chromosome 5, with the B4 clone having a combination of gains and losses and M1 only gains. M1 also exhibited multiple unique deletions in chromosome 8. In 4/64 cells (subclone A), SCS of B4 detected a profile identical to the M1 clone. The MRS profile for B4 best corresponded in the SCS data to the highly prevalent subclone C (13/64 cells). Colored arrowheads denote corresponding regions of gains and losses in the MRS and SCS data, respectively.

(24)

DISCUSSION

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 explosions7_{. 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 analyses clonal ancestry and the evolutionary history of mutations important from a clinical diagnostic and/or prognostic perspective. In tumour 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 towards 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 tumour type. For further validation, we also performed single cell karyotyping based on which phylogenetic analysis can be performed without the additional step of clonal deconvolution24_{. Although only eight cases could be analysed 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 tumours. An alternative approach to our sample number standardization and single-cell validation would be to include strictly the same number of samples for each tumour subtype. However, this would create other biases, such as a failure to include small tumours that provide few samples or, conversely, missing biologically essential phylogenetic branches in large tumours due to relative undersampling.

Understanding what factors drive a high-risk phenotype in pediatric solid tumours is essential for finding new ways to cure these neoplasms for which the greatest threat to patient survival is relapse after chemotherapy. Relapse tumours have, in the vast majority of cases, a genomic profile distinct from that of their corresponding primary tumours1-7_.

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 evolvability per se would be a factor that is strongly associated to the risk of relapse. The present study highlights

(25)

210 | Chapter 8

that very simple phylogenetic parameters, such as branch length and branching frequency, could act as proxy markers for such overall cancer cell evolvability. Our data also enabled a systematic phylogenetic inventory of genetic aberrations characteristic of the tumour types in question. We found that a significant number of such canonical aberrations important for diagnosis or prognosis were confined to branches and thus were often detected only in some tumour regions. This was most pronounced in Wilms tumours where only around one third (27/74) of the detected canonical aberrations were present in all cells in all samples (stem anomalies); the remaining canonical changes were either detected solely in regional branches or modified in branches to create regional variations of a stem anomaly. For Wilms tumours, this high intratumoural variability of the most common genetic markers could explain why relatively few clinically useful molecular prognosticators have been detected in this tumour type despite decades of biomarker research over large datasets25_.

The present study suggests that quantifying genetic variation per se could be a more successful strategy in Wilms tumour than focusing on specific genetic markers.

Even though our data indicated that the investigated high-risk tumours 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 PAX3/7-FOXO1 fusion for the latter. These mutations were concomitant to a large number of other structural chromosomal rearrangements, typical of high-risk tumours of both subtypes. In contrast, the anaplastic Wilms tumours started out largely as other Wilms tumours, often with the characteristic but non-prognostic copy number neutral imbalance of the IGF2/H19 and WT1 loci in 11p. In high-risk tumours, this was subsequently followed by a regional emergence of complex chromosome changes concomitant to the TP53 mutations that signify diffuse anaplasia. This concurrence in time between chromosomal instability and prognostic aberrations in high-risk tumours indicate 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 tumours, but also true for neuroblastoma, in which MYCN amplification has been linked to an attenuated G1 arrest after DNA damage similar to that described for TP53 loss.26_{More specifically, recent studies have indicated that}

(26)

Extensive clonal branching shapes the evolutionary history of high-risk pediatric cancers | 211 to DNA damage27_{. MYCN-amplified neuroblastoma cells have also recently been shown}

to circumvent cell death by upregulating components of an error-prone non-canonical alternative nonhomologous end-joining (NHEJ) DNA repair pathway28,29_{. What triggers}

early-onset chromosomal instability in the neuroblastoma subgroup dominated by structural rearrangements in the absence of amplified MYCN is less clear, but their chromosomal instability phenotype may be related to an early allelic gain of 17q, harbouring the negative p53 modulator PPM1D, which is typically overexpressed in high-risk neuroblastoma30_.

Finally, in alveolar rhabdomyosarcomas, the PAX3/7-FOXO1 gene fusion has been shown to directly upregulate MYCN31,32_{, suggesting a shared mechanism of structural chromosomal}

instability between high-risk rhabdomyosarcoma and neuroblastoma.

In all, our data indicate that the childhood solid tumours most prone to relapse are those that acquire extensive evolvability through oncogenic driver mutations that also attenuate control of genomic integrity. This destabilization of the genome occurs already at initiation of clonal expansion for some tumours, which are thus born bad, while it evolves at a later stage for others. Our findings underscore that a substantial evolutionary plasticity of the neoplastic cells needs to be anticipated when designing novel therapeutic strategies for high-risk childhood cancers.

ACKNOWLEDGEMENTS

This study was supported by grants to D.G. 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, the Gunnar Nilsson Cancer Foundation, and the Medical Faculty of Lund University Sweden. The authors also 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.

(27)

212 | Chapter 8

AUTHOR CONTRIBUTIONS

D.G., J.K., A.V., and F.F. conceived and designed the project. B.B., J.K., A.V., L.H.M., and D.S. coordinated and analyzed whole genome genotyping and sequencing data. N.A. performed phylogenetic analysis. A.V. performed statistical simulations. D.G. performed clinical correlation studies.

Additional information

Supplementary Information accompanies this paper.

Competing interests: The authors declare no competing interests.

MATERIALS AND METHODS URLs

Mitelman Database of Chromosome Aberrations and Gene Fusions in Cancer https://cgap.nci.nih.gov/Chromosomes/Mitelman

European Genome-phenome Archive www.ebi.ac.uk/ega/home

Custom software

http://github.com/avalind/subclonal-enrichment.

Ethics review, patients, and samples

The present study complied with all relevant laws and ethical guidelines regarding human research participants and animal experiments. Analysis of human tumours was approved by the Regional Ethics Board of Southern Sweden (permit no. LU289-11 and LU119-03) with written informed consent obtained from parents or legal guardians of the patients. Data on mutation type and its prevalence in each sample was retrieved from Supplemental Table 1 of Karlsson et al7_{. Corresponding data for WT24 and NB23}

were added as part of consecutive data collection of pediatric cancers at the Lund and Karolinska University Hospitals, Sweden, from 2000-2016. Inclusion criteria for the study was a histopathologically confirmed diagnosis of Wilms tumour, neuroblastoma or rhabdomyosarcoma with prognostic subtyping performed according to present guidelines, availability of at least two primary tumour samples with >50% tumour cells according