ABSTRACT Neuroblastoma (NB) is one of the most lethal childhood cancers due to its propensity to become treatment resistant. By spatial mapping of subclone geographies before and after

chemotherapy across 89 tumor regions from 12 NBs, we find that densely packed territories of closely related subclones present at diagnosis are replaced under effective treatment by islands

of distantly related survivor subclones, originating from a different most recent ancestor compared to lineages dominating before treatment. Conversely, in tumors that progressed under

treatment, ancestors of subclones dominating later in disease are present already at diagnosis. Chemotherapy treated xenografts and cell culture models replicate these two contrasting

CLONAL EVOLUTION DURING METASTATIC SPREAD IN HIGH-RISK NEUROBLASTOMA Article 11 May 2023 EVOLUTIONARY TRAJECTORIES OF SMALL CELL LUNG CANCER UNDER THERAPY Article Open access 13 March 2024

CHEMOTHERAPY INDUCES CANALIZATION OF CELL STATE IN CHILDHOOD B-CELL PRECURSOR ACUTE LYMPHOBLASTIC LEUKEMIA Article 05 July 2021 INTRODUCTION Neuroblastoma (NB) is the most common pediatric

solid malignancy outside the central nervous system. Although the pathogenesis of NB has been thoroughly investigated, the five-year overall survival of high-risk NB remains around 50%1. The

main cause of death in NB is treatment refractory disease2. Treatment resistance occasionally manifests as progressive metastatic disease already under first-line multimodal chemotherapy.

More commonly, a treatment-resistant metastatic relapse emerges after an initial response to therapy and surgical removal of the primary tumor. Relapse is associated with a dismal prognostic

outlook3,4, with few if any standard therapeutic options that reliably promise success. Therefore, molecularly targeted therapy is increasingly being considered and there is now a broad

range of clinical trials, of which many are coupled to the detection of a specific genetic profile in the patient’s tumor5,6. A major challenge when matching treatment to molecular profile

is the fact that the NB genome is neither uniform nor stable within a patient7. Early studies comparing primary NBs to matched relapses, demonstrated significant differences in genomic

profiles over time8,9,10, while more recent studies have shown substantial genetic diversity also in primary tumors11 with ongoing evolution during and following chemotherapy12, and through

the metastatic process13. However, little is known about exactly how NB cells surviving chemotherapy in the primary tumor are related to the lineages dominating prior to treatment. In the

present study we perform high-resolution spatial mapping of subclones having survived chemotherapy in NB and investigate their phylogenetic relationships to clones detected at diagnosis. We

find that subclones surviving effective chemotherapy are typically collaterally related (“cousins”) to those present in the treatment-naïve biopsy, while those found after progression under

failed treatment are usually linear descendents (“sons” and “daughters”) of the subclones found prior to treatment. RESULTS SPATIOTEMPORAL PROFILING UNDER CHEMOTHERAPY REVEALS TWO

CONTRASTING EVOLUTIONARY PATTERNS NB is largely a copy-number aberration (CNA) driven disease14,15,16,17,18,19,20,21. We have reported a software tool (DEVOLUTION) that allows integration of

point mutations and CNA data, where the latter is curated according to constraints of chromosomal evolution22. Please see specification in Methods for Data and Code availability. To perform

high-resolution mapping of subclone territories across anatomic tumor space, we applied DEVOLUTION to mutational data from archived paraffin embedded formalin-fixed NB tissue (Fig. 1a–e;

“Methods”). This allowed a comprehensive comparison between the genomes of NB cells before and after treatment, including both patients with a remaining, extensively viable tumor and cases

with good chemotherapy response where survivor populations were limited to spatially dispersed tumor islets a few millimeters in size. Areas of ≅20 mm2 containing >30% tumor cells were

subjected to whole genome copy number profiling parallel to whole exome sequencing (WES). Mutations detected by WES were validated by targeted resequencing using a custom-made panel based on

the WES findings. From a consecutive patient cohort treated at a single pediatric oncology center over a period of 20 years, we identified 12 patients from whom viable tumor tissue was

available from ≥2 geographically separated tumor regions from ≥2 time points (Fig. 1f, g; Supplementary Data 1a). A similar number of primary tumor samples (median 3.0) per patient from

before treatment as from after chemotherapy were analyzed with both CNA profiling and sequencing; fewer samples were from metastatic sites (median 2.0). Among the 12 patients, two exhbitied

very high-risk genetic profiles including _MDM2_ amplification23 in one and _NF1_ deletion along with _MYCN_ amplification24, in the other. Both patients showed progression of the primary

tumor with metastases under treatment (Patients 1–2). Nine other patients showed a significant therapy response (40–97% primary tumor regression including necrosis, fibrosis and extensive

maturation; Patients 3–11), while one patient had a single small tumor in the neck and was treated with up-front surgery only and later developed a metastatic relapse (Patient 12). All seven

patients whose tumor could be evaluated before and after significant chemotherapy response, showed a strikingly similar phylogenetic relationship between subpopulations detected before and

after treatment (Patients 3–9; Fig. 2a and Supplementary Fig. 1, Supplementary Data 2a): some genetic changes found at diagnosis were no longer detected after treatment while other mutations

emerged, as populations present prior to chemotherapy were replaced by others, typically arising from a different most recent common ancestor (MRCA). In the phylogenetic trees, this

corresponded to one set of collateral branches being replaced by other, parallel branches and hence we denoted this pattern collateral clonal replacement (CCR). CCR occurred irrespective of

whether the pre-treatment sample was obtained from the primary tumor (Patients 3–7) or from a metastasis sampled at presentation of disease (Patients 8–9). That CCR could take place

gradually over time was illustrated by Patient 8 (Fig. 2b). Here, the tumor was biopsied at diagnosis, after ~50 days of therapy (due to slow response) and then sampled again at resection

after ~500 days under treatment. A comparison of these time points revealed a gradual shift towards CCR as populations were replaced. Using an index of genomic diversity (IGD) that was

independent of both sample number and the total number of detected mutations (see Methods), we compared the subclonal heterogeneity among samples obtained before and after chemotherapy. In

the primary tumor, diversity (IGD, Fig. 2c) was found to be higher after than before chemotherapy, even though the mutational burden (phylogenetic branch lengths) did not differ between

these time points (Fig. 2d). CCR in the primary tumor under therapy was thus associated with the emergence of a more diverse clonal landscape. CCR was also observed for the majority of

subclones in the primary tumor when compared to those in a later metastatic relapse (_n_ = 5; Patients 4, 9–12). Phylogenies containing clones detected at metastatic relapse exhibited a

higher degree of irregularity (see “Methods”) than those only based on clones detected in the primary tumor and/or up-front progression (Fig. 2e, f) as the lengths of branches to clones

detected in relapse samples were overall longer than those in the primary tumor. In most relapsed cases, the lineage of relapsing subclones diverged from that of the primary tumor at an

early stage, with CCR encompassing all subclones in the primary versus the relapse. Minor but significant exceptions from this were Patients 4 and 12. In Patient 4, a regional population

corresponding the common ancestor of relapsing subclones was detected in the primary tumor after therapy (Supplementary Fig. 1d:VII, subclone A), while in Patient 12 it was detected prior to

therapy (Supplementary Fig. 1l:VII, subclone H). Notably, these two primary tumors were the ones with least treatment-induced regression; Patient 4 responded largely by tumor cell

ganglionic maturation, while Patient 12 was treated with surgery only. This indicated that inferior killing of cancer cells was associated with an aspect of linear evolution from surviving

ancestral clones when/if the patient progressed. Indeed, in the cases with local and metastatic progression under chemotherapy (Patients 1–2), a subclone ancestral to the genomes of

metastatic populations was ascertained at diagnosis (Fig. 2g–i), supporting features of linear evolution under inferior therapy response. LINEAR EVOLUTION UNDER PROGRESSIVE GROWTH OF

CHEMOTREATED PATIENT-DERIVED XENOGRAFTS Because the number of patients with progressive disease under therapy was limited, we set up a model system to mimic progressive tumor growth under

chemotherapy. We used a well-established NB murine model system of patient derived xenografts (PDXs)25,26. PDXs from a _MYCN_-amplified (MNA) NB (_n_ = 6) were treated by cisplatin

monotherapy, resulting either in stationary disease (_n_ = 3) or progression (_n_ = 3) under treatment (Fig. 3a and Supplementary Fig. 2a:I) as reported27. We also simulated progression

after incomplete surgery by partial excision of non-chemotherapy exposed PDXs followed by regrowth to the same volume (_n_ = 4). Control PDXs (_n_ = 4), grown to the same volume under saline

treatment, were used as a reference. We analyzed two biopsies from each tumor (A and B) and phylogenetic analysis revealed a pattern of linear evolution with the cell culture from which

PDXs were established acting as founder (Fig. 3b, c; Supplementary Fig. 2a:II–III). Most PDXs progressing under chemotherapy or after surgery showed expansion of a clone with partial gain of

1p, all of 1q and an extra copy of 17q, leading up to clonal sweeps (Supplementary Fig. 2b–d). They also had a significantly higher CNA burden and more private mutations than the

cisplatin-stationary and the untreated group, in a fashion that was to some degree time-dependent (Supplementary Fig. 2e); while untreated tumors had the same end-size as tumors progressing

under treatment or after resection, they reached this size in shorter time and had fewer anomalies. Whole exome sequencing revealed clonal evolution in PDXs mimicking clinical tumors,

including parallel evolution of RAS/MAPK pathway mutations in linear radiations from the mother culture (Fig. 3d; Supplementary Data 2c–e). To evaluate evolutionary dynamics also at

progression under multimodal high-dose chemotherapy, a different PDX model, with up-front resistance to multimodel chemotherapy was then analyzed (Fig. 3e; data obtained from Manas et

al.27). Here, PDX tumors were subjected to CNA profiling either after free growth (_n_ = 3) or after progression under treatment (_n_ = 4) with an intense five-drug chemotherapy protocol

mimicking the regimen used to treat high-risk NB patients, i.e. rapid COJEC containing cisplatin (C), vincristine (O), carboplatin (J), etoposide (E), and cyclophosphamide (C). Similar to

the situation with progression under cisplatin monotherapy, we found that both untreated and treated PDX tumors exhibited linear radiation from founder populations, progressing through

clonal sweeps in 1/3 and 3/4 treated and untreated PDXs, respectively (Fig. 3f–g and Supplementary Fig 2f–g). Thus, PDX models with both cisplatin- and COJEC-resistant NB cells independently

confirmed the clinical scenario of a pattern of linear evolution from a common ancestor at upfront progression, radiating in different directions at different anatomic sites (different PDX

mice). SUBCLONE REPLACEMENT AT RELAPSE AFTER EFFECTIVE CHEMOTHERAPY IN PATIENT-DERIVED XENOGRAFTS To compare the results from patient cases with effective chemotherapy treatment to a murine

model, PDX3 tumors were subjected to COJEC (Fig. 4a)27. Mice were randomized into treatment groups when NB PDX tumors reached approximately 500 mm3. After regression under COJEC to <200 

mm3, surgery was performed to remove all visible tumor. The COJEC-treated tumors (9 × 2 samples per tumor) were compared to identically sampled untreated tumors (saline injections; _n_ = 9 ×

 2). Treated and untreated tumors had some common evolutionary features, including expansion in 17/18 tumors of a subclone signified by the same gains in chromosomes 1 and 17 as in the

cisplatin-only monotherapy model detailed above (Fig. 4b, c). A comparison of genetic diversity in treated and untreated tumors showed significant differences between the cohorts (Fig. 4d,

left and center panels; Supplementary Fig. 3a–b; Supplementary Data 2f, g). Despite treated tumors being on average 10 times smaller in volume than those not treated, all treated tumors

(9/9) harbored subclones that were unique to one of the two sampled regions while this was rare (2/9) in untreated tumors (_P_ = 0.0023; two-tailed Fisher’s exact test). In line with this,

the number of subclonal alterations differentiating the two sampled regions as well as the total number of region-unique subclonal alterations was higher in treated than in untreated tumors,

while there was little difference in the total number of aberrations between the COJEC-treated and untreated cohort (Fig. 4d, right panel). The scenario was thus similar to the findings in

the clinical cohort, where there was no overall increase in the total number of aberrations but still an increased genetic intratumor diversity after chemotherapy treatment. Because PDX

tumor sizes were feared to be too small for sampling prior to treatment without inducing populations bottlenecks, the finding of CCR under therapy from the clinical cohort could not be

validated by sampling PDXs before and after therapy as in the patient cohort. However, COJEC-treated PDXs T1 and T5 relapsed locally, enabling comparison of clonal landscapes at two distinct

time points, i.e. just after treatment versus after regrowth (Fig. 4e and Supplementary. 3c). The two PDX relapses showed a strikingly similar route of evolution, with additional 17q

copies, and successive intragenic deletions in the _PTPRD_ and _TCF4_ tumor suppressors27,28,29. While not being linked to treatment resistance per se, deletions of these genes have recently

been shown to be recurrent features of evolving neuroblastomas in PDX systems27. Both PDXs showed genetic aberrations that were unique to each sampling time point, i.e. a pattern of

subclonal CCR (Fig. 4f, g). Similar to relapses in the clinical cohort, there was an overall increase in the total number of CNAs after treatment in the PDXs, with the most complex subclones

being found at relapse. PDX modeling thus supported the association between effective chemotherapy response and CCR at relapse. RECAPITULATION OF COLLATERAL CLONAL REPLACEMENT UNDER

CHEMOTHERAPY IN VITRO To test whether the contrasting scenarios of linear evolution under inferior treatment/progression versus CCR under effective treatment could be reproduced also under

tightly controlled in vitro conditions, we used the cisplatin-sensitive cell line IMR-32, derived from a high-risk MNA NB (Supplementary Fig. 4a). In vitro growth for 42 days (11 passages)

without treatment resulted in expansion of the major subclone up to fixation and concomitant loss of a parallel subclone (Fig. 5a, Samples A1-A3; Supplementary Fig. 4b; Supplementary Data 

2h). Low-dose (0.1 µM) cisplatin treatment initially suppressed cell growth but the cells were proliferating at the end of the experiment. This treatment led to a preserved clonal landscape

(Fig. 5a, Samples B1-B3). In contrast, high-dose (1.0 µM) cisplatin exposure with a 94% reduction in population size, followed by regrowth for 42 days, resulted in CCR, with elimination of a

major subclone and the expansion of clones with novel CNAs in 2/3 parallel experiments (Fig. 5a, Samples C1-C3). Reseeding the same number of IMR-32 cells (6% of original) as the viable

population remaining after cisplatin treatment, followed by expansion to the same final population size as the high-dose cisplatin experiment did not change the clonal landscape (Fig. 5b and

Supplementary Fig. 4c), nor was the landscape largely changed by a series of three mechanical reseeding bottlenecks of around 50% (Fig. 5c, Samples MB1-MB3; Supplementary Fig. 4d), while it

again exhibited CCR after three consequtive cisplatin-bottlenecks each killing around 50% of cells. These results again corroborated that effective chemotherapy treatment is linked to a

pattern of CCR and indicated that this effect is not merely contingent on the population bottlenecks to which effective treatment is associated (Supplementary Fig. 4e). We then turned to the

SK-N-SH line, well known for containing two morphological cell types, of which one is more resistant to chemotherapy than the other30,31. These respective cell types have recently been

shown to harbor expression signatures corresponding to mesenchymal (resistant) and adrenergic (chemosensitive) cell states, respectively31. While the two cell populations have been reported

to have largely similar copy number profiles31, in-depth phylogenetic analysis of SK-N-SH under chemotherapy treatment had not been performed previously. We treated SK-N-SH cells with

cisplatin for 72 h in vitro, followed by 3 days of regrowth, and genetic profiling (Supplementary Fig. 4f:I). Treated cells were reseeded, expanded, and again treated with cisplatin,

followed by profiling for CNAs. Untreated control cells were grown in parallel. As expected, after the first treatment, the cells showed a shift to a more epithelioid morphology,

corresponding to chemoresistant state transition and there was little evidence of cell death after the second treatment, compared to 80% cell death at first treatment. Already after the

first treatment, there was a notable shift in the clonal landscape, with a clonal sweep of a cell population containing gain in 1q (1q+ ) replacing a cell population with only stem

aberrations, present in mother cultures. Concomitant to this, there was an expansion of a nested subclone with imbalances in chromosomes 9 and 11 in treated cells (Supplementary Fig. 4f:II

and g). In contrast, the population with only stem aberrations was retained in untreated controls. Furthermore, the two subpopulations dominating treated cultures were suppressed in control

cells by expansion of a population with 2p+ , which in turn was suppressed close to the detection limit (10%) in treated cells. Hence, even in SK-N-SH cells, with a known capacity to undergo

phenotypic state-transition independent of major genetic alterations31, there were rapid shifts in the genetic landscape with features of CCR at chemotherapy treatment. PHYLOGENETIC

BRANCHES CORRESPOND TO SPATIAL TERRITORIES The fact that tumors in the clinical cohort were derived from well-annotated paraffin blocks made it possible to trace back the approximate spatial

location of each subclone in each tumor. In total 35/122 subclones spanned two or more biopsy locations allowing positioning of shared clones to neighboring areas, while subclones confined

to a single biopsy location could not be mapped spatially beyond that location. Applying these principles to the 12 patients in the clinical cohort, we traced phylogenetic branches back to

anatomic territories (Fig. 6 and Supplementary Fig. 5). In untreated diagnostic biopsies, these territories displayed a dense patchwork (median 0.36 subclones per mm2, median distance

between the center points of clonal territories of 1.7 mm), with no obvious anatomic barriers between subclone territories (Fig. 6a–d, i, k; Supplementary Fig. 5a, b, f, l–m). A similar

situation with densely variegated clones were observed in biopsies from metastatic relapses (Supplementary Fig. 5f, i–k). Subclones in the primary tumor after chemotherapy were more

dispersed (median 0.024 subclones per mm2, median distance between the centerpoints of clonal territories of 6.5 mm; p = 0.0087 pre versus post; p = 0.048 post vs. relapse; Mann–Whitney U

test, two-tailed; Supplementary Fig. 5l–m). Mutations detected after chemotherapy were typically fixed at a clonal level ( >90% of all tumor cells) in vastly dispersed islands of

surviving tumor cells surrounded by necrosis, hemorrhage and other reactive changes (Fig. 6a lower, b right, c right, e–g, j, l). These distinct tumor cell populations typically originated

from multiple phylogenetic branches, indicating a high degree of taxonomic diversity well in accordance to increased IGD after chemotherapy. Subclones identified at diagnosis were not

re-detected after treatment in any of the viable primary tumor areas left after treatment (all areas sampled in Patients 3, 5–8). This indicated that the scenario of CCR under therapy was

not only a result of geographically dispersed sampling but resulted from a true disappearance of subclones under treatment. CLONAL REPLACEMENT UNDER TREATMENT UNCOVERS _MYCN_ AMPLICON

DIVERSITY CCR could potentially be explained by ancestral clones surviving therapy and sprouting new branches under treatment, parallel to those detected in the pre-therapy samples.

Alternatively, NB phylogenies are highly branched already at an early stage, prior to treatment. To explore genetic diversity in early phases of disease we performed an in-depth analysis

(down to 0.05 Mb) of copy number variation around the _MYCN_ oncogene in 2p24 in MNA cases, as gain of _MYCN_ gene is thought to be one of the earliest mutations in NB32,33. In the six MNA

cases, an amplified state ( >5 copies) of _MYCN_ was present across all biopsied areas. However, spatiotemporal heterogeneity in the architecture of amplicons in or around _MYCN_ was also

found in all of them (Supplementary Fig. 6). When chromosomal breakpoints of amplicon cassettes were used to infer phylogenetic relationships, all cases had at least one shared breakpoint

across all samples, indicating that one or a few breakage(s) in 2p24 first arose to create patient-specific amplicon cassettes, followed by further modification later on. In total, only ~30%

of samples (4/14) obtained from the primary tumor before treatment or from metastatic sites growing off-treatment after remission exhibited variant (non-stem) breakpoints around _MYCN_. In

contrast, 76% of samples obtained from chemotherapy treated primary tumors exhibited variant breakpoints (13/17; p = 0.00122; Fisher’s exact test). This was in accordance with the scenario

of increased genomic diversity after treatment found across the entire genome. PHYLOGENETIC BRANCHING IS AN EARLY AND STABLE FEATURE OF NB To further monitor early branching evolution with a

focus on the _MYCN_ amplicon, we performed low-pass single cell whole genome sequencing (scWGS) of single biopsies from nine NB primary tumors (three MNA tumors), resulting in 505 single

cell genomes (Patients S1-S9 in Supplementary Data 1b). A unifying feature among all tumors irrespective of risk group was extensive phylogenetic branching (Supplementary Fig. 7). Branches

typically led up to subsets of cells having identical CNA profiles, with the exception of one case (Patient S7), which had a different CNA profile in every cell (Supplementary Fig. 7d). In

MNA patients S1 and S2, the earliest branching event consisted of structural variation in the _MYCN_ amplicon cassette (Fig. 7a–e and Supplementary Fig. 7a). Both these cases contained cells

having low-level _MYCN_ gain as the first event in their phylogenetic trees (in Patient S2 co-gained with _MAML3_), corresponding to an ancestral population. Clonal evolution from this

stage occurred through structural aberrations characteristic of NB, including deletion of 1p, gain of 17q (both Patients S1 and S2) as well as deletion in 11q (Patient S2 only; Supplementary

Fig. 7a:III). To corroborate that MNA could precede 17q gain (the most common CNA in NB), we further evaluated the concurrence of these imbalances by fluorescence in situ hybridization in a

subset of MNA NBs (Supplementary Data 1c). The presence of tumor cells with _MYCN_ copy number gain while still not having acquired extra copies of 17q was ascertained in 5/6 cases (Fig. 7b

and Supplementary Fig. 7f), the exception being S3 in which neither scWGS nor FISH could identify such cells (Supplementary Fig. 7c:I). Interestingly, in both Patients S1 and S2, clonal

evolution through 1p loss/17q gain appeared to be a requirement for further _MYCN_ amplification, with significantly higher amplicon numbers in cells having acquired 1p loss/17q gain than in

cells with only MNA (Fig. 7f, g and Supplementary Fig. 7a:V). Further branching into distinct subpopulations continued after 1p loss/17q gain in all three MNA cases. To test whether the

tendency of branching evolution was a stable feature of NB cells, we turned to three PDX models from MNA NBs25,26. We used scWGS to monitor evolution from PDX in vivo generations 1–7, and

further through transfer into free-floating three-dimensional tumor organoids analyzed at in vitro passages 7 and 2027. This showed that phylogenetic branching was a consistent feature

across time in vivo, a feature that re-emerged after single-cell bottlenecks occurring at transition to in vitro conditions (Fig. 7h–j). TIMING AND FITNESS OF CHROMOSOMAL CHANGES IMPACT NB

GENOME PROFILES The types of CNA associated with branching evolution in the scWGS data varied across tumors, but a common feature (9/9 cases) was parallel aberrations affecting the same

chromosome in different branches of the same phylogeny. For example, in Patient S4, the first branches emerged through distinct parallel deletions in 3p (Supplementary Fig. 7b:I–III), while

in Patient S8 chromosome 6 was the substrate of the three earliest branching events, followed by parallel alterations of chromosome 17 and subsequently of chromosome 1 (Supplementary Fig. 

7b:IV–V). Similar parallel evolution was detected in Patients S1 (_MYCN_ amplicons; Fig. 7d), S2 (chromosome arms 1p, 9q, 17q, 20q; Supplementary Fig. 7a:III–IV), S3 (chromosomes 1, 2, 4, 6;

Supplementary Fig. 7c:I–II), S5 (chromosome 5; Supplementary Fig. 7c:III–IV), S6 (chromosome arm 17; Supplementary Fig. 7e:I–II), S7 (chromosome arms 1q and 17q; Supplementary Fig. 7d), and

S9 (chromosomes 1 and X; Supplementary Fig. 7e:III–IV). Also, in “numerical only” low-intermediate risk NBs (_n_ = 4), structural rearrangements contributed to parallel rearrangements of

the same chromosome. However, there was a difference between clinical-genetic NB subtypes regarding the time point in evolutionary history when different types of CNA occurred (Supplementary

Fig. 7g:I). In “numerical only” low-intermediate risk NBs (_n_ = 4), early generation branches (immediately following the stem) contained only whole chromosome aberrations, while chromosome

breaks were enriched in subsequent generations. In contrast, in high-risk cases (with and without MNA; _n_ = 5) chromosome breaks were frequent already at the earliest generations and

remained so through successive generations. There was also a difference across clinical-genetic NB subtypes in the fitness cost associated with different chromosomal imbalances. As a proxy

to the survivability following different types CNA, we quantified the proportion of progeny (PoP; average fraction of detected single cells) originating from branches with certain CNA

profiles compared to branches with other profiles (Supplementary Fig. 7g:II–III). In low/intermediate risk NB, the PoP from lineages having undergone whole chromosome gains during preceding

generations was significantly larger than the PoP from lineages with chromosome breaks or losses. In contrast, high-risk NBs showed a similar PoP resulting from branches dominated by whole

chromosome breaks as by gains, while the proportion of progeny was lower for cells having undergone losses. To finally evaluate the timing of phylogenetic branching in NB compared to other

tumors, phylogenies based on scWGS for a reported reference group (_n_ = 4)11 of pediatric tumors less aggressive than NB were reconstructed. NB was found to exhibit extensive branching even

prior to the outgrowth of the major subclone, while this was not the case for the more indolent tumors (Supplementary Fig. 7g:IV–V). In summary, scWGS confirmed that extensive phylogenetic

branching is a consistent feature of NB over time that can emerge already at ancestral stages with low-level _MYCN_ gain. High-risk NBs seem to have developed tolerance to chromosome breaks

already at the time of first branching, resulting in complex patterns of subclone diversity that can form the substrate of CCR, through which new clones emerge to replace those killed off by

chemotherapy. DISCUSSION The present study shows that genetic intratumor diversity in NB emerges by early phylogenetic branching, resulting in distinct subclone territories. This branched

architecture sets the ground for clonal evolution under therapy, a process which we show is extensively determined by whether the subclones dominating the primary tumor are eradicated or

not. If such dominant subclones survive chemotherapy, they seed direct descendants through linear evolution as disease progresses. If eradicated, they can leave room for replacement by a

genetically diverse panorama of collateral relatives through CCR. In the present study, CCR was found with effective treatment irrespective of whether patients had high-risk disease treated

with very intensive chemotherapy (Patients 4–5, 7–11) or had low/intermediate risk disease treated with less intensity (Patients 3 and 6). CCR was also recapitulated when comparing

COJEC-responsive PDX tumors to their local relapses, even though technical restrictions did not allow for an experimental setup perfectly mirroring the clinical patient scenario with

pre-post sampling. Altogether, CCR is consistent with selection under chemotherapy of pre-existing NB cell populations with potential for treatment resistance, as has previously been

described in a range of other malignancies34,35,36,37,38,39,40,41,42,43. The absence of CCR in cell culture by merely imposing a mechanical population bottleneck also supports Darwinian

selection under chemotherapy as one underlying mechanism. On the other hand, linear evolution was observed at progressive growth irrespective of whether chemotherapy was given or not

(patient 12, PDX3 C5, and PDX1 C1-C3 showed extensive linear radiation from detectable ancestral clones but were not exposed to chemotherapy). Previous paired comparisons of tumor genomes

from primary and relapsed NBs have identified a broad repertoire of driver mutations enriched at relapse, in particular RAS-MAPK pathway mutations8,9,10,33, supporting some degree of

selection under effective chemotherapy. However, the diversity among tumor genomes from different relapse patients and the notable absence of bona fide resistance mutations in NB in the

present and previous studies, indicate that heritable factors other than DNA variants could be the main substrate of selection. Recent studies of chemotherapy resistance in NB have indeed

indicated the presence of a spectrum of lineage development states based on super-enhancer-associated transcription factor networks44,45,46,47, where cells with a mesenchymal transcriptional

signature are enriched in post-chemotherapy and relapsed tumors compared to cells with an adrenergic signature. In fact, using the same PDX systems as in the present paper we have shown an

increased expression of early developmental signatures, with features of Schwann cell precursors in treatment resistant and relapsed lineages, and also an increased mesenchymal signature in

relapsed tumors versus an increased adrenal signature in PDXs cured by chemotherapy and surgery27. In this context, our data suggest that different genetic lineages present in the same

patient have a diverse potential for attaining the mesenchymal, chemoresistant state, a hypothesis also supported by the in vitro data on SK-N-SH cells presented here, where state transition

to chemotherapy resistance was accompanied by clear shifts in the clonal landscape. While a certain lineage dependence of such state-transition could be explained by epigenetic

heritability, our finding of distinct subclonal territories in clinical samples suggest that microenvironmental factors may also contribute to priming NB cells towards attaining

chemoresistance. Another alternative is that lineages dominating the primary tumor undergo negative selection at treatment because they are simply slower, albeit in the end not less capable,

than minor collateral populations to reprogram into the mesenchymal state. The present study is not capable of distinguishing these alternatives for several resons. First, it is based

largely on clinical paraffin-embedded samples, small in size, providing very limited opportunities for more comprehensive, phenotypic chatacterization. Second, our study is underpowered and

undersampled (only 12 patients) and should thus be regarded as largely explorative. Future studies should ideally include larger patient cohorts with a broader access to fresh tumor material

for phenotypic studies. An additional improvement for future studies would be a script-based, simulation-validated process that replaces the manual assignment of allele compositions and

clone size in the complex sublone scenarios often present in clinical NB samples. The present study provides important clues on how to assess neuroblastoma patients for the purpose of

precision oncology. Our results show that, at up-front progression new mutations are typically just added to pre-existing major clones, while regression followed by relapse implies a more

radical shift in the clonal landscape through CCR, mandating resampling for additional sequencing. Furthermore, our tracing of lineages back to anatomic territories reveals the critical

importance of taking spatial heterogeneity into account at precision oncology assessment. Subclones were densely packed in diagnostic biopsies and relapse samples, in a fashion where just a

few mm3 of nearby extra tissue in a biopsy would be rich in additional genetic information. In contrast, when evaluating the primary tumor after treatment, phylogenetically distant

populations were found dispersed across small islands of surviving tumor cells. Thus, after therapy, care must be taken to sample several spatially distant tumor regions to capture the full

repertoire of clones. This may prove practically challenging because surviving tumor cell populations are typically nested in vast areas of necrosis and reactive host tissue. A future

alternative for monitoring clonal evolution in NB could be liquid biopsy protocols. Several studies have shown great promise for using cell-free DNA to monitor the temporal dynamics of NB

genomes, particularly in high-risk patients48,49. However, it remains an open question whether these methods will have sufficient sensitivity for regular clinical use. In summary, the

present study shows that evolutionary branching early in NB pathogenesis sets the stage for clonal replacement under effective therapy, which has direct implications for how to sample these

tumors for genomic profiling. METHODS ETHICS STATEMENT The present study was in compliance with all relevant etica regulations: it was approved by the Regionala Etikprövningsnämnden (EPN)

under permit numbers L289-11 genomic analyses; updated as L796-2017 and L605-05 (biobanking; updated as L883-2018). For data sharing, this was complemented by permit from the Swedish Ethical

Review Authority (2023-01550-01). Written informed consent was obtained for all study subjects included from the children’s parents or from legal guardians. For patients ≥15 years of age,

direct consent was also obtained. Information about the study was given by a trained pediatrician and/or pediatric oncologist. There was no compensation to study participants. All animal

procedures were conducted according to the guidelines from the Regional Ethics Committee for Animal Research (permit no. M11-15, 19012-19). PATIENT COHORT To study alterations in

neuroblastoma genotype over the time of disease course, we reviewed pathology files of all patients <18 years of age diagnosed with and treated for a histopathologically verified

neuroblastoma in the Southern Healthcare Region of Sweden from 1998-2018. From this cohort of approximately 50 patients, we selected those from whom there were samples with >50% viable

tumor cells available for DNA-analysis, material from at least two different disease time points, and samples from at least two intratumoral locations from at least one of these time points.

This resulted in totally 12 patients with material available from two or more of the following phases: * (1) At the diagnostic procedure before treatment was started (pre-chemotherapy). For

most patients, core needle biopsies intersecting different tumor regions were obtained. Each biopsy core with viable tumor cells was analyzed as one or more separate samples; DNA from

biopsies were not mixed. * (2) At the surgical resection of tumor after treatment (post-chemotherapy). Here, islands of surviving neuroblastoma cells were identified by light microscopy and

isolated by microdissection of the corresponding paraffin blocks. Each island, or group of islands within the same 5 mm radius, was analyzed as a separate sample. One patient (Patient 8) was

also sampled (by open biopsy) under ongoing chemotherapy. * (3) At the time of metastatic relapse. This material consisted of either core needle biopsies or open biopsy fragments, where

each core needle biopsy/each fragment was analyzed as a separate sample. * (4) In the context of treatment resistant progression. This material consisted of open surgical biopsies at

metastatic sites (Patient 1) or sampling of primary tumor and metastases at autopsy (Patient 2). No specific analyses of sex and gender were performed as this was not within the scope of the

study. Patients were enrolled based on the diagnosis and the availability of biological samples. The cohort contains children of both sexes as expected because the disease under study

affects both sexes almost equally. There is a slight increased risk for boys and there is also a slight male predominance in our small study cohort. Sex was assigned and not self-reported as

most of the study participants were young children. The study cohort is too small for meaningful subdivision according to sex and specifying it at publication may threaten the privacy of

participants. IN VIVO CHEMOTHERAPY EXPERIMENTS The animal studies included correspond to previously published experiments by Mañas et al.27 (Figs. 3, 4) and Braekeveldt et al.25 (Fig. 7),

whereas the phylogenetic analyses are novel. All mice were monitored for weight loss and other signs of toxicity. According to the ethics protocol mice were sacrificed when tumors reached

1800 mm3 or at a humane end point due to overall health deterioration. Exceptions were not made from this protocol. Mice were housed in a controlled environment with ad libitum access to

food and water. Temperature was set to between 20 and 24 °C to ensure optimal thermal comfort for the mice. Relative humidity (RH) was maintained between 45-55%. Lights were turned on from

6:00 to 18:00 (12 h light/dark cycle), including dusk and dawn periods. NMRI nude were purchased from Taconic, while NSG mice were obtained from in-house breeding25,27. The number of animals

used for each experiment is specified in Figs. 3a, e and 4a. For cisplatin experiments (Fig. 3a) female 8–12 weeks old NMRI nude mice were used. For COJEC experiments with progressive

growth (PDX1; Fig. 3e) female 8–12 weeks old mice were used; controls and two COJEC-treated mice were NSG, the remaning COJEC treated were NMRI nude. For COJEC experiments with response

(PDX3; Fig. 4a) female 8–12 weeks old NMRI nude mice were used. PDX cells were cultured as free-floating tumor organoids in serum-free stem cell medium26,50. Mycoplasma infection was tested

for prior to the experiments, using standard methods, and was confirmed to be negative. PDX cells (2 ×106) were suspended in a 100 µl mixture of the medium and matrigel (Corning, Cat

No.354234) (3:1) and injected subcutaneously into the flanks of female mice. Tumor size was measured using a digital caliper and calculated with the formula V = (πls2)/6 mm3 (where l is the

long side and s is the short side). Mice were randomly allocated to control, cisplatin, COJEC or resection groups once their tumor had reached approximately 500 mm3. The cisplatin group was

treated with intraperitoneal injection of 4 mg/kg cisplatin (Selleckchem, Cat No.S1166) each Monday, Wednesday and Friday, and the control group was treated with equal amounts of saline. The

COJEC group was treated with intraperitoneal injections of cisplatin (1 mg/kg; Selleckchem, Cat No.S1166) and vincristine (0.25 mg/kg; Santa Cruz, sc-201434) on Mondays, etoposide (4 mg/kg;

Santa Cruz, sc-3512) and cyclophosphamide (75 mg/kg; Santa Cruz, sc-361165) on Wednesdays, and carboplatin (25 mg/kg; Santa Cruz, sc-202093A) on Fridays. When the tumor size had reduced to

approximately 200 mm3, the tumor was surgically removed and mice were then observed for tumor regrowth. For the resection only group, >90% of the untreated tumor was surgically removed

and the mice were then observed for tumor regrowth. PDXs established through orthotopic implantation of undissociated patient tumor fragments into the paraadrenal space of NSG mice using

Matrigel, were subjected to serial in vivo orthotopic passaging up to eight generations for each model. Tumor samples from in vivo generations 1 and 7 were analyzed here. IN VITRO

CHEMOTHERAPY EXPOSURE Two neuroblastoma cell lines were used; IMR-32 (ATCC CCL-127) was maintained in RPMI-1640 medium (Gibco, Cat No.21875-034) and SK-N-SH (a gift from Dr. Daniel Bexell,

Lund university) was maintained in MEM (Gibco, Cat No.31095029). Both media were supplemented with 10% fetal bovine serum and 1% penicillin-streptomycin and the cells were grown at 37 °C and

5% CO2. The cells were confirmed to be mycoplasma negative prior to the experiments. For passaging and cell counting, cultured cells were washed by phosphate-buffered saline (PBS),

disaggregated by trypsin (HyClone, Cat No.SH30236.01) and centrifuged at 300 × _g_ for 5 min. Cells were collected and saved as a frozen pellet at –80 °C until DNA was extracted. As

reference, the cell cultures from which cells were seeded for each experiment (mother culture) were analyzed. Six different experiments were performed with the IMR-32 cells (Supplementary

Fig. 4), experiments 1–3 (samples A, B and C) were performed in one batch, followed by experiment 4 (samples SB), and lastly experiments 5 and 6 (samples MB cis and MB). * (1) _IMR-32

CONTINUOUS GROWTH WITHOUT BOTTLENECK_ (Supplementary Fig. 4b: III samples A1-A3): One million cells were seeded, and the medium was replaced one day after seeding. One third of the cells

were passaged to new T25 flasks twice per week. In total, the cells were passaged 11 times during 42 days. * (2) _IMR-32 LONG-TERM LOW-DOSE CISPLATIN EXPOSURE_ (Supplementary Fig. 4b: III

samples B1-B3): One million cells were seeded. The next day, the original medium was replaced with medium containing 0.1 μM cisplatin. The medium was changed twice per week to new

cisplatin-containing medium. Cisplatin suppressed cell growth in the beginning but the population started to expand while exposed to cisplatin at the end of the experiment. Overall, the

medium was changed 12 times and the cells were passaged twice during 42 days. * (3) _IMR-32 SINGLE BOTTLENECK INDUCED BY HIGH-DOSE CISPLATIN_ (Supplementary Fig. 4b: III samples C1-C3): One

million cells were seeded, and the cells were treated with 1.0 μM cisplatin for 72 h with start the day after seeding. This reduced the number of cells with 94% compared to cells grown in

absence of cisplatin. Viable cells were cultured until confluent in a 6-well plate. This took 42 days which also became the endpoint for experiments (1) and (2) in order to remove time as

confounder for experiments (1)-(3). * (4) _IMR-32 SINGLE BOTTLENECK WITHOUT CISPLATIN_ (Supplementary Fig. 4c: III samples SB1-SB3): One million cells were seeded and new medium was added

the following day. After 72 h the same number of cells as in group (3) was collected and re-seeded in order to create a bottleneck with similar cell number compared to cells treated with

high-dose cisplatin. The cells were cultured until they reached confluence in a 6-well plate. The total experimental time was 14 days. * (5) _IMR-32 MULTIPLE BOTTLENECKS INDUCED BY

CISPLATIN_ (Supplementary Fig. 4d: III samples MB cis1-MB cis3): One million cells were seeded, and the cells were treated with 0.75 μM cisplatin the day after seeding. 0.75 μM cisplatin was

used because IMR-32 cells could not tolerate three times exposure of 1 μM as used for the single bottleneck experiment. Viable cells were counted after 72 h exposure to cisplatin and

540,000 cells (mean number alive cells of the three cultures) were cultured until confluence of a 6-well plate. This step was regarded as one bottleneck. One million cells were collected

from the 6-well plate and re-seeded to a T25 flask. The same procedure was repeated two more times in order to create three bottlenecks. The number of viable cells increased after the second

and third exposure to cisplatin compared to the first exposure due to cisplatin resistance as in experiment (2), but the same number of cells (540,000) were seeded at each round. Cells were

cultured until they became confluent in a 6-well plate after the third bottleneck. The total experimental time were 53, 46 and 49 days for MBcis1-MBcis3, respectively. * (6) _IMR-32

MULTIPLE BOTTLENECKS WITHOUT CISPLATIN_ (Supplementary Fig. 4d: III samples MB1-MB3): One million cells were seeded and the medium was changed one day after seeding. As in experiment (5),

540,000 cells were reseeded after 72 h in order to create the same size of bottleneck as for MBcis1-MBcis3. These cells were cultured until they reached confluence in a 6-well plate. This

step was regarded as one bottleneck. Two more bottlenecks were created and the cells were collected after reaching confluence in a 6-well plate. The total experimental time was 42 days. To

validate the findings in chemotherapy treated IMR-32 cells, two different experiments were performed with SK-N-SH cells in one batch (Supplementary Fig. 4f). * (1) _SK-N-SH BOTTLENECKS

INDUCED BY CISPLATIN_ (Supplementary Fig. 4f samples cis1.1-5.1 and cis1.2-5.2): An experimental setup involving two consecutive bottlenecks, was performed. Cells to be used as controls of

the original culture were collected at two time points, one at day zero (CC.1) and one three days later (CC.2). At day zero, one million cells were seeded in eight T25 flasks. At day one,

media was changed in all eight cultures and cisplatin was added to five of the flasks to a final concentration of 4 μM. The cells were treated for 72 h resulting in approximately 80% cell

death. Three days later (day 7), the cisplatin treated cells had started to proliferate and their morphology were changed to large substrate adherent cells, making the culture almost

confluent. Cells from each flask (500000) were re-seeded and 4.0 µM cisplatin was added again the following day. The remaining cells from each flask (cis1.1 - cis5.1) were saved as frozen

pellets to be analyzed by SNP-array. After the second hit with cisplatin it took 25 days for cells to become almost confluent and then, at day 33, these cells were pelleted and saved for

analyses (cis1.2 - cis5.2). * (2) _SK-N-SH CONTINUOUS GROWTH WITHOUT BOTTLENECK_ (Supplementary Fig. 5f: samples ctrl 1.1–3.1 and 1.2-3.2): Cells in three untreated cultures were

continuously passaged at 1:3 approximately twice a week (in total eight times) and fractions of the cells were collected as time point controls (ctrl1.1- ctrl3.1 respectively ctrl1.2-ctrl

3.2) in parallel with the cisplatin treated cells. We did not get SNP-array data from cis3.2 due to insufficient DNA amount. DNA EXTRACTIONS From fresh frozen tumor and blood samples, DNA

was extracted using the DNeasy blood and tissue kit (Qiagen, Cat No. 69504). DNA from PDX snap frozen tumor pieces and from frozen IMR-32 and SK-N-SH cell pellets were extracted using the

AllPrep DNA/RNA Mini Kit (Qiagen, Cat No.80204) with standard methods. For formalin fixed paraffin-embedded tissue (FFPE) DNA was extracted using the Allprep DNA/RNA FFPE kit (Qiagen, Cat

No. 80234) with the following modification: during deparaffination, two 10-minute incubations of the samples with 1 ml xylene at RT were performed to ensure complete removal of paraffin. DNA

concentration was measured using Qubit DS HS DNA (Invitrogen Cat No. Q32854). OVERALL APPROACH TO GENOMIC PROFILING AND PHYLOGENETIC ANALYSIS The high proportion of CNAs compared to driver

point mutations complicates evolutionary analyses of NB. Many types of chromosomal rearrangements violate the infinite sites hypothesis on which conventional phylogenetic analyses rest51. In

particular, aneuploidy can affect only the limited number of chromosomes present in a cell and is thus frequently the subject of parallel evolution and back mutation. Applying phylogenetic

methodology to NB cells thus requires adaptation of methods to incorporate features of chromosome-level evolution, for example by allowing parallel evolution of whole chromosome alterations

and back mutation from whole chromosome gains while treating chromosomal breakpoints outside centromeric regions as unique events and loss of heterozygosity as irreversible52. We have

recently reported a software tool (DEVOLUTION, https://github.com/NatalieKAndersson/DEVOLUTION) that allows input of point mutations in parallel to CNA data, where the latter is curated

according to the constraints of chromosomal evolution22. DEVOLUTION integrates sequencing data and CNA profiles by clustering genetic alterations based on their clone sizes (mutated clone

fractions) across samples. These clusters, together with the criterion that the sum of parallel subclones must not exceed 100% in any biopsy is used to deduce the most likely subclone

configurations across biopsies, including if subclones are nested within or parallel to each other. The algorithm generates an event matrix encompassing the detected subclones. This matrix

was used as input for phylogenetic reconstruction to deduce their evolutionary relationship using the maximum parsimony method or maximum likelihood method employing the R package phangorn

(v2.8.1)53, and visualized using the ggplot2 (v.3.3.5) package. Phylogenies were rooted in a cell having no genetic alterations. The manual “How to build MP and ML trees from a segment file

using DEVOLUTION” is available at Zenodo. To perform high-resolution mapping of subclonal territories across anatomic tumor space in NB patients we used archived FFPE tumor tissue. This

allowed us to make a comprehensive investigation of tumor cells left in primary tumors after chemotherapy, even in cases with good chemotherapy response where survivor populations are often

limited to spatially dispersed tumor islets a few millimeters in size, surrounded by necrosis, hemorrhage and fibrosis. Areas of ∼20 mm2 containing >30% tumor cells were subjected to

whole genome copy number profiling parallel to whole exome sequencing (WES). Mutations detected by WES were validated by targeted resequencing using a custom-made panel based on the WES

findings. Analyses were done on FFPE tumor in all cases except three (Patient 6, Patient 9 and Patient 12) from which DNA was extracted from frozen tumor samples. Combining CNA and SNV/indel

data resulted not only in an increased number of mutations being found in each sample, but also the number clones that could be deconvolved from each sample (Supplementary Fig. 1m, n).

Tumor cell fraction had a slight impact on the number of SNVs/indels detected in each sample but no significant impact on the number of detected CNAs. An overview of the workflow from raw

SNP array and sequencing data to phylogenetic trees can be found in Supplementary Fig. 1p. COPY NUMBER ANALYSES BY SNP ARRAY CytoScan HD SNP array (Thermo Fisher/Affymetrix) was used and the

copy number profiles were visualized with Chromosome Analysis Suite (version 3.3.0.139; Thermo Fisher/Affymetrix) to analyze DNA from fresh frozen samples and cell culture experiments. CEL

files were normalized by the R package Rawcopy54, and scatter plots whose x-axis and y-axis denotes log2 ratio and allelic imbalance, respectively, were generated through the R package TAPS

(Tumor Aberration Prediction Suite)55. DNA extracted from FFPE samples were analyzed with the Oncoscan array platform at the Swegene Center for Integrative Biology (SCIBLU) according to the

manufacturer’s recommendations. The CEL-files were converted to OSCHP-files with the OncoScan Console software (version 1.3.0.39; Affymetrix). The Nexus Express software for Oncoscan

(version 3.1) was used for visual inspection of the OSCHP-files. Three text-files containing information regarding probes, segments and SNPs were exported from the Nexus Express software and

used as input to create TAPS plots. Only copy number segments larger than or equal to 0.1 Mb were included in downstream analyses. Constitutional segments were excluded via the Database of

Genomic Variants56. Aberrations identified in the clinical samples were summarized in circos plots, one plot per patient with individual samples illustrated as concentric circles. These

plots were generated with BioCircos (v0.3.4) in R (v4.1). WHOLE EXOME SEQUENCING OF PDX SAMPLES Equal amounts of DNA from two samples from the same PDX were combined and analyzed with whole

exome sequencing at the Center for Translational Genomics at Lund university as detailed in ref. 57. The Sure Select XT HS Library Prep (Agilent) was used and paired-end sequencing (2×150

cycles) was performed on a NextSeq500 (Illumina). bcl2fastq 2.20 (Illumina) was used to convert files to fastq-formatted files and paired end reads were mapped to the human reference genome

(GRCh37 with decoys from the 1000 Genomes’ Project) using BWA-MEM 0.7.1558. Duplicate reads were marked using sambamba 0.6.759. Somatic variant calling was performed using freebayes (with

the –pooled-continuous, --pooled-discrete and -F 0.03 flags; https://arxiv.org/abs/1207.3907v2) and strelka2 (using the --exome flag)60. Contaminating reads from the mouse genome were

removed. A variant was excluded from the study if detected by less than 10 reads, had a variant allele frequency (VAF) below 0.1 in all samples, was present in more than 1% of the reads in

the corresponding normal sample and/or showed inferior quality at visual inspection of the reads in Integrative genomics viewer61. The Polyphen-2 tool62, was used to predict the impact of

mutations on protein function. WHOLE EXOME SEQUENCING AND TARGETED SEQUENCING OF CLINICAL SAMPLES DNA repair was performed on DNA extracted from FFPE-material using the NEB PreCR kit (New

England Biolabs, Cat.No. M0309S) following the manufacturer’s recommendations, using 500 ng of DNA. The repaired DNA was purified using AMPure XP beads (Beckman Coulter, Cat No. A63881) at a

1.8X ratio and eluted in Illumina Truseq Exome kit RSB buffer (Illumina). Library preparation was performed using the Truseq Exome kit (Illumina, Cat No.20020614). DNA fragmentation was

performed on a Covaris S220, and end repair, size selection, 3’ end adenylation, adaptor ligation and DNA enrichment steps were performed according to the protocol recommendations. Libraries

were quantified with the Qubit DS HS DNA kit and pooled in an equimolar manner. Exome capture, amplification (8 cycles) and purification steps were performed as recommended. Samples were

quantified using Qubit dsDNA HS kit, and library quality was assessed on Bioanalyzer, using the Agilent High Sensitivity DNA kit (Agilent Cat No. 5067-4626) before sequencing. DNA samples

were analysed either in-house with the NextSeq 500/550 high output kit v2.5, 300 cycles (Illumina, Cat No. 20024908) on a NextSeq500 (Illumina) or at the Center for Translational Genomics at

Lund University using the NovaSeq 6000 S1 reagent kit v1, 300 cycles (Illumina, Cat No. 20028317) on the NovaSeq6000 (Illumina). Illumina paired-end reads were aligned to the human

reference genome hg19 by BWA-MEM (https://arxiv.org/abs/1303.3997v2). Duplicate reads marking and local realignment were performed by GATK (version 4.0.11.0)63. Mutect (version 1.1.7)64,

GATK Mutect2 (version 4.0.11.0), and MuSE (version v1.0rc)65 were used to identify somatic single nucleotide variants (SNVs) and small insertions/deletions (indels). The SNVs and indels

called by Mutect2 were further filtered with the GATK FilterMutectCalls. The variant vcf files were converted into maf files by the vcf2maf package (https://github.com/mskcc/vcf2maf). DNA

from fresh frozen clinical samples from Patient 6, Patient 9 and Patient 12 were analyzed in the same way as the PDX samples, but the variant calling was performed jointly for all tumor

samples from each patient using Mutect2 (https://doi.org/10.1101/861054) with a panel consisting of 20 unrelated exome samples used as normals, followed by variant filtering with

FilterMutectCalls including the read orientation filter65, according to the best practice guideline for somatic short variant discovery66. Variants were excluded using the same criteria as

for the sequencing of PDX samples but a VAF cut-off at 0.2 was used for the FFPE extracted samples. Very few mutations were detected in samples from Patient 3 and Patient 7, hence a VAF

cut-off value at 0.1 was used for these samples. Targeted re-sequencing of 454 mutation in total that remained after filtering of WES data as specified above was performed with at the Center

for Translational Genomics at Lund University using a Twist Custom probe Capture panel with the Twist library preparation EF Kit (product no. 101058, Twist Bioscience). The samples were

sequenced on a NovaSeq 6000 (Illumina). Raw reads were processed through the Sentieon® unique molecular indices (UMI) aware pipeline (https://support.sentieon.com/appnotes/umi/). Briefly,

UMIs were extracted from raw fastqs and raw reads were aligned to the reference genome (GRCh37) using the sentieon implementation of bwa mem (https://arxiv.org/abs/1303.3997v2), followed by

consensus fastq generation using the sentieon consensus tool (https://www.sentieon.com/products/). The consensus reads were then mapped to the reference genome (again using the sentieon

implementation of bwa mem) and the resulting bam-files were processed through freebayes (https://arxiv.org/abs/1207.3907) with the following settings: --pooled-continuous, --pooled-discrete,

--min-repeat-entropy 1, -F 0.03 on a per patient basis, i.e. all samples from each patient were called jointly. Variants were excluded if detected in >1% in the corresponding normal

sample, had a VAF value below 0.1 in all samples and were covered <100 reads in all samples. The mean total coverage of the genomic location for mutations included in the study was 545x.

All variants were annotated using the VEP tool (https://grch37.ensembl.org/Homo_sapiens/Tools/VEP). To identify relevant mutations to highlight in the phylogenetic trees (Figs. 1–2,

Supplementary Fig. 1) and heat maps (Supplementary Data 1d), missense mutations flagged to have an impact by PolyPhen or SIFT in genes correlated to neuroblastoma, tumors in general or

neuronal differentiation according to Pubmed (https://pubmed.ncbi.nlm.nih.gov) were selected. The transcript with the highest predictions score was shown at protein level. CLONE SIZE

ESTIMATIONS BASED ON SNP ARRAY DATA A log2 ratio (log2R) value was obtained from SNP-array analyses for all identified aberrations. This is the ratio between the total number of alleles at

the analyzed location and a standardized normal copy number which can be described as: $$\log\!2{{{\rm{R}}}}=\frac{{{{\rm{MSF}}}} * {{{\rm{Nt}}}}+\left(1-{{{\rm{MSF}}}}\right) *

{{{\rm{Nb}}}}}{{{{\rm{Np}}}}}$$ (1) The mutated sample fraction (MSF) is the portion of the entire sample, including normal cells, that harbors the aberration, Nt is the copy number of the

aberration, Nb is the copy number of the background cells (either the number of alleles in a parallel clone or the ploidy of the tumor) and Np is the ploidy level the array is normalized

against. The MSF can be calculated with: $${{{\rm{MSF}}}}=\frac{{{{\rm{Np}}}} * \log\!2{{{\rm{R}}}}-{{{\rm{Nb}}}}}{{{{\rm{Nt}}}}-{{{\rm{Nb}}}}}$$ (2) For array data from Cytoscan HD, the

log2 value needs to be divided by a correction factor of 0.53-0.667. The MSF value can also be calculated using the mirrored B-allele frequency (mBAF) where the number of B-alleles is higher

or equal to the number of A-alleles. $${{{\rm{mBAF}}}}=\frac{{{{\rm{Btot}}}}}{{{{\rm{Atot}}}}+{{{\rm{Btot}}}}}$$ (3) The total number of B alleles, Btot, in a sample at a specific location

is the sum of the B-alleles in the tumor (NB) and the B-alleles in the normal cells, which is assumed to be one.

$${{{{\rm{B}}}}}_{{{{\rm{tot}}}}}={{{{\rm{B}}}}}_{{{{\rm{tumor}}}}}+{{{{\rm{B}}}}}_{{{{\rm{normal}}}}}={{{{\rm{N}}}}}_{{{{\rm{B}}}}}\, {{{\rm{x}}}} \, {{{\rm{MSF}}}}+1\, {{{\rm{x}}}}\,

\left(1-{{{\rm{MSF}}}}\right)$$ (4) The number of A-alleles can be described accordingly $${{{{\rm{A}}}}}_{{{{\rm{tot}}}}}={{{{\rm{N}}}}}_{{{{\rm{A}}}}} \, {{{\rm{x\; MSF}}}}+1 \,

{{{\rm{x}}}}\, \left(1-{{{\rm{MSF}}}}\right)$$ (5) If combining (3), (4), and (5), and simplifying the expression, MSF can be calculated with

$${{{\rm{MSF}}}}=\frac{1-2{{{\rm{mBAF}}}}}{{{{\rm{mBAF}}}}\left({{{{\rm{N}}}}}_{{{{\rm{A}}}}}+{{{{\rm{N}}}}}_{{{{\rm{B}}}}}-2\right)-{{{{\rm{N}}}}}_{{{{\rm{B}}}}}+1}$$ (6) The mBAF-value can

be obtained from the TAPS plot since the allelic imbalance (AI) on the y-axis of these plots is defined as $${{{\rm{AI}}}}=\frac{{{{\rm{mBAF}}}}-0.5}{0.5}$$ (7) CLONAL DECONVOLUTION BASED

ON CNAS DETECTED BY SNP ARRAY IN CLINICAL SAMPLES The SNP-array data were visualized in TAPS plots55, where all chromosome segments in one sample are shown in the same plot with the log2R on

the x-axis and the allelic imbalance on the y-axis. Thus, all normal chromosomal segments are displayed in the same cluster with the log2R and allelic imbalance equal to zero, while clonal

CNAs appear in a grid like pattern where each junction represents a specific combination of A and B-alleles. Subclonal CNAs can be identified in the grid since their clusters deviate from

these junctions. All CNAs were annotated with the lowest number of alleles possible, for example a CNA situated between the 1 + 2 and 1 + 3 junctions in the TAPS plot was annotated as a

subclonal 1 + 3. An aberration was called as “mixed” if present with two different allelic compositions in the same sample, as indicated by cluster positions in TAPS plots deviating from the

grid. Events that were too small to be visualized in the TAPS plot were excluded from the study since the allelic imbalance could not be determined with certainty. The median of MSF values

for aberrations that appeared to be clonal based on the TAPS plot were used to establish the tumor cell fraction (TCF) or purity of the sample. This was done both based on the log2R-values

and on the allelic imbalance. Similar results (difference <approx. 0.1) from these calculations for the same aberrations was used as an indication of correct assessment of ploidy-level

and allelic composition. The mutated clone fraction MCF, describing the fraction of the tumor cells that contains an aberration were then calculated as

$${{{\rm{MCF}}}}=\frac{{{{\rm{MSF}}}}}{{{{\rm{TCF}}}}}$$ (8) The MCF values based on the log2R was used for most aberrations, with the exception for CNNIs that don’t affect the log2R -data,

for rare occasions when the log2R was noisy while the data was stable based on mBAF-values or when the visual appearance of the TAPS plots reflected the AI based MCF-values the best.

Amplifications and homozygous deletions were always scored as clonal. The MSF value for a “mixed” aberration was calculated for the largest of the subclones first and the value for the

smaller subclone was calculated as 1-the MSF for the biggest. The MCF values were visualized in a heat map format and variants with similar MCF pattern over samples were grouped together

(Supplementary Data 1d). An average of the MCF values per sample and group were calculated and that value rounded to the nearest 0.1 became the final clone size for all the variants in that

group. When the appearance of the heat map indicated that a consistent phylogenetic tree could not be created, an event could be added as an “inferred stem” or “inferred private”. For

example if all samples derived from one patient demonstrated a gain with the same breakpoints but containing different number of alleles in different samples, the gain with the lowest number

of alleles was inferred as stem in all samples. Another example could be if a whole chromosome aberration was present in several samples, and it was impossible to visualize a logic order of

events for the identified mutations, the whole chromosome events could be annotated as inferred private, i.e. separate events, since they are likely to appear in parallel. Such occasions

are noted in Supplementary Data S1d. The information from the curated heat maps was used as input data for DEVOLUTION to construct phylogenetic trees. Two subclones containing the same

aberration but with different allelic compositions (the scenario we called “mixed”) cannot be nested events since that would indicate that cells in the nested population would have two

different CNA states at the same time. In such cases, we provided the DEVOLUTION algorithm with a matrix indicating that these genetic alterations should not be nested into each other to

avoid illicit evolutionary trajectories. The maximum likelihood and the maximum parsimony trees obtained from DEVOLUTION were identical for all patients with the exception of the trees based

on samples from Patient 6. Here mainly whole chromosome aberrations were detected and the absence of distinct breakpoints indicative of a relationship between subclones made it difficult to

differentiate between parallel private and shared events. In this case, the tree containing the smallest number of back mutations (maximum parsimony) was chosen. CALCULATION OF CLONE SIZES

BASED ON INTEGRATION OF VARIANT ALLELE FREQUENCIES FROM EXOME SEQUENCING WITH COPY NUMBER FROM SNP-ARRAY DATA The variant allele frequency (VAF) is derived from sequencing data and describes

the number of reads detecting a specific mutation divided by all reads originating from tumor cells and diploid normal cells in the sample. $${{{\rm{VAF}}}}=\frac{{{{\rm{M}}}}\times

{{{\rm{MSF}}}}}{{{{\rm{CN}}}}1\times {{{\rm{f}}}}1+{{{\rm{CN}}}}2\times {{{\rm{f}}}}2+2(1-{{{\rm{TCF}}}})}$$ (9) where MSF is the mutated sample fraction, i.e. the proportion of the entire

sample, including the normal cells, that contains the mutations. M is the number of mutated alleles, TCF is an average of the log2R and allelic imbalance values described above, and CN1 is

the total number alleles in the main fraction of the tumor at the location of the mutation, defined by SNP-array, and f1 is the size of fraction 1. CN2 is the number of alleles in a

background clone if present, and f2 is the size of CN2, and $${{{\rm{TCF}}}}={{{\rm{f}}}}1+{{{\rm{f}}}}2$$ (10) If Eqs. (8), (9) and (10) are rearranged and combined, the MCF can be

calculated as $${{{\rm{MCF}}}}=\frac{{{{\rm{VAF}}}}\times ({{{\rm{CN}}}}1\times {{{\rm{f}}}}1+{{{\rm{CN}}}}2\times ({{{\rm{TCF}}}}-{{{\rm{f}}}}1)+2\left(1-{{{\rm{TCF}}}}\right))}{{{{\rm{Mx\;

TCF}}}}}$$ (11) As a starting point, M was assumed to be 1 for diploid samples and 2 for tetraploid samples. For mutations located on a subclonal CNA or on a mixed background (two subclones

with different CNAs), M could take different values. For example, if two different CNAs were present at the location of a specific mutation at for example 2 + 0 at 30% and 1 + 2 at 70%, M

for that variant could be 0.7, 1.4 or 2, depending on the order of events. If several theoretical values for M were possible, the following rules were used to choose the M-value to be used

for the MCF-calculation: * If the MCF value became higher than 1.2, a higher value of M was chosen * The choice of the M-value should not lead to a biologically impossible order of event,

for example a completely lost (nullisomic) allele cannot be gained at a later stage. * The generated MCF values should not violate the pigeon whole principle * The order of events for a CNA

and a SNV should be the same for all samples harboring identical aberrations If more than one solution for M still remained, the value was chosen that introduced the fewest number of novel

subclones so as not to overcall heterogeneity. The value of M chosen for calculation of MCFs for the specific variants are included in Supplementary Data S2b. As benchmarking to our

heuristic approach, we also re-estimated the MCF values with the phylogenetics tool DeCiFer v2.1.468. It can notably quantify the proportion of cells which acquired an SNV or whose ancestors

acquired the SNV inferring corresponding multiplicities across samples. We compared the MCFs estimated with DeCiFer, which extensively confirmed our calculations (see Supplementary Fig. 8

and its legend for details). DeCiFer uses density-based cluster assignment with variable epsilon boundary that can be defined with the n-sigma limit. By default, a fairly strict grouping is

used (1-sigma limit), which means many data points from our dataset were not assigned to any cluster. In addition, we were unable to compare results of Patient 3, a pentaploid case due to

state trees being unavailable in DeciFer for the copy number alterations in question. CLONAL DECONVOLUTION FROM SNP ARRAY DATA OF CELL LINES AND PDX CELLS CNAs and allelic imbalances were

classified as clonal or subclonal based their relative distribution in whole genome TAPS-plots. Subclonal events were subjected to a more precise clonal deconvolution. The MCF values were

calculated for all aberrations using Np=2, and if the MCF value was over or equal to 90 %, the aberration was considered as clonal. The final clonal deconvolution for phylogenetic analysis

was based on the following rules: * Because of an error margin of 10% at clone size estimation by SNP array69, final clone sizes were rounded to the nearest 10%. The sum of subclone sizes

was allowed to reach a maximum of 120%, because two concomitant subclones will yield a 20% error margin. * Tumor cell concentration was normalized to 100% when the total sum of subclones

reached over 100%. * CNNI events (2 + 0) in a sample against a background of trisomy (1 + 2) was assumed to occur by losing an allele from the trisomy if the majority of the other samples

from the same experiment displayed trisomies (1 + 2). For cell line experiments, a private aberration was defined as an aberration unique for each lineage, while for PDX experiments (two

samples per tumor) a private aberration was defined as an aberration unique for each PDX tumor. All other aberrations were defined as shared or stem. CLONAL DECONVOLUTION FROM SEQUENCING

DATA FROM PDX MODELS SNP array data from the two samples from the same PDX tumor were combined and segmental aberrations present in the main clone were used to calculate the mutated clone

fraction (MCF) for variants. Because reads from murine DNA were removed at analyses of PDX samples, the tumor cell content was assumed to be 100 %. The majority ( >90%) of the detected

mutations were located on normal diploid alleles or on segments that were clonally deviant, hence f1 = 1 and the MCF was calculated by: $${{{\rm{MCF}}}}=\frac{{{{\rm{VAF}}}}\times

{{{\rm{CN}}}}}{{{{\rm{M}}}}}$$ (12) Mutations with an MCF > 0.7 were considered to be clonal. A variant detected on an unbalanced segment (1 + 2) was assumed to be clonal if the same

variant was judged as clonal in a balanced setting. This was achieved by choosing the M that produced an MCF close to 1 in the unbalanced setting. The cut-off value for clonality was lowered

to MCF > 0.6 for a mutation detected on a 1 + 1 segment that was shared between several samples and clonal in the majority (90 %) of them. A shared variant present at a low frequency

(VAF < 0.20) in a sample where a clonal sweep had taken place was considered as a technical artifact and removed. DECONVOLUTION OF CLINICAL SAMPLES AFTER TARGETED RE-SEQUENCING MCFs were

calculated by using the average TCF (purity) values obtained by calculations from the log2-ratio and from the allelic imbalance values from the TAPS plots. The number of mutated alleles M

was assumed to be 1 in a diploid setting but changed to 2 if the MCFs > 1.2 in the majority of the samples where the mutation was identified. Our inclusion criteria for variants (see

above) could result in MCF-values above 1.2 in some samples. This could be due to for example low coverage in the specific sample or sublconal background polyploidization. These variants

were kept in the study since they were clearly detected. Their MCF-values were rounded down to 1.0 and thus did not induce false heterogeneity. The final MCF values were obtained by using

the same heat map approach as for the MCF calculations on SNP array data on clinical samples. After quality filtering and deconvolution 359 of the 454 variants from the TWIST design were

retained for phylogeny construction by DEVOLUTION. QUANTITATIVE ANALYSES OF GENETIC DIVERSITY AND PHYLOGENETIC TREES The index of genomic diversity (IGD_)_ was constructed to minimize the

impact of differences in number of aberrations detected per sample and differences in the number of samples per patient, when comparing the diversity between sample types (i.e. samples taken

pre or post treatment or from a metastasis). The lengths of the branches (dSi) starting from the closest common node for all subclones (n) identified for a sample type, to all (n) the

individual subclones (dS1, dS2, … dSn), were annotated together with the distances (DSi) from the start of the stem to the same subclones (DS1, DS2,..DSn). IGD was calculated as the ratio

between the sums of these two different distances:

$${{{\rm{IGD}}}}=\left({{{{\rm{dS}}}}}_{1}+{{{{\rm{dS}}}}}_{2}+..{{{{\rm{dS}}}}}_{{{{\rm{n}}}}}\right)/\left({{{{\rm{DS}}}}}_{1}+{{{{\rm{DS}}}}}_{2}+..{{{{\rm{DS}}}}}_{{{{\rm{n}}}}}\right)$$

(13) IGD varies between 0 and 1, where a higher IGD-score indicates more genetically diverse subclones. Irregularity describes how the phylogenetic tree deviates from a symmetric star

phylogeny. Branch lengths, i.e. the distances from the stem to the subclones in the tree are annotated and irregularity is calculated as the variance of all branch lengths. RECONSTRUCTING

SUBCLONE GEOGRAPHY Leveraging that the paraffin-embedded material was regularly sliced down to around 3 mm depth across biopsies and patients to obtain material for DNA extraction allowed us

to approximate territories across a two-dimensional spaces, with surface areas corresponding to subclone sizes. Relative subclone sizes were calculated as MCFs as described above. While

subclones detected in more than one tissue samples were assumed to be located around the interface between these biopsies, subclones ascertained only in one biopsy was placed in the

remaining space. If found only in a single sample, nested subclones were placed in the geographic center of their mother clones. SINGLE CELL WGS Single cell copy number profiles from shallow

WGS were analyzed for clonal identification used in phylogenetic analyses and fishplots70, of PDX tumors and tumor organoids. Identified copy numbers encompassing less than five consecutive

1 Mb bins in the dataset, were excluded from further analysis, except for high-grade amplifications of the _MYCN_ region in chromosome arm 2p and of the _MAML3_ region in 4q28.3-q31.1, in

which case two consecutive 1 Mb bins were considered true gains11. Single cell copy number profiles from shallow WGS were used as input to a custom-made algorithm, latest version found on

github (https://github.com/NatalieKAndersson/SCEM). The input file is a matrix where each column is a single cell, each row a chromosomal segment, and the matrix elements are the copy

numbers of chromosomal segments in single cells. The algorithm identifies all genetic alterations present in each single cell while applying an optimized cut-off of five consecutive 1 Mb

bins for identifying true imbalances, and combines the information from all analyzed single cells to identify consecutive events encompassing overlapping chromosomal segments. The output was

an event matrix that was used to perform phylogenetic reconstruction to deduce the evolutionary relationship between individual cells using the maximum parsimony method or maximum

likelihood method employing the R package phangorn (v2.8.1)53, and visualized using the ggplot2 (v.3.3.5) package. Phylogenies were rooted in a cell having no genetic alterations. REPORTING

SUMMARY Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article. DATA AVAILABILITY The datasets generated or analysed during the

current study are available in curated raw data format in Supplementary Data 2. Raw data plots from whole genotyping of clinical samples, PDX tumors and cell lines are available at Zenodo

(https://doi.org/10.5281/zenodo.8017767). Input files used for DEVOLUTION and the subclone annotation derived from running DEVOLUTION are available at Zenodo together with single cell WGS

input files for SCEM (https://zenodo.org/records/10727558). Raw data have been deposited the European Genome-phenome Archive (https://ega-archive.org) under accession number EGAS00001007650,

with Dataset IDs as follows: -Whole genome genotyping/targeted sequencing: EGAD00001015012 -Whole exome sequencing: EGAD00001015413 -Single cell whole genome sequencing: EGAD00001015414.

Access to raw genomic data will require assessment by a data access committee of the applicant’s ethics permit as well as technical specifications of how the data should be stored to ensure

compliance with ethics regulations of the study as well as GDPR. A data transfer agreement will have to be made and signed by legal representatives. Source data are provided with this paper.

CODE AVAILABILITY The software DEVOLUTION was used for clustering and phylogenetic analyses of bulk sequencing data (https://doi.org/10.5281/zenodo.13304334)71. The full software including

access links is published in Andersson et al.22 The software SCEM used to construct phylogenies from single cell WGS is available online at https://github.com/NatalieKAndersson/SCEM.

REFERENCES * Maris, J. M., Hogarty, M. D., Bagatell, R. & Cohn, S. L. Neuroblastoma. _Lancet_ 369, 2106–2120 (2007). Article  CAS  PubMed  Google Scholar  * Morgenstern, D. A. et al. The

challenge of defining “ultra-high-risk” neuroblastoma. _Pediatr. Blood Cancer_ 66, e27556 (2019). Article  PubMed  Google Scholar  * London, W. B. et al. Clinical and biologic features

predictive of survival after relapse of neuroblastoma: a report from the International Neuroblastoma Risk Group project. _J. Clin. Oncol._ 29, 3286–3292 (2011). Article  PubMed  PubMed

Central  Google Scholar  * London, W. B. et al. Historical time to disease progression and progression-free survival in patients with recurrent/refractory neuroblastoma treated in the modern

era on Children’s Oncology Group early-phase trials. _Cancer_ 123, 4914–4923 (2017). Article  CAS  PubMed  Google Scholar  * Zage, P. E. Novel therapies for relapsed and refractory

neuroblastoma. _Children (Basel)_ 5, 148 (2018). * Padovan-Merhar, O. M. et al. Enrichment of targetable mutations in the relapsed neuroblastoma genome. _PLoS Genet_ 12, e1006501 (2016).

Article  PubMed  PubMed Central  Google Scholar  * Turajlic, S., Sottoriva, A., Graham, T. & Swanton, C. Resolving genetic heterogeneity in cancer. _Nat. Rev. Genet_ 20, 404–416 (2019).

Article  CAS  PubMed  Google Scholar  * Schramm, A. et al. Mutational dynamics between primary and relapse neuroblastomas. _Nat. Genet_ 47, 872–877 (2015). Article  CAS  PubMed  Google

Scholar  * Eleveld, T. F. et al. Relapsed neuroblastomas show frequent RAS-MAPK pathway mutations. _Nat. Genet_ 47, 864–871 (2015). Article  CAS  PubMed  PubMed Central  Google Scholar  *

Schleiermacher, G. et al. Emergence of new ALK mutations at relapse of neuroblastoma. _J. Clin. Oncol._ 32, 2727–2734 (2014). Article  CAS  PubMed  Google Scholar  * Andersson, N. et al.

Extensive clonal branching shapes the evolutionary history of high-risk pediatric cancers. _Cancer Res_. 80, 1512–1523 (2020). Article  CAS  PubMed  Google Scholar  * Schmelz, K. et al.

Spatial and temporal intratumour heterogeneity has potential consequences for single biopsy-based neuroblastoma treatment decisions. _Nat. Commun._ 12, 6804 (2021). Article  ADS  CAS  PubMed

  PubMed Central  Google Scholar  * Gundem, G. et al. Clonal evolution during metastatic spread in high-risk neuroblastoma. _Bioarxiv_ 2022.08.15 (2022). * Bell, E. et al. The role of MYCN

in the failure of MYCN amplified neuroblastoma cell lines to G1 arrest after DNA damage. _Cell Cycle_ 5, 2639–2647 (2006). Article  CAS  PubMed  Google Scholar  * Bown, N. et al. Gain of

chromosome arm 17q and adverse outcome in patients with neuroblastoma. _N. Engl. J. Med_. 340, 1954–1961 (1999). Article  CAS  PubMed  Google Scholar  * Brodeur, G. M., Seeger, R. C.,

Schwab, M., Varmus, H. E. & Bishop, J. M. Amplification of N-myc in untreated human neuroblastomas correlates with advanced disease stage. _Science_ 224, 1121–1124 (1984). Article  ADS 

CAS  PubMed  Google Scholar  * Caren, H. et al. High-risk neuroblastoma tumors with 11q-deletion display a poor prognostic, chromosome instability phenotype with later onset. _Proc. Natl.

Acad. Sci. USA_ 107, 4323–4328 (2010). Article  ADS  CAS  PubMed  PubMed Central  Google Scholar  * Caron, H. et al. Allelic loss of chromosome 1p as a predictor of unfavorable outcome in

patients with neuroblastoma. _N. Engl. J. Med_. 334, 225–230 (1996). Article  CAS  PubMed  Google Scholar  * Fransson, S. et al. Estimation of copy number aberrations: Comparison of exome

sequencing data with SNP microarrays identifies homozygous deletions of 19q13.2 and CIC in neuroblastoma. _Int J. Oncol._ 48, 1103–1116 (2016). Article  CAS  PubMed  Google Scholar  *

Molenaar, J. J. et al. Sequencing of neuroblastoma identifies chromothripsis and defects in neuritogenesis genes. _Nature_ 483, 589–593 (2012). Article  ADS  CAS  PubMed  Google Scholar  *

Schleiermacher, G. et al. Segmental chromosomal alterations have prognostic impact in neuroblastoma: a report from the INRG project. _Br. J. Cancer_ 107, 1418–1422 (2012). Article  CAS 

PubMed  PubMed Central  Google Scholar  * Andersson, N., Chattopadhyay, S., Valind, A., Karlsson, J. & Gisselsson, D. DEVOLUTION-A method for phylogenetic reconstruction of aneuploid

cancers based on multiregional genotyping data. _Commun. Biol._ 4, 1103 (2021). Article  CAS  PubMed  PubMed Central  Google Scholar  * Martinez-Monleon, A. et al. Amplification of CDK4 and

MDM2: a detailed study of a high-risk neuroblastoma subgroup. _Sci. Rep._ 12, 12420 (2022). Article  ADS  CAS  PubMed  PubMed Central  Google Scholar  * Holzel, M. et al. NF1 is a tumor

suppressor in neuroblastoma that determines retinoic acid response and disease outcome. _Cell_ 142, 218–229 (2010). Article  PubMed  PubMed Central  Google Scholar  * Braekeveldt, N. et al.

Patient-Derived Xenograft Models Reveal Intratumor Heterogeneity and Temporal Stability in Neuroblastoma. _Cancer Res_ 78, 5958–5969 (2018). Article  CAS  PubMed  Google Scholar  *

Braekeveldt, N. et al. Neuroblastoma patient-derived orthotopic xenografts retain metastatic patterns and geno- and phenotypes of patient tumours. _Int J. Cancer_ 136, E252–E261 (2015).

Article  CAS  PubMed  Google Scholar  * Manas, A. et al. Clinically relevant treatment of PDX models reveals patterns of neuroblastoma chemoresistance. _Sci. Adv._ 8, eabq4617 (2022).

Article  ADS  CAS  PubMed  PubMed Central  Google Scholar  * Anwar, M. et al. TCF 4 tumor suppressor: a molecular target in the prognosis of sporadic colorectal cancer in humans. _Cell Mol.

Biol. Lett._ 25, 24 (2020). Article  CAS  PubMed  PubMed Central  Google Scholar  * Veeriah, S. et al. The tyrosine phosphatase PTPRD is a tumor suppressor that is frequently inactivated and

mutated in glioblastoma and other human cancers. _Proc. Natl Acad. Sci. USA_ 106, 9435–9440 (2009). Article  ADS  CAS  PubMed  PubMed Central  Google Scholar  * Biedler, J. L., Helson, L.

& Spengler, B. A. Morphology and growth, tumorigenicity, and cytogenetics of human neuroblastoma cells in continuous culture. _Cancer Res_ 33, 2643–2652 (1973). CAS  PubMed  Google

Scholar  * Thirant, C. et al. Reversible transitions between noradrenergic and mesenchymal tumor identities define cell plasticity in neuroblastoma. _Nat. Commun._ 14, 2575 (2023). Article 

ADS  CAS  PubMed  PubMed Central  Google Scholar  * Otte, J., Dyberg, C., Pepich, A. & Johnsen, J. I. MYCN function in neuroblastoma development. _Front Oncol._ 10, 624079 (2020).

Article  PubMed  Google Scholar  * Karlsson, J. et al. Four evolutionary trajectories underlie genetic intratumoral variation in childhood cancer. _Nat. Genet_ 50, 944–950 (2018). Article 

CAS  PubMed  Google Scholar  * Prieto-Vila, M. et al. Single-cell analysis reveals a preexisting drug-resistant subpopulation in the luminal breast cancer subtype. _Cancer Res_ 79, 4412–4425

(2019). Article  CAS  PubMed  Google Scholar  * Seth, S. et al. Pre-existing functional heterogeneity of tumorigenic compartment as the origin of chemoresistance in pancreatic tumors. _Cell

Rep._ 26, 1518–1532.e9 (2019). Article  CAS  PubMed  Google Scholar  * Brady, S. W. et al. Combating subclonal evolution of resistant cancer phenotypes. _Nat. Commun._ 8, 1231 (2017).

Article  ADS  PubMed  PubMed Central  Google Scholar  * Shlush, L. I. et al. Tracing the origins of relapse in acute myeloid leukaemia to stem cells. _Nature_ 547, 104–108 (2017). Article 

ADS  CAS  PubMed  Google Scholar  * Morgillo, F., Della Corte, C. M., Fasano, M. & Ciardiello, F. Mechanisms of resistance to EGFR-targeted drugs: lung cancer. _ESMO Open_ 1, e000060

(2016). Article  PubMed  PubMed Central  Google Scholar  * Hata, A. N. et al. Tumor cells can follow distinct evolutionary paths to become resistant to epidermal growth factor receptor

inhibition. _Nat. Med_. 22, 262–269 (2016). Article  CAS  PubMed  PubMed Central  Google Scholar  * Bhang, H. E. et al. Studying clonal dynamics in response to cancer therapy using

high-complexity barcoding. _Nat. Med_ 21, 440–448 (2015). Article  CAS  PubMed  Google Scholar  * Morelli, M. P. et al. Characterizing the patterns of clonal selection in circulating tumor

DNA from patients with colorectal cancer refractory to anti-EGFR treatment. _Ann. Oncol._ 26, 731–736 (2015). Article  CAS  PubMed  PubMed Central  Google Scholar  * Kim, J. et al.

Preexisting oncogenic events impact trastuzumab sensitivity in ERBB2-amplified gastroesophageal adenocarcinoma. _J. Clin. Invest_ 124, 5145–5158 (2014). Article  PubMed  PubMed Central 

Google Scholar  * Paiva, B. et al. Phenotypic and genomic analysis of multiple myeloma minimal residual disease tumor cells: a new model to understand chemoresistance. _Blood_ 127, 1896–1906

(2016). Article  CAS  PubMed  Google Scholar  * van Groningen, T. et al. A NOTCH feed-forward loop drives reprogramming from adrenergic to mesenchymal state in neuroblastoma. _Nat. Commun._

10, 1530 (2019). Article  ADS  PubMed  PubMed Central  Google Scholar  * van Groningen, T. et al. Neuroblastoma is composed of two super-enhancer-associated differentiation states. _Nat.

Genet_ 49, 1261–1266 (2017). Article  PubMed  Google Scholar  * Boeva, V. et al. Heterogeneity of neuroblastoma cell identity defined by transcriptional circuitries. _Nat. Genet_ 49,

1408–1413 (2017). Article  CAS  PubMed  Google Scholar  * Gartlgruber, M. et al. Super enhancers define regulatory subtypes and cell identity in neuroblastoma. _Nat. Cancer_ 2, 114–128

(2021). Article  CAS  PubMed  Google Scholar  * Van Roy, N. et al. Shallow whole genome sequencing on circulating cell-free DNA allows reliable noninvasive copy-number profiling in

neuroblastoma patients. _Clin. Cancer Res_ 23, 6305–6314 (2017). Article  PubMed  Google Scholar  * Chicard, M. et al. Whole-exome sequencing of cell-free DNA reveals temporo-spatial

heterogeneity and identifies treatment-resistant clones in neuroblastoma. _Clin. Cancer Res_ 24, 939–949 (2018). Article  CAS  PubMed  Google Scholar  * Persson, C. U. et al. Neuroblastoma

patient-derived xenograft cells cultured in stem-cell promoting medium retain tumorigenic and metastatic capacities but differentiate in serum. _Sci. Rep._ 7, 10274 (2017). Article  ADS 

PubMed  PubMed Central  Google Scholar  * Kimura, M. Theoretical foundation of population genetics at the molecular level. _Theor. Popul Biol._ 2, 174–208 (1971). Article  CAS  PubMed 

Google Scholar  * Heim, S. & Mitelman, F. _Cancer cytogenetics: chromosomal and molecular genetic aberrations of tumor cells_, ix, 632 pages (Wiley Blackwell, Chichester, West Sussex;

Hoboken, NJ, 2015). * Schliep, K. P. phangorn: phylogenetic analysis in R. _Bioinformatics_ 27, 592–593 (2011). Article  CAS  Google Scholar  * Mayrhofer, M., Viklund, B. & Isaksson, A.

Rawcopy: Improved copy number analysis with Affymetrix arrays. _Sci. Rep._ 6, 36158 (2016). Article  ADS  CAS  PubMed  PubMed Central  Google Scholar  * Rasmussen, M. et al. Allele-specific

copy number analysis of tumor samples with aneuploidy and tumor heterogeneity. _Genome Biol._ 12, R108 (2011). Article  CAS  PubMed  PubMed Central  Google Scholar  * MacDonald, J. R.,

Ziman, R., Yuen, R. K., Feuk, L. & Scherer, S. W. The Database of Genomic Variants: a curated collection of structural variation in the human genome. _Nucleic Acids Res_ 42, D986–D992

(2014). Article  CAS  PubMed  Google Scholar  * Yasui, H. et al. A dynamic mutational landscape associated to an inter-regionally diverse immune response in malignant rhabdoid tumour. _J.

Pathol._ 252, 22–28 (2020). Article  CAS  PubMed  Google Scholar  * Li, H. Aligning sequence reads, clone sequences and assembly contigs with BWA-MEM. In _arXiv_ v. 2 edn (arXiv.org, Broad

Institute, 2013). * Tarasov, A., Vilella, A. J., Cuppen, E., Nijman, I. J. & Prins, P. Sambamba: fast processing of NGS alignment formats. _Bioinformatics_ 31, 2032–2034 (2015). Article

  CAS  PubMed  PubMed Central  Google Scholar  * Kim, S. et al. Strelka2: fast and accurate calling of germline and somatic variants. _Nat. Methods_ 15, 591–594 (2018). Article  CAS  PubMed

  Google Scholar  * Robinson, J. T. et al. Integrative genomics viewer. _Nat. Biotechnol._ 29, 24–26 (2011). Article  CAS  PubMed Central  Google Scholar  * Adzhubei, I. A. et al. A method

and server for predicting damaging missense mutations. _Nat. Methods_ 7, 248–249 (2010). Article  CAS  PubMed  PubMed Central  Google Scholar  * McKenna, A. et al. The Genome Analysis

Toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data. _Genome Res_ 20, 1297–1303 (2010). Article  CAS  Google Scholar  * Cibulskis, K. et al. Sensitive detection

of somatic point mutations in impure and heterogeneous cancer samples. _Nat. Biotechnol._ 31, 213–219 (2013). Article  CAS  PubMed  PubMed Central  Google Scholar  * Fan, Y. et al. MuSE:

Accounting for tumor heterogeneity using a sample-specific error model improves sensitivity and specificity in mutation calling from sequencing data. _Genome Biol._ 17, 178 (2016). Article 

PubMed  PubMed Central  Google Scholar  * Van der Auwera, G. A. et al. From FastQ data to high confidence variant calls: The Genome Analysis Toolkit best practices pipeline. _Curr. Protoc.

Bioinforma._ 43, 11 10 1–11 10 33 (2013). Google Scholar  * Yasui, H. Microenvironmental impact on tumour cell phenotype and genotype in adult and paediatric tumours. doctoral thesis, lund

university publications. _Doctoral Dissertation Ser._ 2021, 3 (2021). Google Scholar  * Satas, G., Zaccaria, S., El-Kebir, M. & Raphael, B. J. DeCiFering the elusive cancer cell fraction

in tumor heterogeneity and evolution. _Cell Syst._ 12, 1004–1018 e10 (2021). Article  CAS  PubMed  PubMed Central  Google Scholar  * Mengelbier, L. H. et al. Intratumoral genome diversity

parallels progression and predicts outcome in pediatric cancer. _Nat. Commun._ 6, 6125 (2015). Article  ADS  CAS  PubMed  Google Scholar  * Miller, C. A. et al. Visualizing tumor evolution

with the fishplot package for R. _BMC Genomics_ 17, 880 (2016). Article  PubMed  PubMed Central  Google Scholar  * Karlsson, J. et al. Early evolutionary branching across spatial domains

predisposes to clonal replacement under chemotherapy in neuroblastoma. _Github_ https://doi.org/10.5281/zenodo.13304334 (2024). Download references ACKNOWLEDGEMENTS The present study was

supported by grants from the Swedish Research Council, the Swedish Cancer Society (21 1383 Pj), the Swedish Childhood Cancer Foundation (PR2022-0026), Regional Clinical Research grants (ALF

Projekt0067), the Royal Physiographic Society, and the LMK Foundation. We thank Center for Translational Genomics, Lund University and Clinical Genomics Lund, SciLifeLab for providing

sequencing service. The illustration of a mouse used in Figs. 3a and e and 4a is a slight adaptation of “Vector diagram of a laboratory mouse (black and white)” by Gwilz/Wikimedia Commons/CC

BY-SA 4.0. FUNDING Open access funding provided by Lund University. AUTHOR INFORMATION AUTHORS AND AFFILIATIONS * Division of Clinical Genetics, Department of Laboratory Medicine, Lund

University, Lund, Sweden Jenny Karlsson, Hiroaki Yasui, Natalie Andersson, Caroline Jansson, Geoffroy Durand, Naveen Ravi, Michele Ferro, Minjun Yang, Subhayan Chattopadhyay, Kajsa Paulsson,

 Anders Valind & David Gisselsson * Department of Obstetrics and Gynecology, Nagoya University Graduate School of Medicine, Nagoya, Japan Hiroaki Yasui * Division of Translational Cancer

Research, Department of Laboratory Medicine, Lund University, Lund, Sweden Adriana Mañas, Karin Hansson, Kristina Aaltonen & Daniel Bexell * European Research Institute for the Biology

of Ageing (ERIBA), University of Groningen, University Medical Center Groningen, Groningen, The Netherlands Diana Spierings & Floris Foijer * Department of Pediatrics, Skåne University

Hospital, Lund, Sweden Anders Valind * Department of Pathology, Office of Medical Services, Region Skåne, Lund, Sweden David Gisselsson Authors * Jenny Karlsson View author publications You

can also search for this author inPubMed Google Scholar * Hiroaki Yasui View author publications You can also search for this author inPubMed Google Scholar * Adriana Mañas View author

publications You can also search for this author inPubMed Google Scholar * Natalie Andersson View author publications You can also search for this author inPubMed Google Scholar * Karin

Hansson View author publications You can also search for this author inPubMed Google Scholar * Kristina Aaltonen View author publications You can also search for this author inPubMed Google

Scholar * Caroline Jansson View author publications You can also search for this author inPubMed Google Scholar * Geoffroy Durand View author publications You can also search for this author

inPubMed Google Scholar * Naveen Ravi View author publications You can also search for this author inPubMed Google Scholar * Michele Ferro View author publications You can also search for

this author inPubMed Google Scholar * Minjun Yang View author publications You can also search for this author inPubMed Google Scholar * Subhayan Chattopadhyay View author publications You

can also search for this author inPubMed Google Scholar * Kajsa Paulsson View author publications You can also search for this author inPubMed Google Scholar * Diana Spierings View author

publications You can also search for this author inPubMed Google Scholar * Floris Foijer View author publications You can also search for this author inPubMed Google Scholar * Anders Valind

View author publications You can also search for this author inPubMed Google Scholar * Daniel Bexell View author publications You can also search for this author inPubMed Google Scholar *

David Gisselsson View author publications You can also search for this author inPubMed Google Scholar CONTRIBUTIONS J.K., N.A., S.C., K.P, F.F., D.B., and D.G., designed the study. H.Y.,

A.M., K.H., K.A., C.J., G.D., N.R., J.K., and M.F. performed experimental work. J.K., H.Y., A.M., N.A., K.A., M.Y., N.R., S.C., D.S., F.F., A.V., and D.G. analyzed data. D.G. and J.K.

selected patients from clinical material. CORRESPONDING AUTHOR Correspondence to David Gisselsson. ETHICS DECLARATIONS COMPETING INTERESTS D.B. has received research funding from Healx,

aPODD foundation, and Captor Therapeutics for unrelated work. The other authors declare no competing interests. PEER REVIEW PEER REVIEW INFORMATION _Nature Communications_ thanks Sheila

Singh, John M. Maris, and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. A peer review file is available. ADDITIONAL INFORMATION PUBLISHER’S NOTE

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. SUPPLEMENTARY INFORMATION SUPPLEMENTARY INFORMATION PEER REVIEW FILE

DESCRIPTION OF ADDITIONAL SUPPLEMENTARY FILES SUPPLEMENTARY DATA 1 SUPPLEMENTARY DATA 2 REPORTING SUMMARY SOURCE DATA SOURCE DATA RIGHTS AND PERMISSIONS OPEN ACCESS This article is licensed

under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate

credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article

are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and

your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this

licence, visit http://creativecommons.org/licenses/by/4.0/. Reprints and permissions ABOUT THIS ARTICLE CITE THIS ARTICLE Karlsson, J., Yasui, H., Mañas, A. _et al._ Early evolutionary

branching across spatial domains predisposes to clonal replacement under chemotherapy in neuroblastoma. _Nat Commun_ 15, 8992 (2024). https://doi.org/10.1038/s41467-024-53334-x Download

citation * Received: 03 October 2022 * Accepted: 09 October 2024 * Published: 18 October 2024 * DOI: https://doi.org/10.1038/s41467-024-53334-x SHARE THIS ARTICLE Anyone you share the

following link with will be able to read this content: Get shareable link Sorry, a shareable link is not currently available for this article. Copy to clipboard Provided by the Springer

Nature SharedIt content-sharing initiative