Single-cell analysis of shared signatures and transcriptional diversity during zebrafish development

Abstract

During development, animals generate distinct cell populations with specific identities, functions, and morphologies. We mapped transcriptionally distinct populations across 489,686 cells from 62 stages during wild-type zebrafish embryogenesis and early larval development (3-120 hours post-fertilization). Using these data, we identified the limited catalog of gene expression programs reused across multiple tissues and their cell-type-specific adaptations. We also determined the duration each transcriptional state is present during development and suggest new long-term cycling populations. Focused analyses of non-skeletal muscle and the endoderm identified transcriptional profiles of understudied cell types and subpopulations, including the pneumatic duct, individual intestinal smooth muscle layers, spatially distinct pericyte subpopulations, and homologs of recently discovered human best4+ enterocytes. The transcriptional regulators of these populations remain unknown, so we reconstructed gene expression trajectories to suggest candidates. To enable additional discoveries, we make this comprehensive transcriptional atlas of early zebrafish development available through our website, Daniocell.

Figure 1:. A high temporal resolution single-cell RNAseq timecourse encompassing embryogenesis and early larval development.

(A) Single-cell transcriptomes were collected from whole zebrafish embryos at 50 different developmental stages (colored dots) between 14–120 hpf and then merged with our previous dataset encompassing 3.3–12 hpf (Farrell et al., 2018). Size of dots represents the number of cells recovered from each stage. (B–C) UMAP projection of single-cell transcriptomes, colored by (B) developmental stage (colored as in Fig. A) and (C) curated major tissues. (D) Distribution of the coefficient of variation (CV) of cluster means of expression for each gene (log-transformed). Lower CV indicates relatively similar expression across all clusters, while high CV indicates high variation in expression, either temporally or across cell types. Genes were divided into categories based on thresholds (red dashed lines). (E) Schematic showing our approach for identifying transcriptionally similar cells using an epsilon (ε) neighborhood approach and determining whether each cell was in a ‘short-term’ or ‘long-term’ state (based on the mean of absolute stage difference between the analyzed cell and its ε-neighbors). (F) Timeline bar plots showing the duration of ‘long-term’ cycling cell states identified using a 36-hour threshold. Each bar represents a cell population that was identified as ‘long-term,’ and the length of bar represents the minimum timespan that encompasses 80% of its ε-neighbors. hpf: hours post fertilization, PGCs: primordial germ cells, LL: lateral line.

Figure 2:. Identification of gene expression programs shared by multiple distinct cell types during development.

(A) Binary heatmap showing expression domains of gene expression programs (“GEPs”, x-axis) reused in multiple tissues (y-axis) during development. Annotations for select GEPs are shown. Asterisk indicates the GEPs further investigated in panels B, C, D, and Supplementary Fig. 5. Full GEP annotation can be found in Supplementary Table 2. (B) Dot plot showing expression of Megalin-associated genes and SLC transporters from module GEP-193 that are shared between intestinal lysosome-rich enterocytes (LREs) and pronephros proximal convoluted tubule (PCT) and pronephros proximal straight tubule (PST) cells. (C) Dot plot of top-loaded genes from five epithelial GEPs: one shared universally across all epithelial cell types (GEP-106), two comprising classical epithelial genes (GEP-100, GEP-37), and two tissue-specific (GEP-91, GEP-16). X-axis: cell types with epithelial characteristics across several tissues, and muscle as a non-epithelial control. (D) Dot plot of shared and tissue-specific members of module GEP-94, associated with mucin O-glycosylation. For panels B–D, color: mean expression per cell type; size: percent of cluster cells that express each gene on the Y-axis. B3GNTs: β1,3-N-acetylglucosaminyltransferases; GALNTs: N-acetylgalactosaminyltransferases; RPE: retinal pigmented epithelium; EVL: enveloping layer.

Figure 3:. Sub-clustering of non-skeletal muscle reveals distinct pericyte subtypes.

(A) UMAP projection of 3,866 non-skeletal muscle cells, numbered and color-coded by cluster. (B) Dot plot of selected differentially expressed pericyte-specific markers (x-axis) compared to vascular (vaSMCs) and visceral SMCs (viSMCs). (C) Dot plot of selected differentially expressed genes (y-axis) between the three pericyte clusters (x-axis, P0–P2) and myofibroblasts compared to vascular SMCs (x-axis, vaSMCs). See Supplementary Fig. 6C for additional markers. (D) Expression of pericyte marker genes visualized on the UMAP projection. Color bar shows mean expression level of each gene. (E–G”) RNA in situ hybridization for markers specific to the pericyte-2 population on a flk::mCherry-CAAX background in 5 dpf larvae. Panels E–E” indicate lateral view of the whole zebrafish head for the universal pericyte marker ndufa4l2a and pericyte-2 marker epas1a. Panels F–G” indicate a higher magnification of the brain posterior cerebral vein with epas1a⁺ pericytes. Arrows indicate ndufa4l2a⁺/epas1a⁺ pericytes in contact with as well separated from the posterior cerebral vein. In panels G–G”, arrowheads indicate ndufa4l2a⁺/epas1a⁻ pericytes along other hindbrain vessels while asterisks indicate autofluorescent red blood cells inside blood vessels. (H) Bar graph quantifying the number of ndufa4l2a⁺/epas1a⁺ pericytes that were visible in a similar-sized field of view near the posterior cerebral vein, per animal. (I–K”) RNA in situ hybridization for ndufa4l2a and epas1a across other blood vessels in the zebrafish head including forebrain (I–I”), eye (J–J”), and pharyngeal arches (K–K”). Arrowheads mark cells that are ndufa4l2a⁺; no ndufa4l2a⁺/epas1a⁺ cells were observed in these regions. (L) Jitter plot showing the proportion of ndufa4l2a⁺ pericytes that were also epas1a⁺ in different regions of the zebrafish head. Error bars indicate standard error of mean (S.E.M). PA: pharyngeal arches. Scale bar: 25 μm.

Figure 4:. Distinct gene expression within longitudinal and circular intestinal smooth muscle cells.

(A) Dot plot of top differentially expressed genes (y-axis) between circular and longitudinal SMC clusters (x-axis). (B–D) Feature plots of common (acta2, cald1b) and differentially expressed (il13ra2, tesca, fsta, kcnk18) markers of intestinal SMCs. (C–J) RNA in situ hybridization of markers expressed in longitudinal (il13ra2) and/or circular smooth muscle (fsta) cells lining the zebrafish intestinal tract. (C, E, G, I) anterior intestine view, including the intestinal bulb. (D, F, H, J) Higher magnification images of the intestinal bulb with nuclear (Hoecsht 33322) counterstain showing spatially segregated association of il13ra2 and fsta with differently oriented nuclei. (K) Diagram showing the orientation of circular and longitudinal SMCs along the intestinal tract from a lateral view (top) and cross-section (bottom). (L–M) Expression of longitudinal (il13ra2) and circular (fsta) SMC-specific markers shown from a cross-sectional view. (N) Force-directed layout of an URD-inferred hierarchical tree calculated on foxc1a/b⁻ and prrx1a/b⁻ (putatively non-neural crest derived) pericytes and smooth muscle cells. The trajectory tree is colored by stage, as in Fig. 1A. (O) Expression of intestinal smooth muscle layer-specific TFs and markers visualized on the URD trajectory. Scale bar: 25 μm.

Figure 5:. Subclustering of endodermal derivatives enables molecular characterization of human disease associated cell types.

(A) UMAP projection of 12,592 endodermal cells during zebrafish development, color coded and numbered by cluster. (B) Expression of specific (sftpba, sim1b) and strongly expressed (ihha) pneumatic duct markers visualized on the UMAP projection. Color bar shows expression of each gene. (C) Dot plot of top differentially expressed genes (y-axis) within the pneumatic duct compared to other endodermal derivatives (x-axis). (D–D”) RNA in situ hybridization of two specific markers of the pneumatic duct (sftpba and sim1b). Yellow arrowheads indicate staining in the pneumatic duct and inflated posterior swim bladder, white arrows denote staining in the anterior swim bladder bud primordium that inflates at 21 dpf. (E) Dot plot showing top differentially expressed markers (y-axis) in best4⁺ enterocytes compared to other intestinal cell types. (F) Expression of general intestinal marker (cdx1b) and two best4⁺ enterocyte markers (best4 and otop2) shown on UMAP projection. White arrowhead indicates otop2 expression in the posterior lysosome-rich enterocyte (LRE) cluster, and yellow arrowhead indicates expression within the best4⁺ enterocytes. (G–H”) RNA in situ hybridization of best4⁺ enterocyte marker (best4) against (G–G”) cdx1b and (H–H”) otop2. Yellow arrowheads indicate best4 and otop2 co-expressing cells in the anterior intestine, white arrowheads indicate a posterior patch of otop2 staining likely corresponding to the posterior LREs as shown in (F). Scale bar: 50 μm. EC: enterocyte; prog: progenitors; LREs: lysosome-rich enterocytes; EECs: enteroendocrine cells; PP: pancreatic polypeptide

Figure 6:. Trajectory analysis of zebrafish intestinal cells reveal candidate regulators for best4⁺ enterocyte specification.

(A) Venn diagram showing number of differentially expressed genes shared between human colonic (Smillie et al., 2019), human small intestinal (Burclaff et al., 2022), and zebrafish best4⁺ enterocytes. (B) Average log-fold enrichment of genes in best4⁺ enterocytes compared to best4⁺ enterocyte subtypes in human colon (x-axis, Smillie et al. 2019) and zebrafish intestine (y-axis). Blue: human colon-specific best4⁺ enterocyte markers; red: zebrafish-specific best4⁺ enterocyte markers; black: shared markers. Comparison of zebrafish best4⁺ enterocytes to human small intestinal best4⁺ enterocytes shown in Supplementary Figure 12. (C) Force-directed layout of a URD-inferred hierarchical tree generated with zebrafish intestinal cells between 14–120 hpf, colored by developmental stage as in Fig. 1A. (D) Expression of transcription factors and characteristic markers of individual intestinal cell types visualized on the URD trajectory. (E, F) Temporal dynamics of selected genes along the best4⁺ enterocyte (E) and posterior LRE (F) trajectory; lines show impulse response fits across pseudotime. Y-axis: scaled expression. (G–H”) RNA in situ hybridization of best4 and candidate transcriptional regulator pbx3a in a 5 dpf old zebrafish intestine. (H–H”) Higher magnification of yellow boxes from G–G”. Scale bar – 100 μm.