Patient-Derived Xenograft Models Reveal Intratumor Heterogeneity and Temporal Stability in Neuroblastoma Free

Malignant tumors are often made up of distinct subclonal populations, which can evolve both spatially, that is, in different regions of the tumor, and temporally as a consequence of tumor growth or treatment. This evolving intratumor heterogeneity (ITH) has been implicated as fundamental in cancer progression, metastasis, treatment resistance, and cancer survival (1, 2), particularly in adult cancers. Recent work suggests that pediatric tumors also exhibit genomic diversity.

Neuroblastoma, a childhood cancer originating from the developing sympathetic nervous system, is characterized by a diverse clinical phenotype (3). Patients with high-risk neuroblastomas often have a poor prognosis despite heavy multimodal treatment, and survivors can suffer from long-term, severe side-effects. Emerging data suggest that clonal evolution and genetic ITH also is a feature in neuroblastoma (4–11) and is involved in neuroblastoma treatment response (8, 9, 11). We have previously shown that high genetic diversity is common in childhood cancer after chemotherapy and correlates to poor prognosis (4, 11). Less is known about transcriptional, proteomic, and functional ITH in neuroblastoma; yet, understanding this is important as it most likely influences biomarker discovery, target identification, treatment decisions, and treatment response.

Patient-derived xenografts (PDX) have emerged as promising preclinical cancer models because PDXs can retain many of the molecular and functional features of their ancestral tumors in patients (12). Thousands of PDX models have thus been established in major academic and industrial preclinical drug testing programs. The PDX–Avatar and the coclinical trial concepts depend on the establishment of a single PDX model mirroring a patient tumor. By simultaneously establishing and treating a PDX–Avatar of a patient enrolled in a clinical trial with a new agent, underlying mechanisms can be studied, novel combination strategies can be evaluated, and potential biomarkers identified. The use of well-characterized PDX models thus holds promise to improve the transition of preclinical drug evaluation data to the clinics. However, the functional consequences of ITH for PDX–Avatar and coclinical trial studies have not been addressed. This could be crucial because clonal evolution and ITH are essential factors in the response to anticancer drugs and the development of treatment resistance (2).

We have previously established neuroblastoma patient-derived orthotopic xenografts (PDOX), which retain the histopathologic, stromal, and metastatic features of aggressive patient tumors (13, 14). We prefer to use PDOXs instead of subcutaneous PDXs, because orthotopic tumors have more relevant biology (15) and retain spontaneous metastatic capacity (16). In this study, up to eight in vivo generations of neuroblastoma PDOXs were established through serial passaging in NSG mice. PDOXs were transplanted as undissociated tumor fragments to avoid genetic aberrations and clonal selections associated to culturing procedures (17). Comprehensive molecular analyses of multiple in vivo passages from each PDOX model were performed to examine temporal evolution of human high-risk neuroblastoma. We demonstrate high genetic, transcriptional and phenotypic stability over time, and identify PDOX-derived gene signatures, which correlate with neuroblastoma-associated processes and with clinical characteristics including patient prognosis. These gene signatures are preserved for more than 2 years of in vivo serial passaging in mice. However, using multiple PDOXs derived from a single high-risk patient tumor, we uncovered diverse transcriptional, proteomic, and phosphoproteomic profiles, suggesting a significant spatial ITH of the patient tumor. Our findings highlight functional and molecular spatial ITH as an element of high-risk neuroblastoma. This complicates the interpretation of PDX–Avatar studies (1 mouse, 1 patient) for biomarker discovery as well as target identification and treatment-guiding decisions in personalized oncology.

Animal procedures were performed as described previously (13, 14). Briefly, 4-to-6-week-old female or male NSG mice were housed under pathogen-free conditions and received autoclaved water and food. We have previously established neuroblastoma PDOXs through the implantation of undissociated patient tumor fragments into the paraadrenal space of immunodeficient NSG mice (13, 14). Established PDOXs underwent serial orthotopic transplantation into next generation NSG mice to establish up to eight in vivo generations for each PDOX model. The study was approved by the Regional Ethics Board of Southern Sweden (289-2011). The Malmö–Lund Ethical Committee for the use of laboratory animals approved all animal experiments (M146-13). See SI Appendix for further information on animal procedures.

Xenograft tumors and mice organs were formalin-fixed, embedded in paraffin, and 4 μm tissue sections were cut and analyzed. See SI Appendix for further information.

PDOXs were snap frozen and stored at −80°C for SNP array analysis. Briefly, DNA was extracted from the tissues using the DNeasy Blood and Tissue Kit (Qiagen) according to the manufacturer's instructions. The Affymetrix CytoScan HD platform was used for SNP array analysis as described previously (13).

DNA was extracted from frozen tumors or constitutional blood from neuroblastoma patients and from the PDOXs using standard procedures prior to fluorometric quantitation and DNA integrity assessment on Agilent Tapestation (Agilent). Exome sequencing was performed on tumors from 5 patients and corresponding constitutional DNA for 3 of the patients, and on PDOXs from different in vivo generations (ranging from G1–G8). See SI Appendix for further information.

For RNA sequencing (RNA-seq), 50 ng of RNA was reverse transcribed according to Ion AmpliSeq Transcriptome Human Gene Expression Kit Preparation protocol (Thermo Fisher Scientific). See SI Appendix for further information.

Proteomic and phosphoproteomic analyses were performed on 10 different tumors of PDOX #5 (G1). See SI Appendix for further information.

A publically available RNA-seq patient cohort of 498 neuroblastoma tumors was analyzed using R2: Genomics Analysis and Visualization Platform (http://r2.amc.nl; Tumor Neuroblastoma - SEQC - 498 - RPM - seqcnb1). See SI Appendix for further information.

See SI Appendix.

Neuroblastoma PDOXs were previously established through the implantation of undissociated patient tumor fragments into the paraadrenal space of immunodeficient NSG mice (13, 14). Patient characteristics are presented in Table 1. Here, we performed serial passaging of these PDOXs to establish up to eight in vivo generations (G) of tumor-bearing NSG mice from each established model. PDOXs from G1–G2 were defined as low in vivo G PDOXs (lowG-PDOX) and PDOXs from G5–G8 were defined as high in vivo G PDOXs (highG-PDOX; Fig. 1A). Although the engraftment rate from patient tumor to mouse was around 50%, as reported previously (13), the in vivo propagation rate once a PDOX has been established, is much higher (220 of 233 mice; 94.4%). The in vivo growth time for the generations of various PDOX models is shown in Fig. 1A. PDOXs #1–3 and 5 showed a stable growth pattern during serial in vivo passaging, but, in contrast, PDOX # 4 presented with dynamic growth between the lowG-PDOXs and highG-PDOXs (Fig. 1A). This likely reflects that the corresponding patient tumor sample showed signs of chemotherapy-induced differentiation and slow tumor-cell proliferation. Thus, the growth rates of the PDOXs remained mostly stable upon serial passaging.

Neuroblastoma lowG-PDOXs retain cellular morphology and protein marker expression of the patient tumors (13, 14). When highG-PDOXs were analyzed within this study, they were found to retain the poorly differentiated morphology as well as positive focal-to-regional expression of the neuroblastoma markers synaptophysin, chromogranin A, and tyrosine hydroxylase. There was also retention of strong diffuse expression of the neural cell adhesion molecule (NCAM; CD56; Fig. 1B). HighG-PDOXs also retained spontaneous metastatic capacity to bone/bone marrow, liver, and lungs, as shown by NCAM IHC of tissue from these organs of the mice (Fig. 1C). Thus, the histopathologic and clinical hallmarks of patient tumors and lowG-PDOXs are preserved in highG-PDOXs.

We performed whole-exome sequencing of the original patient tumors and their corresponding lowG-PDOXs and highG-PDOXs. Most of the somatic mutations found in the patient tumors were retained in the corresponding lowG-PDOXs and highG-PDOXs (Fig. 1D). The exception was # 4, where PDOXs displayed several unique mutations as compared with the patient tumor sample. The number of common mutations were 18 (PDOX #1), 27 (PDOX #2), 17 (PDOX #3), 6 (PDOX #4), and 11 (PDOX #5). A number of unique mutations were found in patient tumors (0–6), and in the PDOXs (0–19) as well as mutations unique for individual PDOX generations (0–6; Fig. 1D; Supplementary Table S1). We identified mutations in genes previously implicated in neuroblastoma progression (Fig. 1E). Specifically, NRAS p.Q61K (PDOX #3), ALK p.F1174L (PDOX #4) and p.L1240V (PDOX #5), ARID1B mutations (PDOX #4), and MYCN p.P44L (PDOX #5) mutations have been implicated in aggressive neuroblastoma growth and/or treatment relapse (5, 7, 18–22). MYCN p.P44L is a potentially activating mutation, which could lead to high MYCN levels despite no MYCN amplification (20, 23). A complete list of all mutations is shown in Supplementary Table S1.

SNP array analysis was performed to determine how allelic imbalances evolve in PDOXs following serial in vivo passaging. Comparing patient tumors with their corresponding PDOXs, we found both similarities and PDOX-unique allelic imbalances (Fig. 2A; Supplementary Table S2). However, the latter were a minority of the total number of imbalances and there was a high level of consistency regarding the total number of aberrations between patient tumors and their corresponding lowG and highG PDOXs (Spearman ρ = 0.94, P = 1.3 × 10⁻⁷; Fig. 2B and C). Chromosomal changes typical of neuroblastoma such as 1p36 deletion, MYCN amplification (MNA) and 17q gain were retained in both lowG- and highG-PDOXs of PDOXs #1 to 4 (Fig. 2D). Most structural changes that differentiated the non-MNA (N-MNA) patient tumor #5 and the corresponding PDOXs #5 were located to chromosome regions affected by complex chromothripsis–like changes, a type of chromosomal rearrangement associated with poor prognosis in neuroblastoma (21). Detailed data on the allelic imbalances are presented in Supplementary Table S2.

We performed whole-transcriptome analysis (RNA-seq) on the primary patient tumors (n = 1 for each tumor), and lowG-PDOXs and highG-PDOXs (n = 6–8 for each PDOX model, and n = 2–4 for each G). In total, RNA-seq was performed on 46 tumor samples.

t-SNE clustering was performed on all samples (primary tumors – n = 5; PDOX – n = 41) revealing distinct separation of the different PDOX models based on their overall gene expression (Fig. 3A). The primary tumor samples, however, seemed to form a cluster away from their respective PDOX models. Gene ontology analysis on the differentially expressed genes between the primary tumors and their corresponding PDOXs (ANOVA FDR < 0.01) displayed enrichment of genes involved in primarily stromal-associated processes including immune response, extracellular structure organization, and angiogenesis (Fig. 3B; Supplementary Table S3). This suggests that the reason for the clustering of primary tumors away from the PDOX models was due to the presence of human stromal components, which is expected as PDOX tumors lack human immune and stromal cells after several in vivo passages.

To focus on differences within and between the different PDOXs, we again performed hierarchical and t-SNE clustering of all samples following multi-group ANOVA analysis across the different PDOX models (top 3,000 genes, FDR < 0.001; Fig. 3C and D). Each PDOX model displayed distinct expression patterns of multiple gene clusters. The different expression patterns between the PDOX models remained stable regardless of serial passaging generation, indicating that the inherent differences between models are retained over time. Furthermore, focusing on the varying genes between the PDOX models also resulted in the concurrent clustering of the primary tumor samples with their respective PDOX models (Fig. 3C and D) showing that the temporally stable gene signatures between models are also reflected in their respective tumors of origin.

Across the different tumor models, we identified 6 main clusters of genes (Fig. 3C; Supplementary Table S4). Cluster 2 genes displayed highest expression in PDOXs #1 to 4, all of which are MNA. Ontology analysis of cluster 2 genes indeed confirmed an enrichment of MYC and MYC-MAX target genes as well as genes involved in general transcription and translation (Fig. 3E). In contrast, N-MNA PDOX #5 displayed enrichment of genes involved in cell cycle regulation as well as indications of signaling via WNT and NGF (cluster 6). Interestingly, PDOX #1 displayed a significant enrichment of genes involved in embryonic morphogenesis (cluster 4), which may reflect the young age of the patient (16 months) relative to the others. Furthermore, both PDOX #1 and 4 displayed enrichment of genes involved in neurogenesis and related terms, suggesting potentially varying differentiation states across the different PDOX models. Additional detailed ontology findings are presented in Fig. 3E and Supplementary Table S5.

To investigate the clinical relevance of the identified gene clusters, we analyzed their expression patterns in a publicly available cohort of 498 neuroblastoma patients (Fig. 3F; ref. 24). K-means clustering of patients based on cluster 2 genes, which are associated with the MNA PDOXs #1 to 4, displayed 3 distinct groups of patients (Fig. 3G) with the highest expression of the genes being associated with a poor prognosis and specifically the MNA subset of stage 4 patients (Fig. 3F and G). In contrast, cluster 6 genes that were specific to the N-MNA PDOX #5 displayed their highest expression amongst patients with stage 4 N-MNA and was also associated with a poor prognosis (Fig. 3F and H). Interestingly, patients with MNA displayed the lowest expression of cluster 6 genes, reflecting the mutually exclusive pattern of this signature in the PDOX models. Overall, these data support the clinical relevance of the transcriptional phenotypes observed across the different PDOX models. Expression patterns of all gene clusters in the patient cohort are shown in Supplementary Fig. S1.

Here, we have shown that tumors of patients with neuroblastoma contain higher expression of genes involved in human stromal composition including immune infiltration, ECM processes, and angiogenesis than their corresponding PDOXs. Furthermore, the defining clusters of the different neuroblastoma PDOXs reflect those of their respective primary tumors and these features remain stable and transplantable across multiple generations of serial orthotopic passaging (greater than 2 years of in vivo growth). In addition, these PDOX-derived gene clusters correlate to clinically relevant subgroups of patients with neuroblastoma.

To shed light on functional and molecular spatial ITH, we implanted ten different tumor fragments from a single patient tumor (the high-risk, N-MNA tumor in patient #5) into 10 NSG mice. PDOXs were established in all the 10 mice. Unexpectedly, the time period from implantation until significant tumor burden varied dramatically among the mice. We classified the 10 mice based on the tumor engraftment time into three groups; group 1 (97–149 days in vivo, n = 4), group 2 (198–201 days in vivo, n = 3), and group 3 (317–342 days in vivo, n = 3; Fig. 4A).

SNP array analysis revealed that the various PDOXs retained the typical neuroblastoma-associated chromosomal copy number changes of the patient tumor (e.g., 1p del, chromothripsis-like complex rearrangements of chromosome 5 and 17q gain), and we found no recurrent differences at copy number level between the groups that could explain the differences in growth rate (Fig. 4B). Exome sequencing revealed that all ten PDOXs contained neuroblastoma-associated mutations, ALK p.L1240V and MYCN p.P44L (Fig. 4C). We found no evidence of mutations that were consistently different among the groups 1 to 3.

t-SNE clustering performed on RNA-seq data from the 10 PDOX models displayed intermingled clustering of the 3 groups of samples (Fig. 4D). Unsupervised hierarchical clustering (top 1,000 SD genes) displayed transcriptional heterogeneity between samples, where we could identify 5 gene clusters (Fig. 4E; Supplementary Table S6). Gene ontology analysis of the different gene clusters revealed significant enrichment of processes involved in neurogenesis and neuronal differentiation, mesenchymal transition, and cell cycle regulation as well as stromal composition and angiogenesis (Fig. 4F; Supplementary Table S7). Expression of these clusters identified different subgroups of the PDOXs, suggesting functional and phenotypic heterogeneity amongst the different samples. When these gene clusters were examined between the different groups defined by time-to-tumor burden, although finding trends that may hint at differences between the groups, no significant differences of cluster expression were found (Fig. 4G and H). Furthermore, multigroup ANOVA analysis (FDR < 0.05) between time-to-tumor burden groups revealed no significant differences in gene expression. Taken together, this suggests that transcriptional and phenotypic heterogeneity does exist within biopsies from the same parental tumor, however this heterogeneity does not account for difference in time-to-tumor take. In addition, we examined expression of the 5 gene clusters derived from PDOX #5 biopsies in all primary and PDOX samples (n = 46; Supplementary Fig. S2). Using these clusters, each PDOX model (#1–5) could be distinctly separated from each other, suggesting that differentially regulated processes defined within a given tumor (PDOX #5 in this case) may reflect processes that differ between different tumor models. The fact that the PDOX #5 biopsies did still cluster together, however, supports the notion that ITH is not as prominent as intertumor heterogeneity.

To again address the clinical relevance of the gene clusters defined within the N-MNA stage 4 PDOX #5, we examined their expression in the corresponding subset of patients from the public dataset described above (24). K-means clustering displayed that expression of these gene clusters could identify two subgroups of patients within the N-MNA stage 4 cohort (Fig. 4I–K). These subgroups were primarily defined by differences in cluster 4 and cluster 5 expression (Fig. 4J), with subgroup 2 displaying significantly higher levels of both clusters. In addition, these subgroups differed with significantly with regard to patient prognosis (Fig. 4K).

In addition to genetic and RNA-seq analyses, we performed mass spectrometry–based proteomic and phosphoproteomic analysis to further examine spatial ITH in the ten PDOXs from patient #5. Unsupervised clustering of the top varying peptides in both the proteomic (Fig. 5A; Supplementary Table S8) and phosphoproteomic (Fig. 5B; Supplementary Table S9) analyses revealed heterogeneity amongst the samples. Gene ontology analysis again revealed varying expression of proteins and phosphoproteins involved in stromal contribution as well as neuronal differentiation and axon guidance amongst others (Fig. 5C and D; Supplementary Tables S10 and S11). These results are in line with those obtained from the RNA-seq analysis, suggesting that heterogeneity does exist within PDOX models derived from the same tumor, and that may reflect different states of differentiation as well as varying interaction with and composition of tumor stromal components. Of note, correlation analysis between RNA-seq and proteomic data revealed a lack of correlation between samples derived from the same PDOX, however from different pieces of the tumor (Supplementary Fig. S3). This may be an indication that the heterogeneity observed within the primary tumor also exists within the subsequent PDOXs.

As was observed in the RNA-seq data, varying protein expression patterns were not correlated with time-to-tumor burden, and multigroup ANOVA analysis did not reveal significant differences in protein expression between groups. When grouping the more latent groups together (group 2+3), however, we did identify differentially expressed proteins and phosphoproteins between group 1 and group 2+3 (Supplementary Fig. S3; Supplementary Tables S12 and S13). Notably, protein expression of dihydropyrimidinase-like 3 (DPYSL3) was high in the fast-growing group 1; this has previously been implicated in neuroblastoma progression (25). Additional proteins with higher levels in group 1 included: Ras GTPase–activating protein-binding protein 2 (G3BP2), involved in Ras signal transduction; stathmin (STMN1), involved in neuroblastoma metastasis (26); cofilin-1 (CFL1) implicated in neuroblastoma susceptibility (27), and microtubule-associated protein (MAP2; ref. 28). Among the phospho-proteins, more highly expressed in group 1 were microtubule-binding or -associated proteins (e.g., MAP2, DCX, MAPRE1, DPYSL5, and CLASP1) and proteins linked to neural development processes (e.g., DPYSL3, DCX, DBNL, and DPYSL5). In addition, Ser-9 phosphorylation of GSK-3b is associated to phosphorylation of DPYSL3 at sites Thr-509, Thr-514, and Ser-518 (29), all of which were upregulated in our group 1 dataset (Supplementary Fig. S3; Supplementary Table S13). Notably, Ser-9 phosphorylation of GSK-3b in neuroblastoma leads to inactivation of GSK-3b, which in turn stabilizes the MYCN protein (30). The proteomic differences between the time-to-tumor burden groups were validated with targeted mass spectrometry to confirm these results with an orthogonal method (Supplementary Methods + Supplementary Table S14).

In summary, mass spectrometry–based analysis identified ITH within multiple PDOXs from a single patient tumor. The patterns of ITH found by protein analyses resembled the patterns found across samples using RNA-seq analysis, including variation in genes involved in neuronal differentiation and stromal composition. Although these results did not explain differential times to tumor burden, they do shed light on potential targets in subsets of neuroblastoma cells, all of which are tumorigenic.

The functional consequence of ITH for PDX modeling and its subsequent impact on biomarker discovery and treatment decision–making, is not well understood. Although we and others have previously shown that genetic ITH is a feature of neuroblastoma, studies of transcriptional, proteomic, phosphoproteomic, and functional ITH has been less of a focus. Here, we performed serial orthotopic in vivo passaging, multi-sample implantation, and comprehensive molecular characterization of neuroblastoma PDOXs and their corresponding patient tumors. By establishing ten PDOXs from a single tumor from a patient with high-risk N-MNA, we found distinct growth patterns and diverse molecular profiles suggesting substantial spatial ITH of the patient tumor. Furthermore, we show that neuroblastoma PDOXs can display stable phenotypic, genetic, and transcriptional features despite growing for more than 2 years in NSG mice. Thus, high-risk neuroblastoma can display elements of both intrinsic genetic, transcriptional, and phenotypic stability as well as striking molecular and functional ITH, resulting in diverse tumor phenotypes within a single tumor.

ITH has been highlighted as an important feature of various adult cancers (1) and is implicated in the variable response to treatment. More recently, genetic data have shown that ITH is also found in neuroblastoma (4–11) and is implicated in neuroblastoma treatment response (8, 9, 11). Our study builds on this work to further demonstrate diverse growth patterns and molecular profiles in multiple PDOXs derived from a single patient tumor. It is conceivable that the different main biological processes (neuronal/mesenchymal/stroma and angiogenesis) found in the PDOXs reflect different subpopulations within the patient tumor. Interestingly, time-to-tumor burden did not correlate strongly with genetic or transcriptional variation. Although this may suggest that the strongest representative signals of the tumor's heterogeneity do not directly affect time-to-tumor take, this may also be an indication that different growth rates may be a result of technical variation. Nonetheless, molecular ITH is observed amongst PDOX models derived from a single tumor, which is consistent with findings from colorectal cancer, where diverse tumor cell populations within genetically uniform tumor cell lineages have been observed (31). Further studies examining the differences in response to therapeutics across the different PDOX biopsies are, of course, of extreme interest.

The findings stress the need for multiple biopsy sampling to cover the diversity of high-risk tumors. Considering that ITH is a critical factor involved in the development of treatment resistance (32), and our previous findings suggest that high genomic diversity in childhood cancer correlates with poor prognosis (4, 11), our current results complicate the direct translation of results obtained from personalized PDX–Avatar studies to clinical testing of the corresponding patient. This is particularly important for tumors with high ITH, as these are likely to be the aggressive tumors that we need better treatment for. Given the diversity observed in our studies using multiple samples from a single tumor, one can expect increased geno- and phenotypic diversity when comparing PDXs derived from treatment naïve versus relapsed samples, and from primary versus metastatic sites.

One attempt to handle the issue of high spatial ITH would be to establish multiple PDXs from different regions of the same patient tumor and perform drug testing on each individual PDX. This could be demanding and time-consuming, and sometimes not possible due to shortage of available tissue; however, it could open up possibilities to track and treat subclonal tumor populations involved in chemotherapy resistance. Implantation of heterogeneous patient-derived cell cultures derived from different tumor regions could theoretically also create PDXs with retained ITH. A third option is to establish PDXs through cells derived from liquid biopsies with the assumption that these blood-borne tumor cells would capture the intratumor diversity of the patient tumor.

The results emphasize the need for combination treatment strategies of high-risk neuroblastoma. Furthermore, strategies to target general neuroblastoma oncoproteins, as recently shown using a GPC2-directed antibody-drug conjugate (33), should be further pursued in light of the striking functional and molecular ITH shown in this paper. The phosphoproteomic analysis identified proteins and phosphosites previously implicated in cancer and neuroblastoma progression [e.g., phosphorylation of GSK-3b at Ser9 that stabilizes the MYCN protein in neuroblastoma (30)], as well as novel phosphosites, which provide potential target opportunities for neuroblastoma. Importantly, these findings also highlight the need for mass spectrometry for deciphering ITH and for discovery of novel therapeutic targets. Interestingly, RNA and proteomic data from the same PDOX models did not display significant correlation. The observed varying processes across samples, however, were similar. As extractions for each of these protocols were performed on different pieces of the tumor, this may suggest that the ITH present within the primary tumor may be reestablished within each of the derived PDOX models. Further experiments involving multiple biopsies in both the primary tumor and each of the resulting PDOX tumors may shed light onto this issue.

In contrast to the spatial ITH originating from the ancestral patient tumor, we also found that neuroblastoma can display certain features of temporal stability. We identified stable gene expression signatures for different PDOX models related to neuroblastoma processes, for example, N-MYC signature, neurogenesis, and embryonal morphogenesis, and these signatures were transplantable and correlated to poor clinical outcome in a large cohort of patients with neuroblastoma. These findings reflect an intrinsic temporal stability of a subset of neuroblastoma-associated genes when tumors are grown in a controlled in vivo environment without external selection pressure. Nevertheless, it is still conceivable that neuroblastoma PDX serial passaging also can lead to clonal selection, as recently shown in PDXs derived from adult tumors (34). Technical factors that could influence clonal selection processes in vivo include implantation site (orthotopic vs. subcutaneous), passaging procedure (in vitro cultured cells vs. undissociated tumor fragments), and choice of mouse strain. Furthermore, it is likely that addition of an external selection pressure (e.g., chemotherapy or radiotherapy) would alter the observed temporal stability. For instance, genomic changes following neuroblastoma chemotherapy relapse include mutations predicted to activate the RAS-MAPK pathway and/or induce epithelial–mesenchymal transition (5, 6). The identification of ALK F1174L and ARID1B mutations only in the PDOXs and not in the tumor from patient #4 illustrates the problem of taking only one tumor sample from a patient, as these alterations could have been present in other regions of the original tumor. Other explanations that cannot be ruled out include very low allele frequency of the mutations in the patient samples or de novo appearance of the mutations in the PDOXs.

Our study highlights molecular and functional spatial ITH in high-risk neuroblastoma and demonstrates the importance of assessing ITH, to maximize the usefulness of the data in biomarker discovery, preclinical target identification and drug testing studies (Fig. 5E). A limitation of our study is that only one patient tumor was included in the analysis of ITH. Future studies should extend on our findings to additional neuroblastomas including the different subtypes. Nevertheless, the findings should be important for other pediatric and adult aggressive tumors with high ITH, as ITH has been implicated as a key player in cancer treatment response and resistance. Indeed, recent comprehensive genetic analyses of adult tumors and their corresponding PDXs revealed that PDXs overall resemble the genetic landscapes of their parental tumors, but PDXs can also diverge genetically from the parental tumors, in part due to genomic instability and spatial ITH (34). PDX models can still serve as important tools for tumor assessment and treatment planning, as well as for identification of novel targets. However, our results suggest that PDXs should be used in the context of larger cohorts, including multiple patient tumors, rather than drawing individual conclusions from 1 mouse–1 patient Avatar studies.

No potential conflicts of interest were disclosed.

Conception and design: N. Braekeveldt, K. Stedingk, D. Gisselsson, S. Påhlman, D. Bexell

Development of methodology: N. Braekeveldt, K. Stedingk

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): N. Braekeveldt, K. Hansson, I. Øra, A. Börjesson, R. Noguera, T. Martinsson, D. Bexell

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): N. Braekeveldt, K. Stedingk, S. Fransson, A. Martinez-Monleon, D. Lindgren, H. Axelson, F. Levander, J. Willforss, R. Noguera, J. Koster, T. Martinsson, D. Gisselsson, S. Påhlman, D. Bexell

Writing, review, and/or revision of the manuscript: N. Braekeveldt, K. Stedingk, S. Fransson, D. Lindgren, H. Axelson, F. Levander, J. Willforss, K. Hansson, I. Øra, A.P. Berbegall, R. Noguera, T. Martinsson, D. Bexell

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): N. Braekeveldt, K. Hansson, T. Backman, S. Beckman, J. Esfandyari, A.P. Berbegall, J. Karlsson, T. Martinsson, D. Bexell

Study supervision: N. Braekeveldt, J. Koster, T. Martinsson, D. Bexell

This work was supported by grants from the Swedish Cancer Society (CAN 2017/376 to D.Bexell), Swedish Childhood Cancer Foundation (PR2017-0003 to D. Bexell), the Swedish Research Council (2017-01304 to D. Bexell), Mrs Berta Kamprad Foundation (to S.Påhlman), the SSF Strategic Center for Translational Cancer Research-CREATE Health (to S.Påhlman), the Strategic Cancer Research Program BioCARE (to D. Bexell), Crafoord Foundation (to D. Bexell), Jeanssons stiftelser (to D. Bexell), Mary Béves Stiftelse för Barncancerforskning (to D. Bexell), Ollie och Elof Ericssons stiftelser (to D. Bexell), Berth von Kantzows stiftelse (to D. Bexell), The Royal Physiographic Society in Lund (to D.Bexell), Åke Wibergs stiftelse (to D. Bexell), The Tegger Foundation (to K Stedingk), Gyllenstiernska Krapperupsstiftelsen (to D. Bexell), Skåne University Hospital research funds (to D. Bexell), grants from the Scientific Foundation Spanish Association Against Cancer 2015; FIS contract PI17/01558 and CB16/12/00484 (to R.Noguera), and grants from the ISCIII & FEDER (European Regional Development Fund), Spain (to R.Noguera). The authors would like to acknowledge the support of Science for Life Laboratory/Clinical Genomics, Gothenburg and of the National Genomics Infrastructure (NGI)/Uppsala Genome Center and UPPMAX for providing assistance in massive parallel sequencing and computational infrastructure. Work performed at NGI/Uppsala Genome Center has been funded by RFI/VR and Science for Life Laboratory, Sweden.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.