Neuroblastoma arises in early fetal development and its evolutionary duration predicts outcome

Neuroblastoma, the most frequent solid tumor in infants, shows very diverse outcomes from spontaneous regression to fatal disease. When these different tumors originate and how they evolve are not known. Here we quantify the somatic evolution of neuroblastoma by deep whole-genome sequencing, molecular clock analysis and population-genetic modeling in a comprehensive cohort covering all subtypes. We find that tumors across the entire clinical spectrum begin to develop via aberrant mitoses as early as the first trimester of pregnancy. Neuroblastomas with favorable prognosis expand clonally after short evolution, whereas aggressive neuroblastomas show prolonged evolution during which they acquire telomere maintenance mechanisms. The initial aneuploidization events condition subsequent evolution, with aggressive neuroblastoma exhibiting early genomic instability. We find in the discovery cohort (n = 100), and validate in an independent cohort (n = 86), that the duration of evolution is an accurate predictor of outcome. Thus, insight into neuroblastoma evolution may prospectively guide treatment decisions.

Neuroblastoma, the most frequent solid tumor in infants, shows very diverse outcomes from spontaneous regression to fatal disease. When these different tumors originate and how they evolve are not known. Here we quantify the somatic evolution of neuroblastoma by deep whole-genome sequencing, molecular clock analysis and population-genetic modeling in a comprehensive cohort covering all subtypes. We find that tumors across the entire clinical spectrum begin to develop via aberrant mitoses as early as the first trimester of pregnancy. Neuroblastomas with favorable prognosis expand clonally after short evolution, whereas aggressive neuroblastomas show prolonged evolution during which they acquire telomere maintenance mechanisms. The initial aneuploidization events condition subsequent evolution, with aggressive neuroblastoma exhibiting early genomic instability. We find in the discovery cohort (n = 100), and validate in an independent cohort (n = 86), that the duration of evolution is an accurate predictor of outcome. Thus, insight into neuroblastoma evolution may prospectively guide treatment decisions.
Cancers result from the accumulation of oncogenic mutations 1 . Insights into tumor evolution-and, particularly, the temporal order of driver mutations-are beginning to support diagnosis and treatment 2 . Initially, age-incidence curves were used to estimate the number of rate-limiting mutations in carcinogenesis 3,4 . Recently, the order of driver mutations has been inferred from genome sequencing data [5][6][7][8][9] , and mathematical approaches based on population genetics have been used to reconstruct the clonal evolution of cancers from the statistics of somatic variants 2,10 . Deep whole-genome sequencing (WGS) data are suited for such integrative analyses, because the allele frequencies of neutral somatic variants that hitchhike with driver mutations provide rich information on how drivers shape tumor growth 7,[11][12][13] .
Tumors of early childhood provide a paradigm for cancer evolution in the context of development. A key question is how driver mutations subvert the normal development of the tissue of origin 14,15 . The most common solid tumor in infants, neuroblastoma, arises in the sympathetic nervous system. A striking feature of neuroblastoma is the wide spectrum of clinical outcomes, ranging from low-risk cases requiring light or no treatment to high-risk disease that remains fatal for ~50% of patients 16 . Characteristic mutations, including MYCN amplification, gain of telomere maintenance mechanisms (TMMs) and gains (17q) or losses (1p, 11q) of chromosomal segments, have been associated with high-risk disease 17 . Nevertheless, the prospective stratification of patients into observation and different treatment groups remains a formidable challenge. In this Article, we study how neuroblastomas originate in development and evolve genetically and ask whether this understanding can provide insights into disease severity and outcomes. Article https://doi.org/10.1038/s41588-023-01332-y may erroneously cause a subset of subclonal mutations to be classified as clonal, we analyzed pairs of primary and relapse samples and found that the vast majority (85 ± 5%) of clonal SSNVs were present in both samples, indicating that sampling did not introduce a strong bias. Nevertheless, to estimate MRCA density conservatively we performed all subsequent computations with the measured densities corrected by a factor of 0.85. In the example tumor, the mean SSNV density of the MRCA was one SSNV per 5 million base pairs (bp).
Next, we asked whether clonal chromosomal gains occurred coincident with, or ancestral to, the MRCA. To this end, we quantified SSNVs that were identical and clonal on two copies on trisomic and tetrasomic segments (termed amplified clonal SSNVs 6 ; Fig. 2b). These mutations were acquired before a gain on the respective allele (Fig. 2c, dark green), and hence the number of these mutations correlates positively with the time at which the chromosomal gain occurred 6,29 . To make this timing measure independent of segment length, we use SSNV density. We compared these densities to those at the MRCA based on a negative binomial distribution of SSNVs across the genome (Methods). In the example tumor, nearly all gains had a mean density of one amplified clonal SSNV per 100 Mbp, which is significantly smaller (adjusted P < 0.01) than that of the MRCA (Fig. 2e). Hence, the molecular clock places these gains ancestral to the MRCA in an early common ancestor (ECA) of the tumor. The only exceptions were gains of Chr. 9 and 20q, which had a mutation density consistent with that of the MRCA and hence occurred later than the ECA. On those segments gained ancestral to the MRCA, the densities of amplified clonal SSNVs were statistically indistinguishable (Fig. 2e). Hence, near-triploidization occurred early and as a temporally confined event during the development of this tumor. Thereafter, further genetic evolution occurred before the clonal sweep from the MRCA commenced (Fig. 2f).

Two evolutionary classes of neuroblastoma
We timed the MRCAs for all tumor samples, and the ECAs associated with chromosomal gains. First, we validated these timing approaches by comparison of mutation densities in samples taken following primary diagnosis with those of relapsed neuroblastomas. As expected, timing of the ECA, defining the putative early origin of tumorigenesis, was conserved in all sample types (Extended Data Fig. 2b). By contrast, the timing of the MRCA, defining the origin of a tumor sample, was significantly later in relapsed tumors, consistent with a bottleneck imposed by incomplete tumor resection or cytotoxic therapy (Supplementary  Table 1) and, eventually, regrowth from a small number of surviving cells (Extended Data Fig. 2b,c). Interestingly, the MRCAs of metastases resected after initial diagnosis had SSNV densities indistinguishable from those of the primary tumors (Extended Data Fig. 2b,c), suggesting that metastases had originated around the time when the tumor started to grow from its MRCA.
To infer the early evolution of tumors up to the MRCA, we focused on samples taken at initial diagnosis (primary tumors and metastases, n = 67 in the discovery cohort). Remarkably, mutation densities of the MRCA showed a bimodal distribution (Fig. 3a). This finding suggests that there are two classes of neuroblastoma, one where growth of the resected tumor commenced early (early-MRCA neuroblastoma) and another where it occurred later (late-MRCA neuroblastoma); for mutational signatures, see Extended Data Fig. 2d,e. Of note, although we estimated the mutation density of the MRCA conservatively (Methods), the dichotomy between early and late MRCA emerged robustly.
These data raise the question of whether late-MRCA tumors began to develop later or developed early and evolved for a longer period of time. To address this question in terms of the molecular clock, we analyzed SSNV densities on chromosomal/segmental gains in both tumor classes that defined an ECA. All 20 early-MRCA tumors had such gains (with 70% showing near-triploidization, 35% harboring segmental gains and 25% displaying both features; Fig. 3b). The respective mutation densities timing chromosomal gains and MRCA were statistically

Mutation patterns in a comprehensive neuroblastoma cohort
We assembled a discovery cohort of deep (~80×) WGS data of 100 neuroblastomas (Supplementary Table 1), covering all clinical stages of the disease according to the International Neuroblastoma Staging System (Fig. 1a). Sixty-seven samples were derived from initial diagnoses (including seven from metastases) and 33 from relapsed tumors. Median age at primary diagnosis was 0.7 years for stages 1, 2 and 4S, 3.5 years for stage 3 and 3.9 years for stage 4 (Extended Data Fig. 1a).
Median tumor purity was high (88%), which allowed for reliable estimation of tumor ploidy (Fig. 1b) as confirmed by direct measurement of DNA index (Extended Data Fig. 1b). Tumors had baseline ploidies between two and four, with 55, 33 and 12 samples being near-diploid, near-triploid and near-tetraploid, respectively. Each ploidy class contained tumors of all stages (Extended Data Fig. 1c). Relative to baseline ploidy, we detected in 96% of the tumors characteristic segmental (1q and 17q) or whole-chromosome (2, 7 and 17) gains, as well as segmental (1p, 11q) or whole-chromosome (11) losses (Extended Data Fig. 1d), all of which are probable drivers of neuroblastoma [18][19][20][21][22][23] . Based on tumor cell content and ploidy, the vast majority of these gains and losses were clonal and hence were present in the most recent common ancestor cell (MRCA) of the tumor (Fig. 1c). These data point to aneuploidy as an early feature of neuroblastoma.
In contrast to driver mutations, neutral SSNVs and small indels are continuously accumulated and contain information on tumor evolution. Analysis of mutational signatures assigned the majority of SSNVs to clock-like signatures (SBS1, SBS5, SBS40). Overall the next most abundant signature was SBS18, associated with reactive oxygen species 27 and frequently found in neuroblastoma 28 , followed by SBS3, associated with failure of homologous recombination in dividing cells (Fig. 1e). These signatures overall occurred at similar frequencies clonally and in subclones (Extended Data Fig. 1e).

Neutral SSNVs distinguish sequential driver events
To understand when and how neuroblastomas evolve, we used neutral SSNVs as a molecular clock to time key events: (1) the emergence of the MRCA, from which the clonal sweep that defines the resected tumor emerges, and (2) the acquisition of clonal chromosomal gains (which may have occurred before or coincident with the MRCA). For each tumor, we computed the frequency distribution of somatic variants on all genomic segments of a given copy number. An example tumor, classified as near-triploid (Fig. 2a), had dominant copy numbers 2 and 3 and a smaller portion of segments at copy number 4-a typical configuration in neuroblastoma. On each copy number, we found a large clonal peak (Fig. 2b) comprising mutations occurring on a single chromosomal copy in each tumor cell. In addition, on copy numbers >2, we detected clonal mutations present on the two (and in some cases more) copies of a multiplied chromosome. Collectively, both single-and multiple-copy clonal SSNVs characterize the MRCA (Fig. 2c). Indeed, the density of clonal SSNVs was similar for the different copy numbers (Extended Data Fig. 2a), thus timing the MRCA by means of a molecular clock (Fig. 2d). To evaluate a potentially biasing effect of partial tumor sampling, which  Fig. 3a shows an example), suggesting that aneuploidy was not followed by acquisition of further clonal drivers in these tumors. In the remaining 10% of tumors, separate ECAs and MRCAs were timed but these were very close in regard to mutation density (0.04 SSNVs per Mb between ECA and MRCA; Extended Data Fig. 3b;), implying that they were temporally close events. Indeed, 95% of early-MRCA tumors lacked additional small-scale drivers while a single tumor harbored a MYCN amplification and a mutation in ALK. Hence early-MRCA tumors appear to be driven predominantly by aneuploidization. By contrast, in the majority (55%) of the 47 late-MRCA tumors, we distinguished two well-separated events, ECA and MRCA (Fig. 3c, left), with average distance 0.24 ± 0.1 SSNVs per Mb. In about one-third of these cases, a single early near-triplodization event defined an ECA (Fig. 2); whereas, in the remaining two-thirds, smaller-scale chromosomal gains defined the ECA (Extended Data Fig. 3c). Hence, the majority of late-MRCA tumors showed a signature of early gains followed by further genetic evolution to the MRCA. In the remaining 45% of late-MRCA cases, we could not reliably time a separate ECA (for example, if very short fragments were gained or if gains were found at copy number >4) (Fig. 3c (right) and Extended Data Fig. 3d). Hence, for these tumors we cannot time an early genetic event which, however, leaves open the possibility that small-scale mutations, chromosomal losses or high-level amplifications (four or more copies), neither of which can be timed reliably, preceded the MRCA of the tumor. Consistent with this hypothesis, the mutation densities of the MRCA in tumors without timeable ECA were nearly as high as in those with two distinguishable events (Extended Data Fig. 3e). Moreover, near-tetraploid tumors without timeable ECA showed evidence of sequential events. Here, we found a 2:0 allelic configuration in 11 of 12 near-tetraploid tumors with 1p loss, suggesting that 1p deletion preceded genome doubling. Overall, the ECA of late-MRCA tumors had mutation density indistinguishable from the MRCAs of early-MRCA tumors (Fig. 3d), suggesting that the majority of late-MRCA tumors acquired aneuploidy as early as the early-MRCA tumors.
Finally, we analyzed the timing of gains in genome-wide events: near-triploidization of the genome or genome doubling. In all such tumors, the individual gains involved showed statistically indistinguishable timing, which is consistent with near-triploidization and genome doubling occurring as single catastrophic events (Extended Data Fig. 3f).
In sum, we find that the MRCAs in primary neuroblastoma fell into two evolutionary groups: in early-MRCA tumors, near-triploidization of the genome and/or gains of chromosomal arms or whole chromosomes coincided with the MRCA. Similarly, in more than one-half of late-MRCA tumors, an early ECA was defined by whole-chromosomal  Horizontal lines represent mutation densities at gained segments, as in e, with 95% confidence bounds computed from χ 2 distributions. Segments on which densities of amplified clonal mutations were significantly different from that at MRCA are marked by ** (adjusted P < 0.01); NS, not significant; test and P values as in e. e,f, Color coding for segments is in a.

Long evolution of neuroblastoma predicts unfavorable outcome
The early-MRCA class contained cases of stages 1, 2 and 4 S, typically having a prognosis, but only one case of stage 4 (Fig. 3a). Hence, we asked whether MRCA timing predicts outcome. Remarkably, early MRCA timing clearly identified cases with long event-free survival (Fig. 4a) and long overall survival (Fig. 4b).
To further test this idea, we assembled an independent cohort of primary tumors and metastases, which was enriched for tumors with WGS at the customary coverage of 30× (n = 86; Supplementary Tables 6-10) 22,25 . We evaluated the timing of ECA and MRCA in each sample as described for the discovery cohort. The mutation density of the MRCA showed the same bimodal pattern as in the discovery cohort (Fig. 4c) Fig. 3c, right). Thus, the validation cohort corroborates the scenarios of ECA and MRCA timing found in the discovery cohort.
In the validation cohort, MRCA timing was an accurate predictor of both event-free and overall survival (Fig. 5a,b). Merging the two cohorts (n = 152 primary tumors and metastases, excluding one case lacking survival data), we confirmed this result (Fig. 5c,d). To compare MRCA timing with other predictors of survival, we considered clinical variables used worldwide (stage, age), gain of a TMM (which improves on the clinically used criterion, MYCN amplification 26 ) and a more recently proposed molecular predictor, the mutation status of the RAS/p53 pathway 25 . Early MRCA timing emerged as the most informative predictor of event-free survival (Fig. 5); overall survival was best explained by both MRCA timing and mutations in the RAS/p53 pathway (Extended Data Fig. 4). In sum, survival analyses suggest that extended evolution up to the founding cell of the primary tumor predicts unfavorable outcome.

Genomic instability and telomere maintenance in late-MRCA tumors
Given that both early-and late-MRCA neuroblastomas begin to develop in the same time window (compare with Fig. 3d) but show markedly different durations of evolution and clinical outcome, we asked whether late-MRCA tumors have evolved further than early-MRCA examples simply by chance, or whether there are molecular factors that     Article https://doi.org/10.1038/s41588-023-01332-y (including gains of 17q, 7q and 1q, as well as loss of 11q and 1p) was prevalent in the late-MRCA group. In late-MRCA tumors with timeable ECA, the vast majority of segmental gains that we could time were coincident with the ECA (Figs. 3c and 4e, dark green squares). Taken together, different processes underlie the early acquisition of aneuploidy in the two classes: mis-segregation of entire chromosomes in early-MRCA tumors vis-à-vis genomic instability in late-MRCA tumors (Fig. 6a). Acquired TMMs were strongly enriched in the late-MRCA class (Fig. 6a); 12% of early-MRCA tumors but 81% of late-MRCA tumors gained telomere maintenance-via MYCN amplification, ALT or TERT rearrangement-in the combined discovery and validation cohorts (Fig. 6b). Due to their small size, these structural rearrangements cannot be timed using mutation densities and hence may have been acquired before the emergence of arm-level aneuploidy, together with aneuploidy or subsequently. Early timing of arm-level aneuploidy in more than one-half of the late-MRCA cases (that is, tumors with an ECA; left-hand panels in Figs. 3c and 4e) suggests that telomere maintenance was gained during a secondary event between ECA and MRCA in these cases, which mainly included ALT and TERT rearrangement. In the remaining cases, where arm-level aneuploidy was timed at the MRCA and no ECA was identifiable (right-hand panels in Figs. 3c and 4e), telomere maintenance could also have been acquired in an oncogenic event preceding the MRCA. Interestingly, those cases with MYCN amplification fell predominantly within this latter group (Fisher's exact test, P = 0.006055, odds ratio = 3.9 (1.4, 11.6)), suggesting that MYCN amplification tends to occur earlier than ALT or TERT rearrangement. Indeed, MYCN-amplified tumors generally had an earlier MRCA than tumors with ALT or TERT rearrangement (Fig. 6c).
Collectively, these findings establish a link between genetic evolution of neuroblastoma and the observation that neuroblastomas with extensive arm-level aneuploidy tend to carry a poor prognosis 30 . The typically late-emerging MRCA in these cases indicates that the underlying genomic instability predisposes these tumors to prolonged evolution, including the acquisition of TMMs.

Chromosomal gains occur during sympathetic neurogenesis
Finally, we asked whether our genetic insights into neuroblastoma development could be synthesized into an integrative model of tumor evolution. To this end, we devised population-genetic models and quantified key parameters by fitting the models to our measured data, including the time-dependent incidences of ECA and MRCA and the variant allele frequency (VAF) distribution of somatic variants. In addition, we required the models to reproduce the overall incidence of the disease in the human population. For this reason, we focused on neuroblastomas with poor prognosis ('high-risk', enriched in the late-MRCA class), which occur in around one in 10 5 children (the frequency of low-risk cases, enriched in the early-MRCA class, is not reliably known due to incomplete diagnosis [31][32][33][34]. The models describe a proliferative population of putative cells of origin that are lost by either differentiation or cell death (Fig. 7). On average, µ SSNVs occur per cell division. A rare subset of mutations (including chromosomal gains/losses and smaller-scale events, such as SSNVs or localized MYCN amplification) will be oncogenic drivers and give rise to selected cell clones. As a minimal requirement for the evolution of a high-risk tumor, we accounted for two oncogenic events, defining ECA and MRCA, with mutation frequencies µ 1 and µ 2 per cell division. To determine these parameters, the model was fit to the experimental data using approximate Bayesian computation. Initially, we assumed that cells of origin are generated during fetal development and then remain available for neuroblastoma development. This simple model consistently overestimated the overall incidence of high-risk disease (Fig. 8a). By contrast, a model in which putative cells of origin were available for only a limited time window (Fig. 8b) accurately reproduced both the overall incidence and dynamics of the emergence of ECA and MRCA (Fig. 8c-e and Extended Data Fig. 5a).
The inferred rate of SSNV acquisition (µ) was 3.2 ± 0.4 SSNVs per day, consistent with measurements of somatic mutation rate in the developing central nervous system (5.1 (1.5, 9.0; 95% confidence interval) SSNVs per day) 35 . The model further inferred that oncogenic driver events occurred on average once per 1 million cell divisions (geometric mean of µ 1 and µ 2 ), which is consistent with a global estimate of driver mutation rate of 3.4 × 10 -5 per division in the human genome 13 , given that only a subset of all drivers will cause neuroblastoma. We used SSNV rate to calibrate the molecular clock against real time, which placed the ECA within the first trimester of pregnancy (Fig. 8c), at which time rapidly dividing sympathetic neuroblasts are developing 36 . This first hit (ECA) sustains a subclone in which the MRCA emerges due to a further oncogenic event (Extended Data Fig. 5b) No. of SSNVs (VAF) = ƒ(µ e , no. of SSNVs MRCA )   To obtain an estimate for the population level, these estimates are subsequently averaged across the cohort.
Article https://doi.org/10.1038/s41588-023-01332-y occurred within the first year of life (Fig. 8d). The model fit of Fig. 8 Fig. 5c,d). Hence, confining the analysis to clock-like mutations corroborates the real-time calibration of ECA and MRCA.
A further insight afforded by the population-genetic model is the extent of cell loss in the growing tumor. In general, only a subset of cell lineages will support growth by continued symmetric self-renewal of malignant cells whereas other lineages will terminate by either cell differentiation into nonproliferating states or cell death. We inferred the ratio of self-amplifying tumor cell divisions among all divisions from the subclonal tail of VAF distribution (Methods and ref. 37 ), finding that only ~10% of tumor cell divisions result in growth of late-MRCA neuroblastomas that acquired TMMs (Fig. 8f,g). This inference is consistent with the clinical observation of extensive cell death in such  Article https://doi.org/10.1038/s41588-023-01332-y neuroblastomas. Moreover, the average cell division rate was lower in neuroblastomas with ALT than in those with MYCN amplification or TERT rearrangement (Fig. 8h), again in line with clinical observation; our estimated cell division rates agree quantitatively with those measured in neuroblastoma in vivo 38 . The fraction of self-renewing cell divisions supporting tumor growth was inferred from VAF distribution also for tumors without acquired TMMs, all falling within the early-MRCA class, yielding a higher value of ~30% (Fig. 8g). Hence, cell loss appears to be a stronger selective pressure for late-MRCA neuroblastomas, the majority of which gain TMMs, than for early-MRCA neuroblastomas, which rarely acquire such mechanisms.

Discussion
In this paper we timed genetic events in the evolution of neuroblastoma using the molecular clock of SSNV accumulation and, inferring the rate of SSNV acquisition from the distribution of VAFs, related this clock to real time by factoring in the age at diagnosis. In two-thirds of cases we find that chromosomal gains implicated in the pathogenesis of the disease occurred early, and typically within the first trimester of pregnancy. With respect to further evolution, these cases fall into two distinct classes: in the early-MRCA class the early chromosomal gain event also marked clonal outgrowth of the resected tumor whereas in the late-MRCA class the tumors evolved further before clonal outgrowth commenced. Remarkably, in our cohort MRCA class is an accurate predictor of clinical outcome. This is true regardless of whether we could time an early chromosomal gain, and implies that neuroblastomas with a longer evolutionary history are more aggressive. Because the strong association between MRCA timing and outcome was also present with 30× WGS, the utility of this predictor for patient stratification may be tested in clinical trials.
Our real-time inference shows that neuroblastomas across the entire clinical spectrum acquired aneuploidy within the first trimester of pregnancy, when the adrenal medulla forms from sympathetic neuroblasts. Moreover, matching disease incidence with the population-genetic model suggests that the initial oncogenic event is limited to this time window. The transcriptomes of neuroblastomas most resemble those of sympathetic neuroblasts 15,36 . In this early window, neuroblasts are highly proliferative, which may make them vulnerable to aneuploidy. Finally, the observation that aneuploidy via near-triploidization is a temporally confined event is consistent with the long-standing hypothesis that this karyotype results from endoreduplication of the genome followed by tripolar cell division and selection of the fittest daughter cell 39,40 .
The molecular nature of early aneuploidy is associated with whether the tumor continues to evolve: neuroblastomas with whole-chromosome aneuploidy typically did not evolve further and were overall associated with favorable outcomes; in contrast, most tumors with early genomic instability had unfavorable outcomes. Continued evolution of such tumors has also been noted in a recent study taking multiregion biopsies 41 , emphasizing the potential of spatially resolved genetic and transcriptomic analysis 42,43 . We did not detect specific drivers of genomic instability: in particular, p53 function was not impaired genetically. However, with prevalent 17q gains, reduced expression of TP53 (located on 17p) relative to driver genes expressed on 17q (for example, BIRC5, IGF2-BP1 and BRIP1) may favor genomic instability 44,45 .
Accurate risk stratification in neuroblastoma remains a major concern. Our data suggest a link between diverse criteria-age at diagnosis, segmental versus whole-chromosome gains and losses and acquisition of TMM 25,46 -based on how neuroblastomas evolve. We find that a greater age at diagnosis is often linked to longer evolution of the tumor rather than later origin. Paradoxically, this implies that low-risk tumors reach detectable size earlier than high-risk. Indeed, we infer that low-risk tumors have a substantially lower fraction of tumor cell loss than high-risk tumors and hence should grow faster. Acquired TMMs should, consequently, provide a larger selective advantage in the high-risk, late-MRCA group where they indeed are enriched. Hence, late-MRCA tumors may grow more rapidly only after gaining telomere maintenance (similar to IDH-wild-type glioblastomas 7 ) and hence are diagnosed later. Interestingly, we detected five neuroblastomas with amplified MYCN in the early-MRCA class, and, remarkably, all of these patients survived until the end of the observation periods (1,447-4,177 days), in contrast to the poor survival of patients with MYCN-amplified tumors in the late-MRCA group. Our findings suggest that MRCA timing may be worth considering as a parameter for patient stratification.

Online content
Any methods, additional references, Nature Portfolio reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/s41588-023-01332-y.

Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
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 license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license 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 license, visit http://creativecommons. org/licenses/by/4.0/. Article https://doi.org/10.1038/s41588-023-01332-y

Patient cohorts, tumor samples and ethical approval
Cohorts of primary and relapsed neuroblastoma tumors were retrospectively analyzed. Tumor material was collected as part of the diagnostic workflow of the German Neuroblastoma trial by the Society for Pediatric Oncology and Hematology and collected in the Neuroblastoma tumor bank. All trials were approved by the Ethics Committee of the Medical Faculty, University of Cologne, and collection and use of all tumor tissue material was approved (registry nos. NB97, NB2004 and NB2016). This study includes data of neuroblastoma tumors previously published in Hartlieb et al. 24 and Peifer et al. 26 The study by Hartlieb et al. also contains a subset of tumors from the St. Anna Kinderkrebsforschung at the Children's Cancer Research Institute in Vienna, Austria, as well as tumors analyzed in the registry trial INFORM 47 . The study office of the Neuroblastoma trial in Cologne provided clinical annotations and survival information. All patients or their legal guardian approved the use of tumor material by signed informed consent. For analysis, all resected tumors were divided into four quadrants, all of which were evaluated histologically. MYCN status was assessed as a routine clinical marker for all tumors using fluorescence in situ hybridization. A cross-sectional slice of one quadrant was used for DNA extraction for WGS; the same quadrant was used for ploidy analysis, measuring the DNA index. DNA was isolated using phenol chloroform extraction from fresh-frozen tumor material. Control DNA was isolated from whole blood using the NucleoSpin Blood DNA extraction kit (Macherey-Nagel) according to the manufacturers' instructions. Details of the included samples are provided in Supplementary Tables 1 and 6.
Estimates of tumor cell content were manually adjusted in one case (NBE40) following visual inspection of VAF distribution. Structural variants were excluded in the present study if they had a minimal event score of <5; focal amplifications were defined as regions with copy number ≥10; homozygous deletions were defined as regions with copy number <0.9. Deletions of (parts of) chromosomal arm 1p were defined as p-terminal regions lost relative to 1q and, moreover, with copy number ≤1 in near-diploid, ≤2 in near-triploid or ≤3 in near-tetraploid tumors. In analogy, we annotated deletions of chromosome 11q if the copy number was ≤1 in near-diploid, ≤2 in near-triploid or ≤3 in near-tetraploid tumors, and if 11q was lost relative to 11p. For gains on chromosomes 1, 2, 7 and 17 we required the copy number to be higher than the (rounded) basal ploidy of the tumor. Partial gains on 1q, 2p, 7q and 17q were defined as regions on the respective chromosomal arm that were gained relative to the other chromosomal arm.  50 . Only signatures contributing to ≥5% of the mutations in at least one sample and, in addition, identified in at least 10% of samples, were considered. The contributions of these signatures to each sample were then re-estimated using the R package mmsig v.0.0.0.9000 (ref. 51 ) (setting strandbias=F, bootstrap=F, cos_sim_threshold=0.01, force_include=c("SBS1", "SBS5")). For visualization, signatures SBS1, SBS5 and SBS40 were combined into a single, clock-like mutational signature.
Stratified analysis of clonal and subclonal mutations was performed by classif ying mutations as subclonal if , where ρ is estimated tumor cell content, n var and n ref are the number of variant and reference reads, respectively, and CN is copy number.

Mutation timing. Estimation of numbers of amplified and non-amplified clonal mutations.
We estimated mutation densities at partial and entire chromosomal gains, and at the MRCA of the tumor from the distribution of VAFs among clonal mutations. To this end we counted clonal mutations separately on each autosome, stratified by copy number (CN) state. Regions lacking a CN estimate were assigned to a specific CN state if the measured coverage ratio (CR) and B-allele frequency (BAF) matched the expected ratios within measurement error. Specifically, we required where CN i is the copy number of segment i, ρ the tumor cell content, π the average tumor ploidy and b the number of B-alleles. States with copy number >4 or of size ≤10 7 bp were excluded from the analysis because the statistics become ambiguous for short pieces and high CN states.
On each retained segment we estimated the number of clonal mutations, distinguishing clonal mutations present on a single allele ('non-amplified clonal mutations') from those present on multiple alleles ('amplified clonal mutations'). To this end we applied a statistical framework 52 , distinguishing amplified clonal mutations acquired before a clonal gain and thus present on all A-alleles, or on all B-alleles at a given CN state, from non-amplified clonal mutations acquired either on the non-amplified allele or after clonal gain but before the MRCA and hence present on a single copy; the remaining mutations are subclonal. Accordingly, we expect to find clonal mutations at is the average copy number of a given locus in the sample. Because measured VAFs are randomly distributed around their true values, we fit a binomial mixture model to the clonal mutations. To avoid misclassification of subclonal mutations as clonal, we neglected all variants with VAF < ρ ζ . By symmetry, this cutoff retains half of the non-amplified clonal mutations, introducing a subsequent correction factor of 2 (equation (4)). Defining the weights w = (w 1 , w CN−b , w b ) , for the non-amplified clonal and amplified clonal mutation peak, respectively, we computed the probability of measuring n var, j variant reads and n ref,j reference (nonvariant) reads at the position of the jth SSNV to: where B (n var, j ; n var,i + n ref,j , ρk ζ ) is the binomial probability for drawing n var, j variant reads at sequencing depth n var, j + n ref, j from the peak Article https://doi.org/10.1038/s41588-023-01332-y comprising mutations that are clonal on k chromosomal copies. Defining a uniform prior probability, P(w), for the weights, we computed, up to normalization, the posterior probability as: where N is the total number of SSNVs under consideration. Clonal mutations were then assigned to distinct clonal peaks according to the weights at maximum a posteriori probability (MAP), yielding, on segment l, estimates for the number of clonal mutations on non-amplified chromosomes, n 1,l , on amplified b alleles, n b,l (if b > 1) and on amplified a alleles (with copy number CN − b), n CN−b,l , according to: Mutation timing. Timing MRCA and ECA. Mutation densities (SSNVs per bp) at the MRCA were estimated from the number of non-amplified clonal mutations and total size of the analyzed genome, g = ∑ l g l (ref. 52 ). The number of mutations per copy that were acquired up to the MRCA is If tumor samples were well mixed, or tumors completely sampled, mutation densities at the MRCA could be directly estimated as n l /g. In practice, however, n l may consist of a set of true clonal mutations, acquired before tumor growth, and an additional set of mutations that appear as clonal in the tumor sample due to incomplete sampling. To correct for the latter, false-positive (FP), clonal mutations we compared primary and relapse samples from two such pairs available in our dataset (NBE11/NBE66, NBE51/NBE78). The fraction of conservative clonal mutations in the primary sample that remained undetected in the relapse sample was 15 ± 5% and was taken as the fraction of FP. With this correction, the mutation count, m MRCA , and mutation density, m MRCA , at MRCA were estimated as: respectively. Lower and upper 95% confidence bounds for m MRCA were estimated by bootstrapping the genomic segments 1,000 times. Next, we tested for each gained segment whether amplified clonal mutations were significantly less frequent than expected at the MRCA and, accordingly, assigned the segment to either the MRCA or an earlier time point (in the majority of cases the ECA; here we excluded a small number of segments in seven tumors with a higher density of amplified clonal mutations than the estimated mutation density at the MRCA, because such gains may be subclonal CN alterations erroneously classified as clonal). To this end, we modeled the number of mutations falling on each genomic segment with a negative binomial distribution, which accounts for overdispersion caused by heterogeneous mutation rates along the genome 53 . The probability that the gain of genomic segment l coincided with the MRCA is then computed as: where: .
Here, n k,l;k∈{b,CN−b} is the number of amplified clonal mutations on segment l and g l is the respective segment size. We corrected the P values obtained with equation (7) for multiple testing using Holm correction (false discovery rate ≤ 0.01) and, accordingly, assigned each segment to either the MRCA or an earlier time point.
Finally, we computed the mutation densities at ECAs from the segments with significance level α = 0.01 to: m ECA = ∑ l,p adj,l ≤0.01 n b,l + n CN−b,l ∑ l,p adj,l ≤0.01 g l,b + g l,CN−b . (8) We then tested for each contributing segment whether its mutation load conformed to the ECA according to a negative binomial distribution (in analogy to equation (7)) and computed lower and upper 95% confidence bounds by bootstrapping, as before for the MRCA.
Mutation timing. Translation of mutation densities into estimated weeks p.c. We related mutation densities per haploid genome into weeks post conception (p.c.) by inferring SSNV rates per diploid genome and embryonic day (µλ) and mutation and division rates per day, using the measured VAF distributions and age at diagnosis as outlined in Real-time estimation of cell division rate (with similar results based on all SSNVs or only clock-like SSNVs). Because mutation calling was performed by comparing tumors against a matched blood control, mutation densities correlate with the time post gastrulation (at approximately 2 weeks p.c.). Thus, the mutation density per haploid genome, m, relates to the time p.c., t, according to m (t) = The estimated time of birth was taken as 38 weeks after gastrulation (40 weeks p.c.).

Survival analysis
Survival analysis was performed using the R package survival v.3.1.12 (ref. 54 ). We detected a clear bimodality in the histogram of SSNV density at the MRCA in the discovery cohort and took the upper bin border just before the minimum (0.05 SSNVs per Mb) as threshold to split tumors into groups of early or late MRCA. The same value was found appropriate in the validation cohort.

Modeling neuroblastoma initiation
Modeling emergence of the MRCA. We modeled neuroblastoma initiation by sequential acquisition of driver mutations in two oncogenic events, assuming that both events are infrequent and occurring in early neuroblasts with low probabilities, µ 1 and µ 2 , during cell divisions. We denote the selective advantage conferred by the driver mutations acquired during the two events by r and s. Specifically, we assume that one or more first events generate a precancerous cell population in which a second event creates the MRCA. On this basis, we compute the time-dependent probability of the emergence of the MRCA. Neuroblasts initially expand rapidly, but how this cell population behaves subsequently, because sympathetic neurons differentiate from these precursors, is not known precisely. Hence, two possible scenarios were considered: (1) an initial phase of exponential growth of the neuroblast population is followed by a subsequent phase of differentiation, modeled by exponential decay; and (2) an initial phase of exponential growth is followed by a subsequent phase of precursor homeostasis (for a related model of a two-step process of tumorigenesis, but in a homeostatic tissue, see ref. 55 ). The two phases are associated with distinct rates of cell division, λ 1 and λ 2 , and loss, δ 1 and δ 2 , where λ 1 > δ 1 and λ 2 ≤ δ 2 (with equality in the case of precursor homeostasis).
Thus the population dynamics of neural precursor cells, N(t), are described by: where T denotes the time point at which the cell population peaks. Cells undergoing the first oncogenic event are generated at rates: Article https://doi.org/10.1038/s41588-023-01332-y For simplicity, we take the selective advantage associated with the first oncogenic event into account during the contraction (or homeostasis) phase and neglect it during the initial rapid expansion; this approximation is appropriate when the selective advantage of the first event is comparatively small, which is subsequently borne out by the parameter inference. Thus, M 1 cells harboring the first mutation grow at rate λ 1 − δ 1 during tissue expansion and at rate λ 2 − δ 2 /r during tissue contraction or homeostasis. Hence, depending on the actual value of r > 1, a clone harboring the first event may slowly shrink or expand for t > T. We now ask whether in this clone the second oncogenic event takes place that defines the MRCA of the tumor. The second event occurs at rate: Using the survival probability of the supercritical birth-death process, we have a cell undergoing the second oncogenic event survive with probability: during the expansion phase (E) and, provided that δ 2 sλ 2 < 1, with probability: during the decay or homeostatic phase (D). We are interested in the probability that at least one surviving cell underwent both oncogenic events, P MRCA . There are three possible cases: (1) both oncogenic events occur during precursor expansion, associated with probability P MRCA,1 ; (2) the first oncogenic event occurs during precursor expansion, the second during precursor contraction or homeostasis, associated with probability P MRCA,2 ; and (3) both oncogenic events occur during precursor contraction or homeostasis, associated with probability P MRCA,3 . For each case we assumed that the number of cells with only the first event, M 1 , is small compared with the number of normal cells. We have: We derive expressions for P MRCA,I , P MRCA,2 and P MRCA,3 in Supplementary Note 1a for the case of precursor expansion and decay, and in Supplementary Note 2a for the case of precursor expansion and homeostasis. If precursor expansion is followed by decay, we find: where F = ∫ 1 0 ν 2,I / (ν 2,I z α ) dz and α = δ 1 −sλ 1 . Alternatively, if precursor expansion is followed by homeostasis, P MRCA,2 and P MRCA,3 take the form (Supplementary Note 2a): Equations (14a-e) give the probability that a growing tumor clone has emerged up to time t for the distinct cases.
Modeling the ECA. The ECA is associated with the first oncogenic event, which may have occurred at variable time points before the MRCA. To time the ECA, we computed the conditional probability of having undergone the first oncogenic event before t 1 , given that the second oncogenic event occurred at t 2 , denoted by P(t 1 |t 2 ). The acquisition of the second oncogenic event is proportional to the number of cells with the first event M 1 (equation (11)). As denoted by M 1 (t 2 |τ), the number of cells at t 2 resulting from a first event at τ.P(t 1 |t 2 ) can hence be expressed as: Distinguishing the three cases for the timing of the second event relative to the first, as before, we find for precursor expansion followed by decay for two cases (Supplementary Note 1b). When the second event occurs in the exponential growth phase, then: This corresponds to the classical Luria-Delbrück model. Recalling that t = 0 marks the beginning of neuroblast expansion, the first mutation happens, as expected, with uniform probability during exponential growth. When the second event occurs during the decay phase, the probability for the first event is uniform in the exponential phase and decreases in the decay phase according to: where Θ (⋅) is the Heavyside step function and, of course, 0 ≤ t 1 < t 2 . For precursor expansion followed by homeostasis, P (t 1 |t 2 ) reads (Supplementary Note 2b): and: Relating model and data. To relate the model to the measured data, we translate time into mutation counts. Assuming that SSNVs Article https://doi.org/10.1038/s41588-023-01332-y are accumulated at a constant rate, µ, per cell division, the expected mutation count per cell is Poisson-distributed with rate µλt.
Parameter estimation. We estimated the following parameters: peak size of the neuroblast population (N), relative loss rates in the growth and decay phases (δ 1 /λ 1 and δ 2 /λ 2 , respectively), rate of SSNVs (µ) and rates of first and second oncogenic events (µ 1 and µ 2 ) and their selective advantages (expressed as r and ν 2 ), using equations (9)(10)(11)(12)(13)(14)(15)(16)(17) and Approximate Bayesian Computation with Sequential Monte Carlo sampling (ABC-SMC) as implemented in pyABC 56 . We used a population size of 1,000 parameter samples and prior probabilities as outlined in Supplementary Table 11. Evaluation was abrogated after 25 SMC generations, or if ε ≤ 0.05. We used mutation density estimates at MRCA and ECA of high-risk tumors (primary tumors and metastases) as determined by equations (6) and (8)  For each sampled parameter set we performed the following steps: (1) Sample for each tumor a time point of the MRCA, t 2 . To this end, sample a uniform number x between 0 and 10 -5 , thus accounting for the overall incidence of 10 -5 . Then, from equation (9), determine the time point at which P t = x. To facilitate numerical computation, we approximated the sum in equation (13a) with an integral. To exclude second hits that do not confer a selective advantage during expansion, we required s = max (1, (5) Sample for each sampled t 2 a neutral mutation count from a Poisson distribution with mean µt 2 and divide by the haploid genome length of 3.3 × 10 9 , yielding m ECA,sim . (6) Determine the simulated incidence of the ECA at the experimentally determined mutation loads: I ECA,sim,i = ∑m ECA,sim < i; i ∈m ECA,ev . (7) Determine the simulated incidence of the MRCA at age 10 years, I MRCA,sim,ten years using equation (9). This step was introduced to contrast the incidence at old ages with the clinically observed overall incidence of the order of 10 −5 , which we weighted by assuming an error of 10 −4 .
We computed the cost function, d, as To enforce good fits to the initial phase of the incidence curve for better comparison of the contraction and homeostasis models of neural precursor dynamics, we chose weights Ninety-five per cent posterior probability bounds for the model fits were estimated by simulating the model at each sampled parameter set and cutting off 2.5% at each end of the simulated distribution.
To assess the robustness of the model, we performed an additional parameter estimation on clock-like mutations only. To this end, we multiplied mutation densities at MRCA and ECA by the fraction of mutations generated by clock-like mutational signatures (SBS1, SBS5 and SBS40) and fit the model to the clock-like mutation densities.

Modeling mutation accumulation during tumor growth
Model. Denoting the rate at which tumor cells divide by λ T and the loss rate by δ T , we modeled the number of tumor cells, N T , with an exponential growth model, N T (t) = e (λ T −δ T )t . We assumed that some SSNVs are already present in the founder cell of the tumor and denote their number by n clonal ; these mutations are clonally propagated to the entire tumor and will thus be found at frequency f = 1. Additional SSNVs are acquired during tumor growth, and we denote their number by n subclonal ; these mutations are present in a subset of the tumor only and will thus be found at f < 1. The VAF distribution of a neuroblastoma is hence a superposition of clonal and subclonal mutations accumulated before and during tumor growth, respectively.
To model the number of subclonal mutations, we used a model of neutrally evolving tumors that accounts for the stochastic expansion of neutral subclones while assuming exponential expansion of the tumor mass overall 57 . This model assumes that neutral mutations are acquired at all times at rate μλ T N T (t) and drift stochastically according to a supercritical birth-death process, where 58 : where: Together, this yields the site frequency spectrum, S(i, µ), of subclonal mutations: The total number of subclonal mutations in a tumor is computed as; (21) where N T (t end ) is the number of tumor cells at diagnosis. The number of mutations present in subclones of at least a and at most b cells is, in analogy to equation (21) Because clone sizes a and b are large, we approximate the sum in equation (22) by an integral, yielding: Article https://doi.org/10.1038/s41588-023-01332-y dt. (23) Relating model and data. To relate the model of mutation accumulation to the measured VAF distribution, we modeled mutations on each copy number state separately. To this end, we denote by n f,k the number of mutations at frequency f on genomic segments with copy number state k. Note that f reports the fraction of mutated cells among all tumor cells whereas VAFs report the fraction of mutated alleles in the tissue sample. The two quantities can be converted into each other using the following relation: where ζ is the average copy number of the sample as defined in equation (1). Clonal mutations are associated with frequency f = 1. The number of clonal mutations falling within genomic segments of copy number k, n 1,k , is expected to scale with g k , the fraction of the genome at copy number k, according to To relate true with measured VAFs, we assumed that the latter are binomially distributed around the former according to where C k ∝ Pois(Ĉ k ) is the sequencing coverage on a segment with copy number k, and ρ is tumor cell content. VAF. We excluded mutations with Ṽ AF < 0.1 from the fit and ran the Mobster setting autosetup = "FAST".
This resulted in an additional exclusion of nine cases for which Mobster suggested subclonal selection. Thus we selected 62 tumors with well-resolved subclonal VAF histograms and no evidence of subclonal selection for parameter inference (Supplementary Table 1).
Parameter estimation. We estimated n clonal , µ and δ T /λ T using ABC-SMC as implemented in pyABC 56 . We used a population size of 1,000 parameter samples and prior distributions as outlined in Supplementary  Table 12. Termination criteria were the same as above (Modeling neuroblastoma initiation).
We stratified measured VAFs by copy number, excluding copy numbers >4 or present on, in total, <10 8 bp. For each copy number k, we merged mutations from the clonal VAF peaks constituted by amplified and non-amplified mutations (see Mutation timing for the definition of amplified and non-amplified clonal mutations). To this end, we first assessed the average coverage Ĉ k from all mutations falling on segments with copy number k. Then, we classified mutations at copy number k as amplified clonal if quantile of a binomial distribution with success probability ρl ζ and where l is the B-allele count. These mutations were then merged with those of the non-amplified clonal peak by multiplying their frequencies by l k and adding them l times.
Finally, we computed the cumulative mutation counts of the measured data, F k,ev (f) = ∑VAF k > f , where f runs from 0.05 to 1.00 in steps of size 0.05, and extrapolated the cumulative mutation counts to the whole genome with multiplication by . At a sampled parameter set for n clonal , µ and δ T , the following steps were performed for each copy number state k included in the analysis: (1) Sample for each clonal mutation a sequencing coverage C k according to Pois(Ĉ k ). (2) Sample for each clonal mutation a VAF according to equation (16). (3) Determine n f,k;f≠1 from equation (13), assuming a tumor size of 10 9 cells at diagnosis, and evaluating equation (13) in bins of size 0.05 at the lower limit: S (i, μ).
(4) Sample for each subclonal mutation a sequencing coverage C k according to Pois(Ĉ k ) and a VAF according to equation (16). The cost function is

Real-time estimation of cell division rate
In the previous sections we describe how we inferred the parameters for our models of neuroblastoma initiation and growth from mutation data. This yielded estimates for the rates at which cells are lost, as well as the rates at which neutral and oncogenic mutations are acquired in units of cell divisions. To convert these estimates to real time we estimated the actual cell division rate, using the age distribution at diagnosis and the approximate tumor size at diagnosis. We reasoned that the time span between gastrulation and tumor diagnosis (t D ) consists of two phases: premalignancy up to the formation of the MRCA and malignant growth of the tumor thereafter. Hence, where we assume exponential tumor growth until diagnosis. Thus the cell division rate, λ, can be expressed as To estimate λ with equation (28), we combined the clinical information with parameter estimates from our models as follows: we know t D (the age at diagnosis, A, plus approximately 250 days of embryogenesis after gastrulation) from the clinical data, and estimate that a tumor of a few cubic centimeters consists of the order of N T (t D ) = 10 9 cells. As described below (Mutation timing and Modeling neuroblastoma initiation), we also have estimates for mutation density at the MRCA, m MRCA , and for the SSNV rate per actual cell division and haploid genome, µ. Finally, we obtained an estimate for the effective rate of acquisition of Article https://doi.org/10.1038/s41588-023-01332-y (neutral) SSNVs during tumor growth from the subclonal VAF histograms of 62 tumors (compare with Modeling mutation accumulation during tumor growth), which is related to δ and s as Substituting equation (29) in equation (28), we obtain for each of the 62 tumors (labeled with the index i) an estimate for the division rate with mean, ⟨λ i ⟩ = to the VAF histograms of each of the selected 62 tumors individually hence yields tumor-specific division rates along with relative death rates, which we stratify by TMM.
Finally, one also obtains an estimate for the daily mutation rate during tumor initiation by computing μλ T,i with associated uncertainty μΔλ T,i + λ T,i σ(μ), which relates molecular clock to real time. For this purpose, we average across the inferences from the 44 primary tumors/ metastasis, excluding 18 relapsed tumors among the 62 submitted for analysis.

Statistics and reproducibility
All statistical tests were computed with R (v.3.6.0 and v.4.0.0), and details (statistical tests and whether one-or two-sided, exact sample size, P values and test statistics) are specified in the respective figures and accompanying Source data. This is a retrospective analysis of tumor material that was collected as part of the diagnostic workflow of the German Neuroblastoma trial by the Society for Pediatric Oncology and Hematology and collected in the Neuroblastoma tumor bank. Hence, no statistical method was used to predetermine sample size. All analyzed tumors had a clear copy number profile, tumor cell content ≥25% and no hypermutation genotype; no data were excluded from the analyses. The experiments were not randomized. The investigators were not blinded to allocation during experiments and outcome assessment.

Reporting summary
Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.

Data availability
Subsets of WGS and RNA sequencing data were part of previously published studies 24,26 . Data from these studies are deposited at the European Genome-Phenome Archive (https://www.ebi.ac.uk/ega/) under accession nos. EGAS00001004349 and EGAS00001001308. Additional WGS data generated for this study are available at the European Genome-Phenome Archive under accession nos. EGAS00001004990 and EGAS00001006533. In accordance with the laws of data protection, data are deposited under controlled access. Access can be granted by contacting Frank Westermann (f.westermann@kitz-heidelberg.de) and requires a data access agreement; requests will be replied to within 4 weeks. Variant calls (SNVs, indels, SVs and copy number variations), mutational signatures, model fits and a summary of the mutation profile and modeling results for each tumor can be accessed at Mendeley (https://doi.org/10.17632/ m9pwjbm7c8.1) 70 . All remaining data are available in the Supplementary information. Source data are provided with this paper. Article https://doi.org/10.1038/s41588-023-01332-y Extended Data Fig. 3 | Examples of early evolution. a, Copy number profiles of an example tumor with near-triploid genome and without a timeable ECA. The color panels annotate genomic segments of equal copy number. Next to the copy number profile, the densities of non-amplified (green) and amplified (color-encoded) clonal mutations on segment of length ≥10 7 bp are shown, with the dashed line showing the kernel-density estimate of the distribution of non-amplified clonal mutations. Chromosomal segments of equal copy number were merged and Holm-corrected one-sided P values were computed based on a negative binomial distribution (**, p adj < 0.01, exact P values are provided in Source Data). The bottom panel shows mutation densities at ECA and MRCA with 95% confidence bounds estimated by bootstrapping. The horizontal lines show mutation densities on gained segments with 95% confidence bounds computed from χ 2 distributions. b, Mutation densities at ECA and MRCA of early-MRCA tumors with (left) and without (right) timeable ECA. Solid lines represent maximum likelihood estimates and shaded areas represent 95% confidence intervals, as obtained by bootstrapping. c and d, As in (a) but for a near-diploid tumor with timeable ECA (c) and for a near-tetraploid tumor without timeable ECA (d). Holm-corrected one-sided P values were computed based on a negative binomial distribution (**, p adj < 0.01, exact P values are provided in Source Data). e, Mutation densities at MRCA in early-MRCA (n = 20) and late-MRCA primary tumors (n = 47, thereof n = 26 with ECA and n = 21 without ECA). f, Fraction of polysomic chromosomes in near-triploid (n = 33) and near-tetraploid (n = 12) tumors that were generated in a single oncogenic event. As a control, the fraction of disomic chromosome whose mutation density matched with the MRCA is shown for near-diploid tumors (n = 55). In (e) and (f), boxes show median, 25 and 75% percentiles, whiskers extend to the smallest and largest value within 1.5x interquartile range.
Article https://doi.org/10.1038/s41588-023-01332-y Extended Data Fig. 5 | Neuroblastoma initiation during progenitor expansion. a, Two-dimensional projections of the posterior probability distribution of the model parameters (considering neuroblastoma initiation in a transient population of early neuroblasts). b, Predicted expansion of neuroblasts (grey) and the selected subclone upon acquisition of the first oncogenic event (dark green) in a model of transient neuroblast expansion. Colored areas show the 95% posterior probability bounds (estimated from simulations using 1,000 samples from the posterior probability distribution). Vertical line and shaded area give mean and 95% CI for the estimated time of birth (computed from n = 62 primary tumors). c, Comparison of parameter estimates if fitting the model to all SSNVs or to SSNVs generated by a clock-like mutational process only. For the latter, mutant densities were adjusted by the fraction of SSNVs explained by SBS1, SBS5 or SBS40. For each parameter, median and 80% credible intervals are shown (estimated from n = 1000 samples of the posterior distribution). d, Comparison between the estimated time point at which the MRCA (left, n = 95 primary tumors/metastases with TMM) and the ECA (right, n = 47 primary tumors/metastases with TMM) emerged if using all SSNVs or if using SSNVs that were generated by a clock-like process only. Time is measured in weeks post conceptionem (p.c.).