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

Abstract

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.

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).

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.

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 ISL1^high cells within the medulla at E18.5 developmental stage. Note CART⁺/HUC/D⁺ and ISL1^high 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.

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. 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.

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. 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).

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χ² – p-value of the likelihood ratio test comparing Model 3 and Model 1 (one-tailed chi-squared probability; H₀: 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.

Extended Data Fig. 9. 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.

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).

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).

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 ), showing RNA velocity analysis (for the week 9 embryo 2^nd sample and week 9 embryo 3^rd 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.

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.

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, H₀: 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 χ² probability (Model 1 H₀: clone abundance in sympathoblast cells is independent of its SCP abundance given its chromaffin abundance; Model 2 H₀: clone abundance in chromaffin cells is independent of its SCP abundance given its sympathoblast abundance). Pχ² – likelihood ratio test comparing Model 2 and Model 1 (one-tailed chi-squared probability; H₀: Model 1 is not more likely than Model 2). n=322 independent clones were analyzed. ✗ – model is rejected, ✓ - model is accepted.

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^CreERT2-/+;Rosa26R^YFP/YFP+ 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^YFP+ 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^CreERT2-/+;Rosa26R^TOMATO/+ animals induced at E10.5 and analyzed at E15.5. Filled yellow arrowheads show Ret^TOMATO+ 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^TOM+ 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).

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.. 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. (left) and Dong et al. (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.

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.