Tumor to normal single cell mRNA comparisons reveal a pan-neuroblastoma cancer cell

Neuroblastoma is a childhood cancer that resembles developmental stages of the neural crest. It is not established what developmental processes neuroblastoma cancer cells represent. Here, we sought to reveal the phenotype of neuroblastoma cancer cells by comparing cancer (n=19,723) with normal fetal adrenal single cell transcriptomes (n=57,972). Our principal finding was that the neuroblastoma cancer cell resembled fetal sympathoblasts, but no other fetal adrenal cell type. The sympathoblastic state was a universal feature of neuroblastoma cells, transcending cell cluster diversity, individual patients and clinical phenotypes. We substantiated our findings in 650 neuroblastoma bulk transcriptomes and by integrating canonical features of the neuroblastoma genome with transcriptional signals. Overall, our observations indicate that there exists a pan-neuroblastoma cancer cell state which may be an attractive target for novel therapeutic avenues. One Sentence Summary Comparisons of neuroblastoma cells and relevant normal adrenal cells reveal a uniform cell state that characterizes neuroblastoma across patients and disease phenotypes.


Introduction
Neuroblastoma is a childhood cancer that exhibits a diverse pattern of disease (1), from spontaneously resolving tumors to a highly aggressive cancer. Neuroblastoma arises from aberrant differentiation of the neural crest, mostly within the fetal adrenal medulla, via the sympathetic lineage. Different stages of adrenal medullary cell differentiation have been implicated as the cell of origin of neuroblastoma (2)(3)(4) and proposed to underlie the diversity of clinical phenotypes (2,5).
Advances in single cell transcriptomics have enabled the direct comparison of cancer and normal reference cells. Applied to childhood cancer, such analyses have revealed specific cell types or developmental processes that tumors adopt and adapt (6). Here, we sought to establish, through cancer to normal single cell mRNA comparisons, which developmental processes and cell types neuroblastoma recapitulates and correlate our findings with disease phenotypes.

Building a reference of human fetal medullary cells
In the first instance, we built a normal cell reference for our analyses, by quantitatively defining transcription of cells underpinning human adrenal development. We undertook single cell mRNA sequencing (10x Genomics Chromium platform) of seven adrenal glands obtained from first and second trimester human fetuses (Fig. 1A, table S1). To verify technical and biological replicability we included bilateral adrenal glands in two cases. Following stringent data filtering, including removal of doublets and ambient mRNAs, we obtained count tables of gene expression from a total of 57,972 cells. These cells could broadly be divided into cortical, medullary, mesenchymal cells, and leukocytes (Fig. 1B, fig. S1A), using canonical human adrenal markers ( fig. S1B, table S2). We focused our analyses on medullary cells (n=6,451), from which most cases of neuroblastoma arise (7). The absolute and relative number of medullary cells varied across medullas, likely reflecting technical and biological differences, such as the relative proportion of cortex to medulla (table S1).
Recent single cell analyses of murine fetal adrenals glands have provided transcriptional definitions of four principal medullary cell types (8)(9)(10) : Schwann cell precursor (SCP), bridge cells, chromaffin cells, and sympathoblasts. To search for these cell types within our data, we used two orthogonal unbiased methods (logistic regression; canonical correlation analysis) to quantitatively map murine reference signals onto our human data. Unlike manual annotation of cells by a selection of (biased) marker genes, these approaches avoid the risk of misclassifying cells when interspecies differences on individual markers exist. Accordingly, we identified the human correlates of the four medullary cell types (Fig. 1C-D). Of particular note were the neural crest like SCP cells, identified previously in murine (8), but not human fetal tissues. We validated these cells by single molecule FISH of SCP defining markers (SOX10, ERBB3, MPZ, PLP1) in the fetal medulla and other human fetal organs -fetal gut, eye, skin, and bone ( fig. S2).
Comparing the four medullary cell types across different adrenals did not reveal significant differences in their overall transcriptional "identity" (Fig. 1E), although nuances in transcription may exist across gestational ages (table S3). We therefore merged cell signals from different adrenals to form the reference for cancer-to-normal cell comparisons and for additional analyses on fetal medullary cells, including trajectory analyses and comparison to murine data ( fig. S1C,   fig. S3).

Identification of cancer cells in neuroblastoma
We next sought to establish the relationship between neuroblastoma cells and human fetal medullary development. We generated single cell mRNA readouts from 21 fresh neuroblastoma specimens from two different centres, using two platforms (Chromium 10x, CEL-seq2). We analysed data from these different platforms and centres independently to ascertain technical and biological replicability of our findings. We obtained tissue from untreated patients at diagnosis (n=4), or following treatment with cytotoxic agents at resection (n=17) (table S1). As pre-treated neuroblastomas tend to be largely necrotic at surgical resection, we guided sampling to viable tumor areas through pre-operative metabolic cross-sectional imaging in these specimens (metaiodobenzylguanidine scan), coupled with morphological assessment of frozen sections in 15/17 pretreated cases. Our study cohort represented the three principal prognostic categories of neuroblastoma: low (n=3), intermediate (n=5), and high risk (n=13) cases. In total we obtained 19,723 cells, with variable contribution from each tumor ( Fig. 2A-B; table S1), which segregated into four main cell types: leukocytes (n=10,593), mesenchymal cells (n=4,227), cells bearing markers of Schwannian stroma (n=665), and putative tumor cells exhibiting adrenal medullary-like features (n=3,396).
The first challenge was to identify bonafide cancer cells amongst tumor derived cells. Marker based cell typing alone was of limited value as controversy exists as to which cell types in neuroblastoma represent cancer cells. Although there is a general consensus that neuroblastoma cancer cells exhibit medullary-like features, it has also been suggested that interstitial (mesenchymal) and Schwannian stroma cells, commonly found in neuroblastoma, may be cancerous (11)(12)(13). We therefore used somatic copy number changes to identify cancer cells, by interrogating each cell's mRNA sequence for evidence of the somatic copy number changes underpinning each tumor. Extending a previously developed method (6), we integrated single cell tumor RNA information with single nucleotide polymorphisms (SNP) arrays or whole genome sequencing of the tumor DNA (table S4). This allowed us to assess the allelic imbalance of SNPs in each cell's mRNA sequence across patient specific copy number segments. This analysis revealed that in our cohort only adrenal medullary-like cells were cancerous, but not interstitial or Schwannian stroma cells (Fig. 2C-D).

Comparison of cancer to normal cells at single cell resolution
Next, we investigated which stage of adrenal medullary development cancer cells recapitulate by performing cancer to normal cell comparisons, using two independent methods. We found that neuroblastoma cancer cells did not recapitulate adrenal development but had assumed the state of sympathoblasts ( Fig. 2E-G, fig. S4A). Importantly, although cancer cells formed discrete clusters within, and across patients ( Fig. 2A-B), the cancer to normal cell comparison resolved this diversity into a common sympathoblast-like phenotype (Fig. 2H-I). The sympathoblast state was present in all patients ( Fig. 2H-I).
To validate that neuroblastoma cancer cells principally resemble sympathoblasts, we interrogated neuroblastoma bulk transcriptomes. We curated gene expression profiles of 650 neuroblastomas from two different clinically annotated cohorts that represent the entire clinical spectrum of neuroblastoma (TARGET (n=152) (14) and SEQC (n=498) (15) cohorts). We probed these expression data for transcripts that define the four populations of the normal human fetal medulla: SCPs, bridge cells, sympathoblasts, and chromaffin cells ( fig. S1B). Across cohorts, we found a clear signal of sympathoblast mRNAs (Fig. 3A), thus verifying the sympathoblastic state of neuroblastoma. The sympathoblast signal in bulk transcriptomes was reproduced by using murine reference data ( fig. S4C). Neuroblastomas that arose outside the adrenal gland also exhibited sympathoblast signals (Fig. 3B). The sympathoblast signal was present across the principal risk groups of neuroblastoma, suggesting that it transcends the clinical diversity of neuroblastoma (Fig. 3C).

Features of the neuroblastoma genome corroborate the sympathoblast state
To further corroborate the transcriptional evidence that neuroblastoma cells resemble sympathoblasts, we overlaid features of the neuroblastoma cancer genome with transcriptional developmental pathways. We first analysed how somatic cancer (14) and germline predisposition (16) genes of neuroblastoma are utilised by fetal medullary cells. We found that the majority of these genes were most highly expressed in sympathoblasts ( Fig. 3D; fig. S5A), in particular the two most common neuroblastoma predisposition genes, ALK and PHOX2B, and the cancer gene, MYCN, somatic amplification of which is a key adverse marker used clinically for neuroblastoma risk stratification. Thus, the same genes that operate as cancer genes in neuroblastoma are utilised in normal development predominantly by sympathoblasts. Similarly, genes encoding the synthases of disialoganglioside (GD2), which is a near ubiquitous marker of neuroblastoma that is therapeutically targeted by antibody treatment (17), were predominantly expressed in sympathoblasts (Fig. 3D). We note that the expression of the enzyme catalyzing the final step in GD2 synthesis, B4GALNT1, was confined to sympathoblasts, whereas genes further upstream in the synthesis pathway, ST3GAL5 and ST8SIA1, showed a less restricted pattern of expression within the medulla.
A variety of copy number changes have been described in neuroblastoma, the most recurrent and pertinent of which occur on chromosomes 1, 11, and 17 (18). The presence of these copy number changes carries fundamental prognostic significance and determines treatment intensity in most clinical contexts and treatment protocols. Given the sympathoblastic state of neuroblastoma, it seemed conceivable that somatic copy number changes may impact on the expression of genes that define sympathoblasts. We were able to directly test this hypothesis. For example, we compared gene expression of cancer cells that harbour loss of chromosomes 1 and 11 with cancer cells that do not carry these changes (table S5). We found that the resulting differentially expressed genes were significantly enriched for high confidence sympathoblast marker genes (tf-idf>1, 283 out of 329 markers, p<0.0001, hypergeometric test) (table S6).
Furthermore, the majority of these genes (254/329) exhibited lower expression in cells with chromosome 1 and 11 loss, even if those genes did not reside on chromosomes 1 or 11.
Following on from this, we correlated genomic regions of recurrent chromosome 1 or 11 loss with the genomic localisation of gene expression in fetal medullary cells. The localisation of the boundaries of altered copy number segments is associated with gene expression in human cancer (19). Thus, the genomic position of these boundaries may encode information about the expression profile of the cancer cell of origin. We found regions of sympathoblast specific gene expression that statistically significantly overlapped with the most common boundaries of chromosome 11q loss (p<0.05; permutation test) ( Fig. 3F; fig. S5B) and the most recurrent region of chromosome 1p loss (p<0.05; permutation test) ( Fig. 3F; fig. S5C), thus further corroborating the sympathoblast state of neuroblastoma.

Differences between low and high risk tumors
The measurement of medullary cell signals in bulk transcriptomes indicated that compared to lower risk tumors, the sympathoblast signal was less pronounced in high risk tumors in both cohorts (Fig 3C). This weakened sympathoblast signal may be driven by differences in tissue composition, or by distinct properties of cancer cells themselves. We reanalysed the sympathoblast signal in single cancer cells, this time subdividing neuroblastoma cells by whether they originated from low or from high risk tumors. We found that some high, but not low risk cells, exhibited a significantly decreased sympathoblast signal (Fig 3E). This difference could not be accounted for by technical variation (fig. S4B). Differential gene expression between the two groups of high risk cells showed that the weakened sympathoblast signal was underpinned by, amongst other mRNAs, transcripts defining neuronal and neuroendocrine lineages (table S7).
In aggregate, these findings indicate that a subtype of high risk cell may exist that is characterised by a fainter sympathoblast signal.

Utilisation of fetal transcripts in neuroblastoma
As neuroblastoma cells utilise transcripts of fetal medullary cells, it is conceivable that the expression of some of these genes is relatively restricted to the fetal period and may thus represent formidable therapeutic targets. The precedent for this concept is the disialoganglioside, GD2, which is widely expressed by neuroblastoma cells and some fetal lineages including sympathoblasts (Fig. 4A). In post-natal tissues, however, expression of key genes encoding GD2 synthases is most confined to the central and peripheral nervous systems (20). Therapeutic antibodies targeting GD2 are deployed in most treatment protocols of high risk neuroblastoma.
Although pain, probably mediated by expression in peripheral nerves, is a common adverse effect, significant CNS toxicity is rare. To identify transcripts with a similar distribution to GD2 synthases, we searched for mRNAs utilised by single neuroblastoma and fetal medullary cells (hereinafter referred to as fetal cancer transcript). We then assessed their distribution across panbody post-natal tissues (from GTEx data), whilst validating their expression in neuroblastoma bulk transcriptomes. These analyses revealed two patterns of fetal cancer gene expression across post-natal tissues. One pattern followed the distribution of mRNAs encoding GD2, whereby expression of fetal cancer transcripts were predominantly expressed in neuroblastoma and in CNS tissues. Some genes, including those pertaining to catecholamine synthesis, exhibited an even more restricted expression pattern. Together these analyses reveal a number of fetal cancer genes that may lend themselves as targets in neuroblastoma (table S8).

Differences between neuroblastoma and normal fetal medullary cells
Having established the similarities between neuroblastoma cancer cells and fetal medullary development, we now looked for differences. Previous bulk mRNA analyses have described gene expression changes that characterise neuroblastoma tissues (3) (21) . Our ability to directly compare cancer cells with their normal counterpart, fetal medullary cells, enabled us to distil the transcriptional essence of neuroblastoma cells, which comprises 91 differentially expressed genes (tf-idf > 0.85, table S9), significantly enriched . To verify, and to correlate, these findings with disease phenotypes, we measured mRNA levels of these differentially expressed genes in bulk neuroblastoma transcriptomes and ascertained that a significant number of them had above average gene expression (n=37/89, p<0.001, hypergeometric test; fig. S6). The expression of 10 of these varied significantly by clinical risk group, even after controlling for key determinants of clinical risk, age and MYCN amplification (negative binomial generalised linear model) (Fig.   4B). In addition to known markers of clinical risk (e.g., PRAME), we identified novel markers that have not previously been described and could in theory be used to improve risk stratification in clinical practice in the future. Biologically, the transcription factor gene, SIX3, may be one of the most interesting findings amongst differentially expressed genes, representing an embryonal gene that normally regulates forebrain and eye development (22).

Discussion
The principal finding of our investigation was that the neuroblastoma cancer cell represented an aberrant fetal sympathoblast. Our expectation might have been that neuroblastoma exhibits a cellular hierarchy adopted from adrenal development, similar to, for example, the childhood kidney cancer, Wilms tumor, recapitulating fetal nephrogenesis. However, neither our direct analyses of single cancer cells, nor indirect insights from bulk cancer transcriptomes would corroborate such a hypothesis.
We identified cancer cells by directly assessing SNP ratios in regions of cancer defining copy number changes, rather than inferring copy number from gene expression levels. Accordingly, we did not find evidence that cells other than sympathoblast-like cells were malignant. Therefore, our findings suggest that the sympathoblastic state predominates across the entire spectrum of neuroblastoma. The transcriptional similarity between neuroblastoma cancer cells and sympathoblasts may indicate that neuroblastoma directly derives from sympathoblasts.
However, it is also conceivable that neuroblastoma arises from other cell states of the neural crest lineage that assume the transcriptional state of sympathoblasts upon malignant transformation.
The neuroblastoma transcriptome has been studied extensively from bulk tissues, revealing a plethora of potential targets that have been evaluated over past decades. Nevertheless, leveraging the resolution of single cells to perform precise comparison of cell types, we identified additional potential targets in neuroblastoma including fetal genes that are utilised by neuroblastoma cells with a restricted expression profile in post-natal tissues. Furthermore, direct comparison of cancer and fetal medullary cells revealed neuroblastoma genes that correlate with outcome in two independent validation cohorts, which may prove useful for refined risk stratification.
Overall our findings reveal a surprising homogeneity in the cellular identity of neuroblastoma.
Although future, larger scale analyses may define the nuances of this signal and perhaps identify rare cancer cell types, our study would suggest that the principal normal cell correlate of neuroblastoma is the human fetal sympathoblast.

Ethics statement
Informed consent for research was obtained from participants (or their carers). Studies

Fetal adrenal tissue processing
Tissue was minced with a scalpel and then placed in 15ml falcon tube filled with 5ml Digestion solution (5ml, dilute Liberase TH stock solution (1 vial of 5mg Liberase TH (Sigma 5401135001) powder in 2ml of 1X PBS with final concentration is 2.5 mg/ml.) in 1X PBS (100ul enzyme stock + 4.9ml 1x PBS). Tissue was incubated in 37'C for 30min (water bath) and mixed with 1ml pipette after 15 minutes to facilitate the digestion. 5ml STOP solution (40ml, 2% FBS in 1x PBS (40ml 1x PBS + 800ul FBS)) was added. The mixture was filtered through a 30um strainer and spun down, 750g, 5min at 4'C. If the pellet was red 1ml of 1X RBC (00-4300-54) lysis buffer was added and incubated for 3min at RT. 10ml of STOP solution was added and cells were gently mixed to wash, then centrifuged, 750g, 5min, 4'C. Pellets were resuspended into appropriate amount of 1x PBS or STOP solution, counted and processed on 10x Chromium platform.

Neuroblastoma tissue processing -GOSH pipeline
Surplus tumor tissue obtained at diagnostic biopsy or tumor resection was processed immediately after receipt in the histopathology laboratory (< 1 hour after interventional radiology/surgical procedure). Tissue was minced using a scalpel, and then incubated in RPMI, supplemented with 10% FCS, 1% L-glutamine and 1% Penicillin/Streptomycin, with Collagenase IV (1.6 mg/ml; (patients were treated with debulking type surgical procedure) and preoperative radiation and chemotherapy had been given to the patients (table S1). We extensively optimized the workup protocols and samples were gently dissociated to the single-cell level and sorted by flow cytometry. Single-cell suspensions were subsequently subjected to single-cell RNA-Seq using the CEL-seq2 platform. In brief, freshly resected human neuroblastoma specimens were collected and processed immediately following histological confirmation of the presence of viable tumor tissue. Tumor tissue was minced into 3-4 mm pieces and digested with collagenase type I, II and IV (2 mg/ml) in the optimized medium at 37°C, with agitation for 2 hours. The resulting suspension containing small tumor fragments was passed through a 70μm cells strainer, washed with a cold organoid medium. Half of the tumor fragments were used to generate tumor organoids, while the other half were fragments were further dissociated with NeuroCult kit until the single-cell suspension was obtained. Cells were washed with cold organoid medium and stained with 7AAD or DAPI to enrich for live cells while others were stained with FITC-GD2 or PE-GD2 antibodies to enrich for tumor and PE-CD3 antibodies to enrich for T cells and to exclude erythrocytes. Cells were stored on ice until sorting using a FACS Jazz or FACS AriaII.
Single live cells were sorted into 384-well hard-shell plates with 10 μ l mineral oil, 50 nl of RT primers, dNTPs and synthetic mRNA Spike-Ins and immediately spun down, snap-frozen on dry ice, and stored in -80C to proceed with the total transcriptome amplification, library preparation and sequencing.

10X library preparation and sequencing
The concentration of single cell suspensions was manually counted using a hemocytometer and adjusted to 1000 cells/ul or counted by flow cytometry. Cells were loaded according to the standard protocol of the Chromium single cell 3' kit (v2 and v3 chemistry). All the following steps were performed according to the standard manufacturer's protocol. We used one lane of an Illumina Hiseq 4000 per 10x chip position.

CEL-seq2 library preparation and sequencing
All samples that were processed proceeded with the total transcriptome amplification, library preparation and sequencing into Illumina sequencing libraries as previously described (23).

Copy number detection from DNA data
The ascatNGS algorithm(26) (v4.0.1) was used to estimate tumor purity and ploidy and to construct copy number profiles prior to running the Battenberg algorithm (v2.2.5) (github.com/cancerit/cgpBattenberg) to allow for tumor sub-clonality in bulk DNA sequencing data. For DNA data measured with Illumina SNP arrays, copy number segments were defined by manual inspection of logR and B allele frequency data.

Mapping and quantification of scRNAseq -10X
Single cell RNA-seq data were mapped and counts of molecules per barcode quantified using the 10X software package cellranger (versions 2.0.2 and 3.0.2) to map sequencing data to version 2.1.0 of the build of the GRCh38 reference genome supplied by 10X.

Mapping and quantification of scRNAseq -CEL-seq2
Sharq preprocessing and the QC pipeline have been applied to process the single-cell RNA-seq data as described (27). Read mapping has been done using STAR version 2.6.1 on the hg38 Patch 10 human genome. The featureCounts function of the subread package (version 1.5.2) was used to assign reads based on GENCODE version 26.

Quality control of fetal adrenal single cell data
Ambient mRNA contamination was removed from each channel individually using SoupX (28).
Cleaned count data was then normalised, scaled to have a mean of 0 and standard deviation of 1 for each gene, principal component analysis was performed using highly variable genes, and grouped into clusters using a community detection finding algorithm taking the first 75 principal components as inputs and a resolution parameter of 10. This was performed using the Seurat package in R (29, 30). The clustering parameter was deliberately set to this extremely high value to produce a very large number of small clusters.
Following clustering, we designated each cell as either passing or failing a series of quality control metrics. We marked as failing QC any cell with >30% expression due to mitochondrial genes, fewer than 300 genes detected, fewer than 1000 UMIs detected, or a doublet score greater than 0.2 as determined by Scrublet (31). In addition to excluding these individual cells, we also marked as failing QC any cell belonging to a high resolution cluster containing more than 50% of cells that fail QC due to one of the above filters. Finally, we excluded from all further analysis any cell marked as QC failed for any of the above reasons ( fig. S7A).
Next we calculated a S and G2M phase score for each cell using Seurat (29, 30) and excluded any cell where either score was greater than 0. This removed any cells that showed evidence of being in a phase of the cell cycle other than G0/G1. We excluded these cells from our reference map as cells in the S/G2M phase tend to cluster together based on their phase of cell cycle and not their underlying cell type, reducing the utility of the data as a reference map of cell types present in the developing adrenal gland.
As before, data was normalised by sequencing depth, scaled to 10,000 counts, log transformed, and genes were scaled to have mean 0 and standard deviation 1. Principal component analysis was performed on the scaled data using the most highly variable genes. We selected the number of principal components to use for clustering and visualisation to minimise the molecular cross validated error (32) (https://github.com/constantAmateur/MCVR). Using these principal components, we calculated a uniform manifold approximation and projection (33) two dimensional representation of the data for visualisation and calculated clusters using a k-nearest neighbours algorithm with resolution parameter set to 1.

Quality control of 10x tumor data
We excluded nuclear mitochondrial genes, heat shock proteins and ribosomal genes from our analysis. To filter lower quality cells, we removed any cells that had greater than 20% expression originating from mitochondrial genes, expressed fewer than 300 distinct genes, or contained fewer than 1000 UMIs ( fig. S7B). Data was normalised by sequencing depth, scaled to 10,000 counts, log transformed, and genes were scaled to have mean 0 and standard deviation 1 using the Seurat 'NormalizeData' function. Principal component analysis was performed on the scaled data using the 2,000 most variable genes. Using 75 first principal components, we calculated a uniform manifold approximation and projection (UMAP (33)) two dimensional representation of the data for visualisation and calculated clusters using a community detection algorithm with resolution parameter set to 1. We performed these steps using the Seurat package in R (29, 30).

Quality control of CEL-seq2 tumor data
We removed RNA spike-ins and genes that were not present in the 10X gene list. We excluded nuclear mitochondrial genes, heat shock proteins and ribosomal genes from our analysis. To filter lower quality cells, we removed any cells that had greater than 20% expression originating from mitochondrial genes, expressed fewer than 200 distinct genes, and contained fewer than 500 UMIs (fig. S7C). We used less stringent thresholds for CEL-seq data than for 10X data due to data quality. Data was normalised by sequencing depth, scaled to 10,000 counts, log transformed, and genes were scaled to have mean 0 and standard deviation 1 using the Seurat 'NormalizeData' function. Principal component analysis was performed on the scaled data using 2,000 most variable genes. Using 50 first principal components, we calculated a uniform manifold approximation and projection (UMAP (33) ) two dimensional representation of the data for visualisation and calculated clusters using a community finding algorithm with resolution parameter set to 1. We performed these steps using the Seurat package in R (29, 30).

Cluster annotation
Using well established marker genes of different cell types curated from the literature (table S2), we assigned a cell type to each cluster. Where two clusters were annotated as the same cell type, we merged them together. As further confirmation of our annotation, we next identified marker genes for each population algorithmically (table S6). To do this we used a method which uses the tf-idf metric to identify genes specific to each population, as implemented in the "quickMarkers" function in the SoupX R package (28). We further filtered genes to include only those genes with a p-value less than 0.01 after multiple hypothesis correction (hypergeometric test).

Comparison of mouse and human adrenal medulla
Using an established map of the mouse medulla as our gold standard (8) we mapped this reference to our human data in two orthogonal ways. For this comparison, we restricted both human and mouse transcriptomes to those genes that could be unambiguously mapped between species (1-to-1 orthologs). Using data from mice taken from E12.5 and E13.5 as independent references, we calculated a similarity score using logistic regression (6) for each human medulla cell for each of the four murine medullary populations. Additionally, we used a canonical correlation analysis based method to transfer the murine labels to the human medulla, as implemented in the Seurat R package (29). We identified differentially expressed genes between matching mouse and human medulla using edgeR (34) as described below. To obtain a list of genes with a large effect size difference between mouse and human, we filtered this list to include only genes with FDR less than 0.01 absolute log fold change greater than 1, fraction of cells expressed greater than 80% in the species with highest expression and less than 1% in the species with lowest expression. We further removed any gene that showed a significant difference between species in all four medullary populations (table S10).

Differential expression in early and late medulla
To identify genes that were differentially expressed between early (<21 weeks) and late (>21 weeks) cell types in the adrenal medulla, we performed a differential gene expression analysis between early and late cells using FindMarkers function in Seurat for each medullary cell type (SCPs, Bridge, Sympathoblats and Chromaffin cells). We pre-filtered genes that were detected at <25% frequency in either population. We identified all genes that passed the 0.01 p-value threshold after multiple hypothesis correction (table S3).

Differential expression in high risk cancer cells with low and high sympathoblast signal
We categorised high risk tumor cells as "low sympathoblast signal" and "high sympathoblast signal" by applying k-means clustering with k=2 to the sympathoblast similarity score (logistic regression). To identify genes that were differentially expressed between low and high sympathoblast signal in high risk cancer cells, we performed a differential gene expression analysis between low sympathoblast signal and high sympathoblast signal cells using FindMarkers function in Seurat. We pre-filtered genes that were detected at <25% frequency in either population. We identified all genes that passed the 0.01 p-value threshold after multiple hypothesis correction (table S7).

Pseudotime analysis
We ordered cells into a trajectory based on the similarity of their transcriptomes to create a pseudotime ordering of the medulla using monocle3 (35)(36)(37). We constructed a monocle3 object using data extracted from a Seurat object, rather than re-processing the data in monocle3. We used the "learn_graph" function to build a principal graph of the data. We then used "order_cells" function to place cells on a pseudotime trajectory, selecting a node (as defined in "learn_graph") as a starting point of the trajectory. We used an SCP node based on prior knowledge from murine data that medullary cells originate from SCPs (8). In addition, to further investigate the identity spectrum between bridge and sympathoblasts, we performed the same analysis after subsetting the Seurat object to only include bridge cells and sympathoblasts.

Differential expression of genes along pseudotime trajectory
To identify transcription factors that were differentially expressed along the graph of the medulla, we performed Moran's I test using "graph_test" function in monocle3. We limited our gene set to the 1,665 transcription factors (37), and identified all genes that passed 0.001 p-value threshold after multiple hypothesis correction (table S11). In addition, to further investigate the identity spectrum between bridge and sympathoblasts, we performed the same analysis, this time using all genes, on just the bridge cells and sympathoblasts, and identified all genes that passed 0.001 p-value threshold after multiple hypothesis correction (table S11).

RNA Velocity
RNA velocity was calculated on the foetal adrenal reference map data using velocyto (38) to calculate the number of reads supporting spliced/unspliced isoforms of each gene in each cell.
The velocity field projection to the UMAP embedding was then calculated using the velocyto.R R package, pooling 20 k nearest neighbours and performing a gamma fit to the top/bottom 2% expression quantiles.

Cell similarity calculation
To measure the similarity of a target single cell transcriptome to a reference single cell data-set we used the methodology based on logistic regression outlined in detail in (6). Briefly, we trained a logistic regression model with elastic net regularization (alpha=0.99) on the reference training set. We then used this trained model to infer a similarity score for each cell in the query data set for each cell type in the reference data.
Softmax normalization was not used to allow for the possibility that some cells in the query data set do not resemble any of the cell types in the reference data set. Predicted logits were averaged within each cluster in the query dataset. This approach was implemented using the "glmnet" package in R(39).

Genotyping individual cells in single cell RNA-seq data
The details of the genotyping of individual cells is explained in detail in the supplementary materials with a brief description provided here. For each sample, regions of clonal copy number change and heterozygous SNPs were identified from either whole genome sequencing or copy number arrays (table S4). Taking copy number regions resulting in allelic imbalance, heterozygous SNPs were phased together across the entire copy number aberration. In each single cell transcriptome, the number of counts derived from either the major or minor allele was calculated for each copy number change, by investigating the nucleotide sequence of the mRNA reads allele counter ((40) (https://github.com/cancerit/alleleCount). For each cell, a posterior probability was calculated for both a diploid genome and the tumor genome. Fig. 2 C-D show the posterior probabilities of the tumor genotype for each cell. For downstream analysis, we considered only clusters of cells that had more cells that were definitively tumor (posterior probability > 99%) than definitively normal (posterior probability < 1%).

Identification of cell-type specific genes
To identify genes that were specific to one cell type and no other in the adrenal medulla, we regenerated algorithmically defined markers as above, but with podocyte cells (6) included as a negative control. To use only the most specific markers, we applied a tfidf cut-off of 1 to select genes that were specific to each medullary cell type relative to all other cells in the dataset. Note that a tf-idf cut-off of t implies that the global rate at which a gene occurs in cells is less than exp(-t/rate_local), where rate_local is the rate at which the gene occurs in the cells in the target cluster.
We also removed genes that were expressed in more than 20% of cells in any other single cluster in the dataset. The resulting gene list consisted of genes specific to each cell type in developing adrenal medulla("stringent_markers category in the "dataset" column, table S6).

Adrenal gland cell type signal in bulk RNA-seq data
To identify which adrenal cell type bulk RNA-seq Neuroblastomas resembled most, we measured the expression levels of genes defined in table S6 in two bulk RNA-seq Neuroblastoma datasets, TARGET (14) and SEQC (15). For both data sets, the transcripts per million reads (TPM) value was calculated for each gene in each sample. We defined a gene as being "present" in a sample when its log 2 (TPM) expression value exceeded that of the peak of the log 2 (TPM) distribution across all genes across all samples. For each gene, we then calculated the fraction of samples in which it was present, both globally and stratified by risk group.

Defining neuroblastoma samples outside the adrenal gland
To define which neuroblastoma samples in TARGET (14) arose outside the adrenal gland, we filtered ICD-O(41) description of each tumor by sites of occurrence. We excluded all samples where descriptions were absent or contained words "kidney", "adrenal" , "abdominal", "abdomen", "unknown", "retroperitoneum", "other", identifying 21 samples where tumors occured outside the adrenal gland.

Differential expression between fetal adrenal medulla and tumor cells
To identify the key transcriptomic differences between Neuroblastoma tumor cells and the normal fetal medulla, we performed differential expression analysis between fetal medulla single cells and tumor cells from CEL-seq2 and 10X tumor datasets. For each dataset, we extracted tumor clusters and merged them with fetal adrenal medullary cells. We labelled all tumor cells as "tumor", and used the "quickMarkers" function from SoupX R package (28) to identify tumorspecific genes for each dataset. We used a cut-off of 0.01 p-value after multiple hypothesis correction (hypergeometric test), and genes that were expressed in more than 20% of cells in any other single cluster in the dataset, to define a set of tumor markers.
We merged the two marker lists and calculated an average tfi-df value for each gene present in either or both lists, and used a tf-idf cut off of 0.85 to define a single list of tumor markers. We then filtered the list to exclude genes expressed by more than 25% of leukocytes in each tumor dataset (table S9).

Differential expression between low risk and high risk Neuroblastomas in TARGET
To test if any of the genes defined in table S9 displayed risk group dependence, we compared their expression levels between 4S and high risk bulk Neuroblastomas in the TARGET data (14).
To do this we used a negative binomial model, where we set the square root of the overdispersion value to 0.4 (as genuine biological replicates were not available). We constructed a generalised linear model, with age, MYCN status, and risk group as covariates then used a quasilikelihood F-test to calculate genes for which the risk group coefficient was significantly nonzero (0.05 FDR, edgeR (34) package) (table S12). We did not perform equivalent analysis in SEQC data as the raw counts were not available for this dataset.

Identification of neuroblastoma transcripts not expressed in post-natal tissues
To identify transcripts that were shared between cancer cells and the fetal medulla, we used the "quickMarkers" function in the SoupX R package to find the markers of cancer cells in each tumor dataset (CEL-seq and 10x). We removed genes that were expressed in fewer than 50% of cancer cells and genes that were expressed in more than 20% of cells in any other single cluster each dataset. We overlapped the gene lists for the two tumor datasets, and filtered the list to exclude genes that were expressed in fewer than 10% of adrenal medullary cells in our reference.
We then looked at the average expression (log2(TPM)) of those genes (table S8) in low and high risk neuroblastoma bulk samples and GTEx (42) bulk samples. We sorted the gene list by increasing average expression outside the brain in GTEX data. The top 25 genes (i.e. genes with lowest expression outside the brain) are shown in Fig. 4A.

Association between recurrent CN changes and sympathoblastic expression
To compare regions of recurrent CN changes in neuroblastoma to the genomic pattern of expression in reference cell populations we first defined break-point regions and recurrently changed regions. We defined a recurrently changed region as one which had a gain (or a loss) in 200 or more samples (out of 556) (18). For regions with a well defined transition region (e.g.  Fig. 3F, fig. S8. To calculate the significance of the relative abundance of the sympathoblastic and chromaffin we calculate the average difference in the smoothed expression between the two cell types in the regions of interest. To calculate a p-value for the significance of the observed difference in smoothed expression, we compared the observed value to the value in 1000 randomly selected regions of the genome of the same size. We then defined the p-value to be the fraction of the randomly selected regions with a difference greater than or equal to the observed one ( fig. S5B-C).

Differential expression of cells carrying prognostic CN change
In our 10X data, we identified one sample (GOSH021) that had loss of both chromosome 1p and 11q, which was not present in any of the other tumor samples. To determine the transcriptomic consequences of these prognostically important copy number changes we calculated the differentially expressed genes between those cells with and without the copy number changes.
We limited this analysis to only those cells in the G1 phase of the cell cycle to exclude confounding with differences in proliferation. We then performed a quasi-likelihood F-test using the differential expression testing package edgeR(34) comparing the cells with/without the CN change. This test was performed with default parameters except for estimating the overdispersion parameter of the negative binomial distribution where we set prior.df=0 to prevent information sharing between genes. This was done as there were sufficiently many cells in each group to robustly estimate genewise over-dispersion without sharing information.
To test if this list of differentially expressed genes related to sympathoblastic cells we identified all markers of sympathoblastic genes with a tf-idf cut-off of 1 (28). We then performed a hypergeometric test to measure if this list of sympathoblastic marker genes was over-represented in the genes differentially expressed between tumor cells with and without the prognostic CN changes.

Risk group stratification for SEQC data
As COG risk group information was not available for the SEQC dataset, we stratified samples into risk groups based on MYCN amplification status, age at diagnosis, and disease stage. We defined samples as low risk if age at diagnosis was <18 months and MYCN amplification status was negative, excluding 4s stage samples. We defined samples as high risk if age at diagnosis was >18 months and MYCN amplification status was positive, excluding 4s stage samples.

Code Availability
We have included the source code used to generate the Figures and Tables presented in this analysis as Data S1. The purpose of this code is to provide additional explanation of the analyses described in Methods, such as the precise function call and parameter values used.