1. Developmental Biology
Download icon

An atlas of neural crest lineages along the posterior developing zebrafish at single-cell resolution

Research Article
  • Cited 4
  • Views 2,719
  • Annotations
Cite this article as: eLife 2021;10:e60005 doi: 10.7554/eLife.60005

Abstract

Neural crest cells (NCCs) are vertebrate stem cells that give rise to various cell types throughout the developing body in early life. Here, we utilized single-cell transcriptomic analyses to delineate NCC-derivatives along the posterior developing vertebrate, zebrafish, during the late embryonic to early larval stage, a period when NCCs are actively differentiating into distinct cellular lineages. We identified several major NCC/NCC-derived cell-types including mesenchyme, neural crest, neural, neuronal, glial, and pigment, from which we resolved over three dozen cellular subtypes. We dissected gene expression signatures of pigment progenitors delineating into chromatophore lineages, mesenchyme cells, and enteric NCCs transforming into enteric neurons. Global analysis of NCC derivatives revealed they were demarcated by combinatorial hox gene codes, with distinct profiles within neuronal cells. From these analyses, we present a comprehensive cell-type atlas that can be utilized as a valuable resource for further mechanistic and evolutionary investigations of NCC differentiation.

Introduction

Unique to vertebrates, neural crest cells (NCC) are an embryonic stem cell population characterized as transient, highly migratory, and multipotent. Following their birth from the dorsal neural tube, NCCs migrate extensively, dorsolaterally or ventrally along the main axial levels of the embryo; the cranial, vagal, trunk, and sacral regions (Graham et al., 2004; Le Douarin and Teillet, 1974). Depending on the axial level of their origination, NCCs give rise to different cell types within many critical tissues, such as the cornea, craniofacial cartilage and bone, mesenchyme, pigment cells in the skin, as well as neurons and glia that comprise peripheral ganglia (Hutchins et al., 2018; Epstein et al., 1994; Kuo and Erickson, 2011; Hall and Hörstadius, 1988; Le Douarin and Kalcheim, 1999; Theveneau and Mayor, 2012; Williams and Bohnsack, 2015; Yntema and Hammond, 1954).

During their development, NCCs undergo dramatic transcriptional changes which lead to diverse cellular lineages, making their transcriptomic profiles highly dynamic (Simoes-Costa et al., 2014; Martik and Bronner, 2017; Soldatov et al., 2019; Williams et al., 2019). In support of the model that complex transcriptional programs govern NCC ontogenesis, gene regulatory networks involved in early development of NCCs into broad cell types have been studied at a high level using a combination of transcriptomics, chromatin profiling, and enhancer studies, especially during pre-migratory and early migratory NCC specification along cranial axial regions, across amniotes (Martik and Bronner, 2017; Simoes-Costa and Bronner, 2016; Green et al., 2015; Lumb et al., 2017; Williams et al., 2019; Hockman et al., 2019). For example, during pre-migratory stages the transcription factors FoxD3, Tfap2a, and Sox9 are important for NCC fate specification and in turn regulate the expression of Sox10, a conserved transcription factor that is expressed along all axial levels by early migrating NCCs and within many differentiating lineages (Sauka-Spengler and Bronner-Fraser, 2008; Martik and Bronner, 2017). Gene regulatory networks that are important for select NCC cell fates, like melanocytes and chondrocytes, have been well characterized (reviewed in Martik and Bronner, 2017). Recently, the regulatory circuitry behind glial, neuronal, and mesenchymal fates of vagal NCC was described (Ling and Sauka-Spengler, 2019) where Prrx1 and Twist1 have been described as key differentiation genes for mesenchymal fate. Despite this progress, however, comprehensive knowledge of the genes that are expressed and participate in NCC lineage differentiation programs during later phases of embryogenesis remains to be fully characterized, particularly for posterior tissues (reviewed in Hutchins et al., 2018). Indeed, altered gene expression during NCC differentiation can cause several neurocristopathies, such as DiGeorge syndrome, neuroblastoma, Hirschsprung disease, Auriculo-condylar syndrome, and Klein-Waardenburg syndrome (Barlow, 1984; Bolande, 1997; Brosens et al., 2016; Escot et al., 2016; Vega-Lopez et al., 2017; Wang et al., 2014), further highlighting the need to understand NCC spatiotemporal gene expression patterns during their differentiation into diverse cellular types.

Previous single-cell transcriptomic studies in zebrafish have laid a strong foundation to globally map early lineages of a majority of cell types through early to middle embryonic development (Wagner et al., 2018; Tambalo et al., 2020), and recently this has been extended into the larval stage (Farnsworth et al., 2020). With respect to zebrafish NCC development, the early embryonic window of 11–20 hr post fertilization (hpf) marks the stage of NCC specification and the emergence of their migratory behavior. Further development between 24 and 96 hpf represents the time when NCCs actively differentiate into their many derivatives (Rocha et al., 2020). Concerning the posterior NCC fates, however, many of these cells undergo differentiation programs during the embryonic to larval transition, a developmental stage that emerges between ~48 and 72 hpf. Transcriptomic analysis during this transitional phase would therefore enhance our understanding of the dynamic shifts in cell states that may regulate cellular differentiation programs.

In this study, we leverage the power of single-cell transcriptomics and curate the cellular identities of sox10-expressing and sox10-derived populations along the posterior zebrafish during development. We have utilized the Tg(−4.9sox10:EGFP, hereafter referred to as sox10:GFP) transgenic line to identify NCCs and their recent derivatives (Carney et al., 2006). Using sox10:GFP+48–50 hpf embryos and 68–70 hpf larvae, we identified eight major classes of cells: mesenchyme, NCC, neural, neuronal, glial, pigment, muscle, and otic. Among the major cell types, we annotated over 40 cellular subtypes. By leveraging in depth analysis of each time point separately, we captured the dynamic transition of several NCC fates, most notably we discovered over a dozen transcriptionally distinct mesenchymal subpopulations and captured the progressive differentiation of enteric neural progenitors into maturing enteric neurons. Using Hybridization Chain Reaction (HCR) and in situ hybridization, we validated the spatiotemporal expression patterns of various subtypes. By merging our 48–50 hpf and 68–70 hpf datasets, we generated a comprehensive atlas of sox10+ cell types spanning the embryonic to larval transition, which can also be used as a tool to identify novel genes and mechanistically test their roles in the developmental progression of posterior NCCs. Using the atlas, we characterized a hox signature for each cell type, detecting novel combinatorial expression of hox genes within specific cell types. Our intention is that this careful analysis of posterior NCC fates and resulting atlas will aid the cell and developmental biology communities by advancing our fundamental understanding of the diverging transcriptional landscape during the NCC’s extensive cell fate acquisition.

Results

Single-cell profiling of sox10:GFP+ cells along the posterior zebrafish during the embryonic and larval stage transition

To identify sox10-expressing and sox10-derived cells along the posterior zebrafish during the embryonic to larval transition, we utilized the transgenic line sox10:GFP (Figure 1; Carney et al., 2006; Kwak et al., 2013). Tissue posterior to the otic vesicle, encompassing the vagal and trunk axial region (Figure 1B), was dissected from 100 embryonic zebrafish at 48–50 hpf and 100 larval zebrafishes at 68–70 hpf. Dissected tissues were dissociated and immediately subjected to fluorescence-activated cell sorting (FACS) to isolate sox10:GFP+ cells (Figure 1B; Figure 1—figure supplement 1A,B). Isolated cells were then input into 10X Genomics Chromium scRNA-seq assays and captured at a depth of 2300 cells from the 48–50 hpf time point and 2580 cells from the 68–70 hpf time point (Figure 1C; Figure 1—figure supplement 1C). We performed cell filtering and clustering (Figure 1—figure supplement 1D–I) of the scRNA-seq datasets using Seurat (Butler et al., 2018; Stuart et al., 2019) to computationally identify cell populations based on shared transcriptomes, yielding 1608 cells from the 48–50 hpf time point and 2410 cells from the 68–70 hpf time point, totaling 4018 cells for final analysis (Figure 1—figure supplement 1C). We detected cell population clusters with transcriptionally unique signatures, as shown in heatmap summaries that revealed the top enriched gene signatures per cluster, with 19 clusters (0–18) from the 48–50 hpf time point (Figure 1—figure supplement 2A) and 23 clusters (0–22) from the 68–70 hpf time point (Figure 1—figure supplement 2B), totaling 42 clusters across both time points. Datasets were visualized with the t-Distributed Stochastic Neighbor Embedding (tSNE) method, which spatially grouped cells in each cluster, for both time points examined (Figure 1D,E). The top significantly enriched markers for each cluster at 48–50 and 68–70 hpf are provided in a table in Figure 1—source data 1.

Figure 1 with 5 supplements see all
Single-cell profiling strategy and cell population composition of posterior sox10:GFP+ cells from the posterior zebrafish during the embryonic to larval stage transition.

(A) Confocal image of a sox10:GFP+ embryo at 48 hpf; Hb: Hindbrain; Sc: Spinal cord. A: Anterior, P: Posterior, D: Dorsal, V: Ventral. Scale bar: 50 μM (B) Cartoon illustrations of a zebrafish embryo at 48–50 hpf and an early larval fish at 68–70 hpf depicted laterally to summarize the dissection workflow used to collect posterior sox10:GFP+ cells. (C) Schematic of the 10X Genomics Chromium and data analysis pipeline. (D) tSNE plots showing the arrangement of Clusters 0–18 and where the major cell types identified among sox10:GFP+ cells arrange in the 48–50 hpf dataset. (E) tSNE plots showing the arrangement of Clusters 0–22 and where the major cell types identified among sox10:GFPcells arrange in the 68–70 hpf dataset. (F,G) Dot plots of the identifying gene markers for each major cell type classification in the 48–50 hpf and 68–70 hpf datasets, respectively. Dot size depicts the cell percentage for each marker within the dataset and the color summarizes the average expression levels for each gene.

Figure 1—source data 1

List of marker genes per cluster in the sox10:GFP scRNA-seq datasets.

Tables reporting the Seurat output for genes for each cluster at the 48–50 hpf (Tab 1) and 68–70 hpf (Tab 2) time points, including p-values (<0.01), average log-fold change (≥0.25), adjusted p-values (≤1.0), pct.1 summarizing proportion of cells expressing the individual gene in the cluster, pct.2 showing the proportion of cells expressing the individual gene in all other clusters in the dataset.

https://cdn.elifesciences.org/articles/60005/elife-60005-fig1-data1-v1.xlsx.zip
Figure 1—source data 2

Table summarizing the top identity markers used for major cell type and subtype cellular classifications for each cluster at 48–50 and 68–70 hpf.

https://cdn.elifesciences.org/articles/60005/elife-60005-fig1-data2-v1.pdf

Major classification of sox10:GFP+ cell states

To assess the proliferative state of sox10:GFP+ cells, we determined their G1, S or G2/M phase occupancy, based on expression of proliferative cell cycle marker genes (Figure 1—figure supplement 3I). At 48–50 hpf, 52% of sox10:GFP+ cells were in G1 phase, 31% were in the S phase and 17% in G2/M phase (Figure 1—figure supplement 3G), collectively indicating that 48% of the cells in the 48–50 hpf time point were proliferative. At 68–70 hpf, 64% of cells were in G1 phase, 24% of cells were in the S phase and 12% in G2/M phase (Figure 1—figure supplement 3G), indicating that 36% of the cells were proliferative. The cell cycle occupancy distributions were visualized in tSNE plots, revealing congregations of proliferative and non-proliferative sox10:GFP+ cells (Figure 1—figure supplement 3A,B); aurkb and mcm3 confirmed general occupancy in the G2/M and S phase (Figure 1—figure supplement 3C–F). Together, these data of cell cycle state reflect a general decrease in proliferative cells among sox10:GFP+ populations between 48 and 70 hpf, in agreement with prior observations (Rajan et al., 2018).

Using a combination of gene expression searches of literature and bioinformatics sources, examination of the scRNA-seq transcriptomes indicated that sox10:GFP+ cells exist in several major cell type categories based on the expression of signature marker genes (Figure 1F,G; Figure 1—figure supplement 2E,F). These major cell type categories included: neural, neuronal, glial, mesenchyme, pigment cell, NCC, otic, and muscle; their respective fraction of the datasets was also calculated (Figure 1D–G; Figure 1—figure supplement 2C,D; Figure 1—figure supplement 3H). Neuronal identity refers to cells predominantly expressing neuron markers, such as elavl3/4, while neural cells are defined by a multipotent state with potential towards fates of either glial or neuron identity, and marked by expression of factors such as sox10, dla, and/or ncam1a. Notably, mesenchyme identity represented the largest proportion of the datasets at 61% and 53% of the cells at 48–50 and 68–70 hpf, respectively (Figure 1—figure supplement 3H). Mesenchyme clusters were identified by a combination of mesenchymal gene markers including twist1a/b and prrx1a/b (Soldatov et al., 2019). In addition, cells with an otic vesicle and muscle identity were detected (Figure 1D–G; Figure 1—figure supplement 3H; Figure 1—figure supplement 4), as has previously been described in the sox10:GFP line (Carney et al., 2006; Rajan et al., 2018; Rodrigues et al., 2012; Kwak et al., 2013). Overall, major cell type cluster identities and their top signature marker genes are summarized in Figure 1—source data 2.

Annotation of cellular types among posterior sox10:GFP+ cells

Closer analysis of the 42 cluster gene signatures among the two time points allowed us to annotate cellular identities in more detail (Figure 1—source data 2). Indeed, we identified previously described NCC-derived cell types. For example, the sox10:GFP line has been shown to transiently label sensory dorsal root ganglion (DRG) progenitors between the first and second day of zebrafish development (McGraw et al., 2008; Rajan et al., 2018). We observed sensory neuronal/DRG gene expression in Cluster 17 at 48–50 hpf (Figure 1—source data 2; Figure 1—figure supplement 5) by the markers neurod1, neurod4, neurog1, six1a/b, elavl4 (Carney et al., 2006; Delfino-Machín et al., 2017). Additionally, we identified other NCC-derivatives, including mesenchymal cells (Le Lièvre and Le Douarin, 1975; Kague et al., 2012; Soldatov et al., 2019; Ling and Sauka-Spengler, 2019), pigment cells (Reedy et al., 1998; Higdon et al., 2013), and enteric neurons (Kelsh and Eisen, 2000; Kuo and Erickson, 2011; Lasrado et al., 2017), which we describe in further detail for both time points in Figures 25.

Distinct pigment cell populations are present among sox10:GFP+ cells during embryonic to larval transition.

(A) Cartoon schematic depicting the model for neural crest delineation into pigment cell lineages and the genes that were used to identify each pigment cell population. (B) Dot plot identifying melanophore markers within the 48–50 hpf dataset. Dot size depicts the cell percentage for each marker within the dataset and the color summarizes the average expression levels for each gene. (C) tSNE plots depicting melanophore signature in the 48–50 hpf dataset. Relative expression levels are summarized within the color keys, where color intensity is proportional to expression level of each gene depicted. (D) Dot plot showing distinct pigment chromatophore markers within the 68–70 hpf dataset. Dot size depicts the cell percentage for each marker within the dataset and the color summarizes the average expression levels for each gene. M: melanophore markers; X: xanthophore markers; I: iridophore markers. (E–G) tSNE plots revealing the location of melanophores (E), xanthophores (F), and iridophores (G) in the 68–70 hpf dataset. Relative expression levels are summarized within the color keys, where color intensity is proportional to expression level of each gene depicted. (H) HCR against mitfa and tfec at 48–50 hpf reveals mitfa+ melanophores (white arrowhead) and mitfa+/tfec+ pigment progenitors (red arrowhead). Cropped panels show individual fluorescent channels. (I) HCR against mitfa and tfec at 68–70 hpf presents mitfa+ melanophores (white arrowhead), tfec+ iridophores (blue arrowhead), and mitfa+/tfec+ pigment progenitors (red arrowhead). Cropped panels show individual fluorescent channels. (J) HCR against mitfa and xdh at 68–70 hpf shows mitfa+/xdh+ xanthophores (orange arrowhead). Cropped panels show individual fluorescent channels. Scale bar in H-J: 50 µm.

Figure 2—source data 1

Melanophore populations, shared and unique genes at 68–70 hpf.

Table summarizing the genes (‘elements’) exclusively found in Cluster 4, exclusively found in Cluster 18, and genes found in both Clusters 4 and 18 at the 68–70 hpf time point scRNA-seq dataset.

https://cdn.elifesciences.org/articles/60005/elife-60005-fig2-data1-v1.xlsx.zip

Identification of pigment cell types within sox10:GFP+ scRNA-seq datasets

With robust genetic lineage details published on pigment cell differentiation in zebrafish (Kelsh, 2004; Lister, 2002; Quigley and Parichy, 2002), we sought to validate the scRNA-seq datasets by assessing if we could resolve distinct pigment cell populations. Pigment cell development has been broadly studied in the developing zebrafish, where NCCs give rise to three distinct chromatophore populations: melanophores, xanthophores, and iridophores (Figure 2A). Our annotation analysis of sox10:GFP+ scRNA-seq clusters revealed expression of pigment cell lineage gene markers (Figure 2B–J; Figure 1—figure supplement 2). At 48–50 hpf, melanophores were detected in Cluster eight based on expression of mitfa, dct, tyrp1b, and pmela (Du et al., 2003; Lister et al., 1999; Ludwig et al., 2004; Quigley and Parichy, 2002Figure 1F; Figure 2A), also reflected by the dot and tSNE plots (Figure 2B,C). At 68–70 hpf, we resolved discrete pigment cell populations that included xanthophore, iridophore, and two distinct melanophore clusters (Figure 2D–G, Figure 1—source data 2). The xanthophores mapped to Cluster 15 and were enriched with xdh, aox5, pax7b, mitfa, and gch2 (Nord et al., 2016; Parichy et al., 2000; Saunders et al., 2019; Minchin and Hughes, 2008; Lister et al., 1999; Figure 2A,D,F). Cluster 16 was identified as iridophores, which presented the well-characterized markers: tfec, pnp4a, gpnmb, and atic (Higdon et al., 2013; Lister et al., 2011; Petratou et al., 2018; Petratou et al., 2019Figure 2A,D,G; Figure 1—source data 2). The use of cell cycle markers revealed that two different melanophore clusters at 68–70 hpf (Clusters 4 and 18) were present in different proliferative states (Figure 1—source data 2). While the majority of cells in Cluster four were in G1, Cluster 18 expressed S and G2/M markers, such as pcna and aurkb, suggesting this population to be proliferating melanophores (Figure 1—figure supplement 3B,E,F; Figure 2—source data 1).

At 68–70 hpf, we identified a pigment progenitor population, where iridophore and melanophore markers were co-expressed in Cluster 13 (Figure 2A,D). These undifferentiated pigment progenitor cells expressed tfec in combination with mitfa and have been described recently at 24, 30, and 48 hpf (Petratou et al., 2018). Additionally, Cluster 13 expressed tfap2e, gpx3, and trpm1b (Figure 1—source data 2) whose expression patterns have been previously reported in pigment progenitors (Saunders et al., 2019). Finally, a population of pigmented muscle (Cluster 9) was also found with a weak melanophore signature coupled with expression of the muscle markers ckmb, tpma, tnnc2, and tnnt3b (Figure 1—source data 2; Figure 1—figure supplement 4).

We next performed whole mount HCR to assess the spatial co-expression of mitfa, tfec and xdh. When examining mitfa and tfec at 48–50 hpf (Figure 2H), we detected sox10:GFP+ cells that expressed mitfa, identifying the melanophores (Figure 2H; white arrowhead), and cells that expressed both mitfa and tfec, defining the pigment progenitors (Figure 2H; red arrowhead). At 68–70 hpf, we confirmed the four distinct pigment populations we identified through Seurat (Figure 2B–G):GFP+ melanophores expressing mitfa only (Figure 2I; white arrowhead), iridophores only expressing tfec (Figure 2I; blue arrowhead), and pigment progenitors expressing both mitfa and tfec (Figure 2I; red arrowhead) were detected. When examining xdh and mitfa expression patterns, sox10:GFP+ xanthophores were found to be expressing both markers (Figure 2J; orange arrowhead).

Taken together, the above-described results regarding pigment cell expression patterns validates that the sox10:GFP+scRNA-seq datasets captured discrete NCC-derived populations, and coupled with HCR analysis, shows we are able to validate these cell populations in vivo.

Mesenchyme in the posterior embryo and larvae exists in various transcriptionally-distinct populations

Heatmap analysis of gene expression groups depicted that mesenchyme cells clustered together globally within the datasets (Figure 1—figure supplement 2C,D; Figure 3A,B), with twist1a expression broadly labeling all mesenchyme cells (Figure 1—figure supplement 2E,F). In addition to twist1a, mesenchyme cells also expressed prrx1a/b, twist1b, foxc1a/b, snai1a/b, cdh11, sparc, colec12, meox1, pdgfra (Figure 3A,B), and other known mesenchymal markers such as mmp2 (Figure 1—source data 2Janssens et al., 2013; Theodore et al., 2017). In whole mount embryos at 48 hpf, we observed broad expression of foxc1a and mmp2 along the posterior pharyngeal arches and ventral regions of the embryo via in situ hybridization (Figure 3—figure supplement 1C,D; arrowheads), confirming their expression territories within posterior-ventral mesenchymal tissues.

Figure 3 with 1 supplement see all
Global analysis of mesenchyme cell signatures among sox10:GFP+ cells.

(A,B) A heatmap of signature mesenchyme identity genes within the major cell type classified cells at 48–50 and 68–70 hpf, respectively. Relative expression levels within each cluster is summarized within the color key, where red to blue color indicates high to low gene expression levels. (C) A cluster tree depicting the relationship between general and chondrogenic mesenchyme cellular subtypes. (D) Violin plots summarizing the expression levels for select mesenchyme identity markers within individual clusters at the 48–50 and 68–70 hpf time points, respectively. Data points depicted in each cluster represent single cells expressing each gene shown. (E,L) tSNE plots depicting the co-expression of twist1a (blue) and barx1 (red) or prrx1b (blue) and barx1 (red) in the 48–50 and 68–70 hpf datasets, respectively. Relative expression levels are summarized within the color keys, where color intensity is proportional to expression level of each gene depicted. (F–K) Whole mount HCR analysis reveals the spatiotemporal expression of prrx1b (F), barx1 (G), twist1a (H), sox10:GFP (J) in 48 hpf embryos. (I) A merge of barx1, prrx1b, and twist1a is shown. (K) A merge of barx1, prrx1b, twist1a, and sox10:GFP is shown. White arrowheads denote expression in posterior pharyngeal arch, while yellow arrowheads highlight fin bud expression. (M–R) Whole mount HCR analysis reveals the spatiotemporal expression of prrx1b (M), barx1 (N), twist1a (O), sox10:GFP (Q) in 68 hpf embryos. (P) A merge of barx1, prrx1b, and twist1a is shown. (R) A merge of barx1, prrx1b, twist1a, and sox10:GFP is shown. White arrowheads denote expression in posterior pharyngeal arch, while yellow arrowheads highlight fin bud expression. Ot: otic; Fb: Fin bud. Scale bar: 100 µm.

Further analysis revealed various transcriptionally-distinct populations were present among the sox10:GFP+ cells with a mesenchymal identity. Among these, we detected nine clusters with chondrogenic signatures—identified by expression of mesenchymal signature genes, as well as the chondrogenic markers barx1 and/or dlx2a (Sperber et al., 2008; Sperber and Dawid, 2008; Ding et al., 2013; Barske et al., 2016Figure 3C,D). Feature plot exports revealed distribution of the chondrogenic cells (barx1+) in relation to all other mesenchyme (prrx1b+, twist1a+) cells in the datasets (Figure 3E,L). Within the nine chondrogenic clusters, we discovered gene expression indicative of heterogeneous cell states, ranging from proliferative, progenitor/stem-like, and migratory to differentiating signatures (Figure 1—source data 2). All other mesenchyme clusters (seven in total) were also classified into various progenitor and differentiation categories. Among these categories, the clusters expressed either proliferative progenitor markers, differentiation signatures, or general migratory mesenchymal markers (Figure 1—source data 2). Cluster 14 at 48 hpf and Clusters 1 and 7 at 68–70 hpf exhibited a general mesenchymal signature, but also expressed fin bud marker genes (hand2, tbx5a, hoxd13a, prrx1a, Figure 1—figure supplement 5Yelon et al., 2000; Lu et al., 2019; Nakamura et al., 2016; Feregrino et al., 2019). Additionally, visualization of clusters with general mesenchyme and chondrogenic identities using a cluster tree highlighted potential similar clusters between the time points (Figure 3C). For example, the cluster tree showed proximal location of Cluster 8 at 68–70 hpf and Cluster 2 at 48–50 hpf, which we noted contained clear proliferative chondrogenic gene signatures (Figure 3D; Figure 1—source data 2).

To confirm the spatial co-expression of prrx1b, twist1a, and barx1 within sox10:GFP+ tissues, we utilized HCR analysis (Figure 3F–K,M–R). Corroborating our analysis that mesenchyme-identity populations contained both general and chondrogenic signatures, we found the co-expression of prrx1b, twist1a, and barx1 within sox10:GFP+ domains along the posterior pharyngeal arches (white arrowheads) and fin bud mesenchyme (yellow arrowheads) at both time points (Figure 3F–K,M–R).

Overall, the above-described analyses indicate that sox10:GFP+ mesenchymal cells in the posterior zebrafish exhibit various transcriptional states between the embryonic to larval transition and suggest that posterior mesenchyme exists in various subpopulations during its differentiation.

sox10-derived cells during the embryonic to early larval transition reveal enteric progenitor to enteric neuron progression

At 48–50 hpf, cells with NCC identity were notably detected in Cluster 5, defined by expression of the core NCC markers sox10, foxd3, crestin, and tfap2a (Figure 1F; Figure 4, Figure 1—source data 2; Dutton et al., 2001; Luo et al., 2001; Knight et al., 2003; Stewart et al., 2006). Moreover, Cluster five was found to contain various other genes previously shown to be expressed in zebrafish NCCs; including, vim, snai1b, sox9b, zeb2a, mych, and mmp17b (Figure 4DCerdà et al., 1998; Heffer et al., 2017; Hong et al., 2008; Leigh et al., 2013; Van Otterloo et al., 2012; Wang et al., 2011; Rocha et al., 2020). We reasoned that many of the NCCs had started their respective differentiation programs and were beginning to assume specified lineage profiles. Therefore, we sought to determine if the NCC cluster also contained gene expression profiles of known differentiating NCC types along the posterior body, such as enteric progenitors, also known as enteric neural crest cells (ENCCs).

Enteric neural crest cells and differentiating enteric neurons are present among posterior sox10:GFP+ cell populations.

(A) tSNE feature plots reveal expression of core neural crest cell markers sox10, foxd3, crestin, and tfap2a mapping to the neural crest cell cluster (red arrowhead). (B) tSNE feature plots depict expression of the enteric neural crest cell markers phox2bb, ret, ngfrb and gfra1a within the neural crest cell cluster (red arrowhead). Relative expression levels are summarized within the color keys in (A) and (B), where color intensity is proportional to expression level of each gene depicted. (C) A heatmap reveals expression levels of enteric neural crest cell markers across the eight major cell populations captured in the 48–50 hpf data set (color key denotes cells types represented in color bar on top of heatmap). Neural crest cell cluster highlighted in black rectangle. Relative expression levels within each major cell type cluster is summarized within the color key, where yellow to magenta color indicates high to low gene expression levels. (D) Dot plot of expanded list of neural crest (green line) and enteric neural crest (purple line) cell markers across each major cell type within 48–50 hpf data set. Dot size depicts the cell percentage for each marker within the data set and the color summarizes the average expression levels for each gene. (E,F) Whole mount HCR analysis of 48 hpf embryos reveals co-expression of the enteric neural crest cell markers phox2bb, ngfrb, gfra1a, and crestin in (E), or foxd3, ngfrb, gfra1a and crestin in (F), within the developing gut (dashed outline). Top panels depict merged images of color channels for each HCR probe. Lower panels represent gray-scale images of each separated channel corresponding to the magnified region of foregut (gray rectangle). Arrowheads depict regions where all markers are found to be co-expressed. Hb: Hindbrain, Sc: Spinal cord, pLLg: posterior Lateral Line ganglia, LL: Lateral Line. A: Anterior, P: Posterior, D: Dorsal, V: Ventral. Scale bar: 50 μM. (G) tSNE feature plots reveal expression levels of enteric neuron markers elavl3, phox2bb, gfra1a, nos1, vipb, and ret, within a common region of a neuronal cluster (red arrowhead). Relative expression levels are summarized within the color keys, where color intensity is proportional to expression level of each gene depicted. (H) Dot plot depicts expression levels of pan-neuronal and enteric neuron specific markers across individual clusters generated within the original 68–70 hpf tSNE. Pan-neuronal markers found throughout Clusters 5 and 12, with enteric neuron markers most prominently expressed within Cluster 12. Dot size depicts the cell percentage for each marker within the data set and the color summarizes the average expression levels for each gene. (I) Whole mount HCR analysis depicts differentiating enteric neurons within the foregut region at 69 hpf co-expressing nos1, phox2bb, vipb, and elavl3 (yellow arrowheads). Anterior: Left, Posterior: Right. Scale bar: 50 μM.

ENCCs fated to give rise to the enteric nervous system (ENS), the intrinsic nervous system within the gut, express a combination of NCC and enteric progenitor marker genes over developmental time (reviewed in Nagy and Goldstein, 2017; Rao and Gershon, 2018), which occurs between 32 and 72 hpf in zebrafish (reviewed in Ganz, 2018). Enteric markers in zebrafish include sox10, phox2bb, ret, gfra1a, meis3, and zeb2a (Dutton et al., 2001; Shepherd et al., 2004; Elworthy et al., 2005; Delalande et al., 2008; Heanue and Pachnis, 2008; Uribe and Bronner, 2015). Therefore, we expected to capture a population of ENCCs within our 48–50 hpf dataset. Indeed, within Cluster 5 we observed expression of the enteric markers phox2bb, ret, gfra1a, meis3, sox10, and zeb2a (Figure 4B–D). Using whole mount in situ hybridization, we confirmed the expression of sox10 and phox2bb within ENCCs localized along the foregut at 48 hpf (Figure 3—figure supplement 1A,A’,B,B’; arrowheads). Furthermore, gene orthologs known to be expressed in ENCC in amniotes were detected within Cluster 5, such as ngfrb (orthologue to p75; Anderson et al., 2006; Wilson et al., 2004) and hoxb5b (orthologous to Hoxb5; Kam and Lui, 2015Figure 4B–D).

HCR analysis of 48 hpf embryos validated the co-expression profiles of several ENCC markers along the foregut (Figure 4E–F; foregut in gray box). We observed that a chain of crestin+ cells localized in the foregut contained a subpopulation of cells expressing ngfrb, phox2bb, and gfra1a (Figure 4E; white arrowheads), or expressing foxd3, ngfrb, and gfra1a (Figure 4F; white arrowheads). Together, these HCR data confirm that ENCC markers are co-expressed along the zebrafish gut.

We next asked if we could resolve discrete differentiating enteric neurons over time. Within the 68–70 hpf zebrafish, ENCCs have yet to finish their migratory journey along the gut and have previously been shown to exist in varying stages of neuronal differentiation, where the earliest differentiating neurons are found in the rostral foregut, and the more proliferative undifferentiated ENCCs are continuing to migrate into the caudal hindgut (Elworthy et al., 2005; Olden et al., 2008; Harrison et al., 2014; Uribe and Bronner, 2015; Taylor et al., 2016). During early neuronal differentiation (68–70 hpf), ENCCs display differential enteric progenitor gene expression patterns (Taylor et al., 2016) and neurochemical signatures representative of varying stages of neuronal differentiation and subtype diversification (Poon et al., 2003; Holmqvist et al., 2004; Uyttebroek et al., 2010). Zebrafish early differentiating enteric neurons have been characterized by the mRNA expression of sox10, phox2bb, gfra1a, fgf13b, and ret, as well as the immunoreactivity of Elavl3/4 (Shepherd et al., 2004; Heanue and Pachnis, 2008; Uyttebroek et al., 2010; Taylor et al., 2016). In addition, at this time, enteric neurons express multiple neurochemical markers, with Nos1 being most prominent (Olden et al., 2008; Uyttebroek et al., 2010), a finding consistent with studies performed within the amniote ENS (Hao and Young, 2009; Matini et al., 1995; Qu et al., 2008; Heanue et al., 2016). In light of these previous observations, our 68–70 hpf dataset was expected to contain the transcriptomes of ENCCs captured at various stages of their progressive differentiation into the diverse subtypes of the ENS.

We identified differentiating enteric neurons within the 68–70 hpf dataset based on the combinatorial expression of elavl3, phox2bb, ret, and gfra1a (Figure 4G), which mapped to the neural/neuronal major cell type regions of the dataset (Figure 1E), comprising Clusters 5 and 12 (Figure 1E). Transcripts that encode for the neurochemical marker nos1, and the neuropeptides vip and vipb, a paralogue to vip (Gaudet et al., 2011), were found in a subpopulation of enteric neurons localized to a distal group of the neuronal cluster, likely indicative of a differentiating enteric neuron subtype (Figure 4G; red arrows). We then queried for the presence of a combination of pan-neuronal and enteric neuron markers (Figure 4H). The pan-neuronal markers tuba2, elavl3, stx1b, and gng2/3 (Asakawa and Kawakami, 2010; Kelly et al., 2008; Park et al., 2000; Zheng et al., 2018), as well as the autonomic neuron markers, phox2a and phox2bb (Hans et al., 2013), were present in both Clusters 5 and 12 (Figure 4H). However, the enteric neuron markers, gfra1a, ret, hoxb5b, ngfrb, fgf13b, nos1, vipb, and vip were mostly confined to Cluster 12, suggesting that this cluster contained differentiating enteric neurons (Figure 4H). Indeed, whole mount HCR analysis validated the spatiotemporal co-expression of phox2bb, nos1, vipb, and elavl3 transcripts throughout the foregut of the zebrafish embryo by 69 hpf (Figure 4I; yellow arrowheads). These results suggest that elavl3+/phox2bb+ early differentiating enteric neurons are first seen in the foregut and display an inhibitory neurochemical gene signature, consistent with prior observations in zebrafish and mammalian ENS (Olden et al., 2008; Hao and Young, 2009).

In an effort to examine the enteric neuron populations with finer resolution, Clusters 5 and 12 were subset from the main dataset in Seurat, re-clustered and visualized using a tSNE plot, producing 5 Sub-clusters (Figure 5A). The gene markers from each new Sub-cluster are provided in Figure 5—source data 1. The previously mentioned enteric neuron markers, with the addition of etv1, a recently identified marker of an enteric sensory neuron type, intrinsic primary afferent neurons (IPANs) in mouse (Morarach et al., 2021), were queried and visualized using dot and feature plot allowing the identification of Sub-cluster 3 as a differentiated enteric neuron cluster (Figure 5A–C). nos1, vip, and vipb were enriched in Sub-cluster 3 (Figure 5B,C; Figure 5—figure supplement 1). Interestingly, while expressed at lower average levels than in Sub-cluster 3, the enteric combination markers were also present in Sub-cluster 1 (Figure 5B–C). Sub-cluster 1 formed a central point from which Sub-cluster 3 could be seen emanating as a distal population (Figure 5A). Sub-clusters 1 and 3 likely depict enteric neurons captured at different stages along their progressive differentiation.

Figure 5 with 2 supplements see all
Differentiating enteric neurons captured during key transitional stage of subtype diversification within 68–70 hpf sox10:GFP+ larval cells.

(A) tSNE plot reveals five distinct sub-clusters following the subset analysis and re-clustering of Clusters 5 and 12 from the 68–70 hpf data set. (B) Dot plot depicts expression levels of enteric neuron markers across resulting Sub-clusters. Each marker was expressed at low levels in Sub-cluster 1 and were found to be expressed at higher levels within Sub-cluster 3. (C) tSNE feature plots further depict the expression of enteric neuron markers by illustrating the levels and localization of expression within the Sub-cluster architecture. Feature plots supplement dot plot and demonstrate the prominent expression of enteric neuron markers within Sub-cluster 3, which appears to emanate from Sub-cluster 1. (D,E) Violin and feature plots reveal expression levels of acetylcholine-associated and excitatory neuron markers reported to distinguish enteric IPANs. These markers were found in a discrete pocket of cells forming the distal-most region of Sub-cluster 3 (red arrowhead). Violin data points depicted in each Sub-cluster represent single cells expressing each gene shown. (F) Graphical model summarizes expression patterns observed in 68–70 hpf data set and HCR validation. Common enteric neuroblast capable of diverging into subsequent lineages, IPAN, inhibitory neuron, and interneuron through lineage restricted gene expression. pbx3b promotes assumption of IPAN role through loss of nos1 and vipb and begins expressing calb2a, ache, and slc18a3a. (G) Whole mount HCR analysis reveals co-expression of IPAN marker genes, pbx3b and calb2a, and inhibitory neurochemical marker genes, vipb and nos1 (white arrowheads), within the foregut (dashed white line) at 68 hpf. Vesicular acetylcholine transferase, slc18a3a, was not observed in tandem with pbx3b but was co-expressed with calb2a, vipb, and nos1 (yellow arrowheads). Scale bar: 50 μM. (H) Feature plots reveal expression of opioid receptor genes, oprl1 and oprd1b, within the differentiated enteric neuron Sub-cluster 3. (I–N) Whole mount HCR analysis validates expression of oprl1 in combination with vipb and phox2bb (yellow arrowheads) in enteric neurons localized to the foregut region of a 68 hpf embryo. Scale bar: 10 μM.

Figure 5—source data 1

List of marker genes per Sub-cluster, following subset and re-clustering of enteric Clusters 5 and 12 at 68–70 hpf.

Table reporting the Seurat output for genes for each Sub-cluster (0–4), including p-values (<0.01), average log-fold change (.25), adjusted p-values (1.0), pct.1 summarizing proportion of cells expressing the individual gene in the Sub-cluster, pct.2 showing the proportion of cells expressing the individual gene in all other Sub-clusters in the sub-data set.

https://cdn.elifesciences.org/articles/60005/elife-60005-fig5-data1-v1.xlsx.zip
Figure 5—source data 2

List of enriched pathways within enteric neuron Sub-cluster 3 and genes present in specific opioid proenkephalin pathway identified following PANTHER Overrepresentation Test.

Source Data table contains two sheets. Enriched Gene Pathways sheet lists the enriched pathways identified following the use of the statistical test, Fisher’s exact test that compares the number of genes in a given pathway within the reference Danio rerio list to the number of genes in a given pathway within the enteric neuron Sub-cluster 3 list. This analysis produced 43 statistically overrepresented pathways. Column B: Number of genes within the given pathway present with Danio rerio reference list of 25888 genes. Column C: Number of genes within the given pathway present within Sub-cluster 3 gene list. Column D: Number of pathway genes expected to be present within the Sub-cluster 3 gene list based on the percentage of pathway genes present in the Danio rerio reference list. Column E: Denotes that more pathway genes were present in Sub-cluster 3 gene list than expected. Column F: Fold enrichment of Sub-cluster 3 pathway genes comparative to reference list. Column G: p-values calculated following Fisher’s exact test comparing expected number of pathway genes to number of genes pathway genes in Sub-cluster 3. Column H: p-value following Benjamini-Hochberg false discovery rate (FDR) correction. Opioid proenkephalin pathway sheet lists genes associated with this pathway that are present within the enteric neuron Sub-cluster 3 gene list, notably, opioid receptor oprd1b.

https://cdn.elifesciences.org/articles/60005/elife-60005-fig5-data2-v1.xlsx.zip

Given our hypothesis that the enteric neurons further along a differentiation program were localized to the distal tip of Sub-cluster 3, we asked whether this population of cells contained additional neurochemical or neuron subtype-specific differentiation genes. Within a small pocket of cells in Sub-cluster 3, we detected the expression of calb2a and pbx3b (Figure 5D,E; Figure 5—figure supplement 1), orthologous genes to Calb2 and Pbx3 that have previously been shown to denote adult myenteric IPANs in mammals (Furness et al., 2004; Memic et al., 2018), as well as the two acetylcholine associated genes, acetylcholine esterase (ache, Bertrand et al., 2001; Huang et al., 2019) and vesicular acetylcholine transferase (slc18a3a, Hong et al., 2013; Zoli and Berlin, 2000Figure 5E; red arrowheads). Recently, a scRNA-seq study performed in E15.5 mouse demonstrated the co-expression of Calb2, Pbx3, and Slc18a during ENS development (Morarach et al., 2021). When examining IPAN gene markers in zebrafish larvae, HCR analysis revealed co-expression of pbx3b, calb2a, vipb, and nos1 (Figure 5G; white arrowheads), or co-expression of slc18a3a, calb2a, vipb, and nos1 (Figure 5G; yellow arrowheads), in discrete differentiating enteric neurons within the foregut region of the zebrafish gut at 68 hpf. These data indicate that zebrafish differentiating enteric neurons express IPAN gene signatures during their development. Collectively, our observations suggest that the sox10:GFP+68–70 hpf dataset captured an emerging IPAN population during its transition, where both excitatory and inhibitory neurochemical markers were co-expressed. Therefore, our single-cell analysis in zebrafish suggests that the transcriptional emergence of specific enteric neuron subtypes may be conserved between vertebrate species.

In order to identify novel signaling pathways within the developing enteric neuron population, the significantly enriched gene list from Sub-cluster 3 was processed using gene ontology (GO) pathway enrichment analysis. We found that three opioid signaling pathways were among the top 10 highest fold enriched pathways (Figure 5—source data 2; Mi et al., 2019). These pathways contained the G-protein-coupled receptors, oprl1 and oprd1b, respectively representing nociception/orphan FQ (N/OFQ) peptide (NOP)-receptor and delta-opioid receptor (DOR) subtype members within the opioid receptor superfamily (Figure 5HDonica et al., 2013; Sobczak et al., 2014). Specifically, feature plots showed that expression of the opioid receptors was tightly confined within the pocket of enteric neuron progenitors we identified as undergoing sensory lineage-specification, suggesting the expression of opioid receptor genes within both excitatory and inhibitory neurons at the early stages of enteric neuron differentiation within the zebrafish ENS (Figure 5E and H). Confirming this suggestion using HCR analysis, we observed the combinatorial expression of oprl1, vipb, and phox2bb within migrating enteric neuron progenitors at 68 hpf along the foregut (Figure 5I–N). The presence of opioid receptors within immature enteric neurons undergoing lineage-specification helps us to better understand the complexity of early ENS signaling and highlights an area that requires further investigation.

Based on our observation that the enteric neuron population that comprised Sub-cluster 3 only made up one of five phox2bb+ Sub-clusters (Figure 5A–C), we suspected that the remaining Sub-clusters were made up of closely related autonomic neurons. In order to better visualize specific differences between the Sub-clusters, we viewed them using UMAP (Becht et al., 2019; Mcinnes et al., 2018Figure 5—figure supplement 2). While the identity of the previous tSNE Sub-clusters were maintained, UMAP analysis allowed us to better visualize the separation between the Sub-clusters, which we were able to broadly classify as autonomic neurons based on their shared expression of ascl1a, hand2, phox2a, and phox2bb (Figure 5—figure supplement 2A,B,D). Within the population of autonomic neurons, we were able to distinguish a population of sympathetic neurons within Sub-cluster 2 based on their combinatorial expression of th, dbh, lmo1, and insm1a (Figure 5—figure supplement 2B,E), which could be clearly distinguished from the enteric neuron Sub-cluster 3 (Figure 5—figure supplement 2B,F). Using a cluster tree, we were able to confirm the distinction between enteric neuron Sub-Cluster three and sympathetic neuron Sub-cluster 2 (Figure 5—figure supplement 2C). Taking together the architecture of the UMAP clusters and their respective gene expression signatures, Sub-cluster 0 appears as a common sympatho-enteric neuron pool that lacks the expression of sympathetic and enteric specific neurochemical markers (Figure 5—figure supplement 2A,B). Sub-clusters 1 and 4 emanate as two distinct populations from Sub-cluster 0 and respectively exhibit lower expression levels of sympathetic and enteric markers comparative to the sympathetic and enteric neuron Sub-clusters 2 and 3 (Figure 5—figure supplement 2A,B,E,F). Overall, these results suggest that Sub-cluster 0 cells may represent a pool of immature sympatho-enteric neurons, and that Sub-clusters 1 and 4 both represent further differentiated, yet still immature, pools of enteric and sympathetic neurons captured during the process of lineage specification into their respective terminal enteric and sympathetic neuron populations represented in Sub-clusters 3 and 2.

Atlas of sox10:GFP+ cell types encompassing the embryonic to larval transition

To describe the dynamic transcriptional relationship between sox10:GFP+ cells across both time points, we merged the 48–50 hpf and 68–70 hpf datasets using Seurat’s dataset and Integration and Label Transfer utility (Stuart et al., 2019). The merged datasets were visualized via UMAP, where we detected 27 clusters (Figure 6—figure supplement 1A). We observed that every cluster identified in the 48–50 hpf dataset mapped proximally to clusters at 68–70 hpf (Figure 6—figure supplement 2A). We labeled each cell in the UMAP using the previously described major cell type categories (Figure 1) forming a major cell type atlas (Figure 6—figure supplement 1B). Further refinement of the cell identities based on our previous annotations (Figure 1—source data 2) allowed us to form a higher resolution atlas for each cell type (Figure 6A). The top significantly enriched markers for each major cell type in the atlas are provided in a table in Figure 6—source data 1.

Figure 6 with 2 supplements see all
Integrated atlas of posterior sox10:GFP+ cell types spanning the embryonic to larval transition.

(A) Global UMAP embedding demonstrating the clustering of cell types across 48–50 hpf and 68–70 hpf. Cell labels were transferred from the original curation (Figure 1—source data 2) to the new atlas after its creation, allowing for assessment of cell type organization. (B) Previously identified mesenchyme clusters form a large discernible cluster marked by prrx1b, twist1a, foxc1a, and snai1a, which was separated into both chondrogenic and general mesenchyme, as denoted by its differential expression of barx1 and dlx2a. Importantly, nearly every 48–50 hpf cell type nests with a cluster at 68–70 hpf. (C) Pigment cells clusters reflect differentiation paths described in Figure 4A. Melanophores at 48–50 hpf group near to the 68–70 hpf melanophore cluster, bipotent pigment progenitors bridges both the iridophores and melanophores. Xanthophores cluster separately, reflecting their distinct lineage of origin at this developmental window. (D) Detailed analysis of the larger neural/neuronal cluster shows clear progression of cell fates from progenitor to differentiating glia or neuron. The expression of enteric neuronal markers is distinct from other subtypes at this dataset.

Figure 6—source data 1

List of marker genes per major cell type identity in the sox10:GFP+ merged atlas.

Tables reporting the Seurat output for genes for each major cell type identity including p-values (<0.01), average log-fold change (≥0.25), adjusted p-values (≤1.0), pct.1 summarizing proportion of cells expressing the individual gene in the cluster identity category, pct.2 showing the proportion of cells expressing the individual gene in all cluster identity categories in the atlas dataset.

https://cdn.elifesciences.org/articles/60005/elife-60005-fig6-data1-v1.xlsx.zip

The atlas revealed transcriptionally similar populations among posterior sox10:GFP+ cells across the embryonic to larval transition. Illustrating this point, 48h-Cluster 2 and 68h-Cluster 8 both showed a high degree of similarity, as well as consistent barx1, dlx2a, and twist1a expression (Figure 6B), consistent with our prior analysis (Figure 3C,D). Furthermore, the central node of the pigment region within the atlas was marked by 48h-Cluster 8, which resolved into respective pigment chromatophore clusters at 68–70 hpf (Figure 6C). Specifically, we observed the early specified melanophore population at 48h-Cluster 8 branched into later stage melanophore populations (68h-Clusters 4 and 18). Further, we observed that the common bi-potent pigment progenitor population (68h-Cluster 13) bridged both melanophore clusters and the iridophore 68h-Cluster 16.

Cells within the neural/neuronal clusters assembled such that progenitor cells bridged into differentiating neurons spatially from the top to the bottom of the neural/neuronal region of the atlas (Figure 6D). Six clusters (Clusters 0, 5, 7, 13, 15, 17) were represented from 48 to 50 hpf, and five clusters (Clusters 3, 5, 10, 12, 14) from 68 to 70 hpf. Of interest, the 68–70 hpf neural progenitor populations (Clusters 3 and 10) shared common gene expression with the 48–50 hpf NCC population (48h-Cluster 5), reflected largely by their co-expression of sox10, notch1a, dla, and foxd3 (Figure 6D; Figure 6—figure supplement 1D). We confirmed the spatiotemporal expression domains of notch1a and dla along the hindbrain, spinal cord, and in NCC populations along the post-otic vagal domain at 48 hpf (Figure 3—figure supplement 1E,F; arrowheads), in particular with dla in the ENCCs along the foregut (Figure 3—figure supplement 1F; arrow), a pattern similar to the ENCC makers sox10 and phox2bb (Figure 3—figure supplement 1A,B). Delineated from the neural progenitor cells, we observed a bifurcation in cell states; with one moving toward a Schwann/glial cell fate, while the other branched toward neuronal. The glial arm followed a temporal progression of earlier cell fates at 48–50 hpf (48h-Cluster 15) toward the more mature fates at 68–70 hpf (68h-Cluster 14), both denoted by expression of olig2 and pou3f1, respectively (Figure 6D). From 48h-Cluster 13, we observed the beginning of the neuronal populations, namely 48h-Clusters 0, 7, 13, and 17 and 68h-Clusters 5 and 12. The neuronal progenitor clusters (48h-Clusters 0, 13, and 17; 68h-Cluster 5) formed a spectrum of cell states leading toward the more mature neuronal populations (48h-Cluster 7; 68h-Cluster 12). For example, enteric progenitors culminated into a pool of enteric neurons, with the specific neural signature: vipb, nos1, gfra1a, fgf13b, and etv1 (Figure 6D; Figure 6—figure supplement 1D).

Closer inspection of the pigment, mesenchymal, and neural/neuronal clusters separated by time highlighted both predicted and novel changes in gene expression patterns (Figure 6—figure supplement 2). For example, while the xanthophore differentiation marker xdh demonstrated expected restriction in expression to 68–70 hpf, we identified genes with no known roles in pigment cell development differentially expressed between the two stages, such as rgs16 and SMIM18 (Figure 6—figure supplement 2B). Moreover, within the mesenchyme lineages, both barx1 and snai1b followed expected temporal expression trends, while abracl and id1 both demonstrated novel differential gene expression profiles within the mesenchyme (Figure 6—figure supplement 2C). Lastly, the neural/neuronal lineage showed expected differential gene expression of genes such as etv1 and vipb at 68–70 hpf, while revealing novel expression of nova2 and zgc:162730, which have previously uncharacterized roles in sox10-derived cells (Figure 6—figure supplement 2D). Together, these findings demonstrate dynamic gene expression changes across developmental stages occurring among sox10:GFP+ cells and highlights novel genes for further study.

A hox gene signature within sox10-derived cells in the posterior zebrafish

A common theme examined by many recent and insightful single-cell profile studies of the NCC (Dash and Trainor, 2020; Soldatov et al., 2019) is that the expression of hox genes, which encode for Homeobox transcription factors, display discrete expression patterns between various cell lineages, such as in the cranial NCC. We wondered whether specific hox signatures were expressed within posterior NCC and their recent derivatives. To analyze if we could detect hox gene patterns within the atlas, we queried all the known canonical hox genes within zebrafish as listed on zfin.org (Ruzicka et al., 2019). We detected broad expression of 45 of the 49 zebrafish hox genes across the atlas, with 85% of the cells in the atlas expressing at least one hox gene (Figure 7—figure supplement 1A,J). The four undetected hox genes (hoxc1a, hoxc12b, hoxa11a, and hoxa3a) were not examined further.

A dot plot revealed that specific hox gene expression patterns demarcated distinct tissues, with specific robustness in the neural fated cells (Figure 7A). Common to the neural lineages, we observed a core hox profile which included hoxb1b, hoxc1a, hoxb2a, hoxb3a, hoxc3a, hoxd3a, hoxa4a, hoxd4a, hoxb5a, hoxb5b, hoxc5a, hoxb6a, hoxb6b, and hoxb8a (Figure 7A). One of the top expressed constituents of the core signature, hoxa4a, was also widely expressed in several other lineages (Figure 7A; Figure 7—figure supplement 1). The hox signature applied to the NCC, neural progenitor, enteric progenitor, enteric neuron, glial progenitor, autonomic neuronal progenitor, and CNS lineages described in the atlas. Clustering of all the atlas lineages relying only on hox gene expression highlighted the robustness of the core hox signature to distinguish the neural lineage fates, grouping the differentiating (autonomic neuronal progenitors, enteric neurons, enteric progenitors and CNS neurons) and progenitor lineages (neural progenitors and glial progenitors) into neighboring clades (Figure 7B). Building on the core neural signature unifying the neural fates, slight variations in hox expression between autonomic and enteric lineages distinguished them from one another, which are summarized in Figure 7E. Most notably, considering the lineages in increasing specificity of cell fate, there was a detectable shift in hox expression among the autonomic neural progenitors to the enteric neuronal lineage, demarcated by the increase in hoxb2a, hoxd4a, hoxa5a, hoxb5a, and hoxb5b, accompanied by diminished expression of hoxc3a, hoxc5a, hoxb6a and hoxb8a, which formed a distinctive enteric hox signature (Figure 7—figure supplement 1A,K–N).

Figure 7 with 1 supplement see all
hox genes expressed across cell lineages within the sox10:GFP+ atlas.

(A) Dot plot shows both the mean expression (color) as well as percent of cells (size) per lineage for zebrafish hox genes in the first eight paralogy groups (PG). The full list of hox gene expression profiles per lineage can be found in Figure 7—figure supplement 1. Discrete hox profiles discern specific cell types, which is particularly evident in the enteric neuronal cluster. (B) Clustering of atlas lineages based hox expression profiles groups highlights robust core neural signature, which distinguishes the neural lineages from the remainder of the clades. Neural and glial progenitors formed an intermediate clade between the low-hox expressing lineages and the main neural branch. Additionally, the fin bud mesenchyme, which also has a highly distinctive hox profile, also forms a distinct clade. Subtle variations in hox expression by remaining lineages are further reflected in the remaining portion of the dendrogram; however, these are far less distinct. (C–D) Pairwise comparison of the fraction of cells in either the autonomic neural progenitor lineage (C) or the enteric neurons (D) for the first eight parology groups. Intersection of the gene pairs reflect the fraction of cells with expression for both genes with a log2 Fold change values > 0, with the identical gene intersections along the primary diagonal representing the total number of cells which express that gene in the lineage. Enteric neural hox signature was not only specific to this cell population, but also was abundantly co-expressed. (E) Summary panel describing the specific autonomic and enteric hox signatures detected. A common hox expression profile, referred to as the core signature, was found that is then modified across the specific lineages. (F–M) In situ validation of the chief enteric neural hox signature via HCR. phox2bb (F–J) labels enterically fated neurons at the level of the midgut in larval stage embryos fixed at 70–72hpf. White arrows highlight specific cells of interest. Key hox signature constituents hoxb5a (G) and hoxb5b (H) or hoxd4a (K) and hoxa5a (L) were found to be co-expressed within phox2bb expressing cells (White Arrows). Scale bars in (I,M): 50 μm.

In order to better understand the complexities of hox codes within specific lineages, we performed a pairwise comparison of each hox gene for autonomic and enteric lineages, counting the number of cells which co-expressed each hox pair. Examining the autonomic neuronal progenitors (Figure 7C) and the enteric neuron populations (Figure 7D), both lineages demonstrated pervasive fractions of co-positive cells for combinations of the core hox signature. For example, autonomic neuronal cells were enriched with a high fraction of pairwise combinations for hoxc1a, hoxa4a, hoxb3a, and/or hoxb5b(Figure 7C). The enteric signature was highly enriched in the unique expression of hoxa5a, with co-expression for hoxb5a (36%) and hoxb5b (36%), as well as strong co-expression of hoxd4a or hoxb5a with hoxb5b (Figure 7D).

To confirm that hox core genes were co-expressed within enteric neurons, we sought to validate their expression patterns using HCR probes. As previously described, hoxd4a, hoxa5a, hoxb5a, and hoxb5b all exhibited strong hindbrain expression (Figure 7—figure supplement 1B–IBarsh et al., 2017), confirming the specificity of our probes. At 70–72 hpf, as predicted by our analysis, enteric neurons along the level of the midgut, marked by phox2bb expression (Figure 7F,J), were hoxb5a+/hoxb5b+ (Figure 7G–I) and hoxd4a+/hoxa5a+ (Figure 7K–M). These data confirm that enteric neurons co-express enteric hox code genes during their early development.

With respect to the remaining cluster identities (Figure 7A), many of the populations showed varied hox expression profiles. Both the chondrogenic and general mesenchyme clusters demonstrated hoxa2b expression, as well as weak expression for hoxb2a, hoxb3a, and hoxd4a. Our detection of these hox expression profiles was consistent with prior reports that they are expressed within NCC targets toward the posterior pharyngeal arches, as well as migrating NCC (Minoux and Rijli, 2010; Parker et al., 2018; Parker et al., 2019). We detected the distinct identity of the fin bud mesenchyme (Ahn and Ho, 2008; Nakamura et al., 2016) through the expression of hoxa9b, hoxa10b, hoxa11b, hoxa13b, hoxd9a, and hoxd12a (Figure 7A; Figure 7—figure supplement 1P–Q). The pigment populations, including the pigment progenitors, melanophores, iridophores, and xanthophores, contained generally low levels of hox gene expression. Despite this, we still observed a slight variation of hox expression among the pigment populations. For example, low levels of hoxa4a, hoxb7a, hoxb8a, hoxc3a, and hoxd4a were detected among the iridophore population, while only hoxb7a was detected within a high fraction of xanthophores (Figure 7—figure supplement 1A). Interestingly, these expression profiles are not shared by the melanophore population, which displayed uniformly very low levels of detectable hox expression. Lastly, the muscle, otic, and unidentified cells showed almost no hox expression profile, which serves a foil for the specificity of the signatures outlined. We noted that the ‘pigmented muscle’ cluster weakly mirrored the general neural hox signature, likely a shared signature more reflective of the axial position of the muscle cells rather than a shared genetic profile, as corroborated by their distinct separation of the clusters on the atlas UMAP (Figure 6A).

Overall, these above described hox signatures detected within our scRNA-seq atlas indicates that distinct cell types express unique hox combinations during their delineation. Description of the hox signatures within the sox10 atlas provides further tools to identify these discrete cell populations, as well as exciting new avenues for further mechanistic investigation.

Discussion

We present a single-cell transcriptomic atlas resource capturing diversity of posterior-residing sox10-derived cells during the embryonic (48–50 hpf) to early larval transition (68–70 hpf) in zebrafish. We identified a large number of cell types; including pigment progenitor cells delineating into distinct chromatophores, NCC, glial, neural, neuronal, otic vesicle cells, muscle, as well as transcriptionally distinct mesenchymal cell populations, extending prior whole embryo-based zebrafish single cell studies (Farnsworth et al., 2020) and expanding the resolution at which these cells have been described to date. Our study is the first single-cell transcriptomic analysis covering early ENS development in zebrafish, in addition to analysis of posterior sox10+ mesenchyme and pigment cells present during the late embryonic to larval phase. The developmental window we examined, the embryonic to larval transition, is regarded as an ephemeral phase (Singleman and Holtzman, 2014) and as such is expected to contain the dynamic cell differentiation states that we observed within our atlas. We discovered that distinct hox transcriptional codes demarcate differentiating neural and neuronal populations, highlighting their potential roles during cell subtype specification. We uncovered evolutionarily conserved and novel transcriptional signatures of differentiating enteric neuron cell types, thereby expanding our knowledge of ENS development. Corroborating our transcriptomic characterizations, we validated the spatiotemporal expression of several key cell type markers using HCR. Collectively, this comprehensive cell type atlas can be used by the wider scientific community as a valuable resource for further mechanistic and evolutionary investigation of posterior sox10-expressing cells during development and the ontogenesis of neurocristopathies. The atlas is available via an interactive cell browser (https://zebrafish-neural-crest-atlas.cells.ucsc.edu/).

Collectively, our single-cell datasets captured the transition from enteric neural progenitor to differentiating enteric neuron subtype (Figures 4, 5 and 6). A similar enteric population consisting of Sox10, Ret, Phox2b, and Elavl4 was identified by scRNA-seq in the mouse (Lasrado et al., 2017), indicating zebrafish express conserved enteric programs. As well, a recent scRNA-seq study performed using E15.5 mice, a time point further along in ENS development when compared to our zebrafish study described here, suggests that Nos1+/Vip+ cells represent a post-mitotic immature neuron population capable of branching into excitatory and inhibitory neurons via subsequent differentiation mediated by lineage-restricted gene expression (Morarach et al., 2021). Their model posits that Nos1+/Vip+/Gal+ enteric neurons are capable of assuming an IPAN signature, characterized by the loss of Vip and Nos1, and the gain of Calb, Slc18a2/3, and Ntng1; a process regulated by transcription factors, Pbx3 and Etv1. This model of IPAN formation appears congruent with a previous birth dating study performed in mice, where researchers demonstrated the transient expression of Nos1 in enteric neurons (Bergner et al., 2014). Our observations in zebrafish (Figure 5; Figure 5—figure supplement 1; Figure 6) corroborates the proposed mammalian gene expression model and suggests that we captured a transitional time point where subsequent differentiation is just being initiated (Figure 5F) and suggests an evolutionarily conserved mechanism of ENS formation across vertebrate species.

Within the enteric neuron subpopulation in our dataset, we discovered enrichment of the opioid pathway genes oprl1 and oprd1b, respectively encoding for NOP and DOR class opioid receptors (Figure 5Donica et al., 2013; Holzer, 2004). The expression of opioid receptors has previously been shown in inhibitory interneurons found within the adult ENS where they are known to play functional roles during gut homeostasis (DiCello et al., 2020; Lay et al., 2016; Wood and Galligan, 2004). While the presence and inhibitory effect of opioid receptors is well characterized within the adult ENS, the role of opioid signaling during the earliest stages of enteric neuron maturation and ENS formation has yet to be investigated. Indeed, a recent study performed in zebrafish found that oprl1 was expressed within 7 dpf enteric neurons following bulk RNA sequencing (Roy-Carson et al., 2017). However, these data represent a developmental stage where zebrafish are characterized as free-swimming larvae that display feeding behavior and digestive capability representative of a functioning and more mature ENS (Cassar et al., 2018). As such, our 68–70 hpf dataset, in which we detected the expression of oprl1 and oprd1b, represents the earliest stage in which they have been shown to be expressed within the early developing ENS. The presence of opioid receptor transcripts within the ENS at 68–70 hpf, a time when immature enteric neurons are continuing to migrate and pattern within the developing embryo, highlights an important area of research focusing on the interplay between opioid use and fetal ENS development, a reality that has been on the rise in recent decades as opioid abuse continues to increase.

Finally, we have elucidated a comprehensive, combinatorial code of hox expression which defines specific cell lineages within the context of the sox10 atlas. The fin bud mesenchyme, pigment populations, and neural fates presented the most specific hox codes (Figure 7). Among the neural lineages, we identified a previously undescribed hox signature demarcating the developing enteric neuron population in zebrafish: high relative expression of hoxb2a, hoxa5a, hoxd4a, hoxb5a, and hoxb5b, while also exhibiting low expression of hoxc3a, hoxc5a, hoxb6a, and hoxb8a (Figure 7). While previous bulk microarray studies of differentiating enteric neurons from humans and mice have been shown to express orthologs to the former signature (Heanue and Pachnis, 2008; Memic et al., 2018), our analysis extends knowledge to comprehensively show for the first time specific combinatorial hox expression at a single-cell resolution, thereby providing a heretofore unknown readout of hox heterogeneity among nascent enteric cells. The conservation of an enteric hox signature between zebrafish and other systems points toward a larger conservation of function, which may facilitate translation of future findings between models. These findings imply interesting potential hypotheses wherein combinations of hox codes may represent a molecular address designating the axial site of origination for migrating NCC or indicate dynamic expression profiles which are modified during the NCC to enteric neuron developmental course. Further work is required to test these possible models as well as the functional requirement of the constituent members of the hox code in the developing zebrafish ENS. Considered collectively, these data lend support to a model in which overlapping expression domains of hox genes may facilitate enteric neural subtype differentiation, similar to their function in hindbrain and spinal neurons (Philippidou and Dasen, 2013).

In summary, our study greatly increases foundational understanding of NCC-derived cell fates, as well as other sox10+posterior cell types in zebrafish, thereby complementing ongoing studies in mammalian models and expanding fundamental knowledge of how cells diversify in developing organisms. The spatiotemporal information contained within our zebrafish atlas will serve as a resource for the developmental biology, stem cell, evolutionary biology and organogenesis communities.

Materials and methods

Key resources table
Reagent type (species) or resourceDesignationSource or referenceIdentifiersAdditional information
 Recombinant DNA reagentphox2bbUribe and Bronner, 2015
 Recombinant DNA reagentsox10Dutton et al., 2001
 Recombinant DNA reagentmmp2Mammalian Gene Collection Program Team et al., 2002
 Sequence-based reagentnotch1aUribe LabForward 5'-CAGTGGACTCAGCAGCATC-3' Reverse 5'-CCTTCCGACCAATCAGACAAG-3'
 Sequence-based reagentdlaUribe LabForward 5'-CAGCCAAGTTGCTCAGAG-3' Reverse 5'-GTACAGAGAACCAGCTCATC-3'
 Sequence-based reagentfoxc1aUribe LabForward 5'-ATACGGTGGACTCTGTGG-3' Reverse 5'-CAGCGTCTGTCAGTATCG-3'
 Genetic reagent (Danio rerio)ABZIRCWild-type zebrafish
 Genetic reagent (Danio rerio)Tg(−4.9sox10:egfp)ba2TgCarney et al., 2006GFP Labeled Neural Crest Cells
 Commercial assay, kitmitfaMolecular InstrumentsNM_130923.2
 Commercial assay, kittfecMolecular InstrumentsNM_001030105.2
 Commercial assay, kitxdhMolecular InstrumentsXM_683891.7
 Commercial assay, kitphox2bbMolecular InstrumentsNM_001014818.1
 Commercial assay, kitngfrbMolecular InstrumentsNM_001198660.1
 Commercial assay, kitgfra1aMolecular InstrumentsNM_131730.1
 Commercial assay, kitcrestinMolecular InstrumentsAF195881.1
 Commercial assay, kitfoxd3Molecular InstrumentsNM_131290.2
 Commercial assay, kitvipbMolecular InstrumentsNM_001114555.1
 Commercial assay, kitelavl3Molecular InstrumentsNM_131449
 Commercial assay, kitoprl1Molecular InstrumentsNM_205589.2
 Commercial assay, kitbarx1Molecular InstrumentsNM_001024949.1
 Commercial assay, kitpbx3bMolecular InstrumentsBC131865.1
 Commercial assay, kitprrx1bMolecular InstrumentsNM_200050.1
 Commercial assay, kitslc18a3aMolecular InstrumentsNM_0010775550.1
 Commercial assay, kithoxb5aMolecular InstrumentsNM_131101.2
 Commercial assay, kithoxb5bMolecular Instrumentsbc078285.1
 Commercial assay, kithoxd4aMolecular InstrumentsNM_001126445
 Commercial assay, kithoxa5aMolecular InstrumentsNM_131540.1
 Commercial assay, kittwist1aMolecular InstrumentsNM_130984.2
 Commercial assay, kitSingle Cell 3' v2 Chemistry Kit for 10,000 cells10x GenomicsCG00052https://support.10xgenomics.com/single-cell-gene-expression/library-prep/doc/user-guide-chromium-single-cell-3-reagent-kits-user-guide-v2-chemistry
 Commercial assay, kitNextSeq 500/550 Mid Output Kit v2.5 (150 Cycles)illumina20024904https://www.illumina.com/products/by-type/sequencing-kits/cluster-gen-sequencing-reagents/nextseq-series-kits-v2-5.html
 Chemical compound, drug1-phenyl 2-thiourea (PTU)/E3 solutionKarlsson 741 et al., 2001Sigma-AldrichP7629
 Chemical compound, drugTricaineA5040Sigma
 Chemical compound, drugAccumax bufferA7089Sigma-Aldrich
 Chemical compound, drugHank’s Buffer55021CSigma-Aldrich
 Chemical compound, drugPhusion-HFM0530SNew England Biolabs
 Chemical compound, drugZero Blunt TOPO PCR Cloning Kit451245Invitrogen
 Software, algorithmR v3.6.3R-projectRRID:SCR_001905https://www.r-project.org/
 Software, algorithmSeurat v3.1.1Satija LabRRID:SCR_007322https://github.com/satijalab/seurat
 Software, algorithmFijiPMID:22743772RRID:SCR_002285https://imagej.net/Fiji
 Software, algorithmIMARIS v9.2BitplaneRRID:SCR_007370Bitplane.com
 Software, algorithmPANTHERGENEONTOLOGY Unifying BiologyRRID:SCR_004869http://pantherdb.org
 Software, algorithmCell Ranger v2.1.010x GenomicsRRID:SCR_017344Zheng et al., 2017, https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/what-is-cell-ranger

Animal husbandry, care, and synchronous embryo collection

Request a detailed protocol

Groups of at least 15 adult Tg(−4.9sox10:GFP)ba2Tg (Carney et al., 2006) zebrafish (Danio rerio) males and 15 females from different tank stocks were bred to generate synchronously staged embryos across several clutches. All embryos were cultured in standard E3 media until 24 hr post fertilization (hpf), then transferred to 0.003% 1-phenyl 2-thiourea (PTU)/E3 solution (Karlsson et al., 2001), to arrest melanin formation and enable ease of GFP sorting. While it has been suggested that high concentration (.03%) PTU incubation prior to 22 hpf may cause organism-wide effects,. 003% PTU application after 22 hpf has been shown to have no major effects on NCC survival or pigment cell formation, as described (Bohnsack et al., 2011). Embryos were manually sorted for GFP expression and synchronously staged at 24 hpf. Care was taken such that embryos which exhibited developmental delay or other defects were removed prior to collection. All work was performed under protocols approved by, and in accordance with, the Rice University Institutional Animal Care and Use Committee (IACUC).

Isolation of tissue and preparation of single-cell suspension

Request a detailed protocol

100 embryos between 48- 50 hpf and 100 larvae between 68–70 hpf were dechorionated manually and then transferred to 1X sterile filtered PBS, supplemented with 0.4% Tricaine (Sigma, A5040) to anesthetize. Tissue anterior to the otic vesicle and tissue immediately posterior to the anal vent was manually removed using fine forceps in 48–50 hpf embryos, while tissue anterior to the otic vesicle was removed from 68 to 70 hpf larvae, as schematized in Figure 1. This was to capture as many posterior sox10:GFP+ cells in the later time point as possible. Remaining tissue segments were separated into nuclease-free tubes and kept on ice immediately following dissection. Dissections proceeded over the course of 1 hr. To serve as control for subsequent steps, similarly staged AB WT embryos were euthanized in tricaine and then transferred to sterile 1X PBS. All following steps were conducted rapidly in parallel to minimize damage to cells: Excess PBS was removed and tissue was digested in 30°C 1X Accumax buffer (Sigma-Aldrich, A7089) for 30–45 min to generate a single cell suspension for each sample. At 10 min intervals, tissue was gently manually disrupted with a sterile pipette tip. As soon as the tissue was fully suspended, the cell solutions were then transferred to a fresh chilled sterile conical tube and diluted 1:5 in ice cold Hank’s Buffer (1x HBSS; 2.5 mg/mL BSA; 10 μM pH8 HEPES) to arrest the digestion. Cells were concentrated by centrifugation at 200 rcf for 10 min at 4°C. Supernatant was discarded carefully and cell pellets were resuspended in Hank’s Buffer. Cell solution was passed through a 40 μm sterile cell strainer to remove any remaining undigested tissue and then centrifuged as above. Concentrated cells were resuspended in ice cold sterile 1X PBS and transferred to a tube suitable for FACS kept on ice. The 48–50 hpf and 68–70 hpf experiments were performed on completely separate dates and times using the above described procedures.

Fluorescent cell sorting and single-cell sequencing

Request a detailed protocol

Fluorescent Assisted Cell Sorting (FACS) was performed under the guidance of the Cytometry and Cell Sorting Core at Baylor College of Medicine (Houston, TX) using a BD FACSAria II (BD Biosciences). Zebrafish cells sorted via GFP fluorescence excited by a 488 nm laser, relying on an 85 μm nozzle for cell selection. Detection of GFP+ cells was calibrated against GFP- cells collected from AB wildtype embryos, as well as GFP+ cells collected from the anterior portions of the sox10:GFP embryos. Optimal conditions for dissociated tissue inputs (number of embryos needed, etc.) and FACS gating was determined via pilot experiments prior to collection for subsequent scRNA-seq experiments. Sample preparation for scRNA-seq was performed by Advanced Technology Genomics Core (ATGC) at MD Anderson (Houston, TX). 4905 and 4669 FACS-isolated cells for the 48–50 and 68–70 hpf datasets were prepared on a 10X Genomics Chromium platform using 10X Single Cell 3’ V2 chemistry kit for 10,000 cells. cDNA libraries were amplified and prepared according to the 10X Genomics recommended protocol, with details provided in Figure 1—figure supplement 1C. A 150 cycle Mid-Output flow cell was used for sequencing on an Illumina NextSeq500. Sequencing was aligned at MD Anderson ATGC to the DanioGRCz10 version of the zebrafish genome using the 10X Genomics Cell Ranger software (v2.1.0, Zheng et al., 2017). Gene reads per cell were stored in a matrix format for further analysis.

Data processing and analysis

Request a detailed protocol

The 10x genomics sequencing data was then analyzed using Seurat (Satija et al., 2015; Stuart et al., 2019; Butler et al., 2018) v3.1.1 software package for R, v3.6.3 (R Development Core Team, 2020). The standard recommended workflow was followed for data processing. Briefly, for both the 48–50 hpf and 68–70 hpf datasets, cells which contained low (<200) or high (>2500) genes were removed from analysis. Gene expression was normalized using the NormilizeData command, opting for the LogNormalize method (Scale factor set at 10,000) and further centered using the ScaleData command. Variable features of the dataset were calculated with respect to groups of 2000 genes at a time. Both datasets were evaluated considering the first 20 principal components (PC) as determined by the RunPCA command with a resolution of 1.2 for PCA, tSNE, and UMAP analyses. The resolution value was empirically determined for the FindClusters by iterative assessment of cluster patterns generated values ranging from 0.4 to 2.0. The resolution was set to 1.2 which balanced between over segmentation of cell identities and effective separation of tSNE-mapped clusters. The appropriate PCs were selected based on a Jack Straw analysis with a significance of p<0.01, as generated by the JackStraw command.

Clustering was performed using FindNeighbors and FindClusters in series. We identified 19 clusters in the 48–50 hpf dataset and 23 clusters in the 68–70 hpf dataset. Significantly enriched genes for each cluster were determined via a Wilcoxon Rank Sum test implemented by the FindAllMarkers command. From these expressed gene lists within each cluster, all cluster identities were manually curated via combinatorial expression analysis of published marker genes in the literature, zfin.org and/or bioinformatics GO term analysis via the Panther Database.

Generation of the merged atlas was performed via the FindIntegrationAnchors workflow provided in the Standard Workflow found on the Seurat Integration and Label Transfer vignette. Clustering was performed for the atlas based on the first 20 PCs, consistent with the original datasets. Subsets of the atlas discounted any spuriously sorted cells for clarity. All features plots represent expression values derived from the RNA assay. Sub-clustering of the enteric clusters was performed by sub-setting Clusters 5 and 12 from the 68 to 70 hpf dataset and reinitializing the Seurat workflow, as described above. Clusters were identified based on the first 6 PCs. Detection of cell cycle phase was conducted following the Cell cycle and scoring vignette. Genes used for identification of cell cycle phases can be found in Figure 1—figure supplement 3. Pairwise hox analysis was conducted using tools from the Seurat package in R by assessing the number of cells which had expression for hox gene pairs queried with a log2 fold change greater than 0. Dendrogram Cluster trees were generated using Seurat’s BuildClusterTree function.

Whole mount in situ hybridization

Request a detailed protocol

cDNAs for foxc1a, notch1a, and dla were amplified via high fidelity Phusion-HF PCR (NEB) from 48 hpf AB WT cDNA libraries using primers in the Key resources table. PCR products were cloned using the Zero Blunt TOPO PCR Cloning Kit (Invitrogen), as per manufacturer protocols, and sequenced validated. Plasmids encoding phox2bb, sox10, mmp2 were generously sourced as listed in the key resources table. Antisense digoxigenin (DIG)-labeled riboprobes were produced from cDNA templates of each gene. AB wild type embryos were treated and stained to visualize expression as previously described in Jowett and Lettice, 1994. Following in situ reactions, embryos were post-fixed in 4% Paraformaldehyde (PFA) and mounted in 75% Glycerol for imaging. A Nikon Ni-Eclipse Motorized Fluorescent upright compound microscope with a 4X objective was used in combination with a DS-Fi3 color camera. Images were exported via Nikon Elements Image Analysis software.

Whole mount hybridization chain reaction

Request a detailed protocol

HCR probes were purchased commercially (Molecular Instruments Inc, CA) and were targeted to specific genes based on their RefSeq ID (Key resources table). Whole mount HCR was performed according to the manufacturer's instructions (v3.0, Choi et al., 2016; Choi et al., 2018) on sox10:GFP+ or AB embryos previously fixed at the appropriate stage in 4% PFA.

Confocal imaging and image processing

Request a detailed protocol

Prior to imaging, embryos were embedded in 1% low melt agarose (Sigma) and were then imaged using an Olympus FV3000 Laser Scanning Confocal, with a UCPlanFLN 20×/0.70NA objective. Confocal images were acquired using lambda scanning to separate the Alexafluor 488/Alexafluor 514 or the Alexafluor 546/Alexafluor 594 channels. Final images were combined in the FlowView software and exported for analysis in either Fiji (Rueden et al., 2017; Schneider et al., 2012; Schindelin et al., 2012) or IMARIS image analysis software (Bitplane). Figures were prepared in Adobe Photoshop and Illustrator software programs, with some cartoons created via https://biorender.com/.

Data availability

Request a detailed protocol

The raw sequence read files and processed cellular barcode, gene, and matrix files produced by CellRanger are available in the National Center for Biotechnology Information’s (NCBI) Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/), accession number: GSE152906. The atlas/associated processed Seurat objects are available on the University of California, Santa Cruz (UCSC) Cell Browser (https://zebrafish-neural-crest-atlas.cells.ucsc.edu/). Code written for Seurat data analysis is available on GitHub (https://github.com/UribeLabRice/NeuralCrest_Atlas_2020; UribeLabRice, 2021; copy archived at swh:1:rev:195fe7dd0bf388ad4cd6f01db052f45d8b8fa67e).

Data availability

Sequencing data have been added to GEO under accession GSE152906. The raw sequence read files and processed cellular barcode, gene, and matrix files produced by CellRanger are available in the National Center for Biotechnology Information's (NCBI) Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/), accession number: GSE152906. The atlas/associated processed Seurat objects are available on the University of California, Santa Cruz (UCSC) Cell Browser (https://zebrafish-neural-crest-atlas.cells.ucsc.edu/). Code written for Seurat data analysis is available on GitHub (https://github.com/UribeLabRice/NeuralCrest_Atlas_2020; copy archived at https://archive.softwareheritage.org/swh:1:rev:195fe7dd0bf388ad4cd6f01db052f45d8b8fa67e/).

The following data sets were generated
    1. Howard AA
    2. Baker PA
    3. Ibarra-García-Padilla R
    4. Moore JA
    5. Rivas LJ
    6. Tallman JJ
    7. Corteguera JA
    8. Westheimer JL
    9. Singleton EW
    10. Uribe RA
    (2020) NCBI Gene Expression Omnibus
    ID GSE152906. An atlas of neural crest lineages along the posterior developing zebrafish at single-cell resolution.

References

  1. Book
    1. Barlow AJ
    (1984)
    Neural Crest Cells in Enteric Nervous System Development and Disease
    In: Barlow A. J, editors. In Neural Crest Cells. Elsevier Inc. pp. 101–104.
    1. Delalande JM
    2. Guyote ME
    3. Smith CM
    4. Shepherd IT
    (2008) Zebrafish sip1a and sip1b are essential for normal axial and neural patterning
    Developmental Dynamics: An Official Publication of the American Association of Anatomists 237:1060–1069.
    https://doi.org/10.1002/dvdy.21485
    1. Dutton KA
    2. Pauliny A
    3. Lopes SS
    4. Elworthy S
    5. Carney TJ
    6. Rauch J
    7. Geisler R
    8. Haffter P
    9. Kelsh RN
    (2001)
    Zebrafish colourless encodes sox10 and specifies Non-Ectomesenchymal neural crest fates
    Development 128:4113–4125.
  2. Book
    1. Hall BK
    2. Hörstadius S
    (1988)
    The Neural Crest
    Oxford University Press.
    1. Kelsh RN
    2. Eisen JS
    (2000)
    The zebrafish colourless gene regulates development of non-ectomesenchymal neural crest derivatives
    Development 127:515–525.
  3. Book
    1. Le Douarin N
    2. Kalcheim C
    (1999)
    The Neural Crest
    Cambridge University Press.
    1. Le Lièvre CS
    2. Le Douarin NM
    (1975)
    Mesenchymal derivatives of the neural crest: analysis of chimaeric quail and chick embryos
    Journal of Embryology and Experimental Morphology 34:125–154.
    1. Lister JA
    2. Robertson CP
    3. Lepage T
    4. Johnson SL
    5. Raible DW
    (1999)
    Nacre encodes a zebrafish Microphthalmia-Related protein that regulates Neural-Crest-Derived pigment cell fate
    Development 126:3757–3767.
    1. Lu J-K
    2. Tsai T-C
    3. Lee H
    4. Hsia K
    5. Lin C-H
    6. Lu J-H
    (2019)
    Pectoral fin anomalies in tbx5a knockdown zebrafish embryos related to the cascade effect of N-Cadherin and extracellular matrix formation
    Journal of Developmental Biology 7:15.
    1. Parichy DM
    2. Ransom DG
    3. Paw B
    4. Zon LI
    5. Johnson SL
    (2000)
    An orthologue of the Kit-Related gene fms is required for development of neural Crest-Derived xanthophores and a subpopulation of adult melanocytes in the zebrafish, Danio rerio
    Development 127:3031–3044.
  4. Software
    1. R Development Core Team
    (2020) R: A Language and Environment for Statistical Computing
    R Foundation for Statistical Computing, Vienna, Austria.
  5. Data
    1. Thisse B
    2. Thisse C
    (authors) (2004) Fast release clones: a high throughput expression analysis
    ZFIN Direct Data Submission. ID ZDB-PUB-040907-1.
  6. Data
    1. Thisse C
    2. Thisse B
    (authors) (2005) High throughput expression analysis of ZF-Models consortium clones
    ZFIN Direct Data Submission. ID ZDB-PUB-051025-1.
    1. Yelon D
    2. Ticho B
    3. Halpern ME
    4. Ruvinsky I
    5. Ho RK
    6. Silver LM
    7. Stainier DY
    (2000)
    The bHLH transcription factor Hand2 plays parallel roles in zebrafish heart and pectoral fin development
    Development 127:2573–2582.
  7. Book
    1. Zoli M
    2. Berlin H
    (2000) Distribution of Cholinergic Neurons in the Mammalian Brain with Special Reference to their Relationship with Neuronal Nicotinic Acetylcholine Receptors
    In: Clementi F, Fornasari D, Gotti C, editors. Neuronal Nicotinic Receptors. Handbook of Experimental Pharmacology. Springer. pp. 13–30.
    https://doi.org/10.1007/978-3-642-57079-7_2

Decision letter

  1. Richard M White
    Senior Editor; Memorial Sloan Kettering Cancer Center, United States
  2. Lilianna Solnica-Krezel
    Reviewing Editor; Washington University School of Medicine, United States
  3. Cecilia Lanny Winata
    Reviewer
  4. Kristin Artinger
    Reviewer

In the interests of transparency, eLife publishes the most substantive revision requests and the accompanying author responses.

Acceptance summary:

The manuscript reports comprehensive analysis of single cell sequencing of FACS sorted neural crest cells from the posterior head and trunk of the zebrafish embryo at two developmental time points, generating a ~4,000 cell scRNAseq dataset. By performing standard analysis using Seurat to cluster and annotate the data, the authors identified major cell types within the neural crest at these stages, including pigment cells, craniofacial and enteric neuronal populations, and unique combinatorial hox gene expression among particular neural crest cell types. The results presented here could help other investigators identify populations of neural crest cells as well as how the gene regulatory network functions at later time points in development.

Decision letter after peer review:

Thank you for submitting your article "An atlas of neural crest lineages along the posterior developing zebrafish at single-cell resolution" for consideration by eLife. Your article has been reviewed by three peer reviewers, and the evaluation has been overseen by a Reviewing Editor and Richard White as the Senior Editor. The following individuals involved in review of your submission have agreed to reveal their identity: Cecilia Lanny Winata (Reviewer #2); Kristin Artinger (Reviewer #3).

The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.

As the editors have judged that your manuscript is of interest, but as described below that additional experiments are required before it is published, we would like to draw your attention to changes in our revision policy that we have made in response to COVID-19 (https://elifesciences.org/articles/57162). First, because many researchers have temporarily lost access to the labs, we will give authors as much time as they need to submit revised manuscripts. We are also offering, if you choose, to post the manuscript to bioRxiv (if it is not already there) along with this decision letter and a formal designation that the manuscript is "in revision at eLife". Please let us know if you would like to pursue this option. (If your work is more suitable for medRxiv, you will need to post the preprint yourself, as the mechanisms for us to do so are still in development.)

Summary:

The manuscript reports comprehensive analysis of single cell sequencing of FACS sorted neural crest cells from the posterior head and trunk of the zebrafish embryo at two developmental time points, generating a ~4,000 cell scRNAseq dataset. By performing standard analysis using Seurat to cluster and annotate the data, the authors identified major cell types within the neural crest at these stages, including pigment cells, craniofacial and enteric neuronal populations, and unique combinatorial hox gene expression among particular neural crest cell types. The results presented here could help other investigators identify populations of neural crest cells as well as how the gene regulatory network functions at later time points in development.

Essential revisions:

There was consensus among the reviewers, that while this work presents a useful data set, the current manuscript falls short on conclusions about neural crest drawn from the data. The manuscript discusses at length about gene expression in various NCC populations without much functional analysis or offering new insights. To become suitable for publication in eLife the manuscript would need to be revised to shorten methodological aspects, present new findings and highlight their significance. Based on the data presented, enteric neurons and hox code present such opportunities. For enteric neurons, the current manuscript has missed the opportunity to use the scRNA-seq data in an unbiased way, the authors searched for known signatures. It will be important to demonstrate that this data set can reveal previously undescribed genes and pathways. It is understood that in the current situation of ongoing pandemic and reduced laboratory efforts, functional analyses are not feasible, although validation of expression of some novel genes, would significantly strengthen the manuscript.

1) The Introduction would be improved by a discussion about the known gene regulatory network at these stages and cell types. It will help put the study in context of what is known in the field. This would include expression analysis by RNA in situ and bulk RNA-seq analysis on this population.

2) The use of PTU does not seem to be necessary for the purpose of the experiment as the GFP signal should still be visible even with pigmentation. The reviewers therefore wondered why the authors incorporated this procedure. Although PTU is technically not supposed to affect the early steps of melanophore differentiation, possible implications of this need to be at least clarified as there is evidence that PTU can affect the organism at molecular and physiological levels.

3) There should be a discussion of how the cell type classification markers where chosen for the clustering. For those not in the field, it is not clear. Please add.

4) The Results section should be shortened and condensed significantly. Over half of the Results section (some 12 to 14 pages) is merely describing the authors' process of annotating the data. This could easily be summarized in a full-page figure with an annotated tSNE plot and a dot plot indicating the markers that were used to define the annotations. It would, in fact, be much clearer for readers.

5) Additionally, many other sections of the Results are so verbose that they obscure the authors' meaning. Tightening the prose would significantly improve the manuscript. Paragraph three of subsection “Pigment cell chromatophore lineages resolved” is such an example. Could the authors not just say "Previous work identified melanophores in distinct proliferative and differentiating states at 5 dpf (Saunders et al., 2019). We find these states arise between 50 and 68 hpf, as all melanophores were in the G1 phase at 48-50 hpf, but formed two distinct clusters at 68-70 hpf differentiated by their cell cycle state."

6) The HCR expression analyses are a nice addition to show the expression of some of the genes that are expressed in specific populations. However, the genes shown are already known and thus higher resolution images (For craniofacial figure, pigment I and J are hard to see, enteric figure is good but also would benefit from sectional analysis) to show neural crest cells at the cellular level and potential overall of expression would move the field and be more similar to the single cell sequencing. It would also be helpful for data interpretation. That being said, it would also be interesting to show some of the genes in the clusters that are novel or that have not been shown to be expressed in that cell type as well.

7) "Hierarchical clustering of cells with General Mesenchyme and Chondrogenic identities using a cluster tree": This is not technically correct, assuming the authors used BuildClusterTree (which it seems so from the Materials and methods). This function does not operate on the level of cells, but hierarchically clusters the mean expression signature of the already-calculated clusters.

8) It's not clear why the authors present the data twice, first as each individual time point and then as a combined analysis, when the same clusters are recovered in the combined analysis. It seems as if it would be easier and more informative to just present the combined data.

9) The authors provide an estimated number of cells reported by 10x Cell Ranger pipeline which corresponds to 2300 and 2580 for 48-50hpf and 68-70hpf, respectively. We understood that these numbers correspond to the estimated number of cells that 10x pipeline was able to identify which in the following steps serves as an input for Seurat analysis in R. The total number of usable cells after all QC steps were 1608 (58-60hpf) and 2410 (68-70hpf). However, there was no information on how many cells were loaded on to the 10x chip. This information would help others, especially in the zebrafish community, to better plan single-cell experiments and would provide additional information about cell suspension quality.

10) The reviewers were surprised that there is a large population of neural crest classified cells at these stages. They wondered if those are progenitor cells and do they remain in adults? In addition, it seems like most of the neuronal population in enteric, where are the dorsal root ganglia and sympathetic neuronal populations?

11) In the two stages of neural crest development, it would be interesting to include how the gene expression changes across developmental time. Could you present a combined UMAP with stages colored to see where stages fall and sets of genes that change expression through developmental time in subclusters? Highlighting genes that change over developmental time in each cluster would be very informative.

12) Discussion penultimate paragraph: This is potentially an interesting conclusion that the authors could spend more time on in the paper. Is it worth generating some version of a figure that combines a review of what is known about these signatures from other systems (and which other systems) with what the authors observe in zebrafish? One could envision some form of dot plot with an organismal color code, but obviously the authors may have more creative ideas.

13) (783/784) Mitochondrial genes were regressed out from the data set during the Cell Ranger alignment. Mitochondrial genes, together with ribosomal genes, are useful metrics for the assessment of cell quality. High expression levels of mitochondrial genes may indicate poor sample quality, apoptotic cells, or multiplets. Or if it is limited to a few clusters, it may reveal the nature of some cells (e.g. increased metabolic activity). However, the authors decided to remove mitochondrial genes from the alignment. Why?

14) The paper makes no mention of how the data and code will be made available. It would be standard practice to deposit the FASTQ files as well as the counts table output by CellRanger in GEO. Moreover, it's fairly standard practice (and generally appreciated) to make the processed Seurat objects available (either through lab website, Data Dryad, or some other hosting platform). Additionally, it's generally accepted practice to publish the code that was used to analyze the data and generate the figures.

[Editors' note: further revisions were suggested prior to acceptance, as described below.]

Thank you for submitting your article "An atlas of neural crest lineages along the posterior developing zebrafish at single-cell resolution" for consideration by eLife. Your article has been reviewed by three peer reviewers, and the evaluation has been overseen by a Reviewing Editor and Richard White as the Senior Editor. The following individuals involved in review of your submission have agreed to reveal their identity: Cecilia Lanny Winata (Reviewer #2); Kristin Artinger (Reviewer #3).

The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.

We would like to draw your attention to changes in our revision policy that we have made in response to COVID-19 (https://elifesciences.org/articles/57162). Specifically, we are asking editors to accept without delay manuscripts, like yours, that they judge can stand as eLife papers without additional data, even if they feel that they would make the manuscript stronger. Thus the revisions requested below only address clarity and presentation.

Summary:

The manuscript reports comprehensive analysis of single cell sequencing of FACS sorted neural crest cells from the posterior head and trunk of the zebrafish embryo at two developmental time points, generating a ~4,000 cell scRNAseq dataset. By performing standard analysis using Seurat to cluster and annotate the data, the authors identified major cell types within the neural crest at these stages, including pigment cells, craniofacial and enteric neuronal populations, and unique combinatorial hox gene expression among particular neural crest cell types. The results presented here could help other investigators identify populations of neural crest cells as well as how the gene regulatory network functions at later time points in development.

Revisions:

All reviewers thought that the authors have made significant improvement to the manuscript by removing many of the technical descriptions from their results. However they also thought that more can be done to make the Results section more readable. Many unnecessary details and repetitions remain, which obscures the communication of essential findings. Moreover, the manuscript is perceived as written for a neural crest audience and thus less accessible to a more general (eLife) audience making the authors motivation and essential findings obscured. The specific comments and suggestions of the reviewers that can guide additional revisions are appended below.

Reviewer #1:

My only remaining comment that should not block publication, but perhaps would inspire some additional revisions from the authors, is that overall, I still find that the prose is disappointingly focused on technical aspects of a pretty standard single-cell analysis at the expense of describing discoveries from those analyses. Many sections of the paper end with the main conclusion that the section "demonstrates the power of the sox10 atlas to.…" but I don't know why this is repeatedly considered a result, given that standard and well demonstrated analysis methods were used on data generated using now standard approaches. The new and expanded Discussion section illustrates that the authors have found some interesting results, but I feel that they dilute them significantly by focusing so much on telling us about the power of now-standard single-cell genomics approaches and about the process of annotating the data, when they could instead place the focus on their more interesting findings. For instance, I find this is especially true in the presentation of the mesenchyme - which seems mostly focused on describing existing mesenchyme markers and how they were used to identify the mesenchyme and then the number of clusters found therein, but without significant description of what those clusters are or what defines them.

Reviewer #2:

Despite the new additions, the manuscript remains highly descriptive and a major issue in terms of writing style remains.

1) The authors have made significant improvement by removing many of the technical descriptions from their results. However this did little to make the Results section more readable. Many unnecessary details and repetitions remain, which obscures the communication of essential findings. The language used is often colloquial, with some grammatical errors which makes reading extremely challenging, especially for a non-expert in NCC or single cell biology. Having the manuscript read by a professional scientific editing service might be of help. To cite some examples:

– "We have utilized the Tg(-4.9sox10:EGFP) (hereafter referred to as 117 sox10:GFP) transgenic fish line to identify". Authors were requested to use the term "zebrafish" instead of "fish" since other fishes are used for these types of studies and it is confusing. Please correct this.

– What does "remarkably captured" mean?

– "Taken together, these results show that the scRNA-seq datasets effectively identify discrete subpopulations, and coupled with our HCR analysis, effectively shows we are able to validate these cell populations in vivo". Please clarify which subpopulations are referred to here.

– It is unclear what the authors meant by "early neuronal differentiation". Please specify the exact or range of developmental stage(s) for clarity.

– "Further investigation into these pathways led to the identification of…" – would it be sufficient to simply state for e.g. that "oprl1 and oprd1b, which are members of this pathway, were expressed in subcluster 3"?

– ". patterns of Homeobox transcription factors, known as hox genes…" – how can transcription factors be known as genes? Please rephrase this.

2) The Discussion section is also extremely lengthy, containing many repetitions of the results which could be significantly condensed. There are also a lot of anecdotal information written on opioids – which seems irrelevant and unnecessary.

3) In describing the mesenchymal clusters, the authors used the term "subtypes". It seems more appropriate to call them simply as "clusters" as there were no systematic analyses done (e.g. expression of specific cell type markers) to define whether each of these barx1+ and dlx2a+ clusters represent distinct cell subtypes. Moreover, the number of clusters obtained depends on the threshold set in the scRNA-seq data analysis, which is therefore arbitrary in nature.

Along the same lines, the resolution parameter from Seurat FindClusters function determines the number of obtained clusters. Naturally, the selection of this parameter is arbitrary and needs to be individually tailored for each sample. Looking at the data at different resolutions may help to improve the identification of novel cell types. Is there any particular reason why the authors decided to apply the resolution of 1.2? A clarification on how this threshold was chosen would be helpful.

4) "… pbx3b expression may promote the assumption of an IPAN signature characterized by the presence of calb2a, ache, and slc18a3a and the loss of inhibitor markers nos1 and vipb." – it is unclear how the authors' observations led to this hypothesis as earlier on they stated that "…pbx3b expression was found in combination with calb2a, vipb, and nos1…". Please clarify.

5) "Overall, these results suggest that sub-cluster 0 cells may represent a pool of immature sympatho-enteric neurons, and that sub-clusters 1 and 4 both represent better resolved, yet still immature, pools of enteric and sympathetic neurons." – it is unclear how the authors arrived at this notion, especially since cluster 0, 1, and 4 was not even mentioned at all previously.

6) It would be helpful if the authors also annotate the clusters in Figure 6—figure supplement 2 panel A and Figure 7—figure supplement 1 panel J.

7) "We performed cell filtering and clustering (.) cells which contained low (<200) or high (>2500) genes were removed from analysis".

– As indicated in Figure 1—figure supplement 1 panel C, around 4905 and 4669 cells were loaded for encapsulation into 10x cartridge, resulting in 2300 and 2580 identified cells. After a filtering step that was done on the basis of gene content only, 1608 and 2410 cells were obtained. The expected multiplet rate at such cell suspension concentration should be low (around 2% according to 10x protocol), therefore the great loss of cell numbers is rather surprising. Have the authors checked why the cell loss at the applied threshold level for the 48-50 hpf sample is much higher than expected. Also, did they check the distribution of expressed genes before setting up the thresholds for gene content per cell?

– Additionally, including mitochondrial gene content into the filtering step may help to confidently identify low quality cells.

8) Figure 1—figure supplement 1 panel C – inconsistent number format (comma, dot, space), e.g. mean read per cell (74,017 vs 62109 etc.)

9) “…Major cell type categories were based on the presence of signature marker genes…”. Is the word "presence" an appropriate term in this context? Presence is rather a binary value and takes two options: 1- present, 0 – not present. It seems more accurate to use the term "expression" of these signatures/marker genes.

https://doi.org/10.7554/eLife.60005.sa1

Author response

Essential revisions:

There was consensus among the reviewers, that while this work presents a useful data set, the current manuscript falls short on conclusions about neural crest drawn from the data. The manuscript discusses at length about gene expression in various NCC populations without much functional analysis or offering new insights. To become suitable for publication in eLife the manuscript would need to be revised to shorten methodological aspects, present new findings and highlight their significance. Based on the data presented, enteric neurons and hox code present such opportunities. For enteric neurons, the current manuscript has missed the opportunity to use the scRNA-seq data in an unbiased way, the authors searched for known signatures. It will be important to demonstrate that this data set can reveal previously undescribed genes and pathways. It is understood that in the current situation of ongoing pandemic and reduced laboratory efforts, functional analyses are not feasible, although validation of expression of some novel genes, would significantly strengthen the manuscript.

We have revised the manuscript to balance methodological aspects, while highlighting new findings made possible by leveraging analyses of our single cell RNA-seq datasets. We agree with reviewers that enteric and Hox gene expression patterns are opportunities to focus on novel findings; we have thus expanded analyses of these signatures in the revised manuscript.

First, we present analysis of incipient enteric neuron populations, where we have discovered the presence of cells along a spectrum of differentiation states by using sub-clustering analysis. From this analysis (Figure 5), we found that enteric populations in zebrafish exist in several sub-populations—enteric neuronal progenitors, immature enteric neurons and enteric cells differentiating towards a specific sensory subtype. Leveraging data from the enteric subpopulations we discovered, we performed unbiased pathway analysis and, among many new pathways, detected the presence of opioid pathway receptors within differentiating enteric neurons—a novel finding for the developmental ENS field. Importantly, we validated in situ expression of the opioid receptor genes within enteric neurons along the gut in vivo using HCR. We highlight and discuss the significance of opioid receptor detection within nascent enteric neurons within our revised discussion.

Secondly, we expanded analysis of Hox gene codes within the sox10 atlas presented in the manuscript. By calculating the pairwise, co-expression values of all hox genes expressed within the atlas, we present a data matrix that summarizes combinatorial hox expression within enteric neurons and autonomic neurons (Figure 7). We found for the first time in zebrafish a specific enteric neuron hox code, marked by enrichment of hoxb2a, hoxa5a, hoxd4a, hoxb5a and hoxb5b. We validated in situ HCR expression of four of these hox factors within phox2bb+ enteric cells along the zebrafish midgut in vivo. The pairwise co-expression analysis indicated that enteric neurons present with varying combinatorial states for the enteric code during their differentiation, revealing novel details regarding heterogeneity of early enteric transcriptional states.

Finally, we have demonstrated that the sox10 atlas can be used to showcase novel genes within subpopulations across time. We highlighted known and novel genes for the mesenchymal, pigment and neural populations in Figure 6—figure supplement 2. Extending beyond this and bringing the power of the atlas to the scientific community, we have deposited our merged Seurat atlas dataset to the UCSC Cell Browser, which enables interactive assessment of annotated cell populations and interrogation of genes of interest.

1) The Introduction would be improved by a discussion about the known gene regulatory network at these stages and cell types. It will help put the study in context of what is known in the field. This would include expression analysis by RNA in situ and bulk RNA-seq analysis on this population.

We have added in some examples of recent studies addressing gene regulatory networks of neural crest cells and their early forming derivatives. As highlighted in the Introduction, most of what is known regarding neural crest development using transcriptomics and chromatin profiling studies centers largely on early neural crest specification in pre-migratory and early migratory stages. Filling a gap in knowledge in the field, our study focuses on later phases of neural crest migration and lineage subtype differentiation.

2) The use of PTU does not seem to be necessary for the purpose of the experiment as the GFP signal should still be visible even with pigmentation. The reviewers therefore wondered why the authors incorporated this procedure. Although PTU is technically not supposed to affect the early steps of melanophore differentiation, possible implications of this need to be at least clarified as there is evidence that PTU can affect the organism at molecular and physiological levels.

We have added a sentence in the Materials and methods section to clarify the use of PTU in our studies. We added low dose PTU to embryos after the first day in development to prevent pigmentation formation that enabled ease of GFP sorting prior to dissections. While high dose PTU prior to the first day in development has been shown to adversely affect embryonic development, we have acknowledged this in the Materials and methods and that low dose PTU is not believed to cause defects in neural crest development at the stages we examined in our study.

3) There should be a discussion of how the cell type classification markers where chosen for the clustering. For those not in the field, it is not clear. Please add.

We have added information regarding how clusters were annotated with cell type classifications into the Results and the Materials and methods. We would also like to clarify that Seurat clustering of the single cell RNA-seq transcriptomes was performed based on PC analysis prior to any cell type classifications, standard for the Seurat pipeline, as we describe in the Materials and methods.

4) The Results section should be shortened and condensed significantly. Over half of the Results section (some 12 to 14 pages) is merely describing the authors' process of annotating the data. This could easily be summarized in a full-page figure with an annotated tSNE plot and a dot plot indicating the markers that were used to define the annotations. It would, in fact, be much clearer for readers.

We agree and have shortened and condensed the first half of the Results section. Along with that, we have also summarized/condensed the first original three figures into one major figure with one supplement. New Figure 1 contains annotated tSNE and dot plots summarizing major cell type categories of sox10+ cells for each major time point, as suggested by reviewers.

5) Additionally, many other sections of the results are so verbose that they obscure the authors' meaning. Tightening the prose would significantly improve the manuscript. Paragraph three of subsection “Pigment cell chromatophore lineages resolved” is such an example. Could the authors not just say "Previous work identified melanophores in distinct proliferative and differentiating states at 5 dpf (Saunders et al., 2019). We find these states arise between 50 and 68 hpf, as all melanophores were in the G1 phase at 48-50 hpf, but formed two distinct clusters at 68-70 hpf differentiated by their cell cycle state."

Along with comment 4 above, we have also edited many sections of the results to tighten prose, including the section regarding pigment cells.

6) The HCR expression analyses are a nice addition to show the expression of some of the genes that are expressed in specific populations. However, the genes shown are already known and thus higher resolution images (For craniofacial figure, pigment I and J are hard to see, enteric figure is good but also would benefit from sectional analysis) to show neural crest cells at the cellular level and potential overall of expression would move the field and be more similar to the single cell sequencing. It would also be helpful for data interpretation. That being said, it would also be interesting to show some of the genes in the clusters that are novel or that have not been shown to be expressed in that cell type as well.

We have replaced the HCR images for the mesenchyme craniofacial figure (new figure 3) with cropped (further zoomed in) images to illustrate the distinct expression domains of the mesenchymal markers we wish to highlight. We have replaced the HCR images showing pigment markers in I and J with new images (new Figure 2). We have added new HCR image data to show co-expression of novel genes of interest in the developing enteric neuronal population along the zebrafish gut, including genes of the Opioid family and Hox transcription factors. We thank reviewers for these suggestions.

7) "Hierarchical clustering of cells with General Mesenchyme and Chondrogenic identities using a cluster tree": This is not technically correct, assuming the authors used BuildClusterTree (which it seems so from the Materials and methods). This function does not operate on the level of cells, but hierarchically clusters the mean expression signature of the already-calculated clusters.

Thank you for pointing out this error. We indeed used the BuildClusterTree tool in Seurat to generate our cluster trees in the paper. We have clarified description of this in the text.

8) It's not clear why the authors present the data twice, first as each individual time point and then as a combined analysis, when the same clusters are recovered in the combined analysis. It seems as if it would be easier and more informative to just present the combined data.

We endeavored to show each dataset separately first, as we intended to showcase the sox10 cellular populations present at each time point to then allow for an appreciation of the combined dataset later shown. Indeed, we pulled important data from individual time point data, as evidenced by the analysis of enteric sub-clusters and the populations contained within them at the late time point (Figure 5). On the other hand, we leveraged the combined analysis to ask questions about global hox gene expression trends (Figure 7) and differential expression between time points (Figure 6). This demonstrates that examination of the dataset at multiple scales yields interesting results worth examining.

9) The authors provide an estimated number of cells reported by 10x Cell Ranger pipeline which corresponds to 2300 and 2580 for 48-50hpf and 68-70hpf, respectively. We understood that these numbers correspond to the estimated number of cells that 10x pipeline was able to identify which in the following steps serves as an input for Seurat analysis in R. The total number of usable cells after all QC steps were 1608 (58-60hpf) and 2410 (68-70hpf). However, there was no information on how many cells were loaded on to the 10x chip. This information would help others, especially in the zebrafish community, to better plan single-cell experiments and would provide additional information about cell suspension quality.

We thank the reviewers for this suggestion. We have added this information into the Figure 1—figure supplement 1 and in the Materials and methods.

10) The reviewers were surprised that there is a large population of neural crest classified cells at these stages. They wondered if those are progenitor cells and do they remain in adults? In addition, it seems like most of the neuronal population in enteric, where are the dorsal root ganglia and sympathetic neuronal populations?

Yes, we had wondered the same after initial examination of the datasets. We do not think the cells we classified as neural crest are merely progenitor cells. They express bona fide neural crest markers such as sox10, foxd3, ngfrb and tfap2a, and a zebrafish-specific neural crest marker, crestin. We confirmed crestin co-expression with foxd3 and ngfrb via HCR (Figure 4) and have also seen crestin+ cells in neural crest pathways at these time points before (Uribe et al., 2015, et al., 2018). We believe they are bona fide neural crest based on their signatures. We don’t know if neural crest cells with these signatures are remaining in adult zebrafish, it would certainly be interesting to explore in the future.

In terms of detecting dorsal root ganglia and sympathetic neuronal cells within the sox10 datasets, we have indeed detected dorsal root ganglia progenitor signatures in our first time point, which we described as sensory signatures in the results and shown in Figure 1—figure supplement 6. For sympathetic neuronal populations, upon closer examination of the second time point of our datasets, we have detected sympathetic signatures and created a new figure to describe this, Figure 5—figure supplement 2. We thank reviewers for the suggestions.

11) In the two stages of neural crest development, it would be interesting to include how the gene expression changes across developmental time. Could you present a combined UMAP with stages colored to see where stages fall and sets of genes that change expression through developmental time in subclusters? Highlighting genes that change over developmental time in each cluster would be very informative.

We agree with reviewers. We had originally included a combined UMAP with stages colored to see where stages fall within the atlas, in original Figure 8—figure supplement 1A. We now have expanded this analysis in a new figure, Figure 6—figure supplement 2, to include sets of genes that change expression in sub-clusters over time, as suggested by the reviewers. For this, we focused on subsets of the atlas to examine mesenchyme, pigment and neural/neuronal annotated cellular populations. From this we found expected genes from our prior analysis, as well as novel genes that have not been described before in these sox10+ cell populations and highlight them in the results. Globally, we agree with reviewers that being able to leverage the atlas to examine changing gene expression over time is informative, and thus, we have integrated the dataset into an interactive cell browser, hosted at the UC Santa Cruz Cell browser (described in the paper) where users can interactively view global cell populations, as well as genes of interest. We thank reviewers and hope this resource will help move the field forward.

12) Discussion penultimate paragraph: This is potentially an interesting conclusion that the authors could spend more time on in the paper. Is it worth generating some version of a figure that combines a review of what is known about these signatures from other systems (and which other systems) with what the authors observe in zebrafish? One could envision some form of dot plot with an organismal color code, but obviously the authors may have more creative ideas.

We have edited this statement out of our writing in the Discussion appropriately. While we agree that this is an interesting figure to synthesize on the topic, it is beyond the scope of what we aimed for our atlas and we intend to include this great suggestion in a future review article or future paper.

13) Mitochondrial genes were regressed out from the data set during the Cell Ranger alignment. Mitochondrial genes, together with ribosomal genes, are useful metrics for the assessment of cell quality. High expression levels of mitochondrial genes may indicate poor sample quality, apoptotic cells, or multiplets. Or if it is limited to a few clusters, it may reveal the nature of some cells (e.g. increased metabolic activity). However, the authors decided to remove mitochondrial genes from the alignment. Why?

We thank reviewers for catching this error in our methods description. We clarify that we did not remove mitochondrial genes from the alignment and have removed this statement from the Materials and methods.

14) The paper makes no mention of how the data and code will be made available. It would be standard practice to deposit the FASTQ files as well as the counts table output by CellRanger in GEO. Moreover, it's fairly standard practice (and generally appreciated) to make the processed Seurat objects available (either through lab website, Data Dryad, or some other hosting platform). Additionally, it's generally accepted practice to publish the code that was used to analyze the data and generate the figures.

We thank the reviewers for catching this omission. Indeed, we have deposited the FASTQ datasets to GEO, the details of which are now included at the end of the Materials and methods. Further, we have transformed the Seurat datasets into an interactive cell browser, hosted at the UC Santa Cruz Cell browser (described in the paper) where users can interactively view global cell populations, as well as genes of interest. The code used to analyze datasets are also available on our lab GitHub. We thank reviewers and hope this resource will help move the field forward.

[Editors' note: further revisions were suggested prior to acceptance, as described below.]

Revisions:

All reviewers thought that the authors have made significant improvement to the manuscript by removing many of the technical descriptions from their results. However they also thought that more can be done to make the Results section more readable. Many unnecessary details and repetitions remain, which obscures the communication of essential findings. Moreover, the manuscript is perceived as written for a neural crest audience and thus less accessible to a more general (eLife) audience making the authors motivation and essential findings obscured. The specific comments and suggestions of the reviewers that can guide additional revisions are appended below.

Reviewer #1:

My only remaining comment that should not block publication, but perhaps would inspire some additional revisions from the authors, is that overall, I still find that the prose is disappointingly focused on technical aspects of a pretty standard single-cell analysis at the expense of describing discoveries from those analyses. Many sections of the paper end with the main conclusion that the section "demonstrates the power of the sox10 atlas to.…" but I don't know why this is repeatedly considered a result, given that standard and well demonstrated analysis methods were used on data generated using now standard approaches. The new and expanded Discussion section illustrates that the authors have found some interesting results, but I feel that they dilute them significantly by focusing so much on telling us about the power of now-standard single-cell genomics approaches and about the process of annotating the data, when they could instead place the focus on their more interesting findings. For instance, I find this is especially true in the presentation of the mesenchyme - which seems mostly focused on describing existing mesenchyme markers and how they were used to identify the mesenchyme and then the number of clusters found therein, but without significant description of what those clusters are or what defines them.

We have edited the presentation of the mesenchyme results to focus on what we discovered defines the clusters, driving the annotations. Following annotations, we found various specific gene signatures among the mesenchyme clusters. For example, we found mesenchyme-identity cells expressing chondrogenic markers (barx1), proliferative signatures (pcna), migratory genes (snai2, rac1), stem-like signatures (uhrf1), etc. Description of what the clusters are and what genes defines them are included in the annotation table in Figure 1—figure supplement 5. For clarity, we now refer to the various mesenchyme clusters as transcriptionally-distinct mesenchyme subpopulations and/or clusters and have edited the text to reflect this. We leave further investigation of the mesenchyme populations open to researchers interested in the future.

We thank the reviewer for providing feedback to allow us to further improve the text. Overall, we have made edits throughout and removed a majority of the language that centered upon highlighting the technical methods of the analysis. We have edited the Results and Discussion sections of the manuscript to be more concise, while also ensuring sufficient details were maintained for the reader to understand the datasets. We have re-organized the discussion to highlight novel findings.

Reviewer #2:

Despite the new additions, the manuscript remains highly descriptive and a major issue in terms of writing style remains.

1) The authors have made significant improvement by removing many of the technical descriptions from their results. However this did little to make the Results section more readable. Many unnecessary details and repetitions remain, which obscures the communication of essential findings. The language used is often colloquial, with some grammatical errors which makes reading extremely challenging, especially for a non-expert in NCC or single cell biology.

We thank the reviewer for providing feedback to allow us to further improve the text. Overall, we have made edits throughout and removed a majority of the language that centered upon highlighting the technical methods of the analysis. We have edited the Results and Discussion sections of the manuscript to be more concise, while also ensuring sufficient details were maintained for the reader to understand the datasets. We have re-organized the discussion to highlight novel findings.

Having the manuscript read by a professional scientific editing service might be of help. To cite some examples:

– "We have utilized the Tg(-4.9sox10:EGFP) (hereafter referred to as 117 sox10:GFP) transgenic fish line to identify". Authors were requested to use the term "zebrafish" instead of "fish" since other fishes are used for these types of studies and it is confusing. Please correct this.

We removed the word “fish line” and replaced it with “transgenic line” in this sentence.

– What does "remarkably captured" mean?

We removed the word “remarkably”.

– "Taken together, these results show that the scRNA-seq datasets effectively identify discrete subpopulations, and coupled with our HCR analysis, effectively shows we are able to validate these cell populations in vivo". Please clarify which subpopulations are referred to here.

We have edited this conclusion statement to clarify which populations (pigment cells) are referenced.

– It is unclear what the authors meant by "early neuronal differentiation". Please specify the exact or range of developmental stage(s) for clarity.

We have clarified the range of developmental stages early neuronal differentiation refers to with regard to enteric cells discussed, which is now included in a revised paragraph.

– "Further investigation into these pathways led to the identification of…" – would it be sufficient to simply state for e.g. that "oprl1 and oprd1b, which are members of this pathway, were expressed in subcluster 3"?

We have edited the paragraph describing opioid pathway gene expression within enteric populations for clarity throughout.

– ". patterns of Homeobox transcription factors, known as hox genes…" – how can transcription factors be known as genes? Please rephrase this.

We have edited this phrasing to state, “expression of hox genes, which encode for Homeobox transcription factors…”.

2) The Discussion section is also extremely lengthy, containing many repetitions of the results which could be significantly condensed. There are also a lot of anecdotal information written on opioids which seems irrelevant and unnecessary.

We have performed edits to make the Discussion section more concise than the previous manuscript version. We have removed anecdotal information regarding opioids. Overall, we have made edits throughout the Discussion and removed all the repetitions we could find.

3) In describing the mesenchymal clusters, the authors used the term "subtypes". It seems more appropriate to call them simply as "clusters" as there were no systematic analyses done (e.g. expression of specific cell type markers) to define whether each of these barx1+ and dlx2a+ clusters represent distinct cell subtypes. Moreover, the number of clusters obtained depends on the threshold set in the scRNA-seq data analysis, which is therefore arbitrary in nature.

The reviewer has pointed an area worth further clarification. We originally called the cells “subtypes” due to the fact that they comprised the clusters we had designated with the major cell type identity of mesenchyme, yet also exhibited transcriptionally-distinct signatures which we carefully annotated by various marker genes denoting cell states. The table in Figure 1—figure supplement 5 lists the major cell type designation for each cluster, as well as individual annotations for each cell “subtype”, under each major cell type category. Following all annotations, we found various specific gene signatures among the mesenchyme clusters. For example, we found mesenchyme-identity cells expressing chondrogenic markers, proliferative signatures, migratory genes, stem-like signatures, etc. For clarity, we will refer to the various mesenchyme clusters as transcriptionally-distinct mesenchyme subpopulations and/or clusters and have edited the text to reflect this.

Along the same lines, the resolution parameter from Seurat FindClusters function determines the number of obtained clusters. Naturally, the selection of this parameter is arbitrary and needs to be individually tailored for each sample. Looking at the data at different resolutions may help to improve the identification of novel cell types. Is there any particular reason why the authors decided to apply the resolution of 1.2? A clarification on how this threshold was chosen would be helpful.

As one of the first data analysis methods we enacted when examining our datasets, we indeed looked at the data at different resolutions to help improve identification of novel cell types. We have added a sentence in the Materials and methods to describe this, ”The resolution value was empirically determined for the FindClusters by iterative assessment of cluster patterns generated values ranging from 0.4 to 2.0. The resolution was set to 1.2 which balanced between over segmentation of cell identities and effective separation of tSNE-mapped clusters.”

4) "… pbx3b expression may promote the assumption of an IPAN signature characterized by the presence of calb2a, ache, and slc18a3a and the loss of inhibitor markers nos1 and vipb." – it is unclear how the authors' observations led to this hypothesis as earlier on they stated that "…pbx3b expression was found in combination with calb2a, vipb, and nos1…". Please clarify.

We have edited the paragraph in the Results describing the IPAN signature expression in zebrafish to address this concern and clarify the data. Furthermore, we moved description of the model of IPAN fate acquisition (Figure 5F) to the Discussion, where it is more appropriately described in conjunction with the recent literature in mouse.

5) "Overall, these results suggest that sub-cluster 0 cells may represent a pool of immature sympatho-enteric neurons, and that sub-clusters 1 and 4 both represent better resolved, yet still immature, pools of enteric and sympathetic neurons." – it is unclear how the authors arrived at this notion, especially since cluster 0, 1, and 4 was not even mentioned at all previously.

We agree inclusion of more information about sub-clusters 0, 1 and 4 would strengthen and clarify the sub-cluster results. We have expanded this paragraph in the Results section to describe sub-clusters 0, 1 and 4 in more detail.

6) It would be helpful if the authors also annotate the clusters in Figure 6—figure supplement 2 panel A and Figure 7—figure supplement 1 panel J.

We agree. We have annotated the clusters in Figure 6—figure supplement 2A and Figure 7—figure supplement 1J.

7) "We performed cell filtering and clustering (.) cells which contained low (<200) or high (>2500) genes were removed from analysis".

– As indicated in Figure 1—figure supplement 1 panel C, around 4905 and 4669 cells were loaded for encapsulation into 10x cartridge, resulting in 2300 and 2580 identified cells. After a filtering step that was done on the basis of gene content only, 1608 and 2410 cells were obtained. The expected multiplet rate at such cell suspension concentration should be low (around 2% according to 10x protocol), therefore the great loss of cell numbers is rather surprising. Have the authors checked why the cell loss at the applied threshold level for the 48-50 hpf sample is much higher than expected. Also, did they check the distribution of expressed genes before setting up the thresholds for gene content per cell?

We investigated the distribution of expressed genes prior to setting a final threshold for our analysis, considering appropriate treatment for both time points analyzed. As we would have needed to apply the <2500 gene threshold to the 68-70 hpf dataset, we also accepted this value at the 48-50 hpf dataset to ensure as equal treatment between the two as possible. Selecting a relatively stringent threshold for the number of features (<200 and >2500) allowed us to produce a high-quality dataset. While it is possible that this diminishes discovering possible rare cell populations, the stringent datasets enabled us to confidently investigate the transcriptional landscape of neural crest. Re-evaluation of our 48-50 hpf dataset with a higher threshold (>200 and < 4000) retains a higher fraction of cells (2297 cells), but the inclusion of these cells would not significantly alter our subsequent analysis. At this point we feel that erring on the more stringent side to produce a high-quality dataset is preferrable to the reciprocal, especially when we ultimately intend for our analysis to serve as a basis for the community to use moving forward. Should information contained within those cells be of useful to a user in the future, they have been included in the raw data published on the NCBI’s Gene Expression Omnibus portal.

– Additionally, including mitochondrial gene content into the filtering step may help to confidently identify low quality cells.

We are grateful for the suggestion about using mitochondrial gene content in the filtering step to assist in detecting low quality cells. We have determined that the removal of the cells which express moderate levels of mitochondrial genes minimally impacts the dataset, particularly when paired with the strict threshold described above. To this effect, we have retained cells with appreciable levels of mitochondrial gene expression and feel confident that the resulting dataset is robustly reflective of the diverse transcriptional landscape expected in posterior neural crest at the stages sampled.

8) Figure 1—figure supplement 1 panel C – inconsistent number format (comma, dot, space), e.g. mean read per cell (74,017 vs 62109 etc.)

We have removed all commas in the table within Figure 1—figure supplement 1 to utilize consistent number formatting.

9) “…Major cell type categories were based on the presence of signature marker genes…”. Is the word "presence" an appropriate term in this context? Presence is rather a binary value and takes two options: 1- present, 0 – not present. It seems more accurate to use the term "expression" of these signatures/marker genes.

We have changed the word “presence” and instead use “expression”.

https://doi.org/10.7554/eLife.60005.sa2

Article and author information

Author details

  1. Aubrey GA Howard IV

    Department of BioSciences, Rice University, Houston, United States
    Contribution
    Conceptualization, Data curation, Formal analysis, Validation, Investigation, Visualization, Methodology, Writing - original draft
    Contributed equally with
    Phillip A Baker
    Competing interests
    No competing interests declared
  2. Phillip A Baker

    Department of BioSciences, Rice University, Houston, United States
    Contribution
    Conceptualization, Data curation, Formal analysis, Validation, Investigation, Visualization, Writing - original draft
    Contributed equally with
    Aubrey GA Howard IV
    Competing interests
    No competing interests declared
  3. Rodrigo Ibarra-García-Padilla

    Department of BioSciences, Rice University, Houston, United States
    Contribution
    Conceptualization, Data curation, Formal analysis, Validation, Investigation, Visualization, Methodology, Writing - original draft
    Competing interests
    No competing interests declared
  4. Joshua A Moore

    Department of BioSciences, Rice University, Houston, United States
    Contribution
    Conceptualization, Data curation, Formal analysis, Validation, Investigation, Visualization, Methodology, Writing - original draft
    Competing interests
    No competing interests declared
  5. Lucia J Rivas

    Department of BioSciences, Rice University, Houston, United States
    Contribution
    Validation, Investigation, Visualization, Writing - original draft
    Competing interests
    No competing interests declared
  6. James J Tallman

    Department of BioSciences, Rice University, Houston, United States
    Contribution
    Data curation, Validation, Visualization, Writing - review and editing
    Competing interests
    No competing interests declared
  7. Eileen W Singleton

    Department of BioSciences, Rice University, Houston, United States
    Contribution
    Validation, Investigation, Methodology, Writing - review and editing
    Competing interests
    No competing interests declared
  8. Jessa L Westheimer

    Department of BioSciences, Rice University, Houston, United States
    Contribution
    Validation, Investigation, Visualization, Writing - review and editing
    Competing interests
    No competing interests declared
  9. Julia A Corteguera

    Department of BioSciences, Rice University, Houston, United States
    Contribution
    Validation, Investigation, Visualization, Writing - review and editing
    Competing interests
    No competing interests declared
  10. Rosa A Uribe

    Department of BioSciences, Rice University, Houston, United States
    Contribution
    Conceptualization, Resources, Data curation, Formal analysis, Supervision, Funding acquisition, Validation, Investigation, Visualization, Methodology, Writing - original draft, Project administration, Writing - review and editing
    For correspondence
    rosa.uribe@rice.edu
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-0427-4493

Funding

Cancer Prevention and Research Institute of Texas (RR170062)

  • Rosa A Uribe

National Science Foundation (1942019)

  • Rosa A Uribe

SDB (Choose Development!)

  • Jessa L Westheimer

Houston Livestock Show and Rodeo Research Award

  • Phillip A Baker
  • Joshua A Moore

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

Funding for this project was provided by Rice University, Cancer Prevention and Research Institute of Texas (CPRIT) Recruitment of First-Time Tenure Track Faculty Members (CPRIT-RR170062) and the NSF CAREER Award (1942019) awarded to RAU, a Houston Livestock Show and Rodeo Research Award to JAM and PAB, and a SDB Choose Development! Fellowship award to JLW. We acknowledge the Cytometry and Cell Sorting Core at Baylor College of Medicine, which is funded from the CPRIT Core Facility Support Award (CPRIT-RP180672), the NIH (P30 CA125123 and S10 RR024574), and the expert assistance of Joel M Sederstrom for assistance with flow cytometry. Single cell library preparation, Illumina sequencing, and Cell Ranger alignment was facilitated by Advanced Technology Genomics Core at MD Anderson Cancer Research Center funded by CA016672(ATGC). IMARIS image analysis was performed using Rice University’s Shared Equipment Authority (SEA) IMARIS workstation. We thank George Eisenhoffer and Oscar Ruiz (MD Anderson) for advice regarding flow cytometry and single-cell RNA-seq methodology. We thank Sarah Kucenas (University of Virginia) for helpful advice on glial populations. We thank Robert Naja and Robyn Fenty for technical assistance.

Ethics

Animal experimentation: All animal work was performed under protocols approved by, and in accordance with, the Rice University Institutional Animal Care and Use Committee (IACUC), under protocol 1144724.

Senior Editor

  1. Richard M White, Memorial Sloan Kettering Cancer Center, United States

Reviewing Editor

  1. Lilianna Solnica-Krezel, Washington University School of Medicine, United States

Reviewers

  1. Cecilia Lanny Winata
  2. Kristin Artinger

Publication history

  1. Received: June 14, 2020
  2. Accepted: January 31, 2021
  3. Version of Record published: February 16, 2021 (version 1)

Copyright

© 2021, Howard et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 2,719
    Page views
  • 283
    Downloads
  • 4
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Download citations (links to download the citations from this article in formats compatible with various reference manager tools)

Open citations (links to open the citations from this article in various online reference manager services)

Further reading

    1. Developmental Biology
    2. Physics of Living Systems
    Yonghyun Song, Changbong Hyeon
    Research Article Updated

    Spatial boundaries formed during animal development originate from the pre-patterning of tissues by signaling molecules, called morphogens. The accuracy of boundary location is limited by the fluctuations of morphogen concentration that thresholds the expression level of target gene. Producing more morphogen molecules, which gives rise to smaller relative fluctuations, would better serve to shape more precise target boundaries; however, it incurs more thermodynamic cost. In the classical diffusion-depletion model of morphogen profile formation, the morphogen molecules synthesized from a local source display an exponentially decaying concentration profile with a characteristic length λ. Our theory suggests that in order to attain a precise profile with the minimal cost, λ should be roughly half the distance to the target boundary position from the source. Remarkably, we find that the profiles of morphogens that pattern the Drosophila embryo and wing imaginal disk are formed with nearly optimal λ. Our finding underscores the cost-effectiveness of precise morphogen profile formation in Drosophila development.

    1. Developmental Biology
    2. Neuroscience
    Tania Moreno-Mármol et al.
    Research Article

    The vertebrate eye-primordium consists of a pseudostratified neuroepithelium, the optic vesicle (OV), in which cells acquire neural retina or retinal pigment epithelium (RPE) fates. As these fates arise, the OV assumes a cup-shape, influenced by mechanical forces generated within the neural retina. Whether the RPE passively adapts to retinal changes or actively contributes to OV morphogenesis remains unexplored. We generated a zebrafish Tg(E1-bhlhe40:GFP) line to track RPE morphogenesis and interrogate its participation in OV folding. We show that, in virtual absence of proliferation, RPE cells stretch and flatten, thereby matching the retinal curvature and promoting OV folding. Localized interference with the RPE cytoskeleton disrupts tissue stretching and OV folding. Thus, extreme RPE flattening and accelerated differentiation are efficient solutions adopted by fast-developing species to enable timely optic cup formation. This mechanism differs in amniotes, in which proliferation drives RPE expansion with a much-reduced need of cell flattening.