Spatial and temporal intratumour heterogeneity has potential consequences for single biopsy-based neuroblastoma treatment decisions

Intratumour heterogeneity is a major cause of treatment failure in cancer. We present in-depth analyses combining transcriptomic and genomic profiling with ultra-deep targeted sequencing of multiregional biopsies in 10 patients with neuroblastoma, a devastating childhood tumour. We observe high spatial and temporal heterogeneity in somatic mutations and somatic copy-number alterations which are reflected on the transcriptomic level. Mutations in some druggable target genes including ALK and FGFR1 are heterogeneous at diagnosis and/or relapse, raising the issue whether current target prioritization and molecular risk stratification procedures in single biopsies are sufficiently reliable for therapy decisions. The genetic heterogeneity in gene mutations and chromosome aberrations observed in deep analyses from patient courses suggest clonal evolution before treatment and under treatment pressure, and support early emergence of metastatic clones and ongoing chromosomal instability during disease evolution. We report continuous clonal evolution on mutational and copy number levels in neuroblastoma, and detail its implications for therapy selection, risk stratification and therapy resistance.

N euroblastoma is the most common solid paediatric tumour, accounting for 15% of cancer-related deaths in early childhood. Especially patients directly diagnosed with high-risk neuroblastoma have a poor prognosis 1,2 . High-risk (HR) and ultra-high-risk (UHR) neuroblastoma are characterised by active telomere maintenance and mutations in distinct pathways defining specific clinically relevant molecular risk classifiers 3 . Telomere maintenance as the hallmark of the UHR concept is either conveyed by TERT gene activation due to TERT promoter mutations, TERT re-arrangements or high MYCN expression or by alternative lengthening of telomeres (ALT) 4 . Intriguingly, even high-risk neuroblastomas initially respond to treatment, but eventually relapse as chemotherapyresistant disease. Current target prioritisation efforts in relapsed/refractory patients such as INFORM, MAPPYACTS, NCI-COG MATCH, PROFYLE or ITHER 5-8 rely on molecular profiling based on bulk sequencing approaches in single tumour biopsies. Targeted therapy based on single-biopsy profiles is often inefficient or only exhibits short-term effects, and rarely results in long-term tumour control. Secondary resistance following promising short-term efficacy is clinically observed in particular for activating ALK mutations, the most recurrent actionable single-nucleotide variant (SNV) present in neuroblastoma for which targeted inhibitors of 1st, 2nd and now 3rd generation have entered clinical trials [9][10][11][12] . These clinical observations indicate the existence of therapyresistant subclones in high-risk neuroblastoma. Indeed, neuroblastoma shows extensive genetic intratumour heterogeneity and distinct evolutionary patterns [13][14][15] . Both linear/latebranching as well as parallel/early-branching evolution occurs in neuroblastoma, and impacts clinical behaviour differently, with early branching being associated with adverse outcome 14 . Thus, intratumour heterogeneity is an important feature of neuroblastoma, which deserves in-depth analysis and further interpretation with respect to potential clinical implications.
While neuroblastoma typically has few SNVs, neuroblastoma genomes show pronounced somatic copy-number alterations (SCNA) 16 . Chromosomal instability describes the progressive accumulation of SCNAs over time, equating to changes in the number and structure of chromosomes during tumour evolution. Chromosomal instability results from chromosome segregation errors during cell division, and high levels have been linked to poor prognosis 17,18 and multi-drug resistance 19,20 in many solid tumours. Neuroblastoma is widely considered to be driven by chromosomal instability since high levels of multiple clinically relevant SCNA events occur with a low overall mutational burden. Focal amplifications of the MYCN oncogene, re-arrangements within the TERT locus on chromosome 5p and SCNAs on chromosomes 1p36, 11q and 17q, genetically define a high-risk clinical phenotype 4,21,22 . Despite the clinical relevance of chromosomal instability in neuroblastoma and recent evidence in other cancers that ongoing chromosomal instability keeps shaping cancer genomes during disease progression [23][24][25] , the extent of ongoing chromosomal instability after initial biopsy at diagnosis and the dynamics of SCNA acquisition during neuroblastoma progression and treatment remain unclear.
To address open questions about neuroblastoma intratumour heterogeneity, we combine multi-region transcriptome and whole-exome sequencing, followed by ultra-deep targeted sequencing of 140 spatially and temporally separated tumour samples from 10 clinically heterogeneous patients with neuroblastoma. We describe spatial and temporal intratumour heterogeneity of actionable mutations and distinct patterns of tumour evolution based on SNVs and SCNAs.

Results
Analysis of SNVs in spatially and temporally distinct neuroblastoma samples reveals extensive intratumour heterogeneity. To analyse neuroblastoma heterogeneity within the tumour tissue and across the course of disease in individual patients, we performed genomic and transcriptomic analysis of multiple spatially and temporally distinct samples per patient. We performed whole-exome sequencing (WES, N = 51 biopsies), targeted sequencing (N = 140 biopsies) and RNA sequencing (N = 48 biopsies) from 10 patients at the age of 0.5-14 years (3 girls and 7 boys) treated for neuroblastoma at the Charité Berlin, according to the German NB2004 trial and NB Registry 2016 protocols ( Fig. 1 and Supplementary Fig. 1). Seven of the ten patients had a high-risk, 1 child an intermediate-risk and 2 patients a low-risk disease according to INRG/INSS classification. Between 2 and 10 biopsies per patient were analysed by WES (Fig. 1). For 3 patients, samples from 2 different time points were included. Samples were biopsied from primary tumours at diagnosis (N = 17), after 4-6 cycles of induction chemotherapy at tumour resection (N = 26) or at diagnosis of relapse (N = 8). Primary (N = 45) and metastatic tumour sites (N = 6) were biopsied (Fig. 1). To increase the sensitivity for SNV and SCNA detection, tumour regions with a high tumour purity (>60% for WES and RNA sequencing) were macrodissected from sequential sections. Matched normal blood controls were collected once from each patient at diagnosis. SNVs and SCNAs were called from WES data. To assess the clonality of SNVs detected with WES at higher resolution, we performed ultra-deep targeted sequencing (mean coverage 2500×). To this end, following WES analysis of 51 samples, we included further 89 samples analysing a total of 140 samples by ultra-deep targeted sequencing of 1479 WES-detected SNVs.
Consistent with previous reports 3,4,21,26 , the total number of somatic, non-silent SNVs varied between 0 and 60 (median: 15) per sample (Fig. 1c). As expected, low-risk neuroblastomas had significantly fewer non-silent mutations than high-risk neuroblastomas (p = 9 × 10 −5 , two-sided Wilcoxon test) 3,10,21 . For each patient, SNVs were categorised as "clonal" if detected in all samples from that patient or "subclonal" if detected in a subset of samples from that patient or "subclonal/specific" if exclusively detected in one sample (Fig. 1c). The proportion of clonal SNVs varied from 0 to 87% among samples (average over patients: 37%), while on average 22% of mutations in each sample were categorised as subclonal/specific (Fig. 1c), pointing towards significant intratumour heterogeneity. SNVs were considered to be clonal if present with a variant allele frequency (VAF) greater than 10% in WES or targeted sequencing data in all samples. This threshold was chosen to avoid overestimation of heterogeneity. All other detected SNVs were considered subclonal. Roughly 45% of SNVs were clonal, irrespective of whether the SNVs were silent, non-silent or were detected within cancer-related genes (Supplementary Fig. 1b and Supplementary Data 1). For the large majority of SNVs (87.9%), assignment as clonal or subclonal were consistent between WES and targeted sequencing datasets ( Supplementary Fig. 1c). As expected, the overall proportion of SCNAs markedly differed among patients. In line with previous reports, MYCN amplification, a hallmark of high-risk neuroblastoma, showed little variation between spatially distinct samples, but pronounced differences in patient CB1003 between time points (Fig. 1c). Our data demonstrate extensive spatial and temporal genetic heterogeneity across distinct neuroblastoma biopsies from the same patient in our cohort.
Spatial and temporal heterogeneity of single-nucleotide variants affects therapeutically actionable genes in neuroblastoma.
To determine biologically and potentially clinically relevant genetic intratumour heterogeneity, we assessed non-silent SNVs in known cancer-related genes (Fig. 2). Whereas between two and six non-silent SNVs were detected in 801 selected known cancerrelated genes (Supplementary Data 1 and Supplementary Table 2; defined by the COSMIC database 27 and other studies 8,10,16,28 ) in high-risk tumours, very few or no non-silent SNVs in the selected cancer-related genes were detected in neuroblastomas from other risk groups (Fig. 2a). Most SNVs were detected in all samples (Fig. 2a) from a given tumour, indicating the clonal presence and early evolutionary emergence. Known pathogenic mutations in HRAS were among these clonal mutations, and a clonal ALK mutation was detected in patient CB1003. However, some detected mutations appeared spatially and/or temporally heterogeneous (e.g. ALK R1275Q in patients CB1002 and CB1008, as well as FGFR1 in patient CB1002, Fig. 2a). The persistence of clonal mutations in ATM, PTK2, GLI3 and CLTCL1 in all tumour locations analysed in samples from diagnosis and at relapse in patient CB1002 (Fig. 2b) could indicate the relevance of these SNVs for tumour progression. In contrast, high overall numbers of de novo SNVs and mutations in the druggable genes, ALK, BRAF and FGFR1, were only detected in a subset of relapse samples from patient CB1002 compared to the primary tumour at diagnosis (in line with previous studies reporting an increase in mutational burden and de novo MAPK pathway mutations in relapsed neuroblastoma 26,29   Only SNVs detected in one or more samples are displayed as one column per SNV. Mutations in cancer-related genes are explicitly named, and were indicated as clonal (VAF > 10% in all samples from a tumour, green triangles) in comparison to subclonal oncogenic mutations (red triangles). Row numbers identify samples at distinct time of biopsy collection within the disease course. Columns adjacent to the right of heatmaps indicate further characteristics. Samples denoted with WES were already included in the exome sequencing data. Samples denoted with FFPE were formalin-fixed and paraffin-embedded before analysis was performed, all other samples were fresh-frozen biomaterial. targets clearly presented in the case of patient CB1002 is likely not a rare exception, creating the need to re-evaluate targeted therapy decision-making.

Relapse
In another example, the actionable ALK R1275Q mutation was observed at high variant allele frequency at diagnosis, but not detectable by WES after chemotherapy in patient CB1008 (Fig. 2a). Targeted ultra-deep sequencing, however, was able to recover this mutation in the sample after therapy, albeit at very low allele frequency ( Supplementary Fig. 2f). This was further validated by a droplet digital PCR (ddPCR) assay confirming high VAF of the ALK R1275Q at diagnosis but loss of the mutation in 8 out of 10 samples from the resection (Supplementary Fig. 3). A mutation in the KIF1B tumour suppressor was detected (WES) only after chemotherapy in patient CB1008, which was confirmed by targeted ultra-deep sequencing, suggesting that this mutation occurred de novo during therapy ( Fig. 2a and Supplementary  Fig. 2f). Finally, clonal expansion of a mutation after therapy was also observed in a subset of cases, such as the DPYSL5 30 mutation in patient CB1003 (Fig. 2c). We conclude that actionable and pathogenic SNVs are spatially and temporally heterogeneous in neuroblastoma, revealing clonal evolution prior to and under selective treatment pressure with potential implications for targeted therapy decisions.
High degrees of spatial copy-number heterogeneity and ongoing chromosomal instability detected in a subset of neuroblastomas. To investigate the degree of ongoing chromosomal instability in neuroblastoma, allele-specific SCNAs were called from WES data for all 51 samples from the 10 patients. Quality control revealed a low signal-to-noise ratio in patient CB1006, and this patient was excluded from further analyses, yielding a final set of 48 samples from 9 patients. Allele-specific SCNA profiles underwent reference phasing, phasing-based correction of copy-number states and detection for structural parallel evolution using refphase as previously described 23 . The phased and corrected allele-specific SCNA profiles included information on subclonal SCNAs derived from opposite alleles as mirrored subclonal allelic imbalances (MSAI) and were subjected to tree reconstruction using MEDICC2 31,32 , including reconstruction of the most recent common ancestor to detect early events prior to clonal diversification.
To allow for consistent comparison of SCNA events across patients, SCNAs were aggregated on a per-cytoband level and stratified into clonal (present in all samples from that patient) and subclonal (present in some samples from that patient) events. Strong heterogeneity was detected on the level of SCNAs across the cohort. Across the cohort, between 13 and 92% (median 34%, mean 42%) of the genome of each patient's cancer was affected by at least one SCNA event in at least one sample ( Fig. 3a and Supplementary Fig. 4). We detected subclonal SCNA events in every patient, with a median of 57 segments affected by subclonal SCNA events compared to 132 clonal events. The majority of SCNAs were clonal for most patients (CB1001, CB1002, CB1007, CB1009 and CB1010), indicating early structural evolution followed by relative stasis. Contrary to this, a similar number of clonal and subclonal SCNAs were detected in patients CB1003, CB1004 and CB1005, demonstrating ongoing clonal evolution and chromosomal instability (Table 1). Patient CB1008 showed a disproportionally large number of subclonal, compared to clonal events, indicative of strong clonal diversification and ongoing chromosomal instability between samples obtained at diagnosis and at tumour resection after first cycles of induction therapy. This patient showed repeated copy-number neutral loss of heterozygosity (LOH) events in large areas of the genome (chromosomes 7-9, 12, 14-16, 17 and 19) only in the tumour resection samples, while sharing the clinically relevant SCNA events, LOH of 1p and MYCN amplification, among diagnostic and resection samples. Across the cohort, aggregated SCNA profiles show clear non-uniformity along the genome ( Fig. 3) with hotspots in the clinically relevant regions, 1p36 (clonal LOH or clonal gain events in 8 of 9 patients), 11q (clonal LOH or gains in 5 of 9 patients) and 17q (clonal gain events in 8 of 9 patients). Interestingly, even though MYCN amplifications were detected in only 3 of 9 patients, we observed clonal gains of 2p24 in 6 additional patients, making the MYCN locus the most frequently gained region in our cohort even in low-and intermediate-risk patients. In line with previous reports 33 , subclonal or clonal gains of chromosome 7q were detected in 8 of 9 patients. We also detected repeated MSAI events in 3 of 9 patients (Fig. 4) indicative of parallel evolution, but a low sample size prevented assessment of statistical significance.
The clinical importance of SCNAs in neuroblastoma strongly depends on the size of the events. Segmental SCNAs generally confer a poor prognosis whereas whole-chromosome gains and losses are considered neutral towards patient outcome [34][35][36] . To explore this, we stratified chromosomes into unchanged (diploid), whole-chromosome altered and segmentally altered chromosomes, and separately tallied their fractions in each category for each patient (Fig. 3b). The majority of chromosomes were affected by segmental duplications in most patients. In two patients (CB1007 and CB1009), however, most chromosomes were instead affected by whole-chromosome gains and losses (Fig. 3b). Patients CB1007 and CB1009 were the only low-and intermediate-risk patients in our cohort for whom SCNA data was available, confirming prior observations about the clinical relevance of segmental versus whole-chromosome SCNAs in neuroblastoma 37,38 . In summary, SCNAs are highly spatially and temporally heterogeneous in neuroblastomas, and in some cases are continuously acquired, indicating ongoing chromosomal instability.
Inference of evolutionary trajectories allows relative timing of the emergence of the metastases in neuroblastoma evolution. We next investigated each patient individually for evidence of early clonal diversification by inferring clonal hierarchies from variant allele frequency-based SNV clustering. We observed both linear (CB1003) and branched (CB1002) hierarchies in high-risk patients (Fig. 4a), supporting variable subclonal diversification during high-risk disease. Clonal evolution based on successive mutation acquisition in different spatial and temporal samples from each patient led to between three and eleven genetically distinct subclonal branches and lineages (clone numbers in circles, Fig. 4a). In some cases, samples from a patient (numbers at the lineage ends, Fig. 4) are assigned to different subclonal lineages indicating that the biopsy sample harboured multiple subclones (i.e. samples 14, 15 from CB1002). This clonal intermixing of subpopulations within one sample illustrated the extensive spatial and temporal genetic heterogeneity in our cohort (Fig. 4a). In 7 out of 8 cases for which both SNV and SCNA data were available, tumour phylogenies based on SCNA profiles agreed with the hierarchies inferred from SNVs ( Fig. 4 and Supplementary Fig. 4). In CB1005 the SNV and SCNA topologies did not agree. Here the SCNA phylogeny groups the dominant clones of samples 12 and 17 as well as those of 11, 13 and 14 into a monophyletic clade ( Supplementary Fig. 4). Ancestral reconstruction using MEDICC2 provided additional evidence in support of this topology, including a gain on chromosome 6q shared between 12 and 17 and an additional gain on 7q shared between 11, 13 and 14. The analysis also revealed multiple parallel events, including MSAI on chromosomes 1q and 18p and 18q. We next compared samples from primary and metastatic tumour sites. CB1001 is a patient with high-risk disease that includes a clonal MYCN amplification, and for whom 5 SCNA profiles were available comprising 4 regions from the primary tumour and 1 from a metastatic renal infiltration. This metastatic infiltration probably occurred not by hematogenous dissemination but by aggressive outgrowth of the primary tumour. Limited heterogeneity at both SCNA and SNV levels occurred in the CB1001 primary tumour sites, and most alterations were clonal (Fig. 4b). Both SCNA and SNV data indicated branching of the metastatic clone prior to the emergence of the clones in the primary (Fig. 4b, c and Supplementary Fig. 2a). This "early" branching is clearly evidenced by parallel evolution between the primary and metastatic clades: chromosome 9 harboured an MSAI event in which a single additional copy of 9q was gained in both primary and metastatic samples, but the copies originated from different parental chromosomes (haplotypes) in the metastatic infiltration and primary tumour (Fig. 4e). In addition, all biopsies from the primary tumour also harboured añ 2.5 mb focal LOH deletion event on chromosome 6q, which was absent in the metastatic infiltration (Fig. 4d). Since lost genetic material cannot be regained, this clearly indicates the branching of the metastatic clone before the diversification of the primary samples. The tree topology additionally showed high robustness as measured by jackknife resampling of the underlying SCNA data (Fig. 4c). A similar pattern emerged in patient CB1003, from whom 10 samples were available (5 biopsies at diagnosis, 4 from the resected tumour and 1 synchronously resected lymph node metastasis). SCNA phylogeny clearly separated the two time points (Fig. 4g) and showed a dynamic acquisition of additional MYCN copies during disease progression ( Fig. 4h and Supplementary Fig. 5). Multiple single-copy gain events occurred on 2p24 in the primary tumour at diagnosis that developed into a MYCN amplification (16 copies) at the time of tumour resection ( Fig. 4f, g, h and Supplementary Fig. 5). An intermediate level of MYCN gain, estimated at 6 copies, was detected in the metastasis (Fig. 4h), indicating the initial 2p24 gain occurred early in tumour evolution, followed by additional 2p24 gains in the primary tumour and the metastasis, likely due to increased accumulation of MYCN in circular extrachromosomal DNA elements, as recently described [39][40][41] . Interestingly, chromosome 5q contained a singlecopy gain event clonal in all samples but the metastasis. Two hypotheses are thus supported by the data: (i) emergence of the metastasis before diversification of all primary samples followed by parallel gains on 2p24, or (ii) branching of the metastasis after chr1  chr2  chr3  chr4  chr5  chr6  chr7  chr8  chr9  chr10  chr11  chr12  chr13  chr14  chr15  chr16  chr17  chr18  chr19  chr20  chr21  chr22  chr4    the emergence of the clone in the biopsy samples but before diversification of the resection samples (Fig. 4) implicating a secondary loss of 5q in the metastatic branch. Tree reconstruction with MEDICC2 inferred option (ii) as the more parsimonious (Fig. 4) and jackknife robustness analysis of the tree indicates a 62% support for the metastatic split between biopsy and resection samples. In any case, all samples obtained from the resected tumour showed clear allelic imbalance and a single-copy gain of chromosome arm 9q, which was absent in both the diagnostic biopsies and the metastasis resected synchroneously with the b    (Fig. 4i), confirming additional events and ongoing evolution after biopsy and metastatic dissemination. Thus, SCNA analysis indicates high genomic divergence between primary tumour and metastases in neuroblastoma suggesting that metastases spread early in disease evolution from the primary tumour site, although additional studies are required to further substantiate this.
Gene expression is temporally but not spatially heterogeneous in neuroblastomas. To detect phenotypic heterogeneity in spatially and temporally distinct neuroblastoma samples, highcoverage paired-end messenger RNA sequencing was performed on 48 of 51 samples and gene expression was quantified (Figs. 1b, 5a and Supplementary Fig. 6). Genes known to be specifically expressed in neuroblastoma, such as HAND2, PHOX2B, and DBH, were homogeneously expressed across all samples, indicating that the proportion of mRNA derived from neuroblastoma cells was equally distributed between each sample. RNA expression analysis confirmed the impact of the most frequent copy-number changes over time in the cohort, i.e. the MYCN amplification, and revealed a clear correlation between MYCN copies and gene expression of MYCN and of its target gene TERT. Other SCNAs also clearly correlated with changes in gene expression, such as the gain in chr 7, which was detected in 8 out of 9 patients and which led to an increase of expression of the IL6 gene located on chr. 7p15.3 ( Supplementary Fig. 6a). To further address the consequences of temporal ITH of SCNAs, we performed Gene-Set Enrichment Analysis (GSEA) of known downstream targets of MYCN 42 (Supplementary Data 2) for patients CB1003 and CB1008. In CB1003, we found MYCNregulated genes enriched towards the tumour resection in line with the acquisition of the MYCN amplification at that time point. In CB1008, MYCN-regulated genes were significantly enriched towards diagnosis, in agreement with the decrease of MYCN copy number during tumour progression ( Supplementary  Fig. 6b). Based on bulk RNA-seq data, temporal expression changes between diagnostic versus tumour resection did not include a significant shift in differentiation state-related signatures i.e. mesenchymal versus adrenergic signatures 43 (Wilcoxon test p-values > 0.05) (Supplementary Fig. 6c). A principal component analysis (PCA) utilising expression levels of the most variable genes confirmed that interpatient differences dominated over intratumour heterogeneity in gene expression (Supplementary Fig. 6g and Supplementary Data 4). A different picture emerged when PCA was restricted to genes previously used for expression-based classification between high-and low-risk neuroblastomas 44 (Supplementary Data 3). Low-and high-risk tumours in our cohort were clearly separated, as expected, but two sets of samples collected at different time points in three patients were also clearly separated (Fig. 5a). In patient CB1003, samples from the diagnostic biopsy clustered with low-risk Fig. 4 Primary and metastatic tumour evolution inferred by genetic analysis of clonality. a Single-nucleotide variant-(SNV)-based phylogenetic trees depict variable clonal diversification. Branch lengths correspond to the total numbers of SNVs. Branch colours represent different subclones. Numbers at the end of each branch indicate samples containing all clones (numbers in circles) of the branch and its ancestor clones. Samples were collected at diagnosis (orange), at tumour resection (red) and at diagnosis of relapse (turquoise). Metastatic samples are indicated by asterisks. A deeper data showcase is provided for patient CB1001, showing SNV-(b) and somatic copy-number alteration-(SCNA)-based trees (c). d SCNA plots from sample 15 (rainbow) showing a loss of heterozygosity (LOH) event on chromosome 6q that is absent in the metastastic sample 11 (green). e SCNA plots from the metastatic sample 11 (green shaded) and primary tumour sample 15 (rainbow) showing an mirrored subclonal allelic imbalances (MSAI) event on chromosome 9. SNV tree (f) compared to SCNA phylogeny (g) in CB1003. Time point or spatial position is indicated by shading, none (metastasis), grey (biopsy at diagnosis) and green (tumour resection). h Dynamic acquisition of additional copies of the MYCN locus on chromosome 2p. i Monoallelic gain on chromosome 9q in the tumour resection samples. (green) but not in the metastatic clone.  tumours, while the samples from the resected tumour clustered with high-risk samples, suggesting therapy-induced changes in the expression of risk-associated genes (Fig. 5a). Our data demonstrate that gene expression programmes in neuroblastomas clearly follow temporal heterogeneity in SCNA and can change significantly under therapy.
Temporal and spatial genetic neuroblastoma heterogeneity in a patient affects clinical assessment of ultra-high risk and second-line therapy options. We and others have recently described molecular determinants of clinical UHR stratification 3 . This UHR stratification is based on (i) the presence of telomere maintenance via TERT rearrangement, MYCN-amplification, high TERT expression, or ALT and (ii) mutations in TP53 and RAS/MAPK pathways. We observed high spatial intratumour genetic heterogeneity and subclonal evolution in our patient cohort that affected these UHR determinants (Figs. 1-4) and hypothesise that clinical risk stratification may be affected by this genetic heterogeneity. To test this, we performed molecular-based risk stratification on each of the neuroblastoma biopsies separately using the information on TERT mRNA expression, MYCN copy number and TP53/RAS/MAPK pathway mutations from WES data. Based on these molecular risk factors, almost one third (29%) of biopsies (15 of 51) were classified as UHR (Fig. 5b). In line with our hypothesis, not all biopsies from the same patients were classified homogenously. Whereas biopsies collected at diagnosis in patient CB1003 were not classified as UHR, all biopsies from the primary tumour resected after induction therapy were classified as UHR (Fig. 5b). This suggests the clonal selection of UHR subclones during therapy. In contrast, diagnostic biopsies from patient CB1008 (currently in complete remission for >3 years) were classified as UHR, while biopsies after chemotherapy at tumour resection were not classified UHR. Most prominently, this occurred due to longitudinal changes in MYCN copy number (CB1003) or ALK mutation as well as TERT gene expression (CB1008). In total, 2 out of 10 patients showed a heterogeneous UHR profile. This is important because the neuroblastoma community is currently discussing to base alternative experimental therapy decisions on the presence of druggable mutations and UHR classification at diagnosis in the future, but our data clearly show that a subset of UHR patients still benefits from standard high-risk treatment. We conclude that clonal evolution in neuroblastoma can influence the molecular profile within the tumour, across disease sites and/or during the course of disease in potentially clinically relevant ways. Our findings support the consideration of strategic and methodological changes for target identification and molecular risk classification as well as subsequent treatment decisions in individual patients with neuroblastoma.

Discussion
Intratumour heterogeneity is considered a major cause of cancer progression and treatment resistance 24 . Multi-region sequencing studies have characterised the extent of intratumour heterogeneity in a variety of cancer entities, and linked the evolutionary subtypes with patient outcome 24,35,45,46 . Our present multiregion sequencing study of 10 patients with neuroblastoma shows extensive spatial and temporal intratumour heterogeneity that also affected genes encoding actionable targets, such as ALK and FGFR1, thus, harbouring potential implications for therapeutic target selection. Analysis of phased SCNAs revealed early emergence of the metastatic lineage and ongoing chromosomal instability, including the appearance of a MYCN amplification later in the course of disease. This process of clonal evolution led to clinically relevant changes in molecular risk stratification in some patients. Our in-depth analysis of SNVs and SCNAs assessed a significantly lower number of SNVs, and practically no mutational hits in cancer-related genes in low-and intermediate-risk neuroblastomas. Overall, we detected a low number of non-silent SNVs (range 0-60, median 15), which is in accordance with previous reports 3,4,26 . Less than half of the detected SNVs were clonal, based on our criteria for clonal appearance (present in all samples from a patient with VAF > 10%). To gain insight into the biological and clinical relevance of intratumour heterogeneity, we focused our analysis on non-silent SNVs in cancer-related genes. Non-silent SNVs (2-6 per tumour) in cancer-related genes (HRAS, ALK and ATM) occurred exclusively in high-risk neuroblastomas and were clonal in a subset of patients but also spatially or temporally heterogeneous in others.
Efficacy of a targeted therapy that might be included into second-line treatment for a patient with neuroblastoma depends on precise and reliable diagnostics. We detected spatially heterogeneous SNVs in genes encoding actionable targets (FGFR1 N577K, BRAF V600E and ALK R1275Q) or temporally heterogeneous SNVs in CCND3 or NOTCH2. Spatially heterogeneous mutations suggest that current diagnostic biopsies may miss or overestimate therapeutically actionable alterations. More specifically, application of current target prioritisation algorithms might have obtained different results when analysing different tumour locations, such as relapse tumour samples from our patient CB1002, which harboured different private druggable mutations. Most importantly, even the presence of activating ALK mutations, the most frequent actionable alteration in neuroblastoma 9,35,47 for which targeted therapies are currently being clinically tested [48][49][50][51][52][53] , was spatially and temporally heterogeneous. Such heterogeneity could impact the clinical interpretation of diagnostic biopsies and may obscure reasons for nonresponse to ALK inhibitor therapy in patients in whom the majority of the tumour cells do not contain activating ALK mutations 54 . Our finding that druggable mutations are spatially and temporally heterogeneous demonstrate that clonal neuroblastoma evolution has potential implications for targeted therapy decisions. We postulate that for a subset of patients, bulk sequencing at a single tumour location might sometimes lead to an imprecise molecular diagnosis, potentially leading to the wrong choice for targeted therapy. In the future, the technological advances of single-cell technologies and liquid biopsies might contribute to solve this problem, as repeated multi-region sequencing of tumour tissue is not feasible in the clinical routine.
Allele-specific SCNA profiling and multi-region phasing based on WES data revealed substantial heterogeneity in neuroblastoma chromosomal instability in patients in our cohort. Subclonal SCNAs were observed in every patient and 3 out of 9 patients showed evidence of parallel evolution through detection of MSAI events. Most SCNAs were numerical chromosomal aberrations in the two patients with low-and intermediate-risk disease, correlating with previous reports 37, 38 . In most patients (5 of 9), the majority of the SCNAs were clonal, demonstrating the early appearance of SCNAs in the course of disease. Interestingly, in one patient with mostly clonal SCNAs, extensive intratumour SNV heterogeneity was detected, suggesting that in this patient SCNAs were the tumour-driving force followed by subclonal mutational diversification. In the remaining patients, variable degrees of ongoing clonal evolution and chromosomal instability occurred. Strong temporal heterogeneity was observed in one patient (CB1008) with many subclonal SCNAs and ongoing chromosomal instability between samples collected at initial diagnosis or after therapy, but stable detection of a MYCN amplification and loss in 1p36 both at neuroblastoma diagnosis and following therapy. The heterogeneous profile of cancerrelated gene mutations in the same patient showed the near elimination of an initial ALK R1275Q mutated clone after therapy, further corroborating the extensive temporal heterogeneity. To better understand neuroblastoma evolutionary diversification in each patient, clonal hierarchies were reconstructed from SNVs in biopsies from different time points and spatial distributions of the disease. We observed both linear and branched evolutionary trees, creating three to eleven distinct lineages per patient, and indicating variable clonal diversification across the cohort. Previous studies 13,14 defined four evolutionary trajectories over different anatomic areas for childhood cancers, among them 23 neuroblastomas including two trajectories correlating with poor outcome. We frequently observed simultaneous subclonal SCNAs, clonal segmental SCNAs and clonal or subclonal cancer gene mutations in the same high-risk neuroblastoma, tracing a clonal sweep as was previously described 14 . Due to our small cohort size of 10 patients, correlation with outcome was not conclusive. In contrast to prior reports 14 , our allele-specific SCNA analysis detected clinically relevant 1p and 11q losses as mostly clonal while in one of three cases with MYCN amplifications, the amplification was not clonally persistent but occurred later during the course of disease. Taken together, the high degree of both spatial and temporal neuroblastoma heterogeneity and the multiple evolutionary trajectories observed to date in individual patients with neuroblastoma indicate that many complex routes culminate in severe neuroblastoma courses and will require detailed and precise molecular analyses in further patients to completely understand causal or key points in the process and to base second-line treatment selection on reliable molecular data for every individual patient.
To uncover the evolutionary processes of metastasis, we analysed synchronously resected primary tumour samples and metastatic infiltrations. Both SCNA and SNV data indicated branching of the metastatic clone during disease progression and divergent evolutionary trajectories, with clear evidence of clonal aberrations within the metastasis that were absent in the primary tumour. This was verified by MSAI events on chromosome 9q, allele-specific LOH or MYCN amplifications that differed between the metastases and matched primary tumours. In contrast, Karlsson et al. 14 found that synchronous metastasis mirrored the primary tumour. Branching of the metastatic clone and substantial heterogeneity between a metastasis and primary tumour may explain why neuroblastoma relapses often occur at metastatic sites while primary tumours are successfully eradicated by surgery and first-line chemotherapy. This is also supported by the high number of de novo mutations in relapse cases in our cohort, which indicate a major genotypic shift during (metastatic) recurrence, in line with the previous reports 14,26,29 . Inferring evolutionary trajectories from SNV and SCNA data reveals variable degrees of clonal diversification and an early emergence of metastatic clones in neuroblastoma.
An UHR group has recently been described based on molecular determinants 3,26 for a new concept of clinical risk stratification of patients with neuroblastomas. Notably, UHR molecular determinants were not homogenous throughout biopsies from the same patient in our cohort. This was especially true for patients CB1003 and CB1008, whose samples showed extensive temporal SNV and SCNA heterogeneity. Distinguishing UHR from HR disease may have important clinical consequences, since UHR patients have poor outcomes and may benefit from an earlier switch to alternative treatment approaches 55 . Based on our observation of spatial and temporal genetic heterogeneity of UHR and HR determinants in individual neuroblastoma cases, we believe that precise risk stratification can only be obtained through molecular characterisation of neuroblastomas at diagnosis, resection and recurrence. Our in-depth analysis of neuroblastoma intratumour heterogeneity reveals that extensive genetic heterogeneity and subclonal diversification evolves under therapy and may have important implications for the clinical interpretation of molecular diagnostic results and for the selection of appropriate second-line treatment approaches.

Methods
Sample collection. Tumour specimens were collected from 10 patients enrolled between 2014 and 2018 in the German Society for Pediatric Oncology and Hematology (GPOH) NB2004 trial or NB registry 2016, and treated at the Charité University Medicine Berlin (Germany) according to the trial/registry protocols. All patients and/or their guardians gave written consent for the use of biosamples and clinical data for research in accordance with the local ethics review board of the medical faculty, University of Cologne as the trial sponsor of the NB2004 and the NB registry 2016 (https://clinicaltrials.gov/NCT03042429). Samples were collected by open surgical biopsy either at diagnosis, at tumour resection after 4-6 cycles of chemotherapy (according to neuroblastoma therapy regimen) or at diagnosis of relapse. Fresh samples were immediately snap-frozen in liquid nitrogen and stored at −80°C. Portions of tumour material were formalin-fixed and paraffinembedded (FFPE) in parallel for diagnostics and preservation in the pathology unit. Two to 30 biopsies were taken from geographically separate areas of the tumour body with a minimal distance of 10 mm from each single tumour. From these biopsies, only tumour regions with high tumour cell content (>60% for WES and RNA sequencing, 10% for targeted re-sequencing) were included. A pathologist confirmed the diagnosis and assigned tumour regions for macrodissection with a high content of vital tumour cells on sequential hematoxylin and eosin-stained sections. Peripheral blood collected from each patient was used as a matched germline control for tumour samples. DNA was prepared using the Qiagen DNA Mini kit according to the manufacturer's instructions.
Compiling a list of relevant cancer genes. A list of genes potentially involved in neuroblastoma development was compiled by combining lists of potential cancer driver genes from several sources: the COSMIC database 27 (version 80), a recent review of the mutational landscape in neuroblastoma 10 , pan-cancer studies of paediatric tumours 8,16 and a pan-cancer study of adult cancers 28 . The final list consisted of 801 genes.
Whole-exome sequencing (WES). The SureSelect Human All Exon V6 kit (Agilent) was used to prepare libraries enriched with exonic sequences. The libraries were further prepared for sequencing using the Illumina TruSeq Exome Kit and sequenced on Illumina HiSeq 2500 and Illumina NextSeq sequencers.
Raw data generation and quality control. Read sequences and base quality scores were de-multiplexed and stored in Fastq format using the Illumina bcl2fastq (2.16) software. Adapter remnants and low-quality read ends were trimmed off using custom scripts. The quality of sequence reads was assessed using FastQC (0.11.5) software. All libraries were deemed of high quality and usable for follow-up analyses.
Alignment and coverage evaluation. Reads were aligned to the human genome (assembly GRCh38) using bwa mem (0.7.17) software 56 . Duplicate read alignments were marked using samblaster (0.1.24), and resulting Bam files were stored for later analyses. Overlap with exon target regions was assessed using the BedTools (2.27.0) suite of tools and the target region file provided with the SureSelect Human All Exon V6 kit (Agilent). Mean coverage was computed from these results as the average number of aligned reads per base in the specified exonic target region as were the exact proportions of targeted bases covered by each number of aligned reads. Median coverage was 346× (range: 163×-636×). Aligned reads were further filtered using sambamba (0.6.8) software to remove the secondary and failed alignments and PCR-duplicate reads.
Somatic variant calling. We determined somatic single-nucleotide variants (SNVs) and short insertions/deletions (InDels) using EBCall (2) software 57 . Hereby, each tumour sample is first compared to the matching control sample, and SNVs and InDels supported by a sufficient fraction (≥10%) of aligned reads in the tumour sample and not exceeding a specific fraction of reads in the control sample (<5%) are selected as candidate variants. Based on the complete set of control samples from all patients, a probabilistic model for base mismatches due to technical reasons was constructed, and only candidate variants deemed highly unlikely (p < 0.001) under this model, and thus, more likely to have arisen from biology rather than technical artifacts are considered as somatic variants. A SNV was considered to be clonal if present in all tumour cells with variant allele frequency greater than 10%, based on the assumption that all somatic mutations were heterozygous and tumour cell content in the sample was~60%. The threshold was adjusted accordingly for SNVs in regions with copy-number gains and losses.
Annotation of somatic variants. The effect that each putative somatic variant would have on genes/proteins was annotated using SnpEff (4.3t) software 58 and the knownGene gene annotation from the UCSC genome browser for the hg38 assembly. For each variant, genes potentially affected by the mutation including effects on the protein level were stored together with information including whether the variant had been previously detected and deposited in a variant database such as dbSNP. A smaller list of variants predicted to affect well-known cancer genes was generated that only included non-silent variants in genes contained in the list of cancer-relevant genes described above.
ALK mutation detection using droplet digital PCR (ddPCR). The ALK R1275Q (3824, G>A) hotspot mutation and corresponding wild-type ALK sequence was detected using the QX200 Droplet Digital PCR System (Bio-Rad Laboratories, Munich, Germany). The AluI restriction enzyme (5 U, New England Biolabs, Frankfurt/Main, Germany) was added to each ddPCR reaction to fragment the tumour DNA 59 . Reaction mixtures for duplexed TaqMan ddPCR used 2× ddPCR Supermix for Probes, no dUTP (Bio-Rad Laboratories) and primers/probes adapted from Combaret et al. 60 , and have been previously described 61 . Primers and probes had the following sequences: ALK 1275 -for, 5′-GTCCAGGCCCTGGAAGAG-3′; ALK 1275 -rev, 5′-GGGGTGAGGCAGTCTTTACTC-3′; ALK R1275Q probe, FAM, 5′-FAM-TTCGGGATGGCCCAAGACAT-BHQ1-3′ and ALK 1275 probe, HEX, 5′-HEX-TTCGGGATGGCCCGAGACAT-BHQ1-3′. Manufacturer recommendations for the QX200 Droplet Generator were followed to create droplets, and reaction endpoints were assessed in the QX200 ddPCR Droplet Reader. The fraction of mutant allele was identified using QuantaSoft Analysis software, version 1.7.4.0917. Appropriate non-template, positive and negative controls were included in each ddPCR assay run for software-based generation of specific thresholds. The model of Armbruster and Pry 62 was followed to calculate (using Bio-Rad look-up tables) the false-positive rate (FPR) and detection limit for point mutation analysis. FPR calculation for the point mutation has been previously described 61 , and was based in principle on the minimally required mutant target molecule concentration (copies/µl) and number of false-positive droplets. When the number of droplets detecting the mutation and the mutant target molecule concentration (copies/µl) were both above the thresholds that were set, a sample was scored as positive. Detection limit was determined as previously described 61 .
Determining allele-specific somatic copy-number alterations. WES read counts from all neuroblastoma and matched germline (blood) samples were aggregated at 20,180,418 common germline variants from the NCBI dbSNP database using samtools mpileup (1.9). Variants were removed that (i) were in highly variable regions of the major histocompatibility complex (MHC, chromosome 6: 28,510,120-33,480,577), (ii) were in positions characterised as not uniquely mappable by the UCSC Genome Browser annotation tracks (hoffmannMappability, k100.Unique.Mappability.bb) 63 and (iii) with <20 reads across both neuroblastoma/germline samples. Average genome-wide GC content was calculated in 50 bp windows using sequenza-utils gc_wiggle 64 . Read counts and GC content were forwarded to Sequenza (v3.0.0) 64 for segmentation and calling of allele-specific SCNAs in each biosample. Segmentations were combined for samples from each patient using a joint consistent segmentation approach as previously described 23,25 . The jointly segmented allele-specific SCNAs then underwent multi-region reference phasing using refphase 25 to determine haplotype-specific SCNAs and detect MSAI events indicative of parallel evolution. All haplotype-specific SCNA profiles were manually reviewed and individual mis-segmentations and erroneous sample purity and ploidy estimates were adjusted accordingly. Haplotype-specific SCNA profiles were then used to infer tumour evolution and reconstruct ancestral genomes using MEDICC2 31,32 .
Targeted sequencing. A set of 1479 SNVs were selected from exome sequencing data for targeted sequencing at higher read coverage and in additional samples. This set included all SNVs in cancer-related genes (from our constructed list, Supplementary Data 1) and SNVs present in all patients to ensure that each patient had a minimum number of specific SNVs represented in the targeted sequencing panel. A custom targeted kit (Agilent) was used to enrich DNA from the genomic positions of these 1479 sites. DNA was isolated from 150 samples from the 10 patients, subjected to the targeted DNA enrichment, and sequencing libraries were generated for Illumina HiSeq and NextSeq sequencing (150 nt single-end reads). Reads were aligned to the human genome (assembly GRCh38) using bwa mem software 56 . Variant allele frequencies for each of the 1479 SNVs in each sample were computed using the Rsamtools Bioconductor/R package.
Reconstruction of clonality trees. The sciClone software (version 1.1.) 65 was used to infer clonal composition and possible neuroblastoma phylogeny in each patient. This method uses a Bayesian mixture modelling approach to stratify variants into clusters based on their copy number and variant allele frequency (VAF) in each sample. Sample tumour purity is estimated from the VAF distribution in copy-number neutral regions and used to scale all VAFs. The software setting, binomial-mixture model was used to cluster VAFs, with an upper limit set to 20 potential clusters and the default settings for other parameters.
VAFs of exactly 0 were replaced by a very low number (10 −10 ) prior to running the algorithm, since the clustering algorithm did not converge in extreme cases in which too many VAF entries were exactly 0. Variant clusters from sciClone were assigned to clones, and parental relationships between clones were inferred as follows. Clusters with high cellular fraction (>90%) in all samples were merged into the founding clone to make up the branch of the tree. Other clusters were grouped according to their cellular fraction in subsets of samples and assigned as clones to the branches of the tree, with the requirement that the ancestor clone was present in the same and potentially more samples. Clusters only present in a single sample were merged into clones that made up the leaves of the tree. Clones with a cellular fraction of <5% in a sample were not assigned to the respective sample. Samples were annotated in the tree at the relevant clone furthest away from the tree stem, indicating that the sample contained the respective clone and all of its ancestor clones. Clones were denoted as numbers, branch lengths corresponded to the total numbers of SNVs in each clone and cancer-relevant genes contained in each clone were annotated at the branches. For patients CB1005 and CB1006, the total numbers of SNVs with median VAFs >0 were too low to infer meaningful clonality trees.
RNA sequencing. Sufficiently high-quality RNA was available for 48 of 51 samples used for exome sequencing. RNA was prepared from areas with high tumour cell content macrodissected from cryosections. RNA was isolated using the NucleoSpin RNA kit (Machery&Nagel) according to the manufacturer's protocol. Sequencing libraries were prepared at the Berlin Institute of Health Sequencing Core using the standard TrueSeq stranded poly-A enrichment protocol. Libraries were sequenced on an Illumina HiSeq sequencer as 2 × 75nt paired-end reads. After de-multiplexing, raw read sequences in Fastq format were aligned to the human genome (assembly hg38) and Gencode transcript annotation (version 26) 66 using STAR alignment software 67 . Gene expression was determined as read counts per gene using the featureCounts tool in subread software. All further analyses of read-count data were performed using the DESeq2 Bioconductor/R package 68 . The list of the 104 genes separating neuroblastoma risk groups was generated by analysing publicly available data from 498 primary neuroblastomas 44 for differential gene expression between high-and low-risk tumours using the limma R package. Genes with an absolute log 2 -fold change >1 and a gene-wise false-discovery rate (q) < 0.1 were considered to be significant in this analysis and retained in the final list. Expression of these genes in our patient samples was used for the second principal component analysis.
Detection of MYCN copy number by SNP array and by FISH. For validation of WES-based estimation of MYCN copy numbers, a SNP array (Cytoscan HD, ThermoFisher Scientific, #901835) was performed for selected patient samples as in ref. 15 . CEL files with intensity probe signals were analysed using the Chromosome Analysis Suite (ChAS) software version 4.1.0.90 (r29400) and converted to CYCHP files. Copy numbers are visualised using the allele difference plot, weighted log2 ratio and the smoothed signal. The smoothed signal is used to directly estimate copy numbers.
MYCN copy number was detected on FFPE tissue on single-cell level using FISH analysis. The hybridisation probe XL MYCN amp (MetaSystems) consists of a green-labelled probe hybridising to the MYCN gene region at 2p24 and an orangelabelled probe hybridising to the NMI gene region at 2q23 as reference. Data interpretation was done according to INRG Biology guidelines 69 .
Reporting summary. Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Data availability
Raw sequencing data used in this study are available in the European Genome phenome archive (www.ebi.ac.uk/ega/) under the accession number EGAS00001004735 under the following dataset IDs: EGAD00001008020 (WES), EGAD00001008156 (targeted resequencing), EGAD00001008133 (RNA-Seq) and EGAD00010002225 (SNP array). Download of all datasets are available by contacting the study affiliated Data Access Committee (DAC) via EGA. For WES analysis, NCBI dbSNP (https:// www.ncbi.nlm.nih.gov/snp) and UCSC Genome Browser (https://genome.ucsc.edu/) were used to aggregate somatic mutations. For defining cancer-related genes, public data was downloaded from COSMIC database (https://cancer.sanger.ac.uk). For MES/ADRN differentiation subsets, analyses were done based on datasets downloaded from GEO (www.ncbi.nlm.nih.gov/geo/) under accession number GSE90805. For analysis of differential gene expression between high-and low-risk tumours, datasets were downloaded from GEO with accession GSE49711.

Code availability
Source code in R for running sciClone and for constructing the clonality trees is included as a zipped file in the supplementary online information (Supplementary Software 1).