Clinically relevant treatment of PDX models reveals patterns of neuroblastoma chemoresistance

Chemotherapy resistance and relapses are common in high-risk neuroblastoma (NB). Here, we developed a clinically relevant in vivo treatment protocol mimicking the first-line five-chemotherapy treatment regimen of high-risk NB and applied this protocol to mice with MYCN-amplified NB patient-derived xenografts (PDXs). Genomic and transcriptomic analyses were used to reveal NB chemoresistance mechanisms. Intrinsic resistance was associated with high genetic diversity and an embryonic phenotype. Relapsed NB with acquired resistance showed a decreased adrenergic phenotype and an enhanced immature mesenchymal–like phenotype, resembling multipotent Schwann cell precursors. NBs with a favorable treatment response presented a lineage-committed adrenergic phenotype similar to normal neuroblasts. Novel integrated phenotypic gene signatures reflected treatment response and patient prognosis. NB organoids established from relapsed PDX tumors retained drug resistance, tumorigenicity, and transcriptional cell states. This work sheds light on the mechanisms of NB chemotherapy response and emphasizes the importance of transcriptional cell states in chemoresistance.


INTRODUCTION
Neuroblastoma (NB) is a pediatric tumor that arises from the developing sympathetic nervous system and accounts for over 15% of childhood cancer deaths (1). Patients with high-risk tumors receive a very intense multimodal treatment regime that includes high doses of chemotherapy, surgery, radiation, stem cell transplants, and, in some cases, immunotherapy and targeted therapies. An established regimen to treat high-risk NB is the one set by the International Society of Pediatric Oncology-European Neuroblastoma (SIOPEN) (2,3). Rapid COJEC is the induction chemotherapy step established by SIOPEN (3,4). COJEC consists of high doses of five chemotherapeutic drugs (cisplatin, carboplatin, cyclophosphamide, etoposide, and vincristine) distributed in three courses that are alternatively administered in eight 10-day cycles. However, resistance is common. Some patients with high-risk NB display upfront progression despite therapy; more commonly, patients initially respond but subsequently relapse, developing treatment-resistant disease. Thus, the overall survival for patients with high-risk disease remains lower than 50% (5).
Treatment resistance, which is a major clinical challenge and cause of death, is often linked to various aspects of intratumor heterogeneity. In principle, treatment resistance can arise from clonal evolution and selection of tumor cell subclones or by nongenetic mechanisms, including transcriptional reprogramming (6,7). NB is mainly driven by large chromosomal aberrations and is characterized by a low number of somatic driver mutations, even after relapse, and high inter-and intratumor heterogeneity (8,9). Transcriptional and epigenetic analyses revealed that NB can exhibit at least two phenotypic cellular states, which are commonly referred to as adrenergic (ADR or ADRN) and mesenchymal (MES) (10)(11)(12)(13)(14)(15)(16)(17)(18). Cells in the ADR state are lineage-committed sympathetic noradrenergic cells. Cells in the MES state are immature neural crest cell-like or undifferentiated mesenchymal-like cells. Both states have been detected in established conventional cell lines (10)(11)(12)(13)(14) and in patient samples (13)(14)(15)(16)(17)(18), although some studies have yielded conflicting results regarding the existence of the MES cell state in patients and xenograft models (12,19,20). In vitro, MES cells are implicated in NB treatment resistance (10,12,21,22). However, this association has been difficult to demonstrate beyond cell culture. The effects of treatment pressure on NB cell states are unclear because most patient studies use mainly untreated samples obtained during diagnosis.
The characteristically low number of somatic driver mutations, the small number of patients, and the heterogeneity within and across patients have made it especially challenging to identify the mechanisms of treatment resistance in NB (8,9). A limitation of most preclinical studies of chemotherapy resistance is the use of single agents rather than the combination of multiple drugs as is done for NB patients' treatment. Here, we sought to investigate treatment resistance and relapse in NB in a clinically relevant in vivo setting. We developed and used a COJEC-like treatment protocol that included all five COJEC drugs to treat multiple NB patient-derived xenograft (PDX) models that resemble the genotype and phenotype of NB patient tumors (23,24). With this system, we detailed the transcriptomic and genomic changes occurring during NB treatment and at relapse.

Establishment of a COJEC-like treatment protocol using NB PDX tumors
Using three NB PDX models derived from high-risk NB tumors (23,24), each with 1p deletion, MYCN amplification, and 17q gain (fig. S1A), we modeled NB treatment resistance using a treatment paradigm based on COJEC. PDX1 was derived from a COJECrefractory tumor, whereas PDX2 and PDX3 were derived from COJEC-responsive tumors that subsequently relapsed (fig. S1A). In our COJEC-like protocol, cisplatin, vincristine, etoposide, cyclophosphamide, and carboplatin were administered intraperitoneally in a cycled manner, for six 7-day cycles ( fig. S1B). Two tiers of dosing were designed and tested. The basic protocol (COJEC) was the lower dose and was well tolerated overall, whereas the high-dose protocol (COJEC-HD) induced weight loss that required treatment interruptions to allow mice to recover ( fig. S1C). Both nude and nonobese diabetic scid gamma (NSG) mice were evaluated, and no difference in treatment response was observed between the mouse strains (table S1). Mice were randomized into treatment groups when subcutaneous tumors reached ~500 mm 3 . Mice were treated intraperitoneally for 6 weeks with either saline (control), cisplatin (4 mg/kg, three times/week), or COJEC ( Fig. 1A and fig. S1B).
Consistent with the patient response, PDX1 tumors showed no or limited response to treatment with either cisplatin or the basic COJEC protocol (Fig. 1B). Therefore, we subjected this PDX model to the COJEC-HD protocol (fig. S1B). Both COJEC doses led to smaller tumor volumes than controls on day 13 (median survival day for control group) ( fig. S1D), but only COJEC-HD treatment showed an increase in mouse survival (Fig. 1B). Overall, 19 of 20 PDX1 mice treated with cisplatin, COJEC, or COJEC-HD presented progressive disease (PD), and only one mouse showed stable disease (SD) during treatment (table S1). These findings thus mirror the (lack of) therapy response seen in the corresponding patient.
Mice bearing PDX2 tumors responded well to treatment, presenting significant tumor size reduction (fig. S1D; P < 0.05 for each therapy) and significantly increased survival when treated either with cisplatin or COJEC ( Fig. 1B; P < 0.0001). Most COJEC-treated tumors (five of eight) displayed a partial response (PR) and a significantly smaller tumor volume during treatment (table S1 and fig. S1D; P = 0.0008), whereas four of seven cisplatin-treated tumors showed PD (table S1). One COJEC-treated PDX2 mouse presented with complete response (CR) under treatment and relapsed locally after a period without treatment ( Fig. 1B and table S1).
COJEC-treated PDX3 tumors presented an overall PR, with significantly reduced tumor volume during COJEC treatment (fig. S1D; P < 0.0001) and significantly increased survival ( Fig. 1B; P = 0.0006), consistent with the corresponding patient's clinical response (fig. S1A). Two mice presented a CR when treated with COJEC, and relapsed after treatment was removed (table S1). To resemble the clinical regimen in patients, we performed surgical resection on a subgroup of COJEC-treated mice (n = 9) when tumors had reduced below 200 mm 3 (aimed volume reduction, >60%) (tables S1 and S2). Following surgery, the mice did not receive any further COJEC treatment. At the end of the study, two of nine mice in the COJEC + surgery group had relapsed locally (22% relapse rate), and seven of nine were "cured" (tumors did not relapse). Among the COJECtreated PDX3 mice that did not have surgery, two of seven mice relapsed after presenting a CR, and five of seven had tumors that exhibited a PR or SD but regrew rapidly after treatment cessation ( Fig. 1B and table S1). Among the 17 mice in the COJEC-treated group, one had to be euthanized early because of severe weight loss despite exhibiting SD with regard to the tumor. Because this mouse did not complete the 6-week COJEC protocol, for downstream analyses, this tumor was grouped with the surgically collected samples in a group hereon referred as surgery (Sur, n = 10), composed of cured (n = 7) and not cured (n = 3) PDX3 mice (table S2).
Thus, we developed and evaluated a COJEC-like protocol and showed that high-risk NB PDX models display clinically relevant and unique responses to this treatment, ranging from upfront PD to maintained CR. The results mirror those seen in the corresponding patients.

Histological characterization of COJEC-treated NB PDX tumors
When subjected to immunohistochemical analysis, all tumors were positive for the NB marker paired like homeobox 2B (PHOX2B) (Fig. 1C). Hematoxylin and eosin (H&E) staining and blinded morphological analysis of tumor sections revealed that COJEC-treated PDX3 tumors (surgery group) displayed clear signs of morphological differentiation (enrichment of neurofibrillary matrix), whereas this was less obvious in COJEC-treated PDX2 and absent in PDX1 tumors from mice receiving COJEC or COJEC-HD (Fig. 1, D and E, and fig. S1E). Terminal deoxynucleotidyl transferase-mediated deoxyuridine triphosphate nick end labeling (TUNEL) staining revealed increased cell death in PDX1 only for COJEC-HD treatment and in the PDX3-surgery group (Fig. 1, D and F, and fig. S1F). All tumors displayed intense staining for the proliferation marker Ki67 (fig. S1G). PDX1 and PDX2 tumors exhibited large and collapsed blood vessels (fig. S1G), consistent with a worse prognosis (25). PDX3 tumors had overall smaller and open blood vessels. Within the PDX3-surgery group, no statistically significant difference was seen in morphological differentiation or proliferation between the cured and noncured samples (Fig. 1, G and H).

Clonal dynamics and signs of parallel and convergent evolution
To characterize the genetic changes occurring during and after CO-JEC treatment, we performed whole-genome copy-number analysis of selected tumors within the control and treatment groups, as well as of parental tumor organoids before injection. A high number of subclones were detected in all three PDX models during tumor progression and treatment, and all models clearly showed a complex, branched evolutionary pattern (Fig. 2, figs. S2 and S3, and data file S1). For each model, we defined a "stem" clone as having all copy number aberrations (CNAs) present in all organoids and tumors.
The nonresponsive NB PDX1 had a high number of subclones in tumors from both the control and the treatment groups ( Fig. 2A). One tumor in a mouse that received COJEC-HD treatment exhibited a clonal sweep ( Fig. 2A). A clone tree (Fig. 2B) based on phylogenetic analysis ( fig. S2A) illustrated how whole-genome duplication (WGD) in the parental organoids created an abundance of genomic material and many subsequent whole-chromosome or whole-arm CNAs in all analyzed tumors. Complex branching and parallel evolution were detected, but no evidence of enrichment for specific subclones as a response to COJEC treatment was found.
For each PDX, we summed the number of CNAs in each tumor or parental organoid to calculate the "copy number aberration burden" (CNAB). First, this was done for all aberrations detected in the respective stem as a comparison between the PDX models ( fig. S2B). After that, we calculated the CNAB per tumor in addition to the phylogenetically defined stem clone ( fig. S2C). For PDX1, CNAB was comparable between controls and treated tumors ( fig. S2C), suggesting that COJEC treatment did not increase the number of CNAs.   However, the fraction of private aberrations (that is, clones unique to a single tumor and representing high intertumor heterogeneity) was higher in the treated samples ( Fig. 2A and fig. S2D), suggesting diverse evolutionary trajectories in response to COJEC treatment in the individual tumors despite genetically shared starting material. We analyzed the genetic diversity to describe subclonal variation within each tumor ( fig. S2E). Simpson's index of diversity (Ds; 0 to 1, where 0 = no clonal variation and 1 = all subclones are different) showed no statistical difference between the controls and treated samples in PDX1 ( fig. S2E). Compared to PDX1, PDX2 (Fig. 2, C and D) had lower CNAB ( fig. S2C) and fewer private aberrations ( Fig. 2D and fig. S2D). Treatment-specific copy number losses were detected on chromosomes 22q and 8q (Fig. 2D, green shading, and fig. S2F). The 22q region included a specific loss of a small part of LARGE1 (encoding a glycosyltransferase) that was present in all treated samples. The 8q region included RAB2A (encoding a RAS family member involved in vesicular fusion and trafficking from the endoplasmic reticulum to the Golgi) and CHD7 (encoding a DNA binding protein with a chromodomain and helicase activity). CNAB calculation indicated that the total number of chromosomal aberrations was similar between controls and treated samples ( fig. S2C), as was genetic diversity ( fig. S2E). No clonal sweeps were detected in this model.
Most PDX3 tumors showed a marked response to COJEC treatment. A high number of subclones and several selective sweeps occurred in both control tumors and treated samples (Fig. 2, E and F, and fig. S3, A and B). No specific subclone was selected by COJEC treatment, but we observed continuous clonal dynamics during tumor growth both in the presence and absence of treatment. These dynamics were evident by the trend of increased CNAB in the regrown and relapsed tumors ( fig. S2C), as well as in the linear correlation between high CNAB and time of tumor growth (Fig. 2G). Tumors from the control, cisplatin, and COJEC groups collected around the same time showed comparable CNAB (Fig. 2G). These data suggested that CNA events accumulate over time but that the number of different clones within tumors (genetic diversity) was not increased ( fig. S2E). We found two major CNA events occurring in the majority of the samples: an extra copy of the already gained q region in chromosome 17 (17q++), clone 25 ( Fig. 2, E and F, red), and a new gain spanning part of the p arm and all of the q arm in chromosome 1 (1pq+), clone 21 ( Fig. 2, E and F, orange). These genomic events occurred not only in both controls and treated samples but also in subclones that were phylogenetically not derived from clones 25 and 21 (Fig. 2F, orange text). The accumulation of the same genetic aberrations in tumors that were not clonally related represents a typical example of parallel evolution.
There was evidence of parallel evolution of small deletions of specific genetic regions within each of the PDXs (colored stars, Fig. 2). Many of these genes were found in two or all three PDX models, suggesting convergent evolution in aggressive MYCN-amplified NB. For example, deletions in MACROD2 (20p12.1) and in LSAMP (3q13.31) occurred in controls and treated samples across the three PDXs; thus, they are unlikely related to treatment. MACROD2 microdeletions may cause chromosomal instability, and the gene could function as a tumor suppressor in CNA-driven tumors (26). In all three models, the MACROD2 microdeletion was associated with extensive branching and the accumulation of CNAs (Fig. 2, B, D, and F, and fig. S2, A and F). LSAMP encodes a neuronal surface glycoprotein, and the gene has also been suggested to be a tumor suppressor (27,28). PTPRD (9p23) deletions were detected in PDXs 1 and 3. This gene has been suggested as a tumor suppressor in high-risk NB through destabilization of AURKA and MYCN (29). Partial deletion of the gene ITPR1 (3p26.1) was found in PDX1 (parallel evolution) and in PDX3 (stem). ITPR1 encodes an intracellular calcium channel important for apoptosis in response to endoplasmic reticulum stress (30).
In summary, bulk DNA sequencing and clonal deconvolution revealed high subclonal diversity of CNAs in all PDX models. Nevertheless, the intrinsic resistance model PDX1 presented with extremely high clonal diversity, which was retained after COJEC treatment. Convergent and parallel evolution of NB-associated chromosomal regions (e.g., 1 pq+ and 17q++), as well as microdeletions in specific genes (e.g., PTPRD and MACROD2), were detected in both controls and treated tumors. This suggests that they resulted from tumor-intrinsic evolutionary forces rather than chemotherapy-induced selection. No specific CNA subclone could be identified as responsible for chemotherapy resistance and relapse.

High genetic diversity and clonal evolution over time revealed by single-cell DNA analysis
We performed single-cell DNA sequencing (scDNA-seq) and lowpass single-cell CNA analysis of selected PDX tumors to further dissect the subclonal composition and evolutionary patterns. A representative subset of tumors from all groups in the intrinsic resistant (PDX1) and the responsive (PDX3) models were selected ( Fig. 3 and table S2). scDNA-seq analysis of treatment-resistant PDX1 revealed a 58 to 81% of unique single-cell clones in each tumor (Fig. 3A, gray shades) with no clones detected in more than one tumor ( Fig. 3A  and fig. S4A). No enrichment of distinct subclones after treatment could be found, consistent with the results from bulk DNA analysis. High genetic diversity ( fig. S4B) reflected a high number of cellunique clones in control tumors and in tumors sampled post-COJEC. The WGD event, which was detected in the tumor-initiating organoids ( Fig. 2B and fig. S2A), provided a high genetic diversity and genomic instability, which is seen in some patients (31,32). Although we found no indications that a specific clone could explain treatment resistance in this model, it is possible that a high adaptability caused by genomic instability contributes to an evolutionary mechanism for treatment resistance in PDX1 (Fig. 3B).
In the treatment-responsive PDX3, scDNA analysis revealed two major clones: clone A (1p−, MNA, 5p+, and 17q+) and clone C (1p−, 1pq+, MNA, 5p+, and 17q++) (Fig. 3, C and D). Phylogenetic analysis showed that clone A (Fig. 3E) was more parental and that the aberrations associated with clone C (Fig. 3E, 1pq+ and 17q++) occurred after these changes ( Fig. 3D and fig. S4C). This result is consistent with data from bulk DNA analysis. The parental clone A was partly retained in control tumors and in tumors post-COJEC, but the occurrence of this clone was substantially decreased after regrowth and relapse, which represent tumors subjected to a treatmentinduced genetic bottleneck (Fig. 3C). In contrast, clone C (and daughter clones) was enriched in tumors after treatment or at relapse (Fig. 3, C and D, and fig. S4D) and showed aberrations commonly found in patient tumors (33). However, clone C was common in the controls and was not exclusive of treated tumors. Phylogenetic analysis indicated that, in addition to the dominating clones, unique NB cells can accumulate high numbers of CNAs ( fig. S4, A and C). For example, two cells in PDX3 tumor T6 had gone through WGD, similar to PDX1. Enrichment of specific subclones did not lead to a significant 6 of 20 decrease in genetic diversity in tumors after treatment and at relapse, but a trend was visible ( fig. S4, B and D).
Together, scDNA-seq CNA analyses revealed contrasting patterns of clonal evolution between the treatment-responsive PDX3 and the refractory PDX1. PDX3 showed specific clonal enrichment of the clinically relevant 1pq and 17q regions. In contrast, PDX1 displayed an extremely high level of genetic diversity but no distinct clonal enrichment, consistent with the results from bulk DNA analysis.

Transcriptional signatures of COJEC treatment
We performed RNA sequencing (RNA-seq) to explore the effects of COJEC at a transcriptional level. Unsupervised analysis of all the tumors (number of tumors: n = 30 for PDX1, n = 21 for PDX2, and To establish a transcriptional baseline for each NB PDX, independent of treatment, we analyzed the top 1000 most variable genes expressed in the control tumors ( fig. S5E and data file S2). Six welldefined clusters were identified. Gene ontology (GO) analysis revealed that PDX3 (Ctrl-cluster 1) was mainly characterized by nervous system development genes and negative regulation of the mitogen-activated protein kinase/extracellular signal-regulated kinase (MAPK/ERK) cascade ( fig. S5F), whereas PDX1 (Ctrl-cluster 6) contained mainly embryonic developmental signatures and a signature for regulation of Notch signaling ( fig. S5F). Few significant terms were identified for PDX2 with this analysis (Ctrl-cluster 5). Analysis of differentially expressed genes (DEGs) with a false discovery rate (FDR) < 0.01 between control groups confirmed that PDX1 was characterized by up-regulation of genes associated with embryonic development, the mitotic cell cycle, and metabolic processes Next, we analyzed each PDX model to identify how COJEC treatment leads to transcriptional changes in NB. Unsupervised analysis revealed no clear transcriptional pattern in response to treatment for PDX1 and PDX2 ( Fig. 4, B, C, E, and F). Analysis of the top 1000 most variable genes showed down-regulation of one gene cluster in the COJEC groups for PDX1 and PDX2 (Fig. 4, E and F, and data file S2). In both cases, the down-regulated gene clusters were mainly characterized by protein translation and protein localization ontologies ( fig. S6, A and B). Only PDX3 samples showed a tendency to cluster by treatment, especially the cured samples and those that relapsed (Fig. 4D).
Unsupervised analysis of the top 1000 most variable genes in PDX3 revealed four gene clusters ( Fig. 4G and data file S2). Using GO enrichment analysis, we characterized the clusters as "early development" (cluster 1), "cell cycle" (cluster 2), "mixed" (cluster 3), and "nervous system" (cluster 4) ( Fig. 4H; fig. S6, C to F; and data file S3). Relapsed and regrown tumors were characterized by increased expression of genes in early developmental pathways (cluster 1) and down-regulation of genes in nervous system development (cluster 4) (Fig. 4, G to I). A clear difference in expression was also observed between the surgery-collected samples that were cured (Sur-cured) and those that later relapsed (Sur-not cured), with the latter presenting an overall transcriptional signature similar to that of controls or cisplatin samples, with higher expression of clusters 2 and 3 (Fig. 4G). In contrast, the Sur-cured tumors presented with strong down-regulation of genes involved in the cell cycle (cluster 2) and up-regulation of nervous system development genes (cluster 4) (Fig. 4,  (36)] and by comparing patients who presented tumor recurrence/progression under treatment against those that did not ( fig. S7D). The mixed cluster 3, representative of the Sur-not cured and cisplatin samples, correlated with poor prognosis in patients with NB ( Fig. 4J; P < 0.0001). Surnot cured tumors also showed significantly higher expression of DNA repair and cell cycle genes than the Sur-cured samples (P = 1.77 × 10 −03 and P = 3.4 × 10 −02 , respectively; fig. S7E). In addition, we explored the relevance of the MAPK/ERK and Notch signaling pathways, identified in the controls, and pertinent to many different cancers (14,37,38). The Sur-cured samples presented significantly higher expression of genes involved in inhibition of the MAPK cascade ( fig. S7E; P < 0.01) and down-regulation of Notch signaling ( fig. S7E; P < 0.05).
In summary, treatment-resistant and relapsed NB tumors exhibited features of early embryonic development. Tumors that eventually were cured showed features of nervous system development as well as reduced cell cycle and DNA repair gene expression following COJEC therapy, whereas tumors that eventually relapsed (represented by the Sur-not cured samples) maintained a higher expression of cell cycle genes at the time of surgical removal.

Relapsed NB resembles human Schwann cell precursors
NB is derived from progenitors of the sympathetic nervous system. We analyzed how the gene clusters 1 (relapsed tumors) and 4 (Surcured) relate to the human adrenal medulla during normal development (13). Genes in the cluster 1 signature were highly expressed in immature Schwann cell precursors (SCPs), whereas high expression of cluster 4 genes mainly correlated with neuroblasts and chromaffin cells (Fig. 4K). Thus, COJEC-treated relapsed NB PDX tumors resemble human SCPs in the adrenal medulla, which is consistent with other reports (16,39,40).

Identification of NB phenotypes based on transcriptional analysis
Our analysis revealed significantly enriched expression of nervous system development genes (cluster 4) in Sur-cured tumors ( Fig. 4I; P < 0.01) and early development genes (cluster 1) in relapsed tumors ( Fig. 4I; P < 0.05). Previous studies identified ADR and immature MES-like tumor cell states in NB cell lines and patients tumors (10)(11)(12)(13)(15)(16)(17)(18)(19), but their study using in vivo models has been challenging (12,21), and the implications of chemotherapy are not clear. To explore these cell states, we analyzed the gene expression pattern of six pairs of publicly available ADR and MES-like signatures (10,11,(15)(16)(17)(18) across the PDX3 treatment groups ( fig. S8A). Consistently, Sur-cured tumors showed a higher ADR signature across all individual signatures. All of the immature MES-like signatures were higher in relapsed tumors, and some signatures were also high in Sur-cured samples ( fig. S8A). Given the variety of models and methods used to establish the different ADR/MES signatures, and the diverse definitions of "MESlike" or "immature" in published work, we decided to construct novel merged gene signatures by merging publicly available gene signatures along with ours (Fig. 5A). "Merged ADR" (1443 genes) is a signature composed of our cluster 4 genes and those ADR-associated genes in the six public signatures; "Merged MES" (2125 genes) is a signature composed of our cluster 1 genes with all the MES or immature genes in the six public signatures. We analyzed the merged signatures across treatment groups for PDX3. Consistently, the merged ADR signature was enriched in the Sur-cured tumors, and the merged MES signature was enriched in the relapsed tumors (Fig. 5B). Similar to the association between low expression of cluster 4 with poor patient prognosis (Fig. 4J), low expression of merged ADR correlated with poor prognosis in patients with NB ( fig. S8B). The merged MES signature did not associate with NB patient prognosis ( fig. S8B), and neither did the individual MES signatures (data file S4). In summary, we observed that cured tumors display an enrichment of the ADR cell state, whereas the features of relapsed NB resembled the immature MES-like cell state.

Transcriptional features of relapsed and cured NB
Given the differences observed in the Sur-cured samples across the various MES signatures ( fig. S8A and data file S2), we decided to explore the ADR and MES gene signatures based on therapeutic response. This was performed by plotting the merged ADR and MES lists across the Sur-cured and relapse/regrow groups (Fig. 5C and data file S2). From the merged MES list, we identified two subsignatures: "MES-relapse" (enriched in relapsed and regrown tumors) and "MES-cured" (enriched in the Sur-cured samples). GO analysis showed that the MES-relapse signature is enriched in early developmental pathways, cell proliferation, cell migration, and Notch and Wnt signaling (Fig. 5D). The MES-cured signature is characterized by cell death, cell differentiation, and negative regulation of the MAPK/ERK cascade (Fig. 5D). From the merged ADR list, we identified similar subsignatures (Fig. 5C); the ADR-cured subsignature was characterized by nervous system development and differentiation, and the ADR-relapse subsignature was dominated by cell cycle processes ( fig. S8C).
MES-like signatures and markers obtained from NB patient tumors can be derived from NB cells or from stromal cells (15,16).
Here, we identified human-specific transcripts in the RNA-seq to distinguish human NB cells from mouse cells. Analysis of the fraction of RNA-seq "reads" specific to human and mouse showed that the small contribution of mouse reads was highly similar between PDX treatment groups ( fig. S9, A to C). Only Sur-cured samples in PDX3, which have an overall lower MES signature, had a significantly higher presence of mouse stroma ( fig. S9C; P < 0.05). Relapsed samples had low percentages of mouse reads, indicating that the immature MES-like phenotype detected in relapsed tumors represents the NB cells rather than the stromal cells.
The merged signatures include genes from multiple models and methodologies but are quite extensive, and the MES signature lacks predictive power of prognosis ( fig. S8B and data file S4). In the search for more specific signatures and gene markers that correlate with prognosis both in the PDX model and in patients with NB, we selected the ADR-cured and MES-relapse subsignatures, which are characteristic of different COJEC responses. These subsignatures were filtered through the cohort of 219 stage 3 to 4 patients (overall survival, median cut, FDR < 0.05) to obtain new integrated signatures (Fig. 6A). The integrated ADR signature (294 genes) and integrated MES signature (150 genes) were comprised by representative genes from all the individual signatures (data file S2), hence representing a variety of models and methods, and showed a strong predictive power of prognosis in different cohorts of patients with NB (Fig. 6B, fig. S9D, and data file S4). The integrated ADR and MES signatures were significantly higher in the PDX3 Sur-cured samples or in relapse and regrown samples, respectively (Fig. 6C). The Surnot cured samples presented an intermediate expression of both signatures (Fig. 6C). Consistently, PDX3 samples presented an overall higher expression of the integrated ADR signature, whereas PDX1 showed a higher expression of the integrated MES signature (Fig. 6D and fig. S9E).
Protein production from key genes selected from these lists was confirmed in PDX3 tumors (Fig. 6E). PHOX2B was evenly expressed throughout our models (Fig. 1C) and across the PDX3 samples ( fig. S9F); therefore, it was used as a pan-NB marker, similar to its clinical use. Consistent with the RNA results, Sur-cured tumors showed an overall low protein expression of SRY-Box Transcription Factor 9 (SOX9) and Leucine-rich repeat-containing G-protein coupled receptor 5 (LGR5) (MES markers) and high abundance of Tyrosine hydroxylase (TH) (ADR marker), whereas relapse and regrown tumors showed the opposite pattern (Fig. 6E). Coexpression of PHOX2B with ADR/MES markers confirmed the NB identity of the cells. Sur-not cured samples presented high amounts of both MES and ADR markers at the surgical resection time point, but only the expression of the MES markers was maintained in the paired relapses (Fig. 6E). These findings further point to the importance of the immature MES-like cells and transitional phenotypic states in NB chemoresistance and relapse in vivo.

Correlation of MYCN and other key genetic aberrations with COJEC treatment
The MYCN gene is essential in NB pathogenesis (41,42). On the basis of the scDNA-seq data, we determined the number of MYCN copies for PDX1 and PDX3 (Fig. 7, A and B). The nonresponsive PDX1 had higher MYCN copy numbers compared to the responsive PDX3 (Fig. 7A). Within PDX3, all treatment groups had similar numbers of MYCN copies (Fig. 7B). We next analyzed MYCN RNA expression across PDX models and treatment groups (Fig. 7C). Consistent with CNAs, the amount of MYCN RNA was lower in PDX3 compared to PDX1 (Fig. 7D). Sur-cured tumors presented strong MYCN RNA down-regulation, whereas relapsed and Sur-not cured tumors showed expression in the range of PDX1 and PDX2 (Fig. 7, C and E). Consistent with the PDX1 and PDX2 tumors exhibiting higher MYCN transcripts, these tumors had high expression of a MYCN target gene signature (Fig. 7F) (43). Sur-cured tumors also had lower expression of MYCN downstream target genes (Fig. 7G). At the protein level, we observed that relapse samples had a significantly higher fraction of MYCN-positive cells than Surcured samples (Fig. 7, H and I; P = 3.8 × 10 −03 ). Although the lower fraction of immunostained MYCN-positive cells in the Sur-cured samples could be influenced by a higher presence of stroma in this group, the influence of stromal (mouse cells) content is minimized in the paired RNA data. Overall, NB responding to COJEC showed significant down-regulation of MYCN RNA, protein, and its target genes despite high MYCN copy numbers.
On the basis of bulk DNA-seq data (Fig. 2), we observed parallel and convergent evolution of certain fragment losses that affected eight genes, several of which have been previously associated with NB ( fig. S9G) (28)(29)(30). These specific convergent losses seemed to accumulate across COJEC-treated and cisplatin-treated groups ( fig.  S9G), but they did not correlate with treatment response. Instead, we observed that the transcript abundance of these genes correlated with response in the PDX3 model: Sur-cured samples had significantly higher expression of this eight-gene signature ( fig. S9, G and H; P < 0.05). High expression of this signature also correlated with good prognosis in two sets of patients ( fig. S9I). In summary, gene expression levels of MYCN and the other genetic aberrations analyzed here could be important for NB prognosis and treatment response evaluation.

3D tumor organoids from relapsed NB maintain in vivo phenotype and tumorigenicity
We established free-floating three-dimensional (3D) NB organoids, cultured under serum-free conditions with the addition of epidermal growth factor and fibroblast growth factor 2, from control and CO-JEC-treated tumors for all three PDX models (Fig. 8, A to C). All NB organoids were tested with cisplatin and vincristine as single drugs to assess for acquired resistance to chemotherapy (fig. S10, A to C). NB organoids derived from PDX1 and PDX2 showed no difference in drug responses between control and COJEC treated ( fig.  S10, A and B). However, NB organoids derived from a relapsed PDX3 tumor (tumor T5) showed decreased treatment response as compared to PDX3 control organoids ( fig. S10C).
We developed two ex vivo COJEC treatment protocols and applied these five chemotherapies combined to NB organoids derived from PDX3 tumors. Organoids derived from the relapsed T5 tumor displayed higher cell viability and lower cell death under COJEC treatment compared to control NB organoids (Fig. 8D). This suggested that NB organoids derived from COJEC-treated and relapsed NB PDXs can retain chemotherapy treatment resistance ex vivo.
We performed RNA-seq of the PDX3 control and relapsed organoids (untreated and COJEC-treated ex vivo) to examine whether the organoids maintained transcriptional signatures similar to those observed in vivo. Each sample clustered well within their respective groups (Fig. 8E). We applied our previously identified PDX3 cluster 4 (nervous system) and cluster 1 (early development) gene signatures. Relapsed NB organoids displayed lower expression of the cluster 4 signature compared to controls (Fig. 8F), consistent with our in vivo data (Fig. 4, G and I). COJEC-treated control organoids showed significantly higher cluster 1 expression than their untreated counterparts ( Fig. 8G and fig. S10D; P = 1.94 × 10 −06 ).
We analyzed mRNA expression of key NB genes in all control organoids and relapsed organoids. MYCN expression was significantly higher in organoids derived from relapses compared to controls (Fig. 8H), consistent with the in vivo results. The expression of PHOX2B and the ADR markers DBH, TH, and NCAM1 was lower in relapsed organoids compared to control organoids (Fig. 8H). The MES markers LGR5, MAB21L2, SOX9, and FOXC1, which we identified from the in vivo data, maintained an overall higher expression in relapsed organoids compared to controls (Fig. 8H). Nonetheless, within the control organoids, we observed significantly higher expression of some individual MES markers in the COJEC-treated group ( fig. S10E), consistent with the higher expression of cluster 1 in this group (Fig. 7G and fig. S10D). DEG analysis showed that COJEC-treated control organoids presented with down-regulation of cell cycle genes along with up-regulation of genes in specific morphogenesis and cell motility pathways ( fig. S10, F and G). For relapse organoids, the main DEGs identified between the treated and untreated groups were related to cell cycle ( fig. S10, H and I).
Using the relapse organoid derived from PDX3 T5, we confirmed by scDNA-seq that genetic aberrations were maintained from in vivo to ex vivo (Fig. 8I). The relapsed PDX3 T5 was a clonal sweep of subclone #8 (Fig. 2, E and F). The convergent smaller deletions of genes PTPRD and MACROD2, as well as the deletion of PCDH15 (detected in subclone #8 and several PDX3 tumors after treatment), were present in the COJEC relapse organoids and consistently maintained in the established organoids. Thus, the selective sweep of subclone #8 in the relapse was complete, and no new aberrations were detected in the organoids (Fig. 8I). The subclone with LSAMP deletion (#27, yellow), occurring in approximately 60% of the T5 tumor at the time of surgical resection, was not retained either in the relapse or the organoids, illustrating continuous clonal evolution over time.
Last, we injected the relapsed NB organoids subcutaneously into NSG mice (n = 3) and showed that they maintain tumorigenic capacity in vivo (Fig. 8J). Tumors displayed an undifferentiated morphology, were positive for the NB marker PHOX2B, and maintained high amounts of MES markers SOX9 and LGR5 along with low ADR marker TH (Fig. 8K), in concordance with the parental relapse tumor (Fig. 6E). Thus, NB organoids derived from COJECtreated relapsed PDXs exhibited the following properties: (i) display increased resistance to chemotherapy ex vivo, (ii) maintain transcriptional signatures of their corresponding in vivo tumors, (iii) are genetically similar to the in vivo tumor, (iv) are tumorigenic, and (v) maintain relative expression of ADR-and MES-associated markers.
In summary, in this study, we performed comprehensive analyses of the NB chemoresistance landscape, both in vivo and ex vivo, using a clinically relevant COJEC treatment protocol. We highlight the relevance of high genetic clonal diversity in NB, and we demonstrated that transcriptional phenotypic cell states associate with intrinsic and acquired COJEC treatment resistance (Fig. 9).

DISCUSSION
Treatment resistance and relapse are common for patients with high-risk NB. Chemoresistance is usually studied in vivo by applying a single drug treatment (for example, cisplatin). Here, we developed a clinically relevant COJEC-like treatment protocol including the five chemotherapeutic agents usually given as first-line treatment to children with high-risk NB. This protocol was applied to multiple NB MYCN-amplified PDX models, and genomic and transcriptomic analyses were performed at different time points through the treatment course. Our analyses revealed substantial genetic diversity, as well as high intra-and inter-tumor heterogeneity, at both the DNA and transcriptional levels. NB PDXs in mice with upfront progression during chemotherapy and relapsed NB PDXs after chemotherapy showed an immature transcriptional phenotype that resembles the previously described MES-like cell state (10,11,(15)(16)(17)(18). In contrast, NB tumors from mice that were cured after COJEC and surgical resection showed, at the time of surgical resection, features of nervous system development that resemble the ADR cell state (10,11,(15)(16)(17)(18). The 3D NB organoids derived from relapsed PDXs retained chemoresistance, transcriptional and genetic signatures of the in vivo tumor from the mice, and tumorigenicity, showing that these features can be maintained ex vivo.
Although a mutation-centric view has been dominant for years, recent data also point to transcriptionally plastic oncogenic cell states as mediators of cancer development and treatment resistance (44)(45)(46). Several studies have demonstrated the existence of at least two transcriptional cell states in NB: an immature MES-like cell state and an ADR differentiated cell state (10)(11)(12)(13)(15)(16)(17)(18). There is consensus across studies for the ADR cell state, but inconsistencies exist in the definition of the immature MES-like state, with some data even questioning the existence of this cell state beyond cell culture (12,19,20). The discordances observed between studies could be due to (i) differences in methodology, including technical and bioinformatic methods; (ii) the use of different ADR/MES cell state markers; (iii) the diversity of sample material, such as cell lines, mouse models, PDXs, and patients; or (iv) the use of treatment naïve versus treated tumors.
Our findings showed that these transcriptional cell states occur in patient-derived MYCN-amplified NB models both in vivo and ex vivo. Our results also demonstrated the relevance of the immature Fig. 9. Summary of the genetic and phenotypic traits associated with COJEC treatment resistance in high-risk NB. Intrinsic chemoresistance is characterized by a very high clonal diversity, a lowered ADR signature, and a high immature MES signature. An initial response to COJEC can be characterized, at the surgical resection time point, as favorable (high ADR and low MYCN expression) or incomplete (high cell cycle and MYCN expression). An incomplete response can lead to relapse with acquired resistance, characterized by clonal diversity, enrichment of key genetic aberrations (1pq+ and 17q++), low ADR, and high MES signatures. Created with Biorender.com.
MES-like cell state in intrinsic and acquired NB treatment resistance (Fig. 9). In vivo and ex vivo data suggested that COJEC treatment could select for cells with the immature MES-like phenotype and that these cells might contribute to relapse. In contrast, a predominantly ADR signature right after COJEC treatment correlated with the absence of relapses and a good overall prognosis. Our results also identified multiple cells positive for both ADR and MES markers (Figs. 6E and 8K), indicating the presence of cells in an intermediate state, which has been previously proposed (10,11,17,18). In addition, we showed that the time point of transcriptional analysis is important; for responsive tumors, the immature phenotype was primarily identified in the relapses, and this phenotype was only detected in the absence of treatment in tumors with upfront progression. Our results are consistent with findings in which an NB cellular subtype with immature mesenchymal characteristics was enriched in relapsed patient tumors (10,16). Our data also suggested that transcriptional analysis and histochemical evaluation of various ADR and MES markers, in addition to genomic assays, may be clinically valuable when assessing NB prognosis and treatment course. In particular, at the surgical resection post-COJEC time point, when specific signatures (cell cycle, ADR, and MAPK) and markers (SOX9 and LGR5) could provide valuable information about the COJEC response. In addition, our results pointed to the importance of targeting both cell states from the start of treatment, to address upfront progression and prevent relapses.
We identified examples of convergent evolution between the PDX models, parallel evolution within each PDX model, clonal sweeps, and several examples of smaller deletions of genes that regulate nervous system development and chromosomal stability, both of which may influence chemoresistance. Convergent evolution of small deletions in two or more of the PDX models indicated that Darwinian selection occurred during tumor progression. However, these changes were detected in both controls and treated tumors, suggesting that they are important for NB progression but not necessarily for treatment resistance. We found no recurrent CNAs that could explain resistance or relapse. Thus, although the common NB CNAs lay the oncogenic foundation for NB growth, the development of relapses seems to be primarily mediated by transcriptional [and likely epigenetic (16)] changes. Our findings that COJEC resistance was influenced by the immature MES-like phenotype while tumors lacked specific clonal selection at DNA level suggested that treatment-resistant relapses can develop as a result of Lamarckian induction rather than through pure Darwinian selection. The level of (in)stability of the relapse-associated immature phenotype will be important to understand. The phenotype could, in principle, become stable through acquisition of genetic changes. However, we did not identify specific CNAs that could confer such stability. The identification of the immature phenotype in PDX tumors weeks after cessation of treatment in vivo and in organoids established ex vivo indicated a certain level of intrinsic stability of the phenotype even without treatment pressure.
A limitation in our study is that we analyzed a limited set of NB PDX models, and given the broad heterogeneity of NB, it remains to be shown whether our findings represent a general phenomenon of all MYCN-amplified NBs or only a subset of tumors. We did not investigate MYCN nonamplified high-risk NB because of a lack of stable PDX models for this subtype. Our studies were performed using immunodeficient mice, and the immune system can interact with the NB cell phenotypes (47), which could yield a different outcome in fully immunocompetent individuals. On the other hand, our findings demonstrated that NB can adopt an immature MES-like phenotype even in the absence of a complete immune system.
The in vivo confirmation of two NB phenotypic cell states opens for novel therapeutic opportunities because current treatment protocols are not designed to target different NB cell states. It is conceivable that the cell states are not binary, but NB cells exist along a continuum of states. Our results showed the relevance of the MAPK and ERK cascades, in line with findings that suggest therapeutic targeting of these pathways in chemotherapy-resistant NB (9,16,48).
In addition, other work shows that tumor necrosis factor-related apoptosis-inducing ligand (TRAIL) therapy (22) and immunotherapy (47) could also be of interest. Given the clinical heterogeneity in NB patient tumors, and reflected in our NB PDX treatment responses, treatment resistance can likely develop by multiple mechanisms in individual tumors, as demonstrated in melanoma (7). Nevertheless, successful treatment of NB will likely need to target both, or multiple, cell states. Future studies need to define how and when specific targeting of the NB cell states should be implemented. In this context, the heterogeneous PDX models and 3D organoids that recapitulate NB patient characteristics, and the COJEC-like protocol, can serve as clinically relevant tools to further decipher chemoresistance mechanisms and for the preclinical testing of novel therapeutic strategies against COJEC-resistant NB.

Experimental design
The goal of this work was to uncover the mechanisms involved in COJEC treatment resistance in high-risk NB. We designed and tested a COJEC-like protocol using three NB PDX models. In vivo, mice were injected subcutaneously with dissociated PDX organoids, and upon tumor engraftment, the mice were randomly allocated to treatment groups aiming for a minimum of five mice per group (on the basis of previous experimental experience). In vivo experiments were not blinded. Mice were euthanized on the basis of tumor size, weight loss, overall health deterioration, or end of study time. No tumor-bearing mice were excluded. Only one outlier was identified (PDX3-T10), which was considered as part of the COJEC (only) group for mice survival analysis and as part of the surgery group (not cured) for subsequent tumor analyses (table S2). Samples were blinded for morphological differentiation assessment. Bulk DNA and RNA were collected from the same piece of tumor tissue to allow for direct comparison. Single-nucleotide polymorphism (SNP) and scDNA-seq were performed for samples of interest based on tumor response. RNA-seq was performed for all tumors. All RNA analyses were performed using the embedded analysis tools within the R2 Genomics Analysis Platform and Metascape. For in vitro experiments, three biological replicates with three technical replicates each were performed. For the organoids RNA analysis, two biological replicates with three technical replicates each were analyzed.

Animal experiments
All procedures were conducted according to the guidelines from the regional Ethics Committee for Animal Research (permit nos. M11-15 and 19012-19). NMRI nude mice were purchased from Taconic, and NSG mice were obtained from in-house breeding. NB PDX dissociated organoids (2 × 10 6 cells) were suspended in a 100-l mixture (2:1) of stem cell medium and Matrigel (Corning, catalog no. 354234) and injected subcutaneously into the flanks of female mice. Mouse breed was selected on the basis of the tumor engraftment rate for each PDX; PDX1 was analyzed in both nude (n = 9) and NSG mice (n = 24), PDX2 in NSG mice (n = 22), and PDX3 in nude mice (n = 41). Tumor size was measured using a digital caliper and calculated with the formula V = (ls 2 )/6 mm 3 (l, long side; s, short side of each tumor). Mice were randomly allocated to control, cisplatin, or COJEC group once their tumor had reached approximately 500 mm 3 . The cisplatin group was treated intraperitoneally with cisplatin (4 mg/kg) dissolved in sterile saline (Selleckchem, catalog no. S1166) each Monday, Wednesday, and Friday; the control group was treated with equal volume of saline. The COJEC group was treated with intraperitoneal injections of cisplatin (1 mg/kg; Selleckchem, catalog no. S1166) and vincristine (0.25 mg/kg; Santa Cruz Biotechnology, sc-201434) on Mondays, etoposide (4 mg/kg; Santa Cruz Biotechnology, sc-3512) and cyclophosphamide (75 mg/kg; Santa Cruz Biotechnology, sc-361165) on Wednesdays, and carboplatin (25 mg/kg; Santa Cruz Biotechnology, sc-202093A) on Fridays. A subset of PDX1 mice (n = 7) was treated with a high-dose COJEC protocol: cisplatin (1.5 mg/kg), vincristine (0.5 mg/kg), etoposide (6 mg/kg), cyclophosphamide (100 mg/kg), and carboplatin (35 mg/kg). Treatment was administered for a maximum of 6 weeks, making a maximum of 18 injection days per mouse. Treatment response was defined by the parameters established by the Pediatric Preclinical Testing Program (49) presented in table S1. All mice were monitored for weight loss and other signs of toxicity. Treatment was paused if weight loss is >10% to allow for recovery. When the tumor size had decreased to approximately 200 mm 3 , a subgroup of COJEC-treated PDX3 tumors (n = 9) was surgically removed, treatment was interrupted, and mice were monitored for tumor relapses. Mice were euthanized when tumors reached 1800 mm 3 or at a humane end point if weight loss is >20% from the initial weight or notable health deterioration was observed or 1 year after initial cell injection. Upon collection, tumors were divided into pieces and stored for the corresponding experimental purposes.

Immunohistochemistry and histological analyses
Tumors were fixed in 4% formalin, embedded in paraffin, and cut in 4-m sections for staining. H&E (Histolab Products) staining was performed for histopathological analysis. The slides were blinded for morphological differentiation assessment and evaluated by an independent viewer. TUNEL assay (Abcam, ab206386), using either methyl green or hematoxylin as counterstain, was performed to asses cell death. Quantification was performed using the software QuPath 0.2.3. 3,3'-Diaminobenzidine (DAB) staining was performed using the Autostainer Plus (Dako), and fluorescence staining was performed manually. Heat-induced antigen retrieval was done using sodium citrate buffer (pH 6.

Bulk DNA and RNA extraction
Snap-frozen tumors were added to RNAlater-ICE (Invitrogen, #AM7030) and kept overnight at −20°C before extraction. DNA and RNA were extracted from tumors and organoids using the AllPrep DNA/RNA Mini Kit (Qiagen, #80204).

SNP-array DNA genotyping and analysis
DNA from selected PDX samples was subjected to the Cytoscan HD array (Affymetrix) platform according to the standard methods. The obtained CYCSPH files were imported into Chromosome Analysis Suite (Affymetrix, v.4.1.0.90 r29400). All samples from the same PDX were imported simultaneously, and the chromosomes were assessed manually. Only genetic alterations clearly visible by eye, ≥50 kilo-base pair and with a marker count ≥50, were included in further analysis. Constitutional copy number variants were omitted on the basis of manual comparison against the Database of Genomic Variants with genome GRCh37/hg19. For each genetic alteration, its genomic location, type of aberration (loss/gain), and logarithmic median probe intensity ratio (log 2 R) were extracted. The function Rawcopy (v.1.1) (50) in R was used on the Affymetrix fluorescence array intensity (CEL) files, and its output was used to generate Tumor Aberration Prediction Suite (TAPS; v.2.0) (51) plots of copy number clusters for each chromosome and sample.
Each detected genetic alteration was analyzed in unison with the TAPS plots to determine the allelic composition and the clonality (clonal/subclonal) for the aberrations. The log 2 R, together with the ploidy of the genetic alteration (N t ) as well as the background cells (N p ), was used to compute the mutated sample fraction (MSF), defined as the percentage of cells in that sample that harbor this alteration For copy-number neutral imbalance (cnni:s), the mirrored B-allele frequency (mBAF) together with information about its allelic composition (number of A and B alleles, N A and N B ) was used for MSF computation The MSF was used to compute the mutated clone fraction (MCF) defined as the proportion of cells in the sample having a specific genetic alteration divided by the purity or tumor cell fraction (TCF) of that sample. The TCF was computed for each sample by calculating the mean of the MSF values of all aberrations in a sample defined to be clonal by the visual inspection of TAPS plots. The interval of clonal events was computed using the SD of the MSF values of the alterations used for TCF computation (SD MSF ) as Interval of clonal events = TCF ± 2 · SD MSF ─ TCF defined as the MCF interval in which an alteration is determined as clonal. If a genetic alteration has an MCF in that sample that is within the interval of clonal events, then it was set to 100%; otherwise, the MCF value was kept, and the event was defined as subclonal.

Single-cell low-pass CNA sequencing
Single cells were prepared from viably frozen tumors by mechanical dissociation with sterile scalpels and enzymatic dissociation with Liberase DH (0.15 mg/ml) at 37°C and 250 rpm for 20 to 40 min depending on tumor size. Then, tumor cells were treated with deoxyribonuclease I (2.5 g/ml; Sigma-Aldrich), trypsinized (0.05%), and filtered through a 70-m cell strainer. scDNA-seq was performed as previously described (52)(53)(54). Library preparation was performed from each isolated nucleus. Low-pass (0.01 to 0.02×) whole-genome sequencing produced reads that were aligned to human reference genome (GRCh38). Copy numbers were determined by AneuFinder (version 1.8.0) at 1-Mb bin sizes. Analysis quality was determined as previously described (52). Manual curation of single-cell CNA profiles identified true imbalances based on the previously optimized cutoff of five consecutive 1-Mb bins (33). High-grade amplifications of MYCN (2p24.3-24.2) and of 4q28.3-q31.1 were scored on the basis of two consecutive 1-Mb bins. Specific clones based on unique CNAs were identified for all single cells within the respective tumor.
In total, 18 samples [19 to 85 libraries (cells) per tumor] were successfully analyzed. This resulted in a matrix where each row is a bin, approximately 1 million base pairs of size, each column is a single cell, and the matrix elements are the number of copies of that particular segment. The matrix was used as input to an algorithm written in R. This information was used to construct an event matrix where each row is a genomic event, and each column is a cell or a group of identical cells.

Phylogenetic and genomic analyses
To analyze the evolutionary trajectory of the genomic alterations across the PDX samples, DEVOLUTION (55) was used. The input is a matrix in which each row represents a genetic alteration in a particular sample, along with information about that genetic alteration's genomic location, type of alteration, and MCF. Hence, for bulk genotyping data, we obtain an event matrix using DEVOLUTION, and for single-cell whole genotyping data, we obtain an event matrix using a customized algorithm written in R. The event matrices obtained from DEVOLUTION for bulk data and from the customized algorithm for single-cell data were used to reconstruct phylogenetic trees using the maximum likelihood and maximum parsimony methods using the phangorn (v2.8.1) (56) package and visualized using the ggplot2 (v.3.3.5) package. Genetic diversity was calculated using multiple metrics based on clonal/subclonal information and phylogenetic analysis. First, the fraction of private aberrations (clone/subclone detected in only one tumor within the PDX model) gave a measure of inter-tumor heterogeneity. Second, CNAB was calculated from the number of accumulated CNAs to describe the total burden of CNAs within each tumor where f(clone) is the frequency of each clone/subclone in one tumor, and CNA (clone) is the absolute value of the copy number deviation compared to a diploid cell. The sum is over all genetic alterations and all subclones encompassed by that tumor. To allow for comparison between treatment groups within each PDX, only CNAs that were not present within the stem of that PDX were included in the calculations and graphs. Third, genetic diversity was calculated with inverse Ds as where f equals the frequency of each clone. To visually depict the frequency of clones at certain time points and the emergence of new clones over time, a fish plot based on CNA events at specific time points during treatment was produced using an R script (57).
Comparison between groups and test for trend for the genetic diversity measures was done by one-way analysis of variance (ANOVA). The correlation between CNAB and tumor collection day was investigated with Spearman correlation and linear regression. All statistical analysis was performed by the GraphPad software.
All RNA-seq analyses were performed on the R2 Genomics Analysis and Visualization Platform. The Kocak dataset (35) and the Versteeg dataset (36) were used for patient overall survival analysis (median cut). Visualization of t-SNE plots, heatmaps, and volcano plots was done using R2, FDR < 0.01. Heatmaps present the z score for each gene in each sample, and gene hierarchical clustering was performed using Euclidean distance. Box plots for genetic signatures present the average z score values over the gene set for each sample. Welch's t test correction was used to analyze statistical differences between groups for individual signatures. GO enrichment analysis was performed using the GO Platform (63-65) and Metascape (66). Visualization of GO enrichment and DEG was performed using Revigo (67) and Metascape.
In vitro cultures NB PDX organoids previously established derived from PDXs 1, 2, and 3 (23,24,68) and newly established organoids were cultured as free-floating tumor organoids under serum-free conditions in Dulbecco's modified Eagle's medium and GlutaMAX F-12 in a 3:1 ratio, supplemented with 1% penicillin/streptomycin, 2% B27 without vitamin A, fibroblast growth factor (40 ng/ml), and epidermal growth factor (20 ng/ml). Cells were tested for mycoplasma before experiments. The identity of the NB PDX cells was confirmed using short tandem repeat (STR)/SNP analysis. New PDX-derived tumor organoids were established by performing mechanical dissociation of fresh pieces of tumor using a sterile scalpel, followed 24 hours later by enzymatic dissociation with Liberase DH (0.15 mg/ml; Roche ref. 05401054001).

Cell viability and cell death assays
Cell viability and death were measured using the CytoTox-Glo Cytotoxicity Assay (Promega). Organoids were dissociated using Accutase (Sigma-Aldrich), and single cells were seeded and treated in triplicates in white 96-well plates with clear bottom for 48 hours for three biological replicates. Cells were incubated for 48 hours before treatment to allow the formation of small organoids. Luminescence readouts were performed using a Synergy2 plate reader (BioTek). Drugs used for in vitro treatment were the same used for the in vivo work.

Statistical analyses
All statistical analyses pertinent to DNA and RNA analysis were performed as described in the corresponding method sections. The remaining analyses were performed using GraphPad Prism 9. The log-rank test was used to determine significance for Kaplan-Meier survival curves for in vivo studies. Nonparametric Mann-Whitney test was used for the morphological differentiation analysis, and samples were graded as undifferentiated (0), mixed (1), or differentiated (2). For the analysis of TUNEL cell death quantification, ordinary one-way ANOVA was performed with Dunnett's test correction for multiple comparisons by treatment group compared to controls. For tumor volume comparison, ordinary one-way ANOVA was performed with Welch's t test correction for multiple comparisons. For the organoid drug testing, two-way ANOVA with Šidák correction for multiple comparisons by concentration was performed. Overall, a P < 0.05 was considered significant.