Single-cell transcriptomics of human embryos identifies multiple sympathoblast lineages with potential implications for neuroblastoma origin.

Characterization of the progression of cellular states during human embryogenesis can provide insights into the origin of pediatric diseases. We examined the transcriptional states of neural crest- and mesoderm-derived lineages differentiating into adrenal glands, kidneys, endothelium and hematopoietic tissue between post-conception weeks 6 and 14 of human development. Our results reveal transitions connecting the intermediate mesoderm and progenitors of organ primordia, the hematopoietic system and endothelial subtypes. Unexpectedly, by using a combination of single-cell transcriptomics and lineage tracing, we found that intra-adrenal sympathoblasts at that stage are directly derived from nerve-associated Schwann cell precursors, similarly to local chromaffin cells, whereas the majority of extra-adrenal sympathoblasts arise from the migratory neural crest. In humans, this process persists during several weeks of development within the large intra-adrenal ganglia-like structures, which may also serve as reservoirs of originating cells in neuroblastoma.

Introduction

A number of congenital and cancer-related pathologies originate in human embryos within the sympathoadrenal and aorta-gonad-mesonephros (AGM) regions, where adrenal glands, sympathetic ganglia, hematopoietic cells, kidneys and other tissues develop from the neural crest and mesodermal lineages[1-4]. Although the mouse model provides detailed and experimentally validated knowledge of developmental mechanisms, the differences between mouse and human development might preclude understanding of some disorders. For instance, elucidating the human-specific aspects of development might be critical for future treatments of pediatric cancers including neuroblastoma.

Catecholaminergic cells in sympathetic chain and suprarenal ganglia originate from the neural crest cells migrating ventrally towards dorsal aorta[5,6]. Controlled by a different mechanism, in parallel with this, preganglionic visceral motor nerves bring nerve-associated neural crest-like Schwann cell precursors (SCPs) into the area during the time of adrenocortical development to generate chromaffin cells.[7,8]. Thus, two different, but complementary developmental mechanisms and transitions contribute to the diversity of cell types in the mature sympathoadrenal system. The diversity of progenitors and transitions might be recapitulated in some solid tumors associated with the adrenal medulla or sympathetic ganglia, namely paraganglioma, pheochromocytoma, and neuroblastoma[9-11]. Neuroblastoma is a pediatric cancer, resulting in approximately 15% of total cancer-related deaths of children between the ages of 0-4 years old[12,13]. These tumors commonly emerge near the adrenal gland, and are characterized by high inter- and intra-patient heterogeneity[13,14].

Next to the developing adrenal glands and dorsal aorta, the AGM region generates primordia of adrenal cortex, kidneys, gonads, and liver initiated by the intermediate mesoderm[15-17]. The dorsal aorta hosts hematopoietic progenitors that appear prior to the hematopoietic events in the liver[18-20] (Supplementary Note). Being intermixed within sympathoadrenal and AGM areas, diverse neural crest- and mesoderm-derived cell types influence each other by inducing new cell states and fates.

Knowledge of the transitory states within the context of the entire reconstructed tree of cellular differentiation in a human embryo is currently missing but is critical for fundamental understanding of pediatric tumors, including neuroblastoma, which originates during sympatho-adrenal differentiation.

To respond to these challenges, we examined the cell type heterogeneity and dynamics of the developing human sympathoadrenal region, identifying and characterizing a number of previously unknown cell fates and transitions. We find that during the development of sympathoadrenal system, intra-adrenal sympathetic neuroblasts arise from SCPs. In turn, these immature sympathoblasts can transit into local neuroendocrine chromaffin cells. Such newly established lineage connections might have important implications for understanding of neuroblastoma origin and cell heterogeneity. Finally, our data revealed previously uncharted transitions from mesodermal progenitors into hematopoietic system, endothelial cells, adrenal cortex and early kidney primordium.

Results

Transitions in human and mouse sympathoadrenal fates

We performed single cell transcriptomics analysis (10x Chromium) of isolated individual cells from dissected adrenal glands with surrounding tissue from human embryos at post-conception weeks 6, 8, 9, 11, 12, 14 (Fig. 1a, Extended Data Fig. 1, for technical information see Supplementary Note). Joint analysis of the 74401 cells filtered to meet computational quality control criteria revealed interconnected stage-dependent structure of cell clusters (Fig. 1b-f, Extended Data Fig. 1), reflecting the developmental progression of neural crest-derived sympathoadrenal system and mesoderm-derived mesenchyme, endothelial cell subtypes, hematopoiesis, liver, adrenal cortex and kidneys (Supplementary Figure 1-3) (please see Supplementary Note). The grouping of proliferating cells did not influence the overall topology of the joint embedding (Fig. 1c).

Fig.1

Heterogeneity of human sympathoadrenal region across several developmental stages.

a, The workflow of the single cell RNA sequencing data acquisition. b, A UMAP embedding of the 74401 cells, color-coded by tissue type with the major genes defining each cluster. c, Proliferating cells, defined by MKI67 expression. d, Genes defining each tissue type. e, Cell doublets probability (Scrublet score), and key markers of the major tissue types, shown on the combined UMAP embedding. f, Cells isolated from different developmental stages shown on the combined UMAP embedding. Colors are the same as in (b).

Extended Data Fig. 1

Technical quality controls

a, UMAP embedding of 99553 cells from 11 individual human fetuses (including the 74401 cells shown on Fig.1a and additional samples collected with varying procedure see Methods). b, Numbers of detected genes per cell (upper panel) and numbers of UMIs per cell (lower panel) in each sample. Probability of cell doublets (Scrublet score) per sample. Color legend is the same for all plots. c, Dot plot with key genes for the clusters in (a). d, Probability of cell doublets (Scrublet score) plotted on the general UMAP embedding. e, Numbers of UMIs per cell plotted on a general UMAP embedding. f, Numbers of detected genes per cell plotted on a general UMAP embedding. g, Key genes for cell clusters plotted on a general UMAP embedding h, UMAP embedding of cells from each individual human embryo sample. Note: embryo 2 from week 9 and embryo 2 from week 2 have two biopsies combined (see methods).

We identified the sympathoadrenal differentiation domain based on the expression of FOXD3, SOX10, ELAVL3/4, ISL1, TH, PNMT, and STMN2/4 (Fig. 1b) It included nerve-associated SCPs (SOX10, PLP1, FOXD3), sympathetic neuroblasts (PRPH, STMN2, ISL1, TH) and chromaffin cells (CHGA, PNMT, PENK, TH) of different maturity (Fig. 1d).

To elucidate the sympatho-adrenal cell fate transitions in human at higher resolution, we re-analyzed the 3901 cells of SCPs, chromaffin and sympathetic fates, omitting annotated cell-cycle genes (Fig. 2a-e). The resulting embedding showed that SCPs in the human embryos connect to the STMN2 sympathoblasts and the CHGA chromaffin cells through a “fork-like” transition (arrows 1, Fig. 2a; Extended Data Table 2). We additionally visualized this transition by projecting the cells onto an expression subspace defined by the differential marker genes of the SCPs, chromaffin and sympathetic populations, as identified in Furlan et al.[7] (see Methods, Fig. 2f, Extended Data Fig. 2a,b). The RNA velocity analysis on both embeddings suggested that SCPs transit into sympathoblasts and chromaffin cells (Fig. 2b, f). This “fork-like” transition was observed in all of the investigated stages (Fig. 2d,f; Fig. 3; Extended Data Fig. 2b, Supplementary Note). Such direct and prolonged developmental transition between the SCPs and sympathoblasts in humans would be relevant to the development of neuroblastoma[21].

Fig.2

Heterogeneity and transitions in human sympathoadrenal system.

a, A UMAP embedding resulting from the reanalysis of the SCPs and sympathoadrenal cells subset (3869 cells). Transition from SCPs to sympathoadrenal lineage is indicated by arrows 1, and suggested transitions between sympathoblasts and chromaffin cells are indicated by arrow 2. b, Cells from week 9 embryo 1 belonging to SCPs, chromaffin and sympathetic fates are shown on the UMAP embedding from (a), overlaid with RNA velocity estimates. c, Expression of key marker genes. d, Distribution of the cells from different samples shown on the UMAP embedding. e, Phases of the cell cycle, shown on the UMAP embedding. f, Ternary plot of the SCPs and sympathoadrenal cells subset based on the signatures of cells (from Furlan et al [7]), showing RNA velocity analysis (for the week 9 embryo 2nd sample and week 9 embryo 3rd sample), insets highlight the area of sympathoblast-to-chromaffin transition (for samples W12E2 (int), W12E2 (ext), W9E1 (para) see Methods section). g, Dot plot of cluster-specific marker genes. h, Heatmaps of the gene expression dynamics in the transition from SCPs to sympathoblasts, and from SCPs to chromaffin cells. (from top to bottom) First panels: expression magnitude of SCPs gene expression program in selected cells. Second panels: expression magnitude of sympathoblast or chromaffin cell gene expression programs in selected cells. Third panels: co-expression of the two programs. Forth panels: expression magnitude of the selected genes is visualized as a heatmap (blue – high expression, green – low expression). The middle section in red box highlights co-expression of genes in transition. W – week of gestation, E – sequential number of an embryo.

Extended Data Fig. 2

Quality controls and examples of cells in transition between SCPs and sympathoadrenal fates.

a, Probability of cell doublets (Scrublet score) per sample and sympathoadrenal clustering. b, Fate transitions shown per time point. c, Immunohistochemistry against SOX10 (marker of SCPs), HAND2 (marker of developing sympathoblasts and chromaffin cells), ISL1/2 (marker of sympathoblasts) on the cross-section of adrenal gland at week 6 of human development. Scale bar is 10 μm.

Fig.3

Human-specific transition from sympathoblasts into chromaffin cells.

a, A UMAP embedding of the SCPs, sympathoadrenal and chromaffin populations, same as in Fig. 2a. b, Similar to Fig. 2h, the heatmap shows the gene expression dynamics of the “Transition 2”, between sympathoblasts and chromaffin cells. c, Expression of some genes distinguishing cells in “fork-like” transition and “Transition 2”. d, Immunohistochemistry for SOX10 (marker of SCPs) and HAND2, ISL1 (markers of early sympathoblasts) on cross-sections though week 6 fetal human adrenal gland. Arrowhead points at SOX10+/HAND2+/ISL1+ cell, indicating the presence of transitory cell phenotype in “Transition 1: SCP to sympathetic”. Scale bar is 10 μm for the main panel, 5 μm for the inset. Corrected total cell fluorescence (CTCF) of individual cells for SOX10 and ISL1/2 indicates cells in transition. 60 transiting cells were analyzed in total. e, SOX10 and PENK genes, shown on the UMAP embedding of sympathoadrenal cells. Immunohistochemistry for SOX10 (marker of SCPs) and RNA in situ hybridization for PENK (marker of chromaffin cells) on cross-sections though week 8 fetal human adrenal gland. Arrowhead points at PENK +SOX10+ cells, and PENK + signals in the insets indicating the presence of transitory cell phenotype in “Transition 1: SCP to chromaffin cells”. Scale bar is 10 μm for the main panel, 5 μm for the insets. f, HAND2 and ISL1 genes, shown on the UMAP embedding of sympathoadrenal cells. CTCF of immunostaining for HAND2 and ISL1 at weeks 6 and 11 of human fetal development. Note that at week 11, some HAND2+ cell decrease or lose ISL1/2 immunoreactivity, which indicates their transition towards chromaffin cells fate. (cell number analyzed week 6 medulla region – 60, week 6 region with ganglia-like structures – 60, week 11 medulla region – 45, week 11 region with ganglia-like structures – 42). Immunostaining for HAND2 and ISL1 on cross-sections though week 6 and week 11 of fetal human adrenal gland. Arrowheads point at HAND2+/ISL1+ sympathoblasts in “Transition 2: sympathetic to chromaffin cells”. Scale bar is 10 μm for the main panel and the inset. g, Immunohistochemistry for SOX10 (marker of SCPs) and RNA in situ hybridization for PENK (marker of chromaffin cells) and ISL1 (sympathoblasts) on cross-sections though week 8 and week 12 fetal human adrenal gland, PENK + ISL1 +SOX10- cells are outlined in a sample at week 8 and PENK + ISL1 -SOX10- cells are outlined in a sample at week 12. Scale bar is 10 μm for the main panel and the inset.

The cells within the “fork-like” transition directly co-express key markers of the SCP program with the markers of chromaffin and sympathoblast programs (Fig. 2g,h). Cells in transition showed overlap of the two expression programs and transient expression patterns that could not be explained by cell doublets (Supplementary Note). The experimental validations by immunohistochemistry confirmed the existence of the bioinformatically-predicted SOX10+/ISL1+/HAND2+, SOX10+/ISL1+ and SOX10+/HAND2+ transitory cells, which occupied the intermediate position between the mutually exclusive populations of ISL1+/HAND2+ sympathoblasts and SOX10+ SCPs in human embryonic tissue samples (Fig. 3d, Extended Data Fig. 2c). Cells in transition from SCPs to chromaffin fate were validated by co-localization of PENK RNA and SOX10 protein (Fig. 3e).

We then examined the histological positions of SCPs, sympathoblasts, and chromaffin cells, in the sections of human embryonic adrenal glands. Chromaffin cells in diffuse chain-like conglomerates were mingling with cortical cells and large sympathetic ganglia-like structures in the fetal human adrenal gland (Extended Data Fig. 3a-b). These ganglia-like aggregations of neuronal cells represented distinct large and proliferative structures within the cortical tissue between weeks 6-14 (Extended Data Fig. 3a-c). This appears to be evolutionary conserved, as we detected scattered HuC/D sympathetic neurons in mouse embryonic adrenal glands (Extended Data Fig. 3d) and as reported for different animals[22-29]. Furthermore, within and adjacent to these intra-adrenal ganglia-like structures, we detected SOX10+/ISL1+ SCPs-to-sympathoblast transiting cells (Extended Data Fig. 2c). Within the postnatal human adrenal gland, the smaller aggregations of sympathetic neurons and the solitary sympathetic cells appeared evident by H&E staining, immunohistochemistry for CHGA, VIP, and RNAscope for STMN2, VIP (Extended Data Fig. 3c).

Extended Data Fig. 3

Comparative histological analysis of the developing adrenal glands from human and mouse.

a, The distribution of chromaffin cells as shown by TH immunostaining in the human adrenal gland. Scale bar is 200 μm. Note the diffuse nature of chromaffin cell localization as opposed to the compact medulla in mouse adrenal glands (d). b, c, The ISL1 high/HUC/D+/TH low intra-adrenal ganglia-like structures identified at all investigated stages of the human development are proliferative (according to MKI67 immunoreactivity) and can be discriminated from ISL1 low/HUC/D-/TH high chromaffin cells. Scale bar is 200 μm for the main panel and 20 μm for the insets. Note that chromaffin cells are associated with 2H3+ (Neurofilament+) axons. CART+ extrinsic innervation is found next to HUC/D+ intra-adrenal ganglia-like structures and is not associated with chromaffin cell and sympathetic cells somas at this stage. Later on CART+ expression segregates to chromaffin cell population (for staining see Extended Data Figure 5). Hematoxilin/Eosin staining (H&E) as well as immunohistochemistry against TH, CHGA and VIP of postnatal human adrenal glands, shows intra-adrenal sympathetic neurons in adult human adrenal glands. RNA scope in situ hybridization against CHGA, STMN2, VIP on postnatal human adrenal glands, showing clusters of intra-adrenal sympathetic neurons. Scale bar is 100 μm for the main panel and 20 μm for the insets. d, Analysis of the developing mouse adrenal gland showing HUC/D+ or ISL1high cells within the medulla at E18.5 developmental stage. Note CART+/HUC/D+ and ISL1high cells in the sympathetic ganglion outside of the gland. Scale bar is 100 μm for the main panel and 20 μm for the insets. AM: adrenal medulla, ChCs: chromaffin cells.

We additionally sequenced cells from the sympatho-adrenal region at week 12, where the outer (capsular) layer of adrenal gland was carefully removed for a separate analysis (Extended Data Fig. 4). This enabled us to compare the ganglia-like structures inside the adrenal medulla with the ganglia located on the outside of the gland (for example, suprarenal ganglion (SRG)) and in the sympathetic chain (Extended Data Fig. 4a-f). The comparison revealed no major differences in RNA profiles of intra- and extra-adrenal populations of sympathetic cells (Extended Data Fig. 4a-f). Via immunohistochemistry and FISH, we validated the distribution of predicted sympathetic subtypes in intra- and extra-adrenal populations (Extended Data Fig. 4f). The mature sympathoblasts localized to the central part of the ganglia-like structures as indicated by the activation of PRPH and the loss of EYA4 expression (Extended Data Fig. 4). We detected the expression of CARTPT in the cluster of chromaffin cells in human, and confirmed it by mapping CART-immunopositive cells on the sections of human adrenal glands, where CART co-localized with TH, but not with HUC/D within sympathetic intra- and extra-adrenal populations (Extended Data Fig. 6g-h, Extended Data Fig. 5a-b, Supplementary Note). Contrary to this, Cartpt in mice was restricted to sympathoblasts during development (Extended Data Fig. 5c-d).

Extended Data Fig. 4

Comparison of human intra-adrenal, extra-adrenal sympathoblasts and chromaffin cells.

a, UMAP embedding resulting from the re-clustering of sympathoblasts from 4 embryonic stages (see methods). b, Dot plot of key genes expressed in proliferating and non-proliferating sympathoblasts. c, UMAP embedding with highlighted positions of intra-adrenal and extra-adrenal sympathoblasts. Please note the uniform mixing of both populations. d, UMAP embedding with defined samples contribution. e, Key genes of sympathetic lineage. f, Immunohistochemistry against TH, NPY, VMAT2, NF200, PRPH, HUC/D and RNA in situ hybridization against EYA4, STMN2 on cross-sections of week 8, 9, 11 and 14 human adrenal gland showing that intra- and extra-adrenal sympathoblasts are comprised of the same populations across different developmental stages. Note the EYA4 negative area inside the ganglion at week 14 showing the mature sympathoblasts. g, UMAP embedding showing the positions of SCPs, proliferating sympathoblasts, sympathoblasts, chromaffin cells with corresponding markers. h, Immunohistochemistry against CART, TH, HUC/D and RNA in situ hybridization against PENK on cross-sections of week 14 human adrenal gland showing differences in the marker expression discriminating sympathoblasts and chromaffin cells. Scale bar is 20 μm, scale bar on insets is 10 μm.

Extended Data Fig. 6

Transitions between SCPs and sympathoadrenal cell types

a-c, Trajectories corresponding to Transition 1 and Transition 2 with heat maps show genes with non-linear patterns of mixing.

Extended Data Fig. 5

Expression of CARTPT in the sympathoadrenal domain in human and mouse.

a, UMAP of human sympathoadrenal cells with expression of PRPH (marker of sympathoblasts), PENK (marker of chromaffin cells) and CARTPT showing the specific expression in the chromaffin cells cluster. b, Immunohistochemistry against CART, TH and HUC/D on cross-sections of week 11 human adrenal gland showing external ganglion, internal ganglia-like structure and chromaffin cells. c, UMAP of mouse sympathoadrenal cells with expression of Prph (marker of sympathoblasts), Chga (marker of chromaffin cells), and Cartpt. Note that Cartpt shows the specific expression in sympathoblasts cluster. Scale bar is 20 μm, scale bar on insets is 10 μm. d, Immunohistochemistry against NF200, TH, HUCD, and CART on the cross-sections of E12.5, E13.5, E16.5 mouse adrenal glands showing sympathetic chain ganglia (SG), suprarenal ganglia (SRG) and adrenal medullae (AM). Scale bar is 20 μm for E12.5 and E13.5. Scale bar is 100 μm for E16.5.

The joint embedding of human SCPs, sympathoblasts and chromaffin cells revealed an additional transition (“Transition 2”) between the sympathoblasts and the chromaffin cells (arrow 2, Fig. 2a,b, Fig. 3a). This transition showed direct overlap of the two expression programs (Fig. 3b) and transient expression patterns that could not be explained by cell doublets (Extended Data Fig. 6c). Direct comparison with the “fork-like” transition showed that cells from the “Transition 2” were distinguished by specific markers (SOX4, BEX1, RAMP1, PENK) absent from the “fork-like” transition (Fig. 3c). To ensure that the separation of the two transitions is not due to a potential residual cell cycle signature, we further omitted all genes that showed non-negligible correlation with the key cell cycle genes (see Methods). The results of this re-analysis were consistent with the topology shown in Figs. 2,3 with both transitions remaining stable (Extended Data Fig. 7).

Extended Data Fig. 7

Stability of transitions between SCPs and sympathoadrenal cell types after cell cycle genes removal

UMAP embedding of SCPs and sympathoadrenal populations without cell cycle-associated genes. Note that the exact same cells occupy the same mapped transitions in clustering before and after extensive removal of cell cycle-associated genes (see methods).

We validated the existence of the “Transition 2” state by a quantitative analysis of ISL1 and HAND2 expression and observed a gradually changing mix of predicted ISL1high/HAND2high and ISL1-/HAND2high cells (Fig. 3f). Moreover, the co-localization of PENK RNA (chromaffin cells) and ISL1 RNA (sympathoblasts) in the same cells negative for SOX10 in developing human adrenal glands further supported the independence of “Transition 2” from the SCP-dependent “fork-like” transition (Fig. 3g). The RNA velocity analysis suggested the predominant direction in the “Transition 2” from sympathoblasts to chromaffin cells (Fig. 2b,f), however a reverse direction of transition is not excluded.

In mice, SCPs generate chromaffin cells of adrenal medulla via the intermediate “bridge” cell type[7]. Since the newly observed transitions between cell fates in human vs mouse appear to happen in a different order, we set out to resolve the reasons for this inconsistency between species. In study by Furlan et al [7], lineage-tracing experiments demonstrated that the vast majority of mouse sympathoblasts in suprarenal ganglia and sympathetic chain originate from migratory neural crest cells without transition via the SCP-state. Thus, we hypothesized that there may exist an uncharacterized population of mouse SCP-derived sympathoblasts in the adrenal medulla region. To identify such a population, we performed scRNA-seq of lentivirus-based lineage-traced neural crest derivatives of the mouse adrenal gland region at E13.5. The analysis of the differentiation trajectory and RNA velocity showed the transition from SCPs to chromaffin cells, consistent with earlier studies (transition 1, Fig. 4a-d; Extended Data Fig. 8a-b)[7]. The presence of such transitory SOX10+/HAND2+ cells within mouse adrenal medulla was successfully validated (Fig. 4d,e). The transient cell population connected SCPs to immature chromaffin cells expressing Chga, Chgb, Penk (Extended Data Fig. 8d lower panel) with the bridge cell signature including Htr3a, Thbd and other genes being consistent independently of the sequencing method (Extended Data Fig. 8e).

Fig.4

Immature chromaffin cells give rise to sympathoblasts in a mouse embryo.

a, UMAP embedding and clustering of mouse 2761 SCPs, sympathetic and chromaffin cells, sequenced by 10X Chromium, from the E13.5 mouse embryo trunk, with overlaid RNA velocity analysis. Arrows point at the direction of fate Transitions 1 and 2. b, Heatmaps show gradual gene expression changes during transitions from SCPs into chromaffin (“Transition 1”), and from chromaffin into sympathoblasts (“Transition 2”). The cells are ordered with the direction indicated by the color bar on the top (colors are the same as in (a)). c, tSNE embedding of 356 cells, sequenced by Smart-seq2, from the adrenal gland with associated sympathetic ganglia and color-coded for cell types. Expression of selected sympathetic and chromaffin markers highlight a small group of co-expressing cells (lower panels). d, Genes defining SCPs (Sox10), sympathoblasts (Prph, Isl1), chromaffin cells (Chga), and proliferating cells (Mki67) are shown on the UMAP embedding from (a). Note that Sox10 and Hand2 are co-expressed in transiting cell population from SCPs to chromaffin cells. e, Transition 1: Immunostaining for SOX10 (marker of SCPs), HAND2 and TH (markers of both sympathoblasts and chromaffin cells) on cross-sections though E12.5 mouse adrenal medulla primordium. Scale bar is 100 μm for the main panel and 10 μm for the insets. Arrowheads point at SOX10+/HAND2+ cells in transition from SCPs to chromaffin cells. Transition 2: Immunostaining for CHGB (marker of chromaffin cells) and RNA in situ hybridization for Elavl2 (marker of sympathoblasts) on cross-sections though E13.5 mouse adrenal medulla primordium. Scale bar is 20 μm for the main panel and 10 μm for the insets. CHGB+/Elavl2 + cells in transition from chromaffin cells to sympathoblasts are outlined. f, Scheme illustrating the experimental design of individual cell barcoding by lentivirus library. g, Individual clones with the cells of different fates overlaid on sympatho-adrenal UMAP embedding. h, Distribution of clones with cells of different fates. i, Illustration of potential differentiation models (Model 0-2), with the corresponding likelihood tests. Model 0, P⊥ – probability of independence of observing a clone in chromaffin cells and in sympathoblasts evaluated using Fisher’s exact test (two-tailed, H0: probability of observing a clone in chromaffin cells is independent of whether it was observed in sympathoblasts). Model 1 and 2, P⊥ – probability of conditional independence evaluated by one-tailed χ2 probability (Model 1 H0: clone abundance in sympathoblast cells is independent of its SCP abundance given its chromaffin abundance; Model 2 H0: clone abundance in chromaffin cells is independent of its SCP abundance given its sympathoblast abundance). Pχ2 – likelihood ratio test comparing Model 2 and Model 1 (one-tailed chi-squared probability; H0: Model 1 is not more likely than Model 2). n=322 independent clones were analyzed. ✗ – model is rejected, ✓ - model is accepted.

Extended Data Fig. 8

Integrated analysis of sympatho-adrenal lineage sequenced by different methods across mouse and human species.

a, Probability of cell doublets (Scrublet score) plotted on the UMAP embedding of mouse SCPs and sympathoadrenal cells. b, UMAP embedding of mouse SCPs and sympathoadrenal cells without cell cycle genes. c, tSNE embedding of 356 cells sequenced by Smart-seq2 from the mouse adrenal gland with associated sympathetic ganglia at E13.5. The color-code reflects the identified cell types and shows gene expression aspects of sympathoblasts and chromaffin cell markers (right). d, Joint analysis of sympathoadrenal lineage from mouse embryos at E13.5 as sequenced by 10x Chromium and Smart-seq2 platforms (upper). Early expression of chromaffin cells program in the transient cells (bottom). e, Expression of the Bridge cell marker genes according to the annotation in Smart-seq2 dataset from Furlan et al. Science 2017 in the joint embedding of Smart-seq2 (upper) and 10x Chromium (middle) embedding as integrated by Conos. The lower panels show the Bridge cell markers in the original embedding of mouse scRNA-seq generated by 10x. f, Integrated analysis of sympathoadrenal lineage from human and mouse embryos both sequenced by 10x platform. g, number of clones that are present in different combinations of cell types: SCPs, chromaffin cells and sympathoblasts. h, Illustration of potential differentiation models (Model 1 and 3), with the corresponding likelihood tests. Pχ2 – p-value of the likelihood ratio test comparing Model 3 and Model 1 (one-tailed chi-squared probability; H0: Model 3 is not more likely than Model 1). n=322 independent clones were analyzed. Model 3 that includes the possibility of a direct transition from SCP to sympathoblasts does not show any improvement as compared to Model 1.

However, in this mouse adrenal dataset, we observed a robust transition between the immature chromaffin cells and the population of sympathetic cells (Fig. 4a). While such bridge to sympathoblasts was not apparent in the Smart-seq2 data[7], scoring of cells for expression of sympathetic and chromaffin markers identified several cells that express elements of both signatures (Fig. 4c, Extended Data Fig. 8c,d, Supplementary Data File 1). The bridges reflecting transitions between SCP, chromaffin and sympathoblast populations were consistent in joint embeddings of human and mouse data (Extended Data Fig. 8f).

RNA velocity estimates suggested that the direction of this transition in mouse shows the opposite direction (from chromaffin to sympathoblasts), as compared to the one observed in human adrenal glands (from sympathoblasts to chromaffin) (Fig. 4a). To follow up on this hypothesis, we traced the emerging clonal structures in sympatho-adrenal domain of a mouse using lentiviruses to deliver expressed barcodes to the neural tube before the neural crest immigration (Fig. 4f). Recovering the clonal barcodes from scRNA-seq data (see Methods), we found that the majority of clones covered two or more fates (Fig. 4g,h, Extended Data Fig. 8g, Supplementary Data File 2). Overall, we observed a higher fraction of clones that were shared between the chromaffin and the SCP populations (80%) than between the sympathoblasts and SCPs (64%), consistent with the direction suggested by the RNA velocity. Detailed analysis of the clonal co-occurrence frequencies enables us to quantify the likelihood of various differentiation topologies. A parsimonious model supposes that the SCPs contribute to chromaffin and sympathetic populations directly, through two independent transitions (Model 0, Fig. 4i). However, the frequency of the clones appearing in both chromaffin and sympathetic cells is too high, and Model 0 is clearly rejected with Fisher’s exact test for independence (P ⊥<3.1×10−7). We then examined a model in which, as suggested by the RNA velocity estimates, sympathetic cells contribute to chromaffin population (Model 1, Fig. 4i). This model dictates a conditional independence property. Specifically, Model 1 postulates that the abundance of a given clone in the sympathetic population should be independent of their abundance in the SCPs, given their abundance in the chromaffin population. Indeed, the observed clonal frequencies are consistent with such conditional independence (KCI test P ⊥=0.24, χ 2 test P ⊥=0.60). We note that such conditional independence test does not rule out presence of an additional SCP-to-sympathoblast transition, but instead indicates that such a direct transition is not needed to explain the observed data (Extended Data Fig. 8h). We then evaluated a reverse model, in which sympathetic cells contribute to the chromaffin population (Model 2, Fig. 4i). This model clearly fails the conditional independence tests (KCI test P ⊥=1.7×10−5, χ 2 test P ⊥=2.0×10−19). The clonal data, therefore, argues strongly in favor of Model 1, where chromaffin cells originate from SCPs, and then contribute to the sympathetic population (log likelihood ratio test P<6.9×10−11).

To confirm that the sympathetic neurons in mouse adrenal gland develop from SCPs via an intermediate chromaffin cell state, we performed a lineage tracing experiment by using SCP-specific Plp1 embryos, where genetic recombination was activated at E11.5 – a stage, when the neural crest migration in mouse embryo is over, and all of the SOX10+/FOXD3+ cells are associated with the innervation in the head and trunk regions[7,30,31]. We find that intra-adrenal sympathetic neurons were robustly traced (~30%), with the efficiency of genetic recombination in the nerves under ~80% (Fig. 5a-e). In contrast, only a tiny fraction of sympathetic neurons in the suprarenal and sympathetic ganglia were traced (~2-4%), consistent with our previously published data[7] (Fig. 5b-e). Next we utilized another mouse strain recombining in early sympatho-adrenal lineage – Ret. Activation of genetic recombination in these embryos at E10.5 confirmed the difference in the NCC vs. SCP-derived origin of the sympathoblasts in the extra- and intra-adrenal locations (Fig. 5f-h). Almost 50% of sympathoblasts from the suprarenal and sympathetic chain ganglia were labelled, whereas only ~10% of CART+ sympathoblasts were traced within the adrenal medulla (Fig. 5g-h). Thus, the lineages of SCP-derived and NCC-derived sympathoblasts largely split before E10.5, when the neural crest cells migrate and commit to becoming a sympathetic cell or an SCP along the nerve. As we did not address the origin of intra-adrenal sympathoblasts in previous study[7], new data on the SCP-mediated origin of a significant fraction of intra-adrenal sympathoblasts adds to the overall picture of the developing sympatho-adrenal system in mouse, and reconciles the observations of the current study with Furlan et. al.[7] (Fig. 5i). Taking into account the topology of the transitions apparent from the scRNA-seq data (Fig. 4a), our results show that a fraction of sympathoblasts in mouse, positioned within the adrenal medulla, originates from the immature chromaffin cells, which in turn originate from SCPs via a previously characterized Htr3a bridge. This derivation of the intra-adrenal sympathoblasts from chromaffin population is consistent with the known developmental timing: while the chromaffin cells start to appear starting from E11.5, the CART+ sympathoblasts in adrenal medulla in mouse start to emerge at a notably later time, around E13.5[32] (Extended Data Fig. 7 c-d; Supplementary Note).

Fig.5

Intra-adrenal sympathoblasts originate from Plp1 + SCPs being specified under control of Ret.

a, Immunohistochemistry for yellow fluorescent protein (YFP) (showing Plp1-traced cells), TH (marker of both chromaffin cells and sympathoblasts), and CART (marker of sympathoblasts in mouse) on cross-sections of the developing adrenal medulla (AM), suprarenal ganglion (SRG) and sympathetic ganglion (SG). Genetic tracing in Plp1 animals injected with tamoxifen at E11.5 and analyzed at E15.5. Yellow arrowheads show traced sympathetic neurons Plp1 YFP+/CART+/TH+. Scale bar is 20 μm for the main panels and 10 μm for the insets. b, Efficiency of recombination in the experiment was 77.7±4.5% (mean ± s.e.m., n=4). c, Close view of the outlined area on the cross-section through AM. Filled yellow arrowheads show Plp1 YFP+ sympathoblasts, while empty arrowheads show non-traced sympathoblasts. Scale bar is 10 μm for the insets. d-e, The number of Plp1-traced sympathoblasts per section in SG, SRG and AM. d – actual number and e – percentage of the CART+ cells. Numbers of Plp1 sympathoblasts per section and percentage over the total sympathoblasts for SG: 1.15±0.24 cells (2.72±0.56%), SRG: 2.52±0.56 cells (2.36±0.52%) and AM: 8.03±2.61 cells (30.52±6.48%). P-values of paired t-test: SG vs SRG=0.5165, SG vs AM=0.0197, SRG vs AM=0.0181. Data are presented as mean ± s.e.m., n=3. f, Immunohistochemistry for TOMATO (recapitulating Ret traced cells), TH, and CART on cross-sections of the developing AM, suprarenal ganglion (SRG) and sympathetic ganglion (SG). Genetic tracing in Ret animals induced at E10.5 and analyzed at E15.5. Filled yellow arrowheads show Ret sympathoblasts, while empty arrowheads show non-traced sympathoblasts. Scale bar is 20 μm for the main panels and 10 μm for the insets. g-h, The number of Ret traced sympathoblasts in SG, SRG and AM. g – actual number and h – percentage of the CART+ cells. Numbers of Ret sympathoblasts per section and percentage over the total sympathoblasts for SG: 29.92±3.86 cells (49.39±3.33%), SRG: 46.5±8.14 cells (42.98±1.36%) and AM: 2.08±0.46 cells (8.06±2.11%). P-values of paired t-test : SG vs SRG=0.2891, SG vs AM=0.0144, SRG vs AM=0.0047. Data are presented as mean ± s.e.m., n=3. i. Scheme illustrating the difference in fate transitions in human and mouse. The predominant direction of intra-adrenal sympathoblast differentiation is indicated by the blue arrows. Grey arrows indicate the conventional differentiation of the neural crest cells into the majority of extra-adrenal sympathoblasts and additional ways of sympatho-adrenal differentiation in human embryos. Note: n in all experiments represents the number of biologically independent samples (individual mouse embryos).

This is consistent with our previous findings, where the numerous lineage-traced clones appeared unique for SRG and medulla, supporting that chromaffin cells in medulla and sympathoblast in extra-adrenal ganglia are derived from different neural crest and post neural crest sub-lineages[7]. Importantly, it should be noted that the transition from SCPs to sympathoblasts occurs only to intra-adrenal but not to extra-adrenal cells that represent the majority of captured sympathoblasts in a single cell dataset.

Overall, based on a combination of experimental and single cell transcriptomics-based evidence, we conclude that the large portion of intra-medullary sympathetic neurons in a mouse embryo are derived from SCPs (given the efficiency of the lineage tracing) via immature chromaffin cell state. At the same time, in a human embryo, SCPs can transit both into chromaffin cells and intra-adrenal sympathoblasts. Though it is not currently possible to confirm by clonal tracing, bioinformatics analysis suggests that in humans, these sympathoblasts can, in turn, transit into chromaffin cells (Fig. 5i).

Human sympathoadrenal populations relevant to neuroblastoma

Given the potential relevance of the developmental cell populations within the human AGM and sympatho-adrenal region to neuroblastoma, we directly compared these cell populations with the cell type heterogeneity observed in neuroblastoma. For this, we carried out joint analysis of our data with scRNA-seq of neuroblastoma biopsies from two different patients explicitly described by Olsen et al.[33] (Fig. 6a and Fig. 6c, left) and six patients described by Dong et al.[34] (Fig. 6b and Fig. 6c right). Neuroblastomas appeared highly heterogeneous in their composition, and included malignant cells of both adrenergic and mesenchymal lineages, as reported earlier[33,35,36]. We found that the embryonic abdominal mesenchymal populations aligned with the mesenchymal cells found in tumor samples, the SCP population aligned with Schwann component from the tumor samples, and the majority of the adrenergic tumor cells showed clear similarity to embryonic sympathoblasts with a fraction of tumor cells co-aligning with chromaffin cells (Fig. 6).

Fig.6

Transcriptional signatures of neuroblastoma cells match with transcriptional signatures of fetal human cells.

a, A joint UMAP embedding of the 66597 human fetal cells and two samples of neuroblastoma (5190 and 5288 cells from the samples NB12 and NB26, respectively). The first embedding shows the combined cell populations, with the other three panels showing subsets of cells from fetal data, NB12 and NB26 samples, respectively, using the same embedding. The location of the SCPs, sympathoblasts, and chromaffin clusters is outlined. Lower panels: marker genes of the major cell populations in sympathoadrenal domain. b, A joint UMAP embedding of 66597 human fetal cells and 6 samples of neuroblastoma from Dong et al.[34]. The first embedding shows the combined cell populations, with the other two panels showing subsets of cells from the fetal data and combined neuroblastoma data respectively, using the same embedding. c, The ratio of cell fates annotated based on the fetal data set within the neuroblastoma datasets from Olsen et al.[33] (left) and Dong et al.[34] (right). d, Survival analysis of cluster-specific gene expression signatures in non-MYCN-amplified (n=353 samples) or MYCN-amplified (n=88 samples) neuroblastoma subtypes. P-values are shown for the two-sided log-rank test between Kaplan–Meier curves for the top and the bottom quarter of tumors sorted according to a gene expression signature. For non-MYCN-amplified tumors P-values are: SCPs P=0.16, Proliferating sympathoblasts P<0.0001, Chromaffin cells P=0.015, Sympathoblasts P=0.078. For MYCN-amplified tumors P-values are: SCPs P=0.38, Proliferating sympathoblasts P=0.0011, Chromaffin cells P=0.098, Sympathoblasts P=0.67.

We then examined bulk RNA-seq data from a large cohort of patients with different neuroblastoma subtypes of adrenal localization for potential associations between the patient survival and gene expression signatures of the embryonic SCPs, chromaffin, proliferative and mature sympathoblasts. The worst survival was associated with the transcriptional signature of proliferative sympathoblasts. This association held regardless of the presence or absence of MYCN-amplification in tumors and was independent of cell cycle genes (Fig. 6d, Supplementary Data File 3). The SCPs and mature sympathoblasts signatures were associated with better prognosis in the non-MYCN-amplified tumors, whereas the embryonic chromaffin cells signature was associated with poor prognosis, suggesting possible existence of immature chromaffin-like neuroblastoma cell subtypes. MYCN amplification in tumors worsened the prognosis for the mature sympathoblast signature, although the SCP signature always stayed correlated with positive outcome of the disease (Fig. 6d).

Fetal human adrenal cortex and relation to neuroblastoma

Focused re-analysis of the 37117 cortex and kidney cells, together with the subpopulation of cells connecting them (Fig. 7, Supplementary Figure 5) suggested that the development of human adrenocortical primordium begins with a fate split from an early mesenchymal progenitor state sharing transcriptional signature with the early kidney progenitor cells (SOX18, EGFL7, KDR, CYP26B1, FLT1, CAVIN2, CLDN5) (Supplementary Figure 3, 5). Furthermore, on the kidney side, the development proceeded into the major sub-lineages, including nephron progenitors (SIX2, MEOX1) giving rise to podocytes (PODXL, MAFB) and interstitial progenitors (MEIS1, FOXD1) (Supplementary Figure 3, 5), in agreement with previously published single-cell analysis[37]. In the developing human cortex, we identified a progressive time-dependent shift from a predominant formation of the definitive zone (proliferative and non-proliferative) at weeks 6-9 to the fetal zone and transitional zone at weeks 9-14 (Fig. 7a-c; Supplementary Note).

Fig.7

Heterogeneity of the developing human adrenal cortex, kidney and mesenchyme and their ability to shape microenvironment in neuroblastoma.

a, UMAP embedding showing clusters resulting from reanalysis of 37117 cells of the adrenal cortex, kidney primordium and mesenchyme. b, Cells isolated from different developmental stages are shown on the UMAP embedding (a), with colors representing major sub lineages. c, Expression of cluster marker genes, shown as a dot plot. d, Key genes demarcating different subpopulations of the developing human adrenal cortex, and genes reflecting proliferation, maturation, and the start of hormonal production. e, Survival analysis of organ- or tissue-specific gene expression signatures with survival curves of patients in non-MYCN-amplified (n=353 samples) or MYCN-amplified (n=88 samples) neuroblastoma subtypes. P-values are shown for the two-sided log-rank test between Kaplan–Meier curves for the top and the bottom quarter of tumors sorted according to a gene expression signature. For non-MYCN-amplified tumors P-values are: Adrenal cortex P=0.00055, Kidney P=0.0016, Mesenchyme P=0.11. For MYCN-amplified tumors P-values are: Adrenal cortex P=0.58, Kidney P=0.97, Mesenchyme P=0.55.

Adrenal cortical development proceeds in a tight coordination with cells of the neural crest origin, giving rise to the adrenal medulla. At the same time, the crosstalk between cortex and medulla is not necessary for the formation of chromaffin cells as a cell type[38], but important for a specification of PNMT+ population[39,40]. To analyze the potential role of mesodermal lineage derivatives in formation of neuroblastoma we compared single cell signatures of cells with neuroblastoma patients’ survival prognosis. The presence of the embryonic adrenocortical and kidney signatures in tumors was associated with poor neuroblastoma outcome, for non-MYCN-amplified cases (Fig. 7e, Extended Data Fig. 10). We detected no difference in the outcome of the MYCN-amplified neuroblastomas. Thus, some of the embryonic cortical and kidney-derived signals may modulate the progression of non-MYCN-amplified tumors, whereas the microenvironment dependencies of tumors with MYCN amplification are likely to be distinct.

Extended Data Fig. 10

Correlation of gene expression signature of selected mesoderm-derived populations with survival probabilities of neuroblastoma patients and bootstrap validation for the survival analysis.

a, Correlation of gene expression signature of selected mesoderm-derived populations with survival probabilities of patients with different MYCN-amplified (88 samples) and non-MYCN-amplified (353 samples) neuroblastoma subtypes. P-values of the two-sided log-rank test between Kaplan–Meier curves for the top and the bottom quarter of tumors sorted according to a gene expression signature. Non-MYCN-amplified tumors: immune cells P=0.83, liver P=0.5, erythroid cells P=0.1, intermediate mesoderm P=0.077, endotelial cells P=0.24, Melanocytes P<0.0001. MYCN-amplified tumors: immune cells P=0.41, liver P=0.95, erythroid cells P=0.42, intermediate mesoderm P=0.72, endotelial cells P=0.066, Melanocytes P=0.047. b, Bootstrap validation for the survival analysis of cluster-specific gene expression signatures in non-MYCN-amplified and MYCN-amplified neuroblastoma subtypes. For a set of marker genes, we took 100 random samples of the same size with replacement (see Methods).

Overall, our results reveal previously unrecognized transitions in mesoderm-derived and neural crest-derived human embryonic cell types in the sympatho-adrenal area. This, in turn, may help us to better understand the development of this highly heterogeneous region, including clues mediating local integration. The newly characterized transitions of the sympathoadrenal populations, together with the analysis of local cell type signatures suggest additional hypothesis for the origin of pediatric neurocristopathic cancers, and their connection with the microenvironmental signals.

Discussion

Although the initiation and early development of AGM and sympathoadrenal regions has been previously addressed at histological and anatomical levels in animal models and to some extent in human embryos[23,28,41,42], the inability to resolve human cryptic cell subtypes and to lineage trace human embryos resulted in unresolved human-specific aspects of development. The existing single cell atlases from different animal models did not unambiguously resolve the development of sympathoadrenal system and provided only limited insights into fate transition around dorsal aorta[43-47].

Here we took advantage of single cell transcriptomics and resulting ability to infer the transitions between cell states to investigate the emergence of human sympathoadrenal area and AGM region, including kidneys, hematopoietic, endothelial cell types, and the mesenchymal component (Supplementary Note). Our map of human AGM and sympathoadrenal region provides an important resource for analysis of cancer cell types originating from these regions, including their potential origins, differentiation dynamics, re-enactment of embryonic states or microenvironment[48,49]. Phenotypic plasticity of tumor cells often precludes efficient treatment strategies, and the comparative insights from human developmental studies are particularly valuable for cases of neural crest-related pathologies and cancers[35,36].

Following this logic, our observations are relevant to the possible developmental origin of neuroblastoma. Neuroblastomas do not naturally appear in mice, and the current genetic models do not recapitulate the full spectrum of natural properties of the disease[50]. Current opinion holds that neuroblastoma arises during differentiation of the migratory neural crest cells towards sympathetic neurons. The observed similarity of the adrenergic component of neuroblastoma to healthy sympathoblasts suggests that some immature sympathetic cells or their immediate progenitors might be a point of tumor origin[51]. Surprisingly, almost half of neuroblastomas are found in or adjacent to the adrenal glands[52,53]. Furthermore, this specific localization of tumors points at the relevance of the developmental and cell fate decision processes in the adrenal gland region.

The medulla of adrenal gland contains mostly chromaffin cells[7]. In mice, the embryonic origin of the majority of chromaffin cells differs from the origin of the majority of sympathetic neurons, which are immediately neural crest-derived[7]. Specifically, the neural crest-derived SCPs follow the nerves towards developing adrenal gland medulla, where they separate from the nerve tracks and give rise to chromaffin cells. Only in exceptional cases, SCPs were shown to give rise to significant amounts of sympathetic neuroblasts in specific sympathetic paraganglia[54] (Supplementary Note).

The present study extends our current understanding of SCP potential to give rise to sympathetic neurons. Here we identified a population of intra-adrenal sympathetic neurons that predominantly originate from SCPs in mouse via immature chromaffin state (according to scRNA-seq and lineage tracing data), and in humans rather directly from SCPs (according to scRNA-seq and immunohistochemistry-based validations of cells in transitory state). The intra-adrenal sympathoblasts are born at later developmental time point (from E13.5 in mouse) in comparison with extra-adrenal sympathetic neurons derived predominantly from the neural crest (from E11.5 in mouse). Mouse intra-adrenal neurons are not abundant as compared to numerous chromaffin cells, and we overlooked them in the original interpretation of the lineage tracing data in Furlan et al.[. Similarly, the bridge population connecting chromaffin and sympathoblast cells, has only become apparent in the scRNA-seq data when sequencing an order of magnitude more cells (Fig. 5). Despite the difference in origin, the phenotypic distinction between extra- and intra-adrenal sympathoblasts is not apparent from their transcriptional profiles, in either human or mouse data.

In the human adrenal gland, the transition from SCPs to sympatho-adrenal states persists over several weeks (Fig. 2), unlike a short time window of the neural crest migration, and involves a rich repertoire of transitory cell states. Anatomically, sympatho-adrenal differentiation occurs within the adrenal medulla and large human-specific highly proliferative sympathetic ganglia-like structures inside the developing adrenal cortex. Moreover, the chromaffin and sympathetic subtypes transit are interconnected by additional transitory state.

Thus, the discovered complexity, plasticity, and prolonged differentiation process create unique conditions for neuroblastoma emergence specifically in the human adrenal gland area, and can likely account for inter-patient as well as intra-tumoral heterogeneity of neuroblastomas[35,36,55,56].

Additionally, human sympathoadrenal development differs from the mouse in the order of discovered transitions and at the level of expression of some important markers of sympathetic or chromaffin populations (Cartpt). Such differences might explain why mouse model system does not recapitulate human neuroblastoma disease to a full extent. This developmental perspective is consistent with our results on correlating the populations of individual heterogeneous malignant neuroblastoma cells with cell from human fetal adrenal gland, where the majority of adrenergic tumor cells aligns with embryonic sympathoblasts. In line with this, the high expression level of immature sympathoblast genes was associated with a poor survival prognosis.

Furthermore, heterogeneous neuroblastoma tumors are composed of plastic adrenergic sympathetic and mesenchymal tumor components[35,36], and the interactions between these two tumorigenic cell populations might be related to the interactions between healthy human embryonic sympathetic and mesenchymal cells. In line with this, our analysis of survival curves in different cases of non-MYCN-amplified neuroblastoma suggested a role of local microenvironments organized by mesenchymal, cortical, endothelial and kidney cell types with the least favorable outcomes corresponding to adrenal cortex-derived genes and signals. Altogether, these results may help to explain the origin and human-specific nature of different neuroblastoma subtypes, their internal heterogeneity and inter-patient variation in malignancy.

Online Methods

Ethical aspects of procedures involving human tissue

Human pre-natal tissue was retrieved from clinical routine abortions with oral and written consent from the patient. All procedures regarding the acquisition of human pre-natal tissue for research was approved by the Swedish Ethical Research Authority and the National Board of Health and Welfare (ethical reference is 2018/769-31). The material was donated for a general research purpose especially with the focus on neuronal and nervous system related cell types. Post-conception age was determined using measurement of crown-rump-length (CRL) and anatomical landmarks. After identifying the adrenal glands, this tissue was retrieved and placed in ice-cold PBS pending dissociation.

The post-natal tissue from the 1.5-2.5 years old patients at Princess Máxima Center for Pediatric Oncology, Utrecht, The Netherlands was a part of surgical resections and donated for a general research with a written consent from parents or guardians and after the IRB approval of Princess Máxima Center for Pediatric Oncology, Utrecht, The Netherlands.

The samples of neuroblastoma tumors from Swedish patients were collected after obtaining parents or guardians verbal or written consent. All samples were collected according to permits approved by the Karolinska Institutet/Karolinska University Hospital ethics committees (reference numbers 03-736, 2009/1369-31/1) in agreement with the Declaration of Helsinki. Sampling was performed at clinical routine procedures, treatment and management was according to established national and international protocols and hospital records were used to retrieve the clinical data.

Animal models and strains

The following genetically modified mouse lines were used in the study:

Ret mice were received from D. Ginty laboratory (Harvard University, USA) http://www.informatics.jax.org/allele/MGI:4437245, Plp1 were received from U. Suter laboratory (ETH Zurich, Switzerland) http://www.informatics.jax.org/allele/MGI:2663093, Rosa26R mice were ordered from The Jackson Laboratory (stock number 007914, full strain nameB6.Cg-Gt(ROSA)26Sortm14(CAG-tdTomato)Hze/J), Rosa26R mice were received from The Jackson Laboratory (stock number 006148, full strain name B6.129X1-Gt(ROSA)26Sortm1(EYFP)Cos/J). The animals from non-genetically modified mouse lines CD-1 and C57Bl/6 were obtained from Charles River Laboratories and Janvier.

Ethical aspects of procedures involving animals

Lenti-virus injections were approved by the Swedish Animal Welfare Body (Jordbruksverket and Stockholms Djurförsöksetiska nämnd, permit number Jonas Frisen: N155/16). All other animal experiments were approved by Norra Djurförsöksetiska Nämd, ethical permit Igor Adameyko: N226/15, Jordbruksverket De regionala Djurförsöksetiska nämnderna and Stockholms Djurförsöksetiska nämnd, permit number Igor adameyko: N15907-2019 and by Bundministerium fur Wissenschaft, Forschung und Wirtschaft in Ausrtia, Medical University of Vienna to Igor Adameyko (BMWFW-66.009/0018-WF/V/3b/2017).

Animals were kept in an SPF animal facility with standardized conditions (24°C, 12:12 h light–dark cycle, normal humidity, food and water ad libitum).

2-5-month-old females were placed in time-controlled matings, and the noon of day of plug was defined as embryonic day (E) 0.5. The embryos from non-genetically modified mice were collected at required stage for further analysis. The sex of the harvested embryos was not determined.

Human fetal cell isolation, storage in methanol and rehydration

Tissue was chopped into 1 mm pieces and placed in 500 ul 0.05% Trypsin/0.02% EDTA (Sigma, 59417-C) and incubated at 37°C for 15 min with gentle swirling every 5 min. Then, the tissue was triturated by pipetting until complete dissociation. 500 ul 10% FBS in PBS was added to cell suspension followed by centrifugation at 500 rcf for 5 min at 4°C and two washes with 1000 uL PBS. Cell suspension was pass though the 35 um cell strainer (Falcon, 352235) and spin down at 500 rcf for 5 min at 4°C. Cell pellet was resuspended in 100 ul 0.04% BSA in PBS and 400 ul of ice cold methanol was added for methanol fixation of the cells for storage -80°C. On the day of library preparation the cells were brought to +4°C and centrifuged at 1000 rcf for 10 min at 4°C. Supernatant was removed and 500 uL of rehydration buffer (1X DPBS (Gibco 14190144) containing 1.0% BSA (Sigma, B4287) and 0.5U/μl RNAse Out (ThermoFisher Scientific, 10777019) was added to cell pellet, re-suspended and centrifuged at 1000 rcf for 10 min at 4°C, the wash was repeated two times. If a lot of debris was observed during cell count the cell suspension was subjected for FACS to remove debris. BD FACS Aria Fusion equipped with 100 um nozzle (BD Biosciences, San Jose, CA) instrument was used for sorting. After FACS cells were pelleted again and re-suspended in rehydration buffer to get 700-1200 cells/ul but in no less than 50 uL.

For sample Week 9 embryo 1 together with adrenal gland the tissue between the kidney was processed to obtain cell from sympathetic ganglia and organ of Zuckerkandl. Cells are computationally pulled together (Extended data file 1h), shown as W9E1 (para) on Fig. 2f.

For the sample week 12 embryo 2 the outer layer of epithelium was separated from adrenal cortex and processed to additionally to internal cells obtain cell from sympathetic ganglia outside of adrenal gland (external). Cells are computationally pooled together (Extended data file 1h), shown as W12E2 (int) and W12E2 (ext) on Fig. 2f.

Neuroblastoma cell preparation for sequencing

Fresh tissue from surgery or needle biopsies was placed in Media 199 (Gibco) supplemented with 2% (v/v) FBS on ice. The tumor was cut into pieces (1mm[3]) and treated by Collagenase I-IV (Worthington Biochemical, each at 1m/ml) and Dispase (2mg/ml; Thermo Fisher) supplemented by DNAse I solution (Thermo Fisher Scientific; 1:100) and RNase Out (Thermo Fisher; 1:1000). Dissociation was done at 37°C with rotation for 45 min. Single cells were resuspended in in fetal bovine serum (FBS) with 5% DMSO and cryopreserved in liquid nitrogen. For RNA-seq library preparation cells were thawed and enriched of viable cells by staining with anti-CD235-PE (Biolegend), Calcein AM (Thermo Fisher), and DAPI and/or 7AAD staining followed by FACS sorting of for live and non-erythroid cells. BD FACS Aria III equipped with a 100um nozzle (BD Biosciences, San Jose, CA) instrument was used for sorting.

Lentivirus-based neural plate lineage tracing

A lentiviral transfer plasmid was generated by exchanging the PGK1 promoter from LV-GFP[57] with a stronger EF1a promoter, and inserting a 30 bp variable sequence into the 3’ UTR of GFP (resulting in a library of >5 million unique transcribed barcodes). Lentiviral particles were produced by KI VirusTech at high titer virus (>10^9 TU/ml).

Female CD-1 mice (Mus musculus) age 2-5 month were housed in animal facilities at Karolinska Institute, Sweden. Experimental female mice were placed in time-controlled pregnancy, day of vaginal plug presence was counted as embryonic day (E)0.5. At embryonic day (E) 8.25 (range of injection time 5AM-7:30AM) 46 nl of lentiviral particles was injected into the amniotic fluid using ultrasound-guided nano-injection, as described previously for E9.5 injection[57]. Injections performed by Andersson lab, or the core facility Infinigene at Karolinska Institutet (https://ki.se/en/km/infinigene). Injection at E8.25 targets skin and neural plate before neurulation, allowing tracing of neural crest derivatives from neural plate to final destinations.

Cell isolation from mouse embryo

At E13.5 the embryos were eviscerated. The sex of the harvested embryos was not determined. The trunk region of the embryo was dissected from the skin, liver and intestine. The dissociation of cells was the same as for human fetal cells. Single EGFP+ cells were sorted on a BD Influx FACS equipped with a 140 μm nozzle and a cooling unit with a sample temperature of 4°C and collected into cold HBSS buffer. Collected cells were pelleted and subjected to single cell RNA sequencing library preparation.

RNA seq library preparation and sequencing

Library preparation was done in accordance with 10x Genomics Chromium Single Cell 3’ protocol for Reagent Kits v3.

For human fetal cells we aimed to recover 5000-8000 cells with Illumina NovaSeq 6000 Sequencing System (NovaSeq 6000 S1 Reagent Kit or NovaSeq 6000 S2 Reagent Kit were used). The read setup was the standard recommended from 10X Genomics: Read 1: 28 cycles (Cell barcode and UMI), i7 index: 8 cycles (Sample index), Read 2: 91 cycles. For mouse embryonic sympathoadrenal cells, we aimed to recover 8000 cells with similar parameters.

Analysis of single cell RNA-seq data

Basic process and dataset overview

The sequenced libraries from 13 samples were demultiplexed and mapped to GRCh38 human genome (bundle 3.0.0) with Cell Ranger software (version 3.0.1) with default parameters[58]. Cells with the number of detected genes <200 or >7000, or with the percentage of mitochondrial genes >25% were omitted. The combined dataset contained 99553 cells. The count matrices were analyzed with Seurat (version 3.0.2)[59], regressing out total number of detected molecules and mitochondrial fraction. PCA was estimated on top 2000 variable genes, generating UMAP embedding (top 15 PCs, resolution=0.3, Extended Data Figs. 1,2), and clustering (top 10 PCs).

Subgroups of cells from the combined dataset (e.g. cells from the adrenal medulla or cells from the adrenal cortex, mesenchyme and kidney) were analyzed and reclustered in a similar way as the original dataset. Cell cycle was regressed out by a standard approach: first, for each cell, S phase score and G2M phase score were estimated with CellCycleScoring function from Seurat, next, these scores were regressed out.

The cell-specific doublet scores were estimated with scrublet package (version 0.2.1)[60], separately for each dataset.

Focused analysis of SCP, chromaffin and sympathoblast populations

Cells with the scrublet score <0.2 were used. To attenuate the impact of the cell cycle, genes annotated in the cell-cycle relevant GO categories (GO:0000278, GO:0006260, GO:0007049) were omitted. The samples were aligned with conos package (version 1.3.0)[61] with parameters k=30, k.self=5 and 30 PCA components. Extra-adrenal samples were excluded (week 12 embryo 2 biopsy 1). Only samples with sufficient number of cells in adrenal medulla (>50 cells) were used. UMAP embedding was generated based on the joint graph (using 200 neighbors for each point). To generate a more conservative embedding in the Supplementary Fig. 8d, we additionally omitted all genes showing correlation (Pearson r>0.8) with the core cell cycle markers (MKI67, TOP2A, TUBB4B). We also removed a small potentially contaminating population of cells expressing STAR gene (corresponding to cortex signature and showing the signs of cortical RNA contamination). Conos integration was then performed (with parameters k=30, k.self=5 and 30 PCA components) and UMAP embedding of the graph (using 200 neighbors).

To visualize the cells that co-express transcriptional programs of SCPs, sympathoblasts and chromaffin cells, we defined literature-based set of markers for each of these three cell populations (ternary plot Fig.2f; genes listed on heatmaps in Fig. 2h and 3b). We calculated the aggregate expression score of each signature in each cell, as a total expression magnitude of all genes in a signature. Scores were visualized using Ternary R package (version 1.1.4). To illustrate gene expression changes along the transition from SCPs to sympathoblasts, we selected the cells in which the expression of chromaffin signature was not greater than 70% of the sum of all three signatures (chromaffin, sympathoblasts and SCP). Cells were then ordered according to the ratio of the sympathoblast signature to the sum of the SCP and the sympathoblast signature (Fig. 2h). Co-expression between the two signatures was defined for each cell as a product of the two gene expression signatures. The analysis of transitions from SCPs to chromaffin cells (Fig. 2h) and sympathoblasts to chromaffin cells (Fig. 3b) was performed in the same way.

Cell fate transitions analyses

To analyze the sequence of gene expression changes along the transitions between cell subpopulations, we first fit a principal graph with monocle3 package (version 0.1.3)[62]. For a given bridge, we identified two graph nodes so that the path between them traversed the transition between cell types. With igraph (version 1.1.0)[63], we selected the shortest path between the given two nodes as a subgraph. We filtered the cells whose nearest point of the trajectory belonged to the identified path. To each of these cells, we assigned pseudotime along the selected path with monocle3. The genes that significantly change their expression along the transition were found with graph_test function from monocle3. To visualize the transitions as heatmaps, the expression of each gene of interest was smoothed as a spline (splines::ns in R with 3 degrees of freedom) with monocle.

To find the genes that specifically characterize the transitions (B) but not the cell population before (A) and after (C) the transitions, we implemented a test that checked if the gene expression in the transition (B) can be explained by a linear mixture of expressions of that gene in populations A and C. Gene expression profile of each cell in B was modelled as a linear combination of mean gene expression profiles in A and B with lm function in R without intercept term (lm(B~meanA+meanC-1)). For each gene in each cell in bridge B, we calculated the magnitude of the residuals non explained by the model. We next calculated mean residuals across all cells in bridge B, and normalized them to the standard deviation of the expression of a given gene. The obtained normalized mean residual values were used to prioritize the genes the genes with distinctive expression pattern in the bridge population.

To perform RNA-velocity analysis (Fig.2b,f, 4a), we estimated exonic and intronic gene counts with velocyto (version 0.17.13)[64], masking repeat regions. The resulting loom files were analyzed with scvelo package (version 0.1.24)[65], using dynamical model.

Analysis of bulk data and patient survival

Bulk gene expression in neuroblastomas was obtained from GEO database (GSE49711). We focused on the tumors of stages 3 and 4, stratifying tumors based on MYCN amplification status. Survival curves for tumors with low (lower quartile) and high (upper quartile) expression of a given gene was compared using survminer R package. P-values were calculated with the two-sided log-rank test between Kaplan–Meier curves for the compared groups. For non-MYCN-amplified tumors we analyzed 353 samples (88 samples in each of quartile) and for MYCN-amplified tumors – 88 samples (22 samples in each quartile). We next defined transcriptomic signatures of individual cell types as the total expression of the given cell type’s markers. Separate analyses were performed for the signatures including and excluding the genes associated with cell cycle. Survival analysis for the expressions of gene signatures was performed in the same way as for gene expressions. To estimate how robust the gene signatures were, bootstrap resampling was carried out.

Integration with neuroblastoma

We merged embryonic single cell transcriptomes with neuroblastoma transcriptomes from 2 samples published in Oslen et al., and 6 samples published in Dong. et al. to investigate if the cell types in tumors match those in normal human embryo. The datasets were merged using the conos package (version 1.3.0)[61] with default parameters. We additionally performed label transfer of the cell types from human embryo to the tumor cells.

Extraction of cloneIDs and clone calling for single cell RNA-seq

Raw 10X Genomics Chromium Version 3 sequencing data were preprocessed with Cell Ranger v3.0.1. Cell Ranger was configured to use a custom reference, consisting of the GRCm38 (mm10) genome and an additional sequence representing the H2B-EGFP-N transgene, in which the cloneID region was marked with “N” wildcard characters. The resulting BAM file of aligned sequencing reads was then processed with TREX, a custom Python tool for cloneID extraction and clone calling[66]. TREX uses reads from cells that align to the H2B-EGFP-N transgene. CloneIDs are recovered from those alignments that cover the masked cloneID region. If soft clipping is encountered at one of the bases adjacent to the region, the alignment is assumed to continue ungapped into the region. All cloneIDs with identical UMIs that come from the same cell (have the same cellID) were collapsed to a consensus sequence. To error-correct cloneIDs, they were clustered using a Hamming distance of at most five as linking criterion (with single linkage). For each cluster, all its cloneIDs were replaced with the cloneID occurring most frequently in that cluster. From the resulting final cellID-cloneID combinations, those that are only supported by one UMI and one read were discarded. We also removed the cloneIDs that were supported by only one UMI, and had a high frequency in another cell, as we assumed such cloneIDs to be contaminations. This cleaned data was transformed into a count matrix showing UMI counts for each cloneID in each cell. Clonally related cells were identified as described previously[67]. Briefly, the Jaccard similarity between each pair of cloneID expressing cells was calculated, and a Jaccard score of 0.7 was used as a cut-off to define related cells. Clones were defined as groups of two or more related cells.

Two types of tests for conditional independence were carried out (Fig. 4i): 1) Kernel-based conditional independence (KCI) test was performed to estimate the p-value to observe the given distribution of clones among the three cell types under the null hypothesis that the events ‘a clone contains sympathoblasts’ and ‘the clone contains SCPs’ were independent given that ‘the clone contains chromaffin cells’ (Fig. 4i, Model 1) or the events ‘a clone contains chromaffin cells’ and ‘the clone contains SCPs’ were independent given that ‘the clone contains sympathoblasts’ (Fig. 4i, Model 2). KCI[68] tests were carried out using CondIndTests R package (version 0.1.5); 2) χ[2]-based conditional independence tests were carried out by comparing the deviance of two glm models (Poisson family with log link) using chi squared test (one-tailed p-value). A model that predicted the abundance of a given clone in sympathoblasts based on its abundance in chromaffin cells was compared to a model that predicted the abundance of a given clone in sympathoblasts based on its abundance in chromaffin cells and in SCPs (Model 1). A model that predicted the abundance of a given clone in chromaffin cells based on its abundance in sympathoblasts was compared to a model that predicted the abundance of a given clone in sympathoblasts based on its abundance in chromaffin cells and in SCPs (Model 2). Statistical comparison between Model 1 and Model 2 (Figure 4i) and between Model 1 and Model 3 (Extended Data Figure 8h) was performed similarly using chi-squared likelihood ratio test.

Mouse genetic lineage tracing

The 2-5-month-old females were placed in time-controlled mating, and the day of vaginal plug presence was considered as E0.5. In the case of Ret tracings, heterozygotes for both the Cre and the reporter R26R were used. In the case of Plp1 tracings, embryos were heterozygotes for Cre and homozygous for R26R, which was previously analyzed extensively in another study and with lack of leakage[54] in the sympathoadrenal system. Induction of recombination was initiated by intra peritoneal (i.p.) injection of tamoxifen (Sigma, T5648) dissolved in corn oil (Sigma, 8267) in pregnant females (0.1 mg/g of body weight) at E10.5 and E11.5. The embryos were collected at E15.5. The sex of the harvested embryos was not determined. Statistical analysis was done in GraphPad Prizm 8.1.1(330).

Immunohistochemistry (IHC) and antibodies

Samples were 4% PFA fixed, sucrose cryopreserved, embedded in OCT and transversely sectioned at 14 um. The samples were brought to room temperature, dried and subjected for an antigen retrieval with 1x Dako Target Retrieval solution (Dako, S1699). Slides were washed with PBS containing 0.1% Tween-20 (Sigma, P9416) (PBST) for 10 min x 3 and incubated with primary antibody solution in PBST in humidified chamber at room temperature (RT) overnight. Next, the slides were washed in PBST for 10 min x 3 at RT and incubated with secondary antibody solution and DAPI (Thermo Fisher Scientific, 1:10,000, #D1306) in PBST for 90 min in humidified chamber RT, followed by washes in PBST for 10 min x 3 at RT. Slides were mounted with Mowiol (Sigma, 81381).

The primary antibodies were used in the dilutions: rabbit anti-TH (1:1000, Pel-Freez Biologicals, #P40101-150), sheep anti–TH (1:2000, Novus Biologicals, #NB300-110), goat anti-human SOX10 (1:800, R&D Systems, #AF2864), mouse anti-ISL1 (1:250, DSHB, clone 39.4D5-s), rabbit anti-CART (1:500, Phoenix Pharma., #H-003-62), rabbit anti-KI67 SP6 (1:500, Thermo Scientific, # MA5-14520), mouse anti-HuC/D (1:300, Invitrogen, #A21271), rabbit anti-HAND2 (1:500, Abcam, #EPR19451), mouse anti-neurofilament medium chain (2h3) (1:300, DSHB, #P12839 - NFM_RAT), rat anti-CD31 (1:500, BD, #MEC13.3), rabbit anti-NPY (1:300, Immunostar, 22940), rabbit anti-VMAT2, (1:300, Phoenix Europe GmbH, H-V008), rabbit anti-CHGB (1:500, Synaptic Systems, 259 103), chicken anti-NF200 (1:2000, Abcam, ab4680), rabbit anti-PRPH (1:500, Chemicon, #AB1530), goat anti-GFP (1:500, Abcam, #ab6662), chicken anti- GFP (1:500, Aves Labs Inc., #GFP-1020), chicken anti-mCherry (1:1000, EnCor Biotech, #CPCA-mCherry).

Primary antibodies were detected with Alexa fluor-488, -555, -647 conjugated secondary antibodies raised in donkey against mouse, rabbit or goat (1:1000, Molecular Probes, ThermoFisher Scientific). For immunostaining of formalin-fixed paraffin embedded tissues, sections were stained on a Ventana Benchmark 2 automated stainer using a peroxidase reporter after heat-based antigen retrieval in citrate buffer, ph 5.

RNA scope

RNAscope® on cryosections was performed using the RNAscope® Fluorescent Multiplex Detection Reagent Kit, while sections of FFPE postnatal adrenal glands were hybridized using the RNAscope® Multiplex Fluorescent Detection Kit v2[69]. Probes used against human specific STMN2 (#525211-C3), EYA4 (#510931-C3), PENK (#548301), TH (#441651-C2), CHGA (#311111-C2), VIP (#452751), ISL1 (#478598-C2), and mouse specific Elavl2 (#491961-C2). For combination with IHC, the RNA scope was done first, followed by IHC without antigen retrieval.

Reproducibility

For all immunohistochemical and RNA detection experiments the staining patterns and morphological features of the cells visualized with different combinations of antibodies and probes were independently reproduced at least 3 times on different stages.

Fate convergence of the neural crest cell-derived sympathoblasts and sympathoblasts originating from SCP-dependent pathways in mouse and human.

a, Scheme depicting convergence of sympathoblasts generated from the neural crest cells (NCC, early pathway) with sympathoblasts generated by alternative late pathways in mouse and human on the UMAP embedding. The indicated transitions from NCCs are based on published results and are not captured in the present datasets due to the early onset. b, Ternary plot of extra-medullary and intra-medullary cells from week 12 of human adrenal gland (sequenced separately) indicate convergence of differentiated sympathoadrenal cells independently of their origin.