Spatiotemporal profiling under chemotherapy reveals two contrasting evolutionary patterns
NB is largely a copy-number aberration (CNA) driven disease14-21. We recently reported a software tool (DEVOLUTION) that allows integration of point mutations and CNA data, where the latter is curated according to constraints of chromosomal evolution22. To perform high-resolution mapping of subclone territories across anatomic tumor space, we applied DEVOLUTION to mutational data from archived paraffin embedded formalin-fixed NB tissue (Fig. 1a-e; Online Methods). This allowed a comprehensive comparison between the genomes of NB cells before and after treatment, including both patients with remaining, extensively viable tumor and cases with good chemotherapy response where survivor populations were limited to spatially dispersed tumor islets a few millimeters in size. Areas of @20mm2 containing >30% tumor cells were subjected to whole genome copy number profiling parallel to whole exome sequencing (WES). Mutations detected by WES were validated by targeted resequencing using a custom-made panel based on the WES findings. From a consecutive patient cohort treated at a single pediatric oncology center over a period of 20 years, we identified 12 patients from whom viable tumor tissue was available from ≥2 geographically separated tumor regions from ≥2 time points (Fig. 1f-g; Table S1a). A similar number of primary tumor samples (median 3.0) per patient from before treatment as from after chemotherapy were analyzed with both CNA profile and sequencing; fewer samples were from metastatic sites (median 2.0).
Among the 12 patients, two showed progression of the primary tumor with metastases under treatment (Patients 1-2), nine showed a significant therapy response (40%-97% primary tumor regression including necrosis, fibrosis and extensive maturation; Patients 3-11), while one patient had a single small tumor in the neck and was treated with up-front surgery only (Patient 12). All seven patients whose tumor could be evaluated before and after significant chemotherapy response, showed a strikingly similar ancestral relationship between subpopulations detected before and after treatment (Patients 3-9; Fig. 1h and Fig. S1, Table S2a): some genetic changes are no longer detected while other mutations emerge, as populations present prior to chemotherapy were replaced by others, typically arising from a different most recent common ancestor (MRCA). In the phylogenetic trees, this corresponded to one set of collateral branches being replaced by other, parallel branches and hence we denoted this pattern collateral clonal replacement (CCR). CCR occurred irrespective of whether the pre-treatment sample was obtained from the primary tumor (Patients 3-7) or from a metastasis sampled at presentation of disease (Patients 8-9). That CCR could take place gradually over time was illustrated by Patient 8 (Fig. 1i). Here, the tumor was biopsied after ~50 days of therapy (due to slow response) and then sampled again at resection after ~500 days under treatment. A comparison of these two time points revealed a gradual shift towards CCR as populations were replaced. Using an index of genomic diversity (IGD) that was independent of both sample number and the total number of detected mutations (see Online Methods), we compared the subclonal heterogeneity among samples obtained before and after chemotherapy (Fig. 1j). In the primary tumor, diversity was found to be lower before than after chemotherapy, even though the burden of copy number aberrations (phylogenetic branch lengths) did not differ between these time points (Fig. 1k). CCR in the primary tumor under therapy was thus associated with the emergence of a more diverse clonal landscape.
CCR was also observed for the majority of subclones in the primary tumor when compared to those in a later metastatic relapse (n=5; Patients 4, 9-12). Phylogenies containing clones detected at metastatic relapse exhibited a higher degree of irregularity (see Online Methods) than those only based on clones detected in the primary tumor and/or up-front progression (Fig. 1l-m) as the lengths of branches to clones detected in relapse samples were overall longer than those in the primary tumor. In most relapsed cases, the lineage of relapsing subclones diverged from that of the primary tumor at an early stage, with CCR encompassing all subclones in the primary versus the relapse. Minor but significant exceptions from this were Patients 4 and 12. In Patient 4 (Fig. S1d), a regional population corresponding the common ancestor of relapsing subclones was detected in the primary tumor after therapy, while in Patient 12 (Fig. S1l) it was detected prior to therapy. Notably, these two primary tumors were the ones with least chemotherapy-induced regression; Patient 4 responded largely by tumor cell ganglionic maturation, while Patient 12 was treated with surgery only. This indicated that inferior killing of cancer cells was associated with an aspect of linear evolution from surviving ancestral clones when/if the patient progressed. Indeed, in the cases with local and metastatic progression under chemotherapy (Patients 1-2), a subclone ancestral to the genomes of metastatic populations was ascertained at diagnosis (Fig. 1n-o), supporting features of linear evolution under inferior therapy response.
Linear evolution under progressive growth of chemotreated patient-derived xenografts
Because the number of patients with progressive disease under therapy was limited, we set up a model system to mimic progressive tumor growth under chemotherapy. We used a well-established NB murine model system of patient derived xenografts (PDX)23. PDXs from a stage III NB (n=6) were treated by cisplatin monotherapy, resulting either in stationary disease (n=3) or progression (n=3) under treatment (Fig. 2a and S2a:I). We also simulated progression after incomplete surgery by partial excision of non-chemotherapy exposed PDXs followed by regrowth to the same volume (n=4). Control PDXs (n=4; saline treatment), grown to the same volume, were used as a reference. Whole exome sequencing revealed clonal evolution in PDXs mimicking clinical tumors, including parallel evolution of RAS/MAPK mutations in different lineages (Fig. 2b; Table S2c-e). Phylogenetic analysis revealed a pattern of linear evolution with the cell culture from which PDXs were established acting as founder (S2a:II-III). Most PDXs progressing under chemotherapy or after surgery showed expansion of a clone with partial gain of 1p, all of 1q and an extra copy 17q, leading up to clonal sweeps (Fig. 2c and S2b-d). They also had a significantly higher CNA burden and more private mutations than the cisplatin-stationary and the untreated group, in a fashion that was to some degree time-dependent (Fig. 2d and Fig. S2e); while untreated tumors had the same end-size as tumors progressing under treatment or after resection, they reached this size in shorter time and had fewer anomalies. The PDX model thus confirmed that progressive NB growth under therapy was associated with a pattern of linear evolution from a common ancestor, radiating in different direction at different anatomic sites (different PDX animals).
Subclone replacement under effective chemotherapy in patient-derived xenografts
To compare the result above with effective chemotherapy treatment in the same PDX model, xenografted tumor cells were subjected to an intense five-drug chemotherapy protocol mimicking the regimen used to treat high-risk NB patients (rapid COJEC; Fig. 3a)24. Mice were randomized into treatment groups when NB PDX tumors reached approximately 500 mm3. After regression under COJEC to <200 mm3, surgery was performed to remove all visible tumor. The COJEC-treated tumors (9 x 2 samples per tumor) were compared to identically sampled untreated tumors (n=9 x 2). Treated and untreated tumors had some common evolutionary features, including expansion in 17/18 tumors of a subclone signified by the same gains in chromosomes 1 and 17 as in the cisplatin-only monotherapy model detailed above (Fig. 3b-c). Comparing genetic diversity in treated and untreated tumors showed significant differences between the cohorts (Fig. 3d; Fig. S3a-b; Table S2f-g). Despite treated tumors being on average 10 times smaller in volume than those untreated, all treated tumors (9/9) harbored subclones that were unique to one of the two sampled regions while this was rare (2/9) in untreated tumors (P=0.0023; two-tailed Fisher’s exact test). In line with this, the number of subclonal alterations differentiating the two sampled regions as well as the total number of region-unique subclonal alterations was higher in treated than in untreated tumors, while there was little difference in the total number of aberrations between the COJEC-treated and saline-treated cohort (Fig. 3d, right panel). The scenario was thus similar to the findings in the clinical cohort, where there was no overall increase in the total number of aberrations but still an increased genetic intratumor diversity after treatment.
Because PDX sizes were feared to be too small for sampling prior to treatment without inducing populations bottlenecks, the finding of CCR under therapy from the clinical cohort could not be validated by sampling PDXs before and after therapy. However, COJEC-treated PDXs T1 and T5 relapsed locally, enabling comparison of clonal landscapes at two distinct time points, i.e. just after treatment versus after regrowth (Fig. 3e and Fig. S3c). The two PDX relapses showed a strikingly similar route of evolution, with additional 17q copies, and successive intragenic deletions in the PTPRD and TCF4 tumor suppressors24-26. Both PDXs showed clonal and subclonal aberrations that were unique to each time point, i.e. a pattern of subclonal CCR (Fig. 3f-g). Similar to relapses in the clinical cohort, there was an overall increase in total number of CNAs after treatment in the PDXs, with the most complex subclone being found at relapse. The PDX model thus supported the association between effective chemotherapy response and CCR.
Recapitulation of collateral clonal replacement under chemotherapy in vitro
To test whether the contrasting scenarios of linear evolution under inferior treatment/progression versus CCR under effective treatment could be reproduced also under highly controlled conditions, we used the cisplatin-sensitive cell line IMR-32, derived from a high-risk MNA NB (Fig. S4a). In vitro growth for 42 days (11 passages) without treatment resulted in expansion of the major subclone up to fixation and concomitant loss of a parallel subclone (Fig. 4a, Samples A1-A3; Fig. S4b; Table S2h). Low-dose (0.1 µM) cisplatin treatment led to a preserved clonal landscape for the same period of time (Fig. 4a, Samples B1-B3). In contrast, high-dose (1.0 µM) cisplatin exposure with a 94% reduction in population size, followed regrowth for 42 days, resulted in CCR, with elimination of a major subclone and the expansion of clones with novel CNAs in 2/3 parallel experiments (Fig. 4a, Samples C1-C3). Reseeding the same number of IMR-32 cells (6% of original) as the viable population remaining after cisplatin treatment, followed by expansion to the same final population size as the high-dose cisplatin experiment did not change the clonal landscape (Fig. 4b and S4c), nor was the landscape largely changed by a series of three mechanical reseeding bottlenecks of around 50% (Fig. 4c and S4d), while it again exhibited CCR at cisplatin-bottlenecks each killing around 50% of cells. These results again corroborated that effective chemotherapy treatment is linked to a pattern of CCR and indicated that this effect is not merely contingent on the population bottlenecks to which effective treatment is associated (Fig. S4e).
Phylogenetic branches correspond to spatial territories
The fact that tumors in the clinical cohort derived from well-annotated paraffin blocks made it possible to trace back the approximate spatial location of each clone in each tumor. In total 35/122 subclones spanned two or more biopsy locations allowing positioning of shared clones to neighboring areas, while subclones confined to a single biopsy location could not be mapped spatially beyond that location. Applying these principles to the 12 patients in the clinical cohort, we traced phylogenetic branches back to anatomic territories (Fig. 5 and Fig. S5). In untreated diagnostic biopsies, these territories displayed a dense patchwork (median 0.36 subclones per mm2, median distance between the center points of clonal territories of 1.7 mm), with no obvious anatomic barriers between subclone territories (Fig. 5a-d, i, k; Fig. S5a-b, f, l-m). A similar situation with densely variegated clones were observed in biopsies from metastatic relapses (Fig. S5f, i-k). Subclones in the primary tumor after chemotherapy were more dispersed (median 0.024 subclones per mm2, median distance between the centerpoints of clonal territories of 6.5 mm; p=0.0087 pre versus post; p=0.048 post vs. relapse; Mann-Whitney U test, two-tailed; Fig. S5l). Mutations detected after chemotherapy were typically fixed at a clonal level (>90% of all tumor cells) in vastly dispersed islands of surviving tumor cells surrounded by necrosis, hemorrhage and other reactive changes (Fig. S5a lower, b right, e, g, j, m). These distinct tumor cell populations typically always originated from multiple phylogenetic branches, indicating a high degree of taxonomic diversity well in accordance to increased IGD after chemotherapy. Subclones identified at diagnosis were not re-detected after treatment in any of the viable primary tumor areas left after treatment (all areas sampled in Patients 3, 5-8). This indicated that the scenario of CCR under therapy was not only a result of geographically dispersed sampling but resulted from a true disappearance of dominant subclones under treatment.
Clonal replacement under treatment uncovers MYCN amplicon diversity
The extensive genetic diversity across territories described above, suggested that branching evolution emerged early in NB pathogenesis. To explore genetic diversity in early phases of disease we performed an in-depth analysis (down to 0.05 Mb) of copy number variation around the MYCN oncogene in 2p24 in MNA cases, as gain of MYCN gene is thought to be one of the earliest mutations in a large subset of NBs27,28. In the six MNA cases, an amplified state of (>5 copies) of MYCN was present across all biopsied areas, However, spatiotemporal heterogeneity in the architecture of amplicons in or around MYCN was also found in all of them (Fig. S6). When chromosomal breakpoints of amplicon cassettes were used to infer ancestral relationships, all cases had at least one shared breakpoint across all samples, indicating that one or a few breakage(s) in 2p24 first arose to create patient-specific amplicon cassettes, followed by further modification later on. In total, only ~30% of samples (4/14) obtained from the primary tumor before treatment or from metastatic sites growing off-treatment after remission exhibited variant (non-stem) breakpoints around MYCN. In contrast, 76% of samples obtained from chemotherapy treated primary tumors exhibited variant breakpoints (13/17; p=0.00122; Fisher’s exact test). This was in accordance with the scenario of increased genomic diversity after treatment found across the entire genome and supported that phylogenetic branching was an early event in NB pathogenesis.
Phylogenetic branching is an early and stable feature of NB
To further monitor early branching evolution with a focus on the MYCN amplicon, we performed low-pass single cell whole genome sequencing (scWGS) of single biopsies from nine NB primary tumors (three MNA tumors), resulting in 505 single cell genomes (Patients S1-S9 in Table S1B). A unifying feature among all tumors irrespective of risk group was extensive phylogenetic branching (Fig. 6 and Fig. S7). Branches typically led up to subsets of cell having identical CNA profiles, with the exception of one case (Patient S7), which had a different CNA profile in every cell (Fig. S7d). In MNA patients S1 and S2, the earliest branching event consisted of structural variation in the MYCN amplicon cassette (Fig. 6a-e and S7a). Both these cases contained cells having low-level MYCN gain as the first event in their phylogenetic trees (in S2 co-gained with MAML3), corresponding to an ancestral population. Clonal evolution from this stage occurred through structural aberrations characteristic of NB, including deletion of 1p, gain of 17q (both S1 and S2) as well as deletion in 11q (S2 only; Fig. S7a:III). To corroborate that MNA could precede 17q gain (the most common CNA in NB), we further evaluated the concurrence of these imbalances by fluorescence in situ hybridization in a subset of MNA NBs (Table S1c). The presence of tumor cells with MYCN copy number gain while still not having acquired extra copies of 17q was ascertained in 5/6 cases (Fig. 6b and Fig. S7f), the exception being S3 in which neither scWGS could identify such cells (Fig. S7c:I).
Interestingly, in both Patients S1 and S2, clonal evolution through 1p loss/17q gain appeared to be a requirement for further MYCN amplification, with significantly higher amplicon numbers in cells having acquired 1p loss/17q gain than in cells with only MNA (Fig. 6f-g and Fig. S7a:V). Further branching into distinct subpopulations continued after 1p loss/17q gain in all three MNA cases. To test whether the tendency of branching evolution was a stable feature of NB cells, we turned to three PDX models from MNA NBs23,29. We used scWGS to monitor evolution from PDX in vivo generations 1-7, and further through transfer into tumor organoids analyzed at in vitro passages 7 and 20. This showed that phylogenetic branching was a consistent feature across time in vivo, including passing through single-cell bottlenecks (Fig. 6h-j).
Timing and fitness of chromosomal changes impacts NB genome profiles
The types of CNA causing branching evolution in the scWGS data varied across tumors, but a common feature (9/9 cases) was parallel aberrations affecting the same chromosome in different branches of the same phylogeny. For example, in Patient S4, the first branches emerged through distinct parallel deletions in 3p (Fig. S7b:I-III), while in Patient S8 chromosome 6 was the substrate of the three earliest branching events, followed by parallel rearrangements of chromosome 17 and subsequently of chromosome 1 (Fig. S7b:IV-V). Similar parallel evolution was detected in patients S1 (MYCN amplicons), S2 (chromosome arms 1p, 9q, 17q, 20q), S3 (chromosomes 1, 2, 4, 6; Fig. S7c:I-II), S5 (chromosome 5; Fig. S7c:III-IV), S6 (chromosome arm 17; Fig. S7e:I-II), S7 (chromosome arms 1q and 17q; Fig. S7d), and S9 (chromosomes 1 and X; Fig. S7e:III-IV). Also, in “numerical only” low-intermediate risk NBs (n=4), structural rearrangements contributed to parallel rearrangements of the same chromosome. However, there was a difference between clinical-genetic NB subtypes regarding the time point in evolutionary history when different types of CNA occurred (Fig. S7g:I). In “numerical only” low-intermediate risk NBs (n=4), early generation branches (immediately following the stem) contained only whole chromosome aberrations, while chromosome breaks were enriched in subsequent generations. In contrast, in high-risk cases (with and without MNA; n=5) chromosome breaks were frequent already at the earliest generations and remained so through successive generations.
There was also a difference across clinical-genetic NB subtypes in the fitness cost associated with different chromosomal imbalances (Fig. S7g:II-III). As a proxy to the survivability following different types of CNA, we quantified the proportion of progeny (PoP; average fraction of detected single cells) originating from branches with certain CNA profiles compared to branches with other profiles. In low/intermediate risk NB, the PoP from lineages having undergone whole chromosome gains during preceding generations was significantly larger than the PoP from lineages with chromosome breaks or losses. In contrast, high-risk NBs showed a similar PoP resulting from branches dominated by whole chromosome breaks as by gains, while the proportion of progeny was lower for cells having undergone losses. To finally evaluate the timing of phylogenetic branching in NB compared to other tumors, phylogenies based on scWGS for a reported reference group (n=4)11 of pediatric tumors less aggressive than NB were reconstructed. NB was found to exhibit extensive branching even prior to the outgrowth of the major subclone, while this was not the case for the more indolent tumors (Fig. S7g:IV-V). In summary, scWGS confirmed that extensive phylogenetic branching is a consistent feature of NB over time that can emerge already at ancestral stages with low-level MYCN gain. High-risk NBs seem to have developed tolerance to chromosome breaks already at the time of first branching, resulting in complex patterns of subclone diversity that can form the substrate of CCR when new clones emerge to replace those killed off by chemotherapy.