Single-Cell Transcriptomics Reveals Evolutionary Reconfiguration of Embryonic Cell Fate Specification in the Sea Urchin Heliocidaris erythrogramma

Abstract

Altered regulatory interactions during development likely underlie a large fraction of phenotypic diversity within and between species, yet identifying specific evolutionary changes remains challenging. Analysis of single-cell developmental transcriptomes from multiple species provides a powerful framework for unbiased identification of evolutionary changes in developmental mechanisms. Here, we leverage a "natural experiment" in developmental evolution in sea urchins, where a major life history switch recently evolved in the lineage leading to Heliocidaris erythrogramma, precipitating extensive changes in early development. Comparative analyses of single-cell transcriptome analysis (scRNA-seq) developmental time courses from H. erythrogramma and Lytechinus variegatus (representing the derived and ancestral states, respectively) reveal numerous evolutionary changes in embryonic patterning. The earliest cell fate specification events and the primary signaling center are co-localized in the ancestral developmental gene regulatory network; remarkably, in H. erythrogramma, they are spatially and temporally separate. Fate specification and differentiation are delayed in most embryonic cell lineages, although in some cases, these processes are conserved or even accelerated. Comparative analysis of regulator-target gene co-expression is consistent with many specific interactions being preserved but delayed in H. erythrogramma, while some otherwise widely conserved interactions have likely been lost. Finally, specific patterning events are directly correlated with evolutionary changes in larval morphology, suggesting that they are directly tied to the life history shift. Together, these findings demonstrate that comparative scRNA-seq developmental time courses can reveal a diverse set of evolutionary changes in embryonic patterning and provide an efficient way to identify likely candidate regulatory interactions for subsequent experimental validation.

Keywords: cell fate specification; evolution of development; life history evolution; scRNA-seq; transcriptional regulation.

Fig. 1.

Comparison of single-cell developmental transcriptomes. a) Time tree of sea urchin species with high-quality reference genomes. Egg and larva sizes are approximately to scale; the ancestral life history is characterized by small eggs and feeding larvae (planktotrophy) and the derived life history by large eggs and nonfeeding larvae (lecithotrophy). Egg and larva sizes from Mortensen (1921), Emlet et al. (1987), Williams and Anderson (1975); topology from Láruson (2017); divergence times from Zigler et al. (2003) and Láruson (2017). b) UMAP plots of scRNA-seq developmental time course for L. variegatus. The large plot shows cells color coded by cluster and labeled according to inferred cell types; the two smaller plots show cells color coded by time point (upper) and the centroids of the six time points common to both species (lower plot). The earliest time point (2 hpf) cells are labeled “totipotent” based on blastomere separation experiments (Hörstadius 1973; Davidson et al. 1998). c) UMAPs of scRNA-seq developmental time course for H. erythrogramma. Organization parallels (b). For individual marker gene expression, see supplementary fig. S2, Supplementary Material online. Clusters in (b) and (c) are colored with the same encoding to facilitate comparison between species (some cell types are present only in one species or the other). Early blastomeres are not totipotent in this species (Henry and Raff 1990; Henry et al. 1990), but a distinct cluster of cells (labeled pluripotent) retains pluripotency into the larva (McDonald et al. 2024). d) Comparison of cell-type proportions and larval morphology. Proportions of four cell types in 24 hpf larvae (see supplementary tables S1 and S2, Supplementary Material online, for cell counts at all stages). Simplified diagrams of larvae are not to scale; colors match bar plot.

Fig. 2.

Integrated single-cell transcriptomes. Developmental transcriptomes of 1:1 orthologs in the two species were integrated using CCA (Butler et al. 2018), followed by dimensionality reduction. a) UMAP color coded by species. Relevant cell clusters are indicated. b) UMAP color coded by developmental time point. The earliest stages sampled are at the upper left, and development generally progresses down and to the right. c) UMAP color coded by degree of differentiation based on CytoTRACE (Gulati et al. 2020). Differentiation tracks developmental time points, with the lowest scores corresponding to the earliest time points and the highest scores corresponding to clusters of differentiated cells. d–g) Cells from 6 and 9 hpf H. erythrogramma do not overlap L. variegatus cells from the same time points (top UMAPs), and instead overlap earlier time points (bottom UMAPs). h, i). Cells from 20 and 24 hpf H. erythrogramma converge on L. variegatus cells from the same time points.

Fig. 3.

Temporal shifts in transcriptomes. a) Heat maps showing degree of similarity among scRNA-seq transcriptomes for four embryonic cell lineages in the two species. Assignment of cells to lineages is based on optimal transport (see Materials and Methods). Boxes highlight the most similar time points. N-sk. Mes., nonskeletogenic mesenchyme; Skel. Mes., skeletogenic mesenchyme. b) Line plots of the most similar time points in (a) reveal an overall delay in H. erythrogramma, with most points above the line defined by a slope of 1. c) Heat map showing degree of similarity for the entire scRNA-seq transcriptome at each stage based on CIDER (Hu et al. 2021). The first few time points in H. erythrogramma are most similar to earlier time points in L. variegatus. d) Plot of CytoTRACE scores (Gulati et al. 2020) during development show an initial delay in differentiation in H. erythrogramma, followed by convergence after ∼16 hpf. e) Line plot showing developmental time of morphogenetic events for comparison. Again, there is an overall delay in H. erythrogramma.

Fig. 4.

Evolutionary changes in timing of differentiation. Optimal transport was used to predict the likely fate for each cell at five stages, based on transcriptomes at 24 hpf (see Materials and Methods). Triangle plots show transcriptomes predictive of blastocoelar cell, skeletogetogenic cell, or any other cell fate near the top, right, and bottom vertices, respectively; cells with undifferentiated transcriptomes occupy the center. Corresponding UMAPs are shown below. Note the much earlier differentiation of skeletogenic cells in L. variegatus and the slightly earlier differentiation of blastocoelar cells in H. erythrogramma. See text for additional interpretation.

Fig. 5.

Measuring transcript co-expression. a) If gene A encodes a transcription factor that regulates the expression of gene B, both must be transcribed in the same cell (with some rare exceptions). Cluster-level co-expression measures whether both genes are expressed within the same cell cluster (pseudobulk). However, cluster-level co-expression can occur even if the two genes are never expressed in the same cells. This can happen when a cluster contains cells with diverging transcriptional states, a situation that arises during every cell fate specification event throughout development. Cell-level co-expression is a more stringent criterion that requires transcripts from both genes A and B to be present in the same cell, and is the definition applied in the present study. b) Example of an experimentally validated regulatory interaction that is reflected in the distribution of cell-level co-expression. alx1 encodes a transcription factor that activates transcription of dri within the skeletogenic cell lineage (Oliveri et al. 2008). Co-expression of alx1 and dri is detectable but low or very low and in a minority of cells in endoderm and ectoderm. Circles indicate regions shown at 2× magnification to the right of each UMAP. Robust and nearly universal co-expression occurs only within the skeletogenic cell lineage, precisely where it is predicted to occur. These results are consistent with an experimentally tested interaction and further imply that alx1 does not influence dri expression outside the skeletogenic cell lineage. In this and all subsequent co-expression UMAPs, light dots represent very low co-expression (only 1 read detected from either or both genes), and dark dots represent moderate to high co-expression (at least 2 reads detected from both genes).

Fig. 6.

Inference of evolutionary changes in regulatory interactions based on proportion of co-expressing cells. a) Simplified version of the skeletogenic portion of the ancestral dGRN present in Camarodonta sea urchins with feeding larvae (adapted from Figs. 6 and 7 of Rafiq et al. (2012); built on data from Kurokawa et al. (1999), Davidson et al. (2002a), Oliveri et al. (2002, 2008), Ettensohn et al. (2003), Oliveri and Davidson (2004), Sharma and Ettensohn (2010), Rafiq et al. (2012)). The three primary activators of skeletogenic-specific transcription (top) feed directly or indirectly into a large set of effector genes, some of which are illustrated (bottom). b) Co-expression analysis of gene pairs involved in 11 experimentally validated regulatory interactions in the ancestral state (McDonald et al. 2024), represented by L. variegatus and compared with expression in H. erythrogramma. Numbers correspond to interactions in (a). The top two plots for each interaction show expression of regulator and target based on bulk RNA-seq (Israel et al. 2016), with a log₂  y axis. The dashed line indicates very low expression (an average of 5 counts per million reads across time points, averaged across 3 biological replicates), which is effectively the lower limit of reliable detection (Israel et al. 2016). Supplementary fig. S11, Supplementary Material online shows larger plots with values. The plot directly below shows the proportion of cells that co-express both genes based on scRNA-seq, with a linear y axis; these time points begin at 6 hpf, the first time point common to both data sets. Note that y axes are not equivalently scaled because genes have a wide range of expression and co-expression levels. Most gene pairs show a strong peak of co-expression at 9 hpf in L. variegatus, which then drops as skeletogenic cells stop dividing while other cell lineages continue to proliferate. In contrast, this peak is notably absent in H. erythrogramma; instead, co-expression is initially zero or very low at 9 hpf and rises modestly 16 to 24 hpf. These results point to a general delay in co-expression of regulatory and target in H. erythrogramma relative to L. variegatus.

Fig. 7.

Inference of evolution changes in regulatory interactions based on the distribution of co-expressing cells. Co-expression analysis of experimentally validated regulatory interactions. UMAPs show the location of cells with co-expression of indicated regulator and target. L. variegatus = green dots and H. erythrogramma = orange dots; dark colors indicate cells with >2 UMIs for both regulator and target gene; pale colors indicate low co-expressing cells, where one or both genes have 1 or 2 UMIs. Boxes indicate areas shown at 2× in the right-hand column and arrows indicate skeletogenic cells in H. erythrogramma. See supplementary fig. S6, Supplementary Material online for a visual guide to interpreting evolutionary differences.

Fig. 8.

Inference of evolutionary loss of a regulatory interaction. a) Density plots showing expression of regulator (alx1) and target (foxB) genes in both species. Note that both alx1 and foxB transcripts are readily detected in both species. b) Co-expression plots. The complete absence of co-expression in H. erythrogramma suggests that the ancestral alx1-foxB regulatory interaction has been lost in this species. See supplementary fig. S6, Supplementary Material online for a visual guide to interpreting evolutionary differences.