Single-nuclei transcriptomes from human adrenal gland reveal distinct cellular identities of low and

Childhood neuroblastoma has a remarkable variability in outcome. Age at diagnosis is one of the most important prognostic factors, with children less than 1 year old having favorable outcomes. Here we study single-cell and single-nuclei transcriptomes of neuroblastoma with different clinical risk groups and stages, including healthy adrenal gland. We compare tumor cell populations with embryonic mouse sympatho-adrenal derivatives, and post-natal human adrenal gland. We provide evidence that low and high-risk neuroblastoma have different cell identities, representing two disease entities. Low-risk neuroblastoma presents a transcriptome that resembles sympatho- and chromaffin cells, whereas malignant cells enriched in high-risk neuroblastoma resembles a subtype of TRKB+ cholinergic progenitor population identified in human post-natal gland. Analyses of these populations reveal different gene expression programs for worst and better survival in correlation with age at diagnosis. Our findings reveal two cellular identities and a composition of human neuroblastoma tumors reflecting clinical heterogeneity and outcome.

Neuroblastoma (NB) is a pediatric cancer arising from the sympathoadrenal cell lineage frequently originating in the adrenal glands (AG)1. This malignancy represents 8–10% of all childhood cancer cases, and is responsible for 15% of all pediatric oncology deaths worldwide2. A clinical hallmark of neuroblastoma is heterogeneity, featuring outcomes ranging from lethal progression to spontaneous regression. The risk classification predicting the clinical behavior of the malignancy and its response to treatment, utilizes the INRGSS criteria (i.e., International Neuroblastoma Risk Grouping Staging System)1,3. One of the most significant and clinically relevant factors for this risk classification is age. Children younger than 18 months at the time of diagnosis display better prognosis (i.e., low-risk) than children diagnosed at a later age, and aging is in turn associated with poorer outcome (i.e., high-risk)4,5. Other prognostic markers are used to assign patients to specific risk groups, for example, ploidy, chromosomal alterations, MYCN amplification, and expression of neurotrophin receptors, such as TRKB (encoded by NTRK2) associated with high-risk and poor outcome. In contrast, neurotrophin receptor TRKA expression (encoded by NTRK1) is associated with low-risk and favorable outcome2. The reason why the age of the patient at the time of diagnosis is one of the strongest predictor of risk and outcome is not understood.

Previously, it has been reported that the majority of mouse chromaffin cells forming the adrenal medulla originate from an embryonic neural-crest progeny, specifically, from multipotent Schwann cell precursors (SCPs). SCPs are nerve associated cells that migrate along the visceral motor nerves to the vicinity of the developing adrenal gland, and form ~80% of the chromaffin cells6. The remaining 20% is directly derived from a migratory stream of neural-crest cells (NCCs) that commit to a common sympathoadrenal lineage in proximity to the dorsal aorta7,8, and is considered as the source of neuroblastoma9,10. SCPs also give rise to paraganglia during mouse embryonic development, such as chromaffin cells in Zuckerkandl’s organ (ZO) and to some sympathetic neurons11. In mouse, the ZO reaches a maximum cell number shortly before birth, in contrast to human where the peak of cell number is reached around the third year of life, indicating species-specific developmental differences12. However, SCPs are retained for a rather short time during mouse embryonic development and disappear around E156. Thus, it remains unknown how human chromaffin cells are produced and regenerated after birth and if any distinct population of cells serves as their progenitor pool. We hypothesize that postnatal chromaffin cells are derived from a different source than their developmental origin (namely SCPs), and aberrations in this cell source could explain the age-dependent risk stratification of neuroblastoma.

Here we deep-sequence full-length coverage RNA from single nuclei of tumors (n = 11) across different risk groups and cross-compare cell clusters transcriptomes with healthy postnatal adrenal gland (human n = 3, mouse n = 5). We further include recently published single-cell sequencing datasets from 10X single-sequenced NB tumors (n = 8)13, E12–E13 embryonic mouse adrenal anlagen6, human fetal adrenal gland14, and the transcriptional profiles of neuroblastoma mesenchymal-/NCC-like and (nor-)adrenergic cell lineages15,16. We identify a cluster of TRKB+ cholinergic cells unique to human postnatal adrenal gland, that differ from previously described embryonic Schwann cell precursors (SCP). This TRKB+ population of cells shares a specific gene signature with a cluster of undifferentiated cells of mesenchymal nature enriched in high-risk neuroblastoma. The gene signature of these mesenchymal cells is in turn coupled with lower patient survival probability and older age-at-diagnosis when tested in a larger cohort of 498 neuroblastoma patients17. Conversely, more differentiated noradrenergic cells are over-represented in low-risk cases, and share specific gene signatures with adrenal human and mouse transcriptomes that resemble the signature of sympatho- and chromaffin cells. Our results suggest that high-risk neuroblastomas are characterized by a population of progenitor cells that resemble a cell type in postnatal adrenal gland with migratory and mesenchymal signatures, while the low-risk neuroblastoma resembles postnatal and developing chromaffin cells and sympathoblasts.

To understand why high-risk neuroblastomas arise in children older than 18 months, we first cataloged normal cell populations in postnatal adrenal gland which is a common location for this pediatric malignancy. We define the identity of normal cell populations in the adrenal gland of both postnatal human (n = 3, 1536 single nuclei) and mouse (n = 5, 1920 single whole-cells) by single- nuclei/cell RNA-sequencing (SmartSeq2, see the “Methods” section) to an average depth of 485,000- and 669,000 reads per nuclei/cell, respectively (Supplementary Fig. 1a, b, Supplementary Data 1). Technical and biological features of the transcription profile acquired from nuclei (human) or from whole cells (mouse) are summarized in Supplementary Fig. 1c.

Postnatal human and mouse adrenal glands share cell populations but exhibit differences in chromaffin cells

The cell populations in human and mouse adrenal glands were annotated under the expectation of recovering both adrenal cortex- and medulla-associated cells (Fig. 1a for human, 1b for mouse, Supplementary Data 2). A reference guide of normal adrenal cell populations was generated by assigning an identity to each cluster, by cross-referencing significantly upregulated transcripts with canonical markers curated from the literature (Supplementary Data 2). Adrenal medulla cells were identified by the expression of a panel of nor- and adrenergic markers, including PNMT, TH, DBH, CHGA, and CHGB (Fig. 1d, Supplementary Data 2). In human, hC4 was identified as the chromaffin cell cluster (“NOR” panel were significantly upregulated, FDR < 0.01, Welch’s t-test, Fig. 1a, d, Supplementary Data 3, whereas in mouse two chromaffin cell clusters (i.e., mC11 and mC15) were identified (noradrenergic markers exemplified in “NOR” panel were significantly upregulated, FDR < 0.01, Welch’s t-test, Fig. 1b, d, Supplementary Data 2, Supplementary Data 4). Nevertheless, the mouse chromaffin population mC15 shared a more significant specific gene signature with the human chromaffin cluster hC4 (FDR < 0.01, Fisher’s exact test, Fig. 1c, Supplementary Data 5). To understand the differences between the two mouse chromaffin populations (mC11 and mC15), we investigated the expression of genes that differed significantly between them. The expression of CHGA, CHGB, PHOX2A, and PHOX2B were significantly higher in mC15 than in mC11 (FDR < 0.01, Welch’s t-test), while the expressions of PNMT was higher in mC11 than in mC15 (FDR < 0.01, Welch’s t-test, Fig. 1d, Supplementary Data 2). In addition, mC15 cluster exhibited a significantly higher expression of a different repertoire of cholinergic muscarinic and nicotinic receptors (mAChR and nAChR) than mC11 including CHRM1, CHRNA3, CHRNA7, and CHRNB4 (FDR < 0.01, Welch’s t-test, Supplementary Fig. 2a, b, Supplementary Data 2). In contrast to mouse chromaffin cells, human postnatal PNMT+ chromaffin cells (cluster hC4) showed a significant expression of the sympathoblast marker PRPH (FDR < 0.01, Welch’s t-test, Fig. 2a).

a 1536 single nuclei of human AG from three different patients were sequenced with Smart-seq2 to an average depth of 485,000 reads per cell. Cells with high-quality (n = 1322) were selected and further processed with PAGODA. Cells were grouped into ten different clusters, including cortex (in brown and gray colors), chromaffin (blue hC4), mesenchymal (purple hC7), endothelial (light blue hC6), and immune cells (i.e., T cells hC10 and macrophages hC2 in green colors). b 1920 single cell of mouse AG from five different samples were sequenced with Smart-seq2 to an average depth of 670,000 reads per cell. Cells with high-quality (n = 1763) were selected and further processed with PAGODA. Cells were grouped into nineteen different clusters, including cortex (in brown and gray colors), chromaffin (blue mC15 and light blue mC11), mesenchymal (red mC6), capsule (purple mC13), endothelial (light blue mC2 and mC4), glial (orange mC10), and immune cells (i.e., T cells mC17 and mC18, and macrophages mC3, mC8, and mC14 in green colors). Human and mouse adrenal glands zonation49 is illustrated in (c). c, d A comparison of the specific gene signature between human and mouse revealed similar transcription signatures for mesenchymal, endothelial, and immune clusters, nevertheless, the two different postnatal chromaffin cells could only be differentiated in mouse. d Hierarchical clustering of cells and expression of a panel of markers including noradrenergic (NOR), mesenchymal (MES), endothelial (Endot.), and migratory (Migrat.). e–g A population of cells with a significant high expression of progenitor markers (i.e., SOX6, ERBB3, RTTN, FDR<0.01) and high differentiation potential found uniquely in human AG, is sourcing the chromaffin cells as indicated by velocity, entropy, and pseudotime analyses. h In this process, gene expression elapses from an undifferentiated stem-like- (i.e., RTTN+) to an adrenergic signature (PNMT+), passing by a noradrenergic stage (i.e., DBH+), as indicated by the pseudotime of the underlying cellular process. Results in (c), (e), and the top inserts in (d) and (h) represent cell clusters by color in (a).

Full size image

a tSNEs representing the expression of indicated genes in human cholinergic progenitor (pink) and chromaffin (blue) populations. The bars next to the tSNEs illustrate the expression measured as the logarithm of the read counts per 10,000. b–e Overview of tile-scanned images (×20) of postnatal human adrenal glands (AG) at indicated age. Scalebar of overview: 200 μm, zoom of boxed image: 10 μm. b–d RNAscope in situ hybridization (ISH) for TH (green), CLDN11 (red), and NTRK2 (white) mRNA and counter stained with DAPI (blue). NTRK2+ CLDN11+ double positive cells were found in adrenal capsule and medulla exclusive from TH positive cells. e RNAscope ISH of a 4-year-old AG labeled with for TH (green), CHRNA7 (red), and CLDN11 (white) mRNA and nuclear counter-stain (DAPI) as indicated. For all RNAscope experiments, the signal distribution patterns and cell morphological features were shown by the different combination of probes and independently reproduced three times on different samples.

Full size image

Identification of a cholinergic progenitor population in the human postnatal adrenal gland

Unexpectedly, we found a population of cells unique to the human postnatal adrenal gland (hC1) which significantly expressed various progenitor and migratory markers, including SOX6, BCL11A, ERBB3, NTRK2 (encoding for TRKB), RTTN, PTPRZ1, TP63, ASXL3, POU6F2, SEMA3E, LAMA3, DOCK7, and CLDN11 (FDR < 0.01, Welch’s t-test, “Progenitors” and “Migrat.”, highlighted by a dotted line, Figs. 1d and 2a, Supplementary Data 2 and Supplementary Data 3). Furthermore, hC1 progenitor population presented a proliferating nature and significantly expressed cell cycling genes such as MKI67, ASPM, and BUB1 (“Cell Cycle” panel was significantly upregulated FDR < 0.01, Welch’s t-test, Figs. 1d and 2a and Supplementary Data 3). These progenitors, however, did not express previously described multipotent Schwann cell precursors (SCPs) markers such as SOX10, FOXD3, and S100B (Fig. 1d, Supplementary Data 2). In support of this finding, previously described SCPs from mouse adrenal anlagen at embryonic day E12/E136, and from human fetal adrenal glands at 8–14 post-conception weeks (PCW)14 shared no significant specific gene signature with the human postnatal progenitor cluster hC1 (Fisher’s exact test, Supplementary Fig. 2c, d, Supplementary Data 5). No progenitor cells other than SCP, chromaffin, and sympathoblast populations have been described for fetal adrenal gland13,14. In this regard, the only gene signatures shared with human fetal adrenal gland and mouse embryonic anlagen belonged to cell-cycle genes expressed in cycling sympathoblast (mouse E13 and human 8–14 PCW) and cycling chromaffin (human 8–14 PCW, FDR < 0.05, Fisher’s exact test, Supplementary Data 5). An estimation of cell velocities (i.e., computational reconstructions of cells trajectory and fates that uses transcript splicing to calculate the direction and speed of differentiation)18,19 suggest that this type of cells (i.e., hC1) repopulates chromaffin cells in postnatal human adrenal gland (Fig. 1e). A velocity-driven gene trajectory analysis using pseudotime indicates that the progenitor cells transit from precursors cells with high differentiation potential to differentiated chromaffin cells (Fig. 1f–h). Furthermore, both progenitor (hC1) and chromaffin population (hC4), express the nicotinic acetylcholine receptor nAChRα7 (CHRNA7), suggesting that progenitor cells are cholinergic in nature (FDR < 0.01, Welch’s t-test, Supplementary Data 2).

To validate the expression and spatial context of the human cholinergic progenitor cells (hC1), a series of RNAscope in situ hybridizations (ISH) was performed in postnatal adrenal glands from 0- and 4-year-old children and adults (Fig. 2b–e, Supplementary Figs. 2e and 3). We first elucidated the anatomy of the glands in each section for medulla (TH), cortex (CYP11B2), and capsule (RSPO3) (Supplementary Fig. 2e). To identify cells belonging to the hC1 progenitor population, we tested markers identified as significantly upregulated in hC1 population: NTRK2 and CLDN11 (Fig. 2b–d and Supplementary Fig. 3a–c). NTRK2+CLDN11+double positive cells were found in human adrenal capsule and medulla exclusive from TH+cells at all ages, but most abundantly at 0 year of age, suggesting that these cells decline over age. To confirm that the postnatal progenitor CLDN11+ cells are of cholinergic nature, we performed in situ RNA-hybridizations with the cholinergic nicotinic receptor CHRNA7(i.e., nAChR α7) together with TH and CLDN11. Both cell types, TH+ chromaffin and CLDN11+progenitor cells, express nAChR α7 (Fig. 2e). To further confirm that the postnatal progenitor cells in human postnatal gland are different from the previously described embryonic multipotent SCPs, we performed in situ RNA-hybridizations together with SCP/glia marker SOX10 (Supplementary Fig. 3d–f). NTRK2+ cells were exclusive from SOX10+ cells in human postnatal adrenal gland at all ages.

Different neuroblastoma risk groups present differences in cell population composition

Next, we used single-nuclei transcriptomics to characterize eleven neuroblastoma samples across different clinical risk groups, and genetic subsets (Fig. 3a, b, Supplementary Data 1). Deep-frozen samples obtained from surgical resections were used for nuc-Seq (SmartSeq2), yielding a total of 3212 high-quality nuclei with an average of 709,676 high-quality reads per nuclei (see “Methods” and Supplementary Fig. 1). Cluster analysis using PAGODA identified ten cell populations, classified as undifferentiated (nC2 and nC3), noradrenergic (NOR clusters nC5, nC7, nC8, and nC9), and stroma clusters: mesenchymal stroma (MSC nC1), endothelial (nC4), macrophages (nC6), and T cells (nC10) (Fig. 3a, Supplementary Data 6). The identity of each cluster was assigned by cross-referencing cluster-defining transcripts with canonical markers curated from the literature (Supplementary Data 2). Clusters nC2 and nC3 (referred as “undifferentiated”), presented a significant high expression of progenitor markers PROM1, RTTN, ERBB3, POU6F2, and the migratory marker CLDN11. Remarkably, the undifferentiated nC3 clusters presented a significantly high expression of MYCN, ALK, BRCA1, and BRCA2 genes and progenitor markers BCL11A, NTRK2, SOX5, SOX6, TP63, LGR5, and USH2A (FDR < 0.01, Welch’s t-test, Fig. 3c and Supplementary Data 2, Supplementary Data 6). Oppositely, the noradrenergic clusters (NOR, nC5, nC7, nC8, and nC9) expressed significantly high NTRK1 (TRKA) and noradrenergic markers TH, DBH, PHOX2A, PHOX2B, and ISL1 (FDR < 0.01, Welch’s t-test, Fig. 3c, d and Supplementary Data 2, Supplementary Data 6).

4224 single nuclei of human NB tumors were sequenced to an average depth of 676,000 reads per cell. a Cells with high-quality (n = 3212) were grouped into ten different clusters with PAGODA. NB risk groups, INSS stages, and samples are exemplified in (b). c Representation of the hierarchical clustering of cells: undifferentiated (light pink nC2 and pink nC3), mesenchymal stroma (MSC red nC1), endothelial (light blue nC4), T-cells (light green nC10), macrophages (green nC6), and noradrenergic (“NOR”, blue nC5, nC7, nC8, and nC9) and various markers of interest (also in insert d). The bars next to the tSNEs in (d) illustrate the expression measured as the logarithm of the read counts per 10,000. e, f Using a gene enrichment-based approach, the noradrenergic clusters of NB (nC5, nC7, nC8, and nC9) proved to share a significant number of upregulated genes with sympathetic noradrenergic15 and “ADR” (i.e., adrenergic)16 transcriptional signatures. Oppositely, the undifferentiated cluster nC3 shared a significant number of upregulated genes with neural-crest cell-like signature15 and MES (i.e., mesenchymal)16 signatures. g A comparison of the gene-specific signature of NB clusters with markers defining clusters in eight tumors previously sequenced with 10X13, suggests a more significant similarity of the noradrenergic clusters (nC5, nC7, nC8, and nC9) with tumor clusters 1, 2, and 3; and of the undifferentiated nC3 cluster with the tumor cluster 313. Only clusters with the most significant similarities are displayed. h An analysis of the CNVs indicates a significant number of cells with rearrangements (indicated by red, Fisher exact tests) in the undifferentiated cluster nC3 for most samples, while a remarkable number of rearrangements was observed in fewer samples for the NOR (nC5), MSC (nC1), and undifferentiated nC2 clusters. i A subset of the predicted CNVs was validated for each sample, and illustrated in the tSNE for cells with an expected rearrangement (i.e., expected state) overlapping the validated region. Colors displayed in (e), (f), and (g) represent cell clusters by color in (a). An arrow next to the bar indicates a FDR < 1 × 10−50.

Full size image

Neuroblastoma samples of different clinical risk groups and stages contained different contribution of cells within these clusters reflecting the clinical heterogeneity of these clinical groups. High-risk cases had a significantly higher proportion of cells in stroma and undifferentiated clusters, whereas the low-risk ones including spontaneous regressing (4S) cases, consisted mostly of cells belonging to the noradrenergic (NOR) clusters (FDR < 0.01, Chi-square test, Fig. 3b, Supplementary Fig. 4). In agreement with these observations, a deconvolution of 172 neuroblastoma bulk-sequenced samples (NB172 NCI TARGET project) indicated that a significantly higher proportion of cells in the stroma and undifferentiated nC3 clusters was present in high-risk patients (n = 139), and a higher proportion of cells in the noradrenergic clusters nC7, nC8, and nC9 in low-risk cases (n = 14, FDR < 0.01, Chi-square test, Supplementary Fig. 4). Some differences observed between single-nuclei and bulk-seq for nC2 and nC5 clusters could be a consequence of the clinical heterogeneity within patients of the same risk groups (Supplementary Fig. 4).

Mesenchymal markers PRRX1, YAP1, and PDGFRA were significantly expressed in both clusters, undifferentiated nC3 and MSC nC1 clusters (FDR < 0.01, Welch’s t-test, Fig. 3c, d and Supplementary Data 2), in contrast to PDGFRB whose upregulation was significantly higher in MSC cells (nC1, FDR < 0.01, Welch’s t-test, Fig. 3c, d and Supplementary Data 2). Consistently with this observation, the undifferentiated nC3 cluster shared a significant number of upregulated genes with the neural-crest cells-like signature (group II) and the mesenchymal (MES) gene expression signatures previously described in neuroblastoma cell lines15,16, while the noradrenergic clusters nC5, nC7, nC8, and nC9 presented a significantly high expression of genes characterizing the sympathetic noradrenergic (group I) and ADR (adrenergic) signatures (FDR < 0.01, Fisher’s exact test, Fig. 3e, f, Supplementary Data 7). In addition, the specific gene signatures from each neuroblastoma cluster were compared with the recently reported markers for six cell clusters in eight GOSH-cohort tumors sequenced with 10X13. The noradrenergic clusters nC5, nC7, nC8, and nC9 resembled more significantly the GOSH sympathoblast-like tumors clusters 1, 2, and 3, whereas the undifferentiated cluster nC3 showed a higher significant resemblance with the less differentiated GOSH-tumor cluster 3 (FDR = 2.87 × 10−36, Fisher’s exact test, Fig. 3g, Supplementary Data 7).

To differentiate malignant- from stroma cell clusters, an analysis of the genome rearrangements was conducted, and further experimentally validated. Significant copy number changes were computed for each sample using as background the changes observed in immune cells (i.e., T cells nC10 and macrophages nC6, only in the seven samples with immune cells as detailed in “Methods”). For most of samples, a significantly higher number of cells presented gains in 17q, losses in 1p or in 11q, in the undifferentiated nC3 cluster (K87, K10, 23, K55, and K3). Fewer samples presented a remarkable number of cells with rearrangements in the MSC nC1 (K87), undifferentiated nC2 (23, K3, K6, and K87), and NOR nC5 clusters (K6 and K87, FDR < 0.05, Fisher’s exact test, Fig. 3h). 1p deletion was confirmed by microarrays for samples 23, K10, and K55, 11q loss for samples 23 and K87, and 17q gains for samples K10, K3, and K87. Cells with predicted significant gains and losses from these samples are illustrated in Fig. 3i.

Next, we validated the expression of significantly enriched markers identified in neuroblastoma clusters (Fig. 4, Supplementary Fig. 5). RNAscope in situ hybridizations in a MYCN-amplified neuroblastoma case (K10, Supplementary Data 1) as independently confirmed by comparative genomic hybridization (CGH)20, revealed a high intratumoral heterogeneity for TH, MYCN, and ALK expression, with tumor regions with intense staining of MYCN and ALK, negative for TH (Fig. 4a, zoom 2). Other tumor regions revealed some diffuse stainings of few TH+ positive cells with enlarged nuclei positive for MYCN, but negative for ALK (Fig. 4a, zoom 1), whereas other regions did not exhibit expression of TH, MYCN, and ALK (zoom 4). Interestingly, in situ hybridizations for NTRK1 and NTRK2(encoding TRKA and TRKB, respectively) revealed that the identified TH+ cells with enlarged nuclei are NTRK1+ but NTRK2−, in contrast to all surrounding NTRK1− NTRK2+ cells with small nuclei (Fig. 4b). It is possible that the identified TH+ cells with enlarged nuclei (TH+, MYCN+, NTRK1+, NTRK2−, ALK−) might account for cells undergoing differentiation21. Expression of neurotrophin receptor NTRK1 is characterized as a marker for favorable outcome and low-risk, while NTRK2 expression is associated with poor outcome and high-risk2. In agreement with these findings, cells for the favorable low-risk neuroblastoma case (K6, stage 4S, Supplementary Data 1) were homogeneously NTRK1+, TH+, NTRK2− (Supplementary Fig. 5a).

Overview of tile-scanned images (20x) high-risk neuroblastoma (K10, MYCN amplified) using RNAscope in situ hybridization. Scalebar of overview: 500 μm; zoom of boxed image: 10 μm. aTumor labeled with RNAscope ISH for TH (green), MYCN (red), and ALK (white) mRNA, and counter stained with DAPI (blue). Dashed circles indicate cells with large nuclei in a region (box #1), that are TH and MYCN positive, but negative for ALK. Box #2 indicates part of the tumor with majority of cells double positive for MYCN and ALK but negative for TH. Box #4 indicates cells negative for all probes: TH, MYCN, and ALK. b Adjacent section from (a) labeled for TH (green), NTRK1 (red), and NTRK2 (white). Box #1 indicate cells with large nuclei in a region that are TH+ and NTRK1+, and negative for NTRK2. Surrounding cells with small nuclei are positive for NTRK2 only. Box #2,3 visualizes tumor region with majority of cells positive for NTRK2 that are negative for TH and NTRK1. c Adjacent section stained for PDGFRA (green), LGR5 (red), and NTRK2 (white) mRNA is highlighting cells double positive for LGR5 and NTRK2 (box #1) or double positive for PDGRFA, NTRK2 (box #2). Box #3 indicates a region of the tumor that is positive for NTRK2 only. d Adjacent section stained for mesenchymal markers PDGFRA (green), CLDN11 (red), and LGR5 (white). Similar to c: some tumor regions (box #1) highlight cells double positive for CLDN11 (red) and LGR5 (white) that are negative for PDGFRA (box #1), as were other regions show cells that are double positive for PDGFRA (green) and CLDN11 (red) (box #2, 3, and 4). For all RNAscope experiments, the signal distribution patterns and cell morphological features were shown by the different combination of probes and independently reproduced three times on different samples in (a) and (b), and independently reproduced four times on different samples in (c) and (d).

Full size image

In addition, we validated the expression of the mesenchymal marker PDGFRA that was identified to be significantly highly expressed in high-risk neuroblastoma cluster (nC3). PDGFRA expression in the high-risk MYCN-amplified case (i.e., K10) was observed in both NTRK2+ (Fig. 4c zoom 2), and also in CLDN11+ cells (Fig. 4d, zoom 2, 3, and 4), whereas in LGR5+ cells showed no PDGFRA signal (Fig. 4c zoom 1 and Fig. 4d zoom 1). In contrast, the low-risk stage 4S case was entirely negative for PDGFRA and PRRX1 expression, however, homogeneously expressed in all cells noradrenergic markers DBH and PHOX2B(Supplementary Fig. 5b, c). Stage 2B neuroblastoma was more heterogeneous in DBH and PHOX2B expression (Supplementary Fig. 5d), with all cells lacking expression of PRRX1 and very few cells positive for PDGFRA (Supplementary Fig. 5e).

A cell population enriched in high-risk neuroblastoma resembles the human postnatal adrenal progenitor population

High-risk neuroblastoma commonly arises in children older than 18 months and frequently originates in the adrenal glands1. To understand how the observed neuroblastoma populations in high and low-risk cases are related to normal postnatal and embryonic adrenal tissues, we performed a comparative analysis of their transcriptional profiles (Figs. 5, 6, and 7). Significant expression of (nor-)adrenergic (“NOR”) markers (i.e., CHGA, CHGB, DBH, TH, PHOX2A, and PHOX2B) is common for healthy sympatho-adrenal cells and tumor NOR populations identified in low-risk neuroblastoma cases (FDR < 0.01, Welch’s t-test, Fig. 5a, Supplementary Data 2). In agreement with this observation, specific gene signatures from neuroblastoma NOR clusters were found to be significantly shared with both, mouse embryonic (E12–E13) and human (8–14 PCW) chromaffin and sympathoblast clusters, and also with human and mouse postnatal chromaffin clusters (FDR < 0.01, Fisher’s exact test, Fig. 5b–e, Supplementary Data 8, Supplementary Fig. 6a). A more significant resemblance was observed between the NOR nC7 and nC9 clusters with fetal (8–14 PCW) sympathoblast cells (FDR = 4.47 × 10−33 and FDR = 5.38 × 10−32, respectively, Fisher’s exact test), while the NOR nC8 cluster presented a closer similitude to fetal (8–14 PCW) chromaffin and sympathoblast clusters (FDR = 5.45 × 10−18 and FDR = 1.08 × 10−18, respectively, Fisher’s exact test, Fig. 5e, Supplementary Data 8). A comparative analysis with the mouse embryonic transcriptional profiles showed significant gene signatures shared between NOR nC7 and nC9 clusters with both, mouse embryonic chromaffin and sympathoblast populations (FDR < 0.01, Fisher’s exact test, Fig. 5b, c, Supplementary Fig. 6a, Supplementary Data 8). In support of these findings, neuroblastoma NOR clusters presented a high signature score (measuring transcriptional resemblance, see “Methods”) for bridge, sympathoblast and chromaffin cells in mouse and human developing adrenal glands, and for human postnatal chromaffin cells (Fig. 6a–c). Furthermore, the significant specific gene signatures of NOR nC7 and nC8 clusters were significantly associated with a better outcome in a large neuroblastoma cohort (i.e., Bonferroni-corrected p-value < 0.01, 498 cases in SEQC17, see “Methods”, Fig. 6d).

a Heatmap illustrating the normalized expression magnitude for selected genes organized following the hierarchical clustering (PAGODA) as shown in Figs. 1 and 3. b–d Venn diagrams illustrating the significantly shared specific gene signature. Shared genes are listed in Supplementary Data 8. Indicated FDRs were calculated with a Benjamini–Hochberg correction on p-values obtained with Welch’s t-tests as detailed in the “Methods”. b Venn diagrams of human neuroblastoma clusters compared to adult mouse postnatal clusters, c embryonic mouse clusters previously described6, and d human postnatal adrenal gland clusters. e Comparison of the specific signature from the human neuroblastoma clusters and the reported transcriptional signal for human fetal adrenal glands cell clusters (Dong et al.14). An arrow next to the bar indicates a FDR < 1 × 10−50. +The original cluster annotations are currently debated and the included labels here correspond to those given by Kildisiute et al.45, and Bedoya-Reina and Schlisio46. Colors displayed in (b), (c), (d), and (e) represent cell clusters by color in (a).

Full size image

tSNE of neuroblastoma clusters illustrating the signature score of genes (as detailed in “Methods”) characterizing a human postnatal, b embryonic14 adrenal glands, and c mouse embryonic adrenal anlagen (E13) cell clusters6. The bars next to the tSNEs in a, b, and c illustrate the scale of the signature score. A greater signature score indicates a larger transcriptional resemblance to the cluster of interest. SCP refers to the cluster of mouse multipotent Schwann cell precursors at E13; and Bridge refers to the cell cluster that defines transiting cells from SCP toward the chromaffin population at E13. d Kaplan–Meier curves for neuroblastoma clusters with significant differences (Bonferroni-corrected [bonf. p], log-rank tests in survival for 498 SEQC neuroblastoma patients with a low (red) and high (blue) signature gene expression17. e–g Using a gene enrichment-based approach (Benjamini–Hochberg corrected, Fisher’s exact tests), the specific signature genes for the noradrenergic nC7, nC8, and nC9, endothelial nC4, macrophages nC6, and undifferentiated nC3 clusters were found to be significantly upregulated in low- and intermediate (inter.)-risk patients and also in individuals with better survival. In contrast, the signature genes of noradrenergic nC9 and undifferentiated nC3 clusters presented a significantly enrichment in patients with high-risk and poor survival. The signature genes of the undifferentiated nC3 cluster presented also a remarkable significant correlation with age at diagnosis. h Signature genes associated with poor survival for noradrenergic nC9 and undifferentiated nC3 clusters are significantly correlated with age at diagnosis, and oppositely those associated with better survival are inversely correlated with age. Genes are included in Supplementary Data 9. A solid arrow next to the bar indicates a FDR < 1 × 10−50. Colors of bars displayed in (e), (f), (g), and (h) represent the neuroblastoma cell clusters by color in (a).

Full size image

a tSNE illustrating the projection of postnatal human adrenal gland cell clusters (i.e., states, reference top right) onto neuroblastoma (query top left), based on the detection of mutual nearest neighbor cell anchors (as detailed in “Methods”). Joint visualization of the neuroblastoma cell clustering (bottom) after cell classification with colors corresponding to the transferred adrenal gland cell states. b Heatmaps of the average probabilities of best matching reference states for each query cluster (i.e., average cell state probability, top), including number of cells (bottom). c tSNE using Euclidean distances in the quartile-normalized matrix of the average normalized expression from single-cell/nuclei sequenced cell populations from (A) human (i) neuroblastoma and (ii) adrenal gland, and (B) mouse (i) adrenal gland and (ii) derived adrenal anlagen at E12 and E136. Mesenchymal, endothelial, and neural-crest derived populations from human and mouse, adrenal gland, and neuroblastoma group accordingly to gene expression in (i) adrenergic and noradrenergic cells, (ii) undifferentiated cells, (iii) mesenchymal, and (iv) endothelial cells. dHierarchical clustering of transcriptional similarities between cell clusters supported with unbiased p-values41. Green numbers in the branches indicate a high support (i.e., unbiased p-value > 95) that has a lowest support at 0 and a highest at 100.

Full size image

In contrast to the neuroblastoma NOR cells, the undifferentiated neuroblastoma clusters (nC3, nC2) did not share specific gene signatures with mouse embryonic and postnatal sympathoblast nor chromaffin populations (Fig. 5b–e). The nC3 cluster, however, shared a marginally significant number of specific-signature genes with cycling populations in the human fetal gland (i.e., cycling sympathoblast and chromaffin, FDR < 0.05, Fisher’s exact test, Fig. 5e). These were nonetheless only cell-cycle related genes (i.e., ASPM, BUB1, CENPE, CLSPN, and ESCO2, Supplementary Data 8).

Importantly, and in contrast to the fetal adrenal gland comparison, the neuroblastoma nC3 cluster shared a specific gene signature with the human cholinergic progenitor population (hC1) identified in postnatal human adrenal gland (FDR < 0.01, Fisher’s exact test, Fig. 5d, Supplementary Data 8). In addition, nC3 and hC1 shared a significantly high expression (FDR<0.01, Welch’s t-test) of NTRK2, mesenchymal (i.e., COL1A2, COL6A3, COL12A1), migratory (i.e., CLDN11, DOCK7), and progenitor genes (i.e., BCL11A, ERBB3, RTTN, TP63, ASXL3, POU6F2, and SOX6, Supplementary Data 2, Supplementary Data 3 and 6). Furthermore, the undifferentiated neuroblastoma cluster (nC3) constituted a larger proportion of cells in high-risk samples in both the 11 nuc-Seq-samples and in a deconvolved cohort of 172 patients (FDR < 0.01, Chi-square test, Supplementary Fig. 4), and its specific-signature genes have a marginally higher expression average in patients with a lower survival probability in a larger cohort (i.e., Bonferroni-corrected p-value < 0.05, 498 cases in SEQC17, see “Methods”, Fig. 6d).

Recently, embryonic multipotent Schwann cell precursors (SCP) have been suggested to be a potential source of neuroblastoma origin6. However, embryonic SCPs in both, human fetal and mouse, shared only specific gene signatures with neuroblastoma stroma clusters MSC nC1 and endothelial nC4 (Fig. 5c, d, Supplementary Fig. 6a, Supplementary Data 8). A detailed look into their shared gene signatures suggests a high expression of genes related to cell adhesion and motility (i.e., GO term, FDR < 0.05, Fisher’s exact test) and did not include SCP lineage markers (i.e., SOX10, S100B, FABP7, or PLP1, Supplementary Data 8).

We further analyzed in detail how the transcriptional profiles of the different neuroblastoma clusters were associated with the survival and age at diagnosis in a larger cohort (i.e., 498 cases in SEQC17, Fig. 6e–h, Supplementary Data 9). We compared the specific gene signatures from each neuroblastoma cluster with genes in 498 patient samples with (1) differential expression by risk groups (Fig. 6e), (2) significant correlation with age-at-diagnosis (Fig. 6f), and (3) differential expression by survival outcome (Fig. 6g as detail in “Methods”). Genes in the specific signatures from endothelial nC4, macrophages nC6, undifferentiated nC3, and NOR clusters (nC7, nC8, and nC9) were significantly enriched in both low/intermediate risks and better survival cases (FDR < 0.01, Fisher’s exact test, Fig. 6e–g, Supplementary Data 9). Nevertheless, genes in the specific signatures from undifferentiated nC3 and NOR nC9 clusters were also significantly enriched in both high-risk and worst survival patients (FDR < 0.01, Fisher’s exact test, Fig. 6e–g, Supplementary Data 9). Remarkably, a more significant directly correlation was observed for the specific gene signature of the undifferentiated nC3 cluster with age at diagnosis (i.e., they are more likely to be expressed in later ages, FDR = 8.12 × 10−31, Fisher’s exact test), while a more significant inversely correlation was found for the NOR nC9 cluster (FDR < 1 × 10−50, Fisher’s exact test, Fig. 6f, Supplementary Data 9).

A detailed analysis of the specific gene signatures for the undifferentiated nC3 and NOR nC9 clusters revealed that the two cell populations feature different gene programs that can be either (1) directly correlated with age at diagnosis and associated with poor survival (more significant in nC3 than in nC9), or otherwise (2) inversely correlated with age at diagnosis and resulting in better survival (more significant in nC9 than in nC3, FDR < 0.01, Fisher’s exact test, Fig. 6h, Supplementary Data 9). The genes inversely correlated with age at diagnosis, and also those resulting in better survival, are significantly associated with cell adhesion and differentiation for both nC3 and nC9 clusters (GO term, FDR < 0.01, Fisher’s exact test, Supplementary Fig. 6d–g, Supplementary Data 10). Oppositely, genes in NOR nC9 associated with a later age at diagnosis and worst survival are significantly enriched in translation, RNA processing, and splicing, while genes in the undifferentiated nC3 cluster associated with worst survival are enriched in DNA damage and repair and cell cycle, and those associate with a later age of diagnosis are enriched in cell motility (GO term FDR < 0.01, Fisher’s exact test, Supplementary Fig. 6d–g, Supplementary Data 10).

We further investigated the commonalities of human postnatal adrenal gland and neuroblastoma clusters, and performed a comparative analysis based on the detection of matching mutual nearest neighbors22,23 (Fig. 7a, b). We used human adrenal gland for the projection of cell states onto neuroblastoma data (i.e., query). In agreement with the identified significant shared gene signatures (shown in Fig. 5d), human postnatal chromaffin (hC4) population matched to the malignant NOR clusters (nC5, nC7, nC8, and nC9), whereas the adrenal gland progenitor population (hC1) matched the malignant undifferentiated neuroblastoma clusters nC2 and nC3 (Fig. 7a, b). Populations of endothelial, macrophages and mesenchymal cells (i.e., stroma cells) in neuroblastoma have a large number of cells with a high state probability of (i.e., best transcriptional match using joined labeling) to similar populations in human adrenal gland.

In addition, we interrogated the relationships between embryonic mouse adrenal anlagen, neuroblastoma, and healthy human postnatal adrenal cell populations. We performed an Euclidean distance-based analysis of the transcriptional similarities between all populations, excluding cortex and immune cells (Fig. 7c, d). This analysis suggested that these populations could be clustered into four groups, namely (1) (nor-) adrenergic, (2) undifferentiated, (3) mesenchymal, and (4) endothelial cells (Fig. 7c). Consistently with the identified significant shared gene signatures (shown in Fig. 5b–d), a statistical analysis of these Euclidean distances (Fig. 7d) indicated that (1) SCPs, glial, endothelial, and mesenchymal cells from human and mouse adult and embryonic E12–E13 populations share a similar transcriptional profile, (2) human adrenal progenitor cluster hC1 resemble that of the neuroblastoma undifferentiated populations nC2 and nC3, and (3) chromaffin cells from human and mouse (postnatal and embryonic) have a similar transcriptional profile, that resemble the neuroblastoma noradrenergic clusters nC5, nC7, nC8, and nC9.

With the aim of understanding the transcriptional basis of the clinical heterogeneity in neuroblastoma, we analyzed 3212 high-quality single cell/nuclei from 11 tumors from patients across different risk groups. We sequenced to an average of ~710,000 reads per nuclei/cell, and compared them to 3456 human and mouse developing and postnatal AG cells sequenced to similar depths. While a larger number of cells could provide a more comprehensive view of the different cell types within the tumors, the large sequencing depth combined with full-length coverage offers higher sensitivity. A recent study compared the results obtained with Smart-Seq2 and 10X Genomics (10X) and concluded that the detection sensitivity was higher in Smart-seq2 than 10X. In particular, the gene capture rate of 20 cells sequenced with Smart-Seq2 is comparable to that of 1000 cells with 10X24. Taking advantage of the enhanced sensitivity offered by Smart-Seq2, we consistently recovered the same populations of cells in neuroblastoma tumors of the same risk groups, and compared them with postnatal and developing mouse and human adrenal glands.

The way chromaffin cells are repopulated in the postnatal adrenal gland remains unknown. While during mouse embryonic development, multipotent Schwann cell precursors (i.e., SCPs)6,11 and sympathoadrenal progenitors7 can give rise to both chromaffin cells and sympathoblasts, their postnatal existence has not been demonstrated. In this study, we identified a population of cells unique to postnatal human adrenal gland that presented precursor features of chromaffin cells, namely (1) the expression of progenitor and migratory markers (i.e., NTRK2, ERBB3, BCL11A, ASXL3, RTTN, SOX6, DOCK7, LAMA3, SEMA3E, and CLDN11), (2) a high differentiation potential, and (3) a computed (by RNA velocity) directionality of the in vivo transitions from progenitor to chromaffin cells. Both human progenitor and chromaffin cells express the cholinergic receptor nAChRs α7 (encoded by the nicotinic acetylcholine receptor gene CHRNA7), suggesting that progenitor cells are of cholinergic nature. nAChRα7 is a critical component for the cholinergic system mediating ligand-gated ion channel activation of calcium-dependent signaling pathways25, and is involved in the response to cortisol and oxidative stress26,27. As this postnatal progenitor population of chromaffin cells does not express the SCP markers SOX10, S100B, nor FOXD3and is not present in the mouse adrenal gland, it is possible that it is unique to humans. RNAscope ISH for markers characteristic of this progenitor population at different ages in human postnatal adrenal glands validated their existence and suggested that this population declines with age.

Comparative analysis of the transcriptional profiles of postnatal human adrenal glands and neuroblastoma tumors of different clinical risk groups revealed that the identified postnatal cholinergic progenitor population (hC1) had a remarkable transcriptional resemblance to the undifferentiated cluster (nC3) characteristic of high-risk neuroblastoma. Both populations, high-risk undifferentiated cluster (nC3) and progenitor population in postnatal gland (hC1) shared a significantly high expression of neurotrophic receptor NTRK2 (encoding TRKB), mesenchymal (i.e., COL1A2, COL6A3, and COL12A1), migratory (i.e., LAMA3, CLDN11, and DOCK7), and progenitor genes (i.e., BCL11A, ERBB3, RTTN, TP63, ASXL3, POU6F2, and SOX6). In contrast, the transcriptional profile of the clusters identified in low-risk neuroblastoma did not resemble the human progenitor population identified in postnatal gland, instead they showed a noradrenergic (NOR) transcriptional signature (i.e., NTRK1, CHGA, CHGB, DBH, TH, and PHOX2B) matching the transcriptome of human and mouse postnatal chromaffin cells, as well as fetal human (8–14 PCW) and mouse embryonic (E13) sympathoblast and chromaffin populations.

Furthermore, our results indicate that the high-risk-associated cell cluster (nC3) has a transcriptional program that changes with age-at-diagnosis and is correlated with lower survival probabilities of patients. This program favors the expression of genes associated with cell motility and metastasis in older patients and with replication stress (i.e., DNA repair and cell cycle) in worst outcome cases. Oppositely, the low-risk-enriched NOR cell clusters presented transcriptional programs associated with greater survival probabilities and younger ages at diagnosis. Within these clusters, the NOR nC9 was an exception as it also showed a transcriptional program that presented a worst outcome in correlation with age. Remarkably, this program is associated with RNA-(mis-)splicing and not with replication stress. Our results suggest that at a younger age (at time of diagnosis) the expression of cell adhesion and differentiation genes from NOR- (and to less extent undifferentiated-) cells signals a better survival probability. Oppositely, a higher expression in older patients of RNA splicing genes from NOR cells, and to a larger extent of replication-stress genes in undifferentiated cells is likely to signal a worst outcome.

In this regard, low-risk neuroblastoma is common in younger children with <18 months of age at the time of diagnosis1,5 and NTRK1 expression (encoding TRKA) is a strong prognostic histological marker for favorable cases, whereas TRKB (encoded by NTRK2) is associated with poor outcome2. One possibility is that different progenitor populations source neuroblastomas at different ages, thereby accounting for the remarkable clinical heterogeneity that is coupled with age. Thus, age-at-diagnosis is a strong outcome predictor. Specifically, favorable neuroblastomas might originate from embryonic developmental errors when SCPs differentiate to chromaffin and neuroblast, whereas unfavorable neuroblastomas could arise from errors during postnatal development when TRKB+ cholinergic progenitors repopulate chromaffin cells post-natally.