Dissecting mechanisms of chamber-specific cardiac differentiation and its perturbation following retinoic acid exposure

The process of early heart development is complex, involving multiple progenitors that contribute differentiated progeny to distinct regions of the heart (Kelly et al., 2014; Meilhac and Buckingham, 2018; Meilhac et al., 2014; Yao et al., 2021; Zhang et al., 2021). Shortly after gastrulation, Mesp1+ cells from the lateral plate mesoderm (LPM) form two distinct progenitor populations termed the first heart field (FHF) and the second heart field (SHF) (Chabab et al., 2016; Chiapparo et al., 2016; Devine et al., 2014; Lescroart et al., 2014; Saga et al., 1996). Together, these form the cardiac crescent (CC), which later balloons anteriorly into the primitive heart tube (PHT) and undergoes looping (Meilhac et al., 2014). The Isl1+ SHF forms posteriorly to the FHF and persists as a progenitor population, continuing to proliferate and contribute cells to the developing heart (Cai et al., 2003; Laugwitz et al., 2005). Derivatives from the FHF and SHF are distributed asymmetrically within the four chambers of the heart (Cai et al., 2003; Liang et al., 2013; Meilhac et al., 2014; Mjaatvedt et al., 2001; Waldo et al., 2001; Zaffran et al., 2004) and dysregulation in these progenitors is thought to underlie a number of congenital heart defects (CHDs) (Bruneau, 2008, 2020; Kloesel et al., 2016; Pierpont et al., 2018).

Defects such as atrial/ventricular chamber size imbalances are not adequately described by the existing FHF/SHF model, raising questions on whether specific subpopulations of progenitors are responsible for these phenotypes and about the precise mechanisms governing atrial and ventricular specification. Retinoic acid (RA) acts as a morphogen driving atrial differentiation (Bernheim and Meilhac, 2020; Niederreither et al., 2001; Perl and Waxman, 2020; Stefanovic and Zaffran, 2017). However, the transcriptional mechanisms downstream of this signaling pathway remain elusive. A number of efforts aimed at understanding the mechanisms of atrial/ventricular lineage segregation have performed bulk and single cell RNA-sequencing (scRNA-seq) as early as embryonic day (E) 9 (DeLaughter et al., 2016; Li et al., 2016). Although informative, these studies did not capture the key transcriptional processes that occur before chamber morphogenesis. More recent studies have profiled the transcriptome of the developing heart shortly after gastrulation (Ivanovitch et al., 2021; Lescroart et al., 2018; Tyser et al., 2021) or at the onset of chamber morphogenesis (Mantri et al., 2021; de Soysa et al., 2019; Xiong et al., 2019; Zhang et al., 2021). These studies highlight the transcriptional heterogeneity of cardiac progenitors and have identified markers that segregate distinct cell populations in a spatio-temporal manner. However, the specification and differentiation of atrial or ventricular cells, particularly before acquisition of chamber-specific marker expression, remains to be explored further. We and others have recently expanded on work demonstrating that atrial/ventricular lineage segregation occurs early during mammalian development, perhaps before FHF/SHF formation (Bardot et al., 2017a; Chabab et al., 2016; Garcia-Martinez and Schoenwolf, 1993; Keegan et al., 2004; Lescroart et al., 2014; Yutzey and Bader, 1995). We have shown that transient expression of Foxa2 in the primitive streak labels a population of cardiac mesoderm that migrates at the leading edge of the LPM. Foxa2 lineage-tracing experiments show that these early migrating cells contribute to the majority of the ventricular, but not atrial, myocardium (Bardot et al., 2017a).

To interrogate the molecular identity of atrial and ventricular progenitors, as well as lineage relationships between progenitors and their differentiated progeny, we performed scRNA-seq of Foxa2Cre;mTmG embryos at the CC (E8.25), PHT (E8.75) and heart tube (HT) (E9.25) stages. RNA velocity revealed relationships between distinct progenitors and atrial/ventricular cells and identified putative driver genes of differentiation that segregated over inferred transcriptional time. Through differential expression analysis we identified a number of candidate markers that correlate with atrial or ventricular differentiation over time. Lastly, we found that exogenous RA negatively impacts ventricular development by preventing differentiation from separate progenitors towards a ventricular fate. This led to defects in ventricular size and shifts in differentiation trajectories towards other lineages. In summary, we provide a transcriptional atlas across cardiac progenitor specification and differentiation at single cell level, illustrating deployment of distinct mechanisms along progenitor differentiation trajectories that are disrupted when cardiac development is perturbed by exogenous RA.

We sub-dissected cardiac regions as well as the surrounding areas comprising the gut tube and pharyngeal structures of Foxa2Cre;mTmG embryos, which allowed us to identify ventricular-fated progenitors at early stages (Bardot et al., 2017a) (Fig. 1A and Fig. S1A). Three age-matched embryos from the CC, PHT and HT stages each were pooled and a minimum of 10,000 cells sequenced (CC: 14,355; PHT: 13,036; HT: 21,318). We performed clustering and differential gene expression (DGE) analysis to identify multiple mesoderm and cardiac cell types, as well as cells from ectoderm and endodermal lineages (Fig. S1B). Importantly, we found that we can identify Foxa2 lineage-traced cells through EGFP expression without cell-sorting, which might have resulted in loss of rare cell types. We merged the data from the three time points to create a time-resolved atlas of the developing heart and surrounding tissues, identifying 27 clusters that separated broadly according to germ layer identity (Fig. S2A). DGE analysis determined signature genes defining populations of differentiated cardiomyocytes and associated progenitors, pharyngeal endoderm and mesoderm, as well as various definitive ectoderm and endoderm cell types among others (Fig. S2B-D).

To focus on the cardiovascular lineages, we subset and re-clustered cardiomyocytes, cardiac progenitors and mesoderm cells along with neural crest, endocardial, endothelial and epicardial cells at each stage (Fig. S3). We merged these data to resolve 27 clusters including heart field progenitors, cardiomyocytes at various stages of differentiation, as well as other supporting cell types of the heart, which were identified through known markers and candidate genes identified by DGE (Fig. 1B-D). Progenitors and differentiated cells segregated in part by cell cycle state, with Pdgfra+ mesoderm populations more commonly found in S phase and differentiated Nkx2-5+ progeny more commonly found in G2 M/G1 phase (Fig. 1E). This was also reflected by the contribution to each cluster from the different samples, with differentiated cardiomyocytes [cluster (C) 4, C2, C17] derived primarily from the PHT and HT stages, whereas the LPM (C26) primarily derived from the CC stage (Fig. 1F). Pdgfra+/Isl1+ SHF progenitors segregated into Tbx1+ anterior and Aldh1a2+/Hoxb1+/Foxf1+ posterior domains (Hoffmann et al., 2014; Rana et al., 2014; Stefanovic et al., 2020) (Fig. 1C,D). Intriguingly, although we did identify an FHF population expressing Nkx2-5/Hcn4/Tbx5 when examining the CC stage alone (Figs S3 and S4A), when merged with later time points these cells did not form a discrete population. Instead FHF cells clustered with multiple differentiating cardiomyocyte populations (C7, C12, C14), underscoring the identity of the FHF as a population of cells that most closely resembles differentiating cardiomyocytes. Differentiated cardiomyocytes segregated into several subpopulations including Kcna5+ atrial and Myl2+ ventricular cells of left/right identity (Fig. 1C; Fig. S4B,C). Atrial cells did not segregate discretely into left/right populations, though did show markers of left/right identity such as Pitx2 (Fig. S4D). We also identified populations corresponding to the Rspo3+ atrioventricular canal (AVC) and Hcn4+/Shox2+ sinus venosus (SV), which clustered near atrial cardiomyocytes.

In addition to cardiomyocytes, we identified Sox2+ neural crest cells, as well as Wt1+ epicardial cells (Fig. 1C,D). Localized distantly to the cardiomyocytes we identified several Pecam1+ endothelial populations (C9, C13, C16), which clustered closely with Nfatc1+ endocardial cells (C23). Both populations are proximal to Kdr+ (C3) and Etv2+ (C24) clusters that likely represent hematoendothelial cells. Curiously, we also identified a population of cardiomyocytes (C0) and mesoderm cells (C15) which clustered separately from the other heart field progenitors/cardiomyocytes and associated more closely with endothelial cells. We were unable to determine the exact source of this observation as these cells expressed similar markers to other mesoderm and cardiomyocytes. Positioning of these clusters could suggest that they represent a distinct subset of mesoderm/cardiomyocytes, or alternatively it could be related to comparatively lower transcript expression in C0/15 and the adjacent endothelial populations.

We were intrigued to find multiple clusters of differentiated Tnnt2+ cardiomyocytes (C0, C2, C4, C5, C7, C12, C17), only some of which had clear atrial/ventricular or left/right identity. To determine if these represented functional subtypes we performed DGE and GO term analysis to identify differentially regulated processes in each cluster (Fig. 1G; Fig. S4C). We found that processes involved in ‘mitochondrial ATP synthesis coupled electron transport’ or ‘striated muscle hypertrophy’ segregated myocardial cells into two broad categories of more or less mature cells. Clusters that did not show enrichment in energetic processes were instead characterized by developmental processes of migration and/or differentiation. Others, such as those in C0, were enriched for terms involved in cardiac muscle contraction and ribosome biogenesis, but few other more advanced processes, leading us to conclude that these represented early myocardial precursors. Cells of the SV (C12) were enriched for genes involved in ‘SA nodal action potential’, suggesting that they contained a subpopulation of pacemaker cells. When subclustered, the SV cells broadly segregate into an Nkx2-5 population expressing sinoatrial nodal markers such as Smoc2, Hcn4 and Shox2 and an Nkx2-5+ population that instead expressed markers of atrial cardiomyocytes (Fig. S5A). We found processes involved in aortic and pulmonary valve development to be enriched in one of the Tnnt2+ clusters (C5). Given the proximity to Isl1/Tbx1-expressing anterior second heart field (aSHF), we hypothesized that this cluster might represent aSHF-derived cells of both the outflow tract (OFT) and right ventricle (RV). We subclustered C5 and found that the cells separated into three distinct populations consisting of two populations of Myl2+ cells and a separate population enriched for Crabp1, Rgs5, Flna and 3632451O06Rik (also known as Armh4), markers of the OFT (Fig. S5B). The clustering of developmental intermediates with divergent fates (such as the OFT and RV cardiomyocytes) suggests that, at early stages of differentiation, transcriptional similarity due to a shared origin may predominate over cell type-specific gene expression.

We next tested the ability of our data to provide more granular detail on cell-type annotation from equivalent populations in other model systems. Recently, Rossi et al. reported the generation of gastruloids from mouse embryonic stem cells that mimic development of cardiac cells and other lineages (Rossi et al., 2021). Using scRNA-seq data from these gastruloids, we subset equivalent populations to those profiled in our data and performed supervised machine learning using a random forest classifier, using our embryo scRNA-seq data as a training dataset (Tan and Cahan, 2019). This allowed us to map the cell types present in cardiac organoids to their most equivalent cell type in the embryo (Fig. 2). We found that pharyngeal and somitic mesoderm, as well as hematoendothelial progenitors in the cardiac organoid, mapped with high accuracy to their equivalent in vivo populations. Cardiomyocytes within the organoids, which had previously been annotated as early and late cardiomyocytes, mapped primarily to aSHF and posterior second heart field (pSHF) cells, as well as the developing SV and AVC. Few, if any, cardiomyocytes mapped to differentiated atrial or ventricular cardiomyocytes, consistent with Rossi et al.'s description of the cardiac regions representing early stages of cardiac development. We applied a similar approach to data obtained from human pluripotent stem cell-derived heart forming organoids (HFOs) (Fig. 2) (Drakhlis et al., 2021). Mesenchymal cell types identified previously by Drakhlis et al. were sub-categorized between the LPM and a/pSHF population, indicating that HFOs recapitulate aspects of heart field progenitor formation. Furthermore, populations identified by Drakhlis et al. as cardiomyocytes mapped mostly to myocardial/OFT intermediates derived from the aSHF (C5), and showed little overlap with bona-fide atrial or ventricular cardiomyocytes (Drakhlis et al., 2021). This underscores that myocardial populations formed in current organoid models represent early stage myocardial intermediates rather than chamber-specific cardiomyocytes, as expected, and demonstrates the potential of our in vivo transcriptional atlas to provide additional detail on cellular heterogeneity and differentiation state of cells in organoid models.

Having identified a number of intermediary cell types in our data, we next aimed to define lineage relationships between cardiac progenitors and their progeny, as well as dynamic changes in genes that inform these relationships. Toward this end, we applied scVelo to estimate RNA velocity, which quantifies the relative abundance of spliced and unspliced transcripts for any given gene in each cell (Bergen et al., 2020). For each gene, the ratio of pre-mature (unspliced) and mature (spliced) mRNA counts is fitted to a model of splicing kinetics constituting a constant transcriptional state, where positive velocity indicates upregulation of a given gene (i.e. a higher abundance of unspliced mRNA than expected in a steady state). From this, we calculated a velocity vector representing the direction and speed of movement of individual cells, revealing several trajectories of differentiation (Fig. 3A). We observed two parallel differentiation streams which began from the pSHF and aSHF (C1 and C21/22, respectively) and terminated in differentiated cardiomyocytes. In addition, we observed diverging differentiation of the aSHF toward the pharyngeal mesoderm and cardiac lineages. Derivatives from the aSHF (e.g. RV and OFT) were clustered near populations of differentiated ventricular cells and Sox2-expressing neural crest cells. These observations are consistent with important interactions between the neural crest and OFT during septation of the aortic and pulmonary arteries (Kirby et al., 1983). Accordingly, RNA velocity also identified a stream of differentiation within the neural crest cells, the terminal position of which diverged and rested with the AVC. This association is consistent with the role of the neural crest in driving endocardial cushion development (Phillips et al., 2013). To understand the temporal ordering of cells along this trajectory, we calculated a latent time score (a measure of the internal clock of the cell based on transcriptional dynamics) and the velocity length and confidence, measurements of the speed of differentiation and congruency of vector fields, respectively (Fig. 3B-E). This captured the transition from cycling low velocity progenitors towards differentiating intermediates, which increased their velocity length as they express differentiated markers, indicating that heterogeneity observed in SHF progenitors is related to differentiation state (Figs 1B and 3D).

To identify the genes underlying the vector field and inferred lineages, we determined the genes with the highest differential velocity as these may represent putative drivers of cell identity or differentiation. We generated a ranked list of the top 100 dynamic genes for each cluster (Table S1). Furthermore, we subset each parallel differentiation trajectory initiated from the pSHF and aSHF and identified the top 300 most dynamically regulated genes across the differentiation stream (Fig. 3F,I). Not unexpectedly, about half of the top dynamically regulated genes were shared across the LPM/pSHF and aSHF differentiation stream; however, 136 genes were unique to each differentiation stream (Tables S2 and S3). In order to visualize the expression patterns of these rapidly changing genes across the pSHF and aSHF, we plotted the expression patterns of the top 300 across latent time, as well as a number of other candidates of interest (Fig. 3F-K; Fig. S6A,B). This allowed us to determine the ordered cascade of these genes across differentiation, and identify at which point these dynamically regulated genes were being expressed. By plotting the splicing dynamics and expression over time of individual candidates from this list of 300 genes, we were able to better understand how these patterns varied for individual populations. We found that genes such as Vsnl1 were most strongly upregulated in cell types found at the earlier/intermediary stages of differentiation (C12, C14) and were later downregulated in atrial/AVC cells (Fig. 3H). In contrast, Ezh2, which was also highly dynamically regulated was highly expressed at earlier stages but became downregulated later during differentiation. This is consistent with previous studies demonstrating a role for Ezh2 in repressing the cardiac progenitor genes in differentiated cardiomyocytes (Delgado-Olguín et al., 2012).

In addition to these examples we identified numerous additional candidates with previously unrevealed cell type-specific regulation across the LPM/pSHF and aSHF differentiation streams. A number of the top dynamically regulated genes identified were specific to one differentiation stream, and had no previously identified role in heart development. These included candidates such as Prtg, which plays roles in neurons and the developing tooth germ (Fig. 3K) (Takahashi et al., 2010; Wong et al., 2010). Prtg was strongly upregulated in differentiating aSHF cells but had low expression in differentiated cardiomyocytes, similar to another dynamically regulated gene Meis2. Meis2 has been implicated in cardiac and facial defects and is essential for cranial and cardiac neural crest development (Louw et al., 2015; Machon et al., 2015), but its expression pattern within the differentiating aSHF has not previously been recognized, and raises the question of whether observed defects in patients are purely due to a neural crest etiology (Fig. 3I). Also dynamically regulated with a cell type-specific expression was Smyd1, which was strongly upregulated in a subset of the differentiated aSHF-derived cells (C5) (Fig. S6D). Smyd1 plays a well-established role in orchestrating early heart development, particularly in the OFT/SHF (Franklin et al., 2016; Rasmussen et al., 2015; Wang et al., 2021), and the selective upregulation of this gene within distinct aSHF derivatives upregulation of this gene within the aSHF derivatives raises interesting questions about whether this occurs in a more granular fashion than previously appreciated. Further studies of these and other candidates identified here may shed new light on the mechanistic role of these genes in a cell type-specific and temporally defined context.

We next studied the transition from pSHF and aSHF cells toward their progeny at the earliest stage of differentiation, which was captured by C14 and C8, respectively, based on RNA velocity length (Fig. 4A). We unexpectedly found that gene set enrichment analysis uncovered that processes involved in ‘left ventricle development’ were differentially regulated in C14 (Fig. 4A), despite its presumed origin from the pSHF. When examining distribution of Foxa2-lineage traced cells we found that EGFP+ cells formed a distinct line that segregated from the atrial cell types and extended from the LPM population to ventricular cardiomyocytes (Fig. 4B). This expression pattern is reminiscent of the distribution of Tbx5+/Hcn4+ cells from the CC that correlated with the FHF, which contributes to the left ventricle (LV) (Fig. 4A). We interrogated the expression of Stard10 and Lbh, markers of early atrial and ventricular identity, respectively, and found Stard10 expression restricted to atrial and SV cells (C4, C12), whereas Lbh was expressed across both differentiation streams, and its expression closely overlapped with EGFP+ cells (Fig. 4B). These data indicate that the stream of rapid differentiation captured by cluster 14 does not simply describe a population of differentiating atrial cells derived from the pSHF. Rather, it likely represents a mix of pSHF derivatives as well as ventricular-fated cells from the LPM that differentiate through similar mechanisms as those employed by the pSHF. This might suggest that intermediary cells at this stage of development are defined more so by their progenitor of origin and its transcriptomic signature during differentiation rather than eventual atrial/ventricular identity.

Pairwise DGE analysis of the parallel trajectories identified differences in various pathways including Wnt, Hh, TGF-β and FGF signaling (Fig. 4C and Table S4). We observed higher expression of canonical Wnt genes such as Wnt2, Sfrp1 and Sfrp5 in pSHF cells (C14), whereas the noncanonical Wnt ligand Wnt5a was highly expressed in differentiating aSHF cells (C8). Differentiation of aSHF progenitors further involved upregulation of Bmp4, as well as Fgf8, which was expressed exclusively in the differentiating aSHF, consistent with previously published data (De Zoysa et al., 2020; High et al., 2009; Pradhan et al., 2017). Although some key components of the cardiac gene regulatory network (GRN) such as Mef2c and Gata6 displayed shared expression between C8 and C14, the differentiating LPM/pSHF shows increased expression of Tbx5 and Tbx20, as well as Gata4 in comparison with differentiation of the aSHF. In contrast, aSHF cells demonstrated increased expression of a number of Irx and Etv family members (Fig. 4C). In this way, our data uncovered both known and new candidates differentially expressed in a/pSHF cells, confirming that differentiation of aSHF and pSHF occurs through activation of shared as well as discrete pathways characteristic for each path.

We next plotted the splicing dynamics and expression over time of the most highly dynamically regulated genes for C8 and C14, reasoning that these might represent genes governing early transition from the LPM/heart field progenitors towards cardiomyocytes. Cap2, a gene implicated in cardiac conduction and sudden cardiac death, was among the top dynamically regulated candidates in C14 (Field et al., 2015; Peche et al., 2013). Cap2 showed rapid upregulation compared with cycling pSHF cells, suggesting that Cap2 may play an important early role in development, consistent with its role as a regulator of actin dynamics (Fig. 4D). Ccdc141 was also among the most highly dynamically regulated genes in C14 and is implicated in heart rate regulation (van den Berg et al., 2017), but no functional studies during heart development have been reported to date. Nebl and Nexn, two genes encoding components of the Z-disc, were among the top five most dynamically regulated genes for C8 and have been associated with familial cardiomyopathies (Fig. 4E). Similar to many candidates identified here, their role during development has not yet been studied; their upregulation at the earliest stage of the differentiation stream suggests that these may play a key role early on. We also identified dynamic regulation of Mest, which is known to be expressed in the trabeculae during development but otherwise has no known function in heart development (King et al., 2002). Mest was downregulated in C14 but rapidly upregulated in the aSHF derivative population (C5), perhaps indicating a transient requirement for this gene, or context-specific roles at different stages of development. These examples highlight interesting patterns not only in gene expression levels, but also in the dynamic changes in the regulation of their expression within individual populations and along differentiation trajectories. This is true both for highly dynamic genes identified using scVelo, as well as well-studied genes involved in cardiac differentiation such as Nkx2-5, the expression of which was positively regulated during differentiation but was negatively regulated following acquisition of myocardial identity later in the differentiation process (Fig. 4F). We observed that Myl7 (an atrial marker in adult cardiomyocytes) was widely expressed in differentiating progenitors irrespective of atrial/ventricular identity but became downregulated at later stages of differentiation in ventricular cells. This underscores a previously established notion that some canonical postnatal atrial/ventricular markers are expressed in a non chamber-specific manner early during development, potentially indicating that they perform different functions during early development compared with later stages or alternatively, highlighting processes that are only transiently required in some lineages but remain key throughout life in others. A more thorough understanding of candidate genes with shifting dynamics, in particular those which are unique to LPM/pSHF and aSHF differentiation trajectories may help to further delineate these populations and their distinct differentiation potential.

Our data indicate that the differentiation streams captured by the clustering and RNA velocity analysis did not directly capture segregation of atrial/ventricular lineages. To understand the mechanisms governing formation of atrial/ventricular lineages, we performed pairwise DGE analysis of differentiated atrial and ventricular cardiomyocytes (C4 and C2/C17, respectively) (Fig. 5A). We identified 3695 differentially expressed genes, including expected markers such as Irx4, Myl2, Mpped2 and Stard10 and markers not previously associated with atrial/ventricular identity (Fig. 5B; Table S5). To predict potential functional consequences of these atrial/ventricular candidate genes, we performed GO and KEGG pathway enrichment analysis (Fig. 5C). We observed upregulation of processes involved in ‘regulation of calcium ion transport’ and ‘mitochondrial respiratory chain complex assembly’, and ‘cardiac myofibril assembly' within ventricular cells, consistent with the increased load on ventricular compared with atrial cells. We also observed differences in a number of pathways including Hh, JAK/STAT and Ras signaling, as well as a number of metabolic processes such as ‘thyroid hormone synthesis’ and ‘steroid biosynthesis’.

To investigate mechanisms underlying atrial/ventricular progenitor segregation at the earliest stages, before expression of chamber-specific genes, we examined the CC stage alone, which was composed of multiple progenitors including the LPM, aSHF and the pSHF, as well as a transcriptionally distinct Hcn4+/Tbx5+ population representing the FHF (Fig. 5D-F). DGE and GO enrichment analysis between EGFP+ and EGFP– cells identified 127 differentially expressed genes (Fig. 5G; Table S6) and uncovered processes involved in ‘lipoprotein biosynthesis and remodeling’ (Fig. 5H). A number of lipoproteins such as Apoe and Apom were upregulated in the EGFP+ cells. We had also observed genes such as Apobec2, a member of the apoB mRNA editing enzyme, to be among the top differentially expressed genes between ventricular and atrial cells (Fig. 5B; Table S5). The spatiotemporal distribution of these lipoproteins at early stages of cardiac development has not been characterized extensively, opening the question as to whether metabolic differences may be implicated in the earliest segregation of atrial and ventricular cell types during development, potentially as early as the CC.

In summary, we identified candidate genes of atrial/ventricular identity within differentiated cardiomyocytes and early progenitors, identifying a putative class of molecules that may characterize the ventricular lineage before chamber formation. Alternatively, candidates enriched in Foxa2 lineage-traced cells at the CC stage may represent components relevant in early differentiating cells, as opposed to later differentiating cells, equally offering interesting opportunities to further dissect mechanisms of early cardiac specification. Experimental validation of the expression pattern and loss-of-function experiments of such candidates will determine their spatiotemporal distribution within the developing heart and their relevance for early heart specification.

Previous studies have shown that an anterior-posterior gradient of RA directs differentiation of cardiac precursors towards an atrial fate (Bernheim and Meilhac, 2020; Perl and Waxman, 2019, 2020; Waxman and Yelon, 2009). Accordingly, loss of the RA synthesis enzyme Raldh2 (Aldh1a2) in the heart leads to loss of the atrial chambers (Niederreither et al., 2001). Exposure to exogenous RA in utero in turn is associated with a number of cardiac and facial defects, indicating that retinoids can act as teratogens (Piersma et al., 2017). Despite the well-known importance for RA signaling during development, the comprehensive transcriptional effects induced by exogenous exposure to RA are incompletely understood in the context of chamber specification.

We and others have shown that manipulations in RA signaling during cardiac development lead to defects in FHF/SHF morphology as well as dysregulation of key components of the cardiac GRN (Bardot et al., 2017a; Bertrand et al., 2011; De Bono et al., 2018; Hochgreb et al., 2003; Lin et al., 2010; Rana et al., 2014; Stefanovic and Zaffran, 2017; Xavier-Neto et al., 1999). It is well established that the effects of RA are dose-dependent, and that different concentrations can even have opposite effects (Perl and Waxman, 2019). To test the effect of RA levels we exposed embryos to different concentrations of RA at E7.25 (Fig. 6A,B). As expected, exogenous RA caused patterning defects across the entire embryo in a dose-dependent manner (Fig. 6B). Intraperitoneal injection of 65 mg/kg RA caused severe defects in the head and pharyngeal structures, along with stunted formation of caudal structures. Exposure to 16.25 mg/kg RA caused milder defects, with grossly intact development of the tail and limb bud and milder, yet nonetheless hypomorphic, effects on head and pharyngeal structures (Fig. 6B). We also observed dysregulation in heart chamber size and overall reduced embryo size. Low concentrations of RA (3.25 mg/kg) led to decreased ventricular size at E10.5 and mild craniofacial abnormalities, leading us to conclude this was more suitable for understanding selective effects of exogenous RA on early heart development.

Embryos from pregnant dams injected with 3.25 mg/kg RA intraperitoneally at E7.25 were collected at equivalent stages as in control embryos (Fig. 6C). We merged samples from CC-, PHT- and HT-stage embryos across control and RA-injected embryos and annotated cluster identity and sample contribution as previously, establishing similar UMAP projections compared with the control samples of the earlier analysis, consistent with the mild morphological changes observed (Fig. 6D-G; Fig. S7A,B).

Consistent with our morphological observations the frequency of cells contributing to ventricular cardiomyocytes (C3, C24) was decreased in RA-injected embryos relative to control (Fig. 6E). We thus first asked whether exogenous RA affected the contribution of cells to the aSHF or pSHF populations. We observed no differences in the relative proportion of cells from RA-injected embryos to aSHF, pSHF and LPM populations (C8/C17, C5/C7 and C12, respectively) (Fig. 6E). We next asked whether exogenous RA influenced the allocation of ventricular-fated progenitors at early stages of cardiogenesis. Using a quantitative whole-mount immunofluorescence approach (Bardot et al., 2017b), we quantified the percentage of the Nkx2-5+ CC that was composed of Foxa2 lineage-traced cells in control and RA-injected embryos at E8.25 (Fig. 6H; Fig. S7E). We found no significant differences between both groups, indicating that changes in ventricular chamber size were not due to a shift in chamber-specific fate of progenitors and instead suggested that the observed defects might originate in aberrant differentiation of progenitors.

In order to identify trajectories of differentiation, we visualized RNA velocity across the various populations from both control and RA-injected embryos (Fig. 7A), noting dysconnectivity between the aSHF cells (C8, C17) and the differentiating aSHF derivatives (C10), which cluster adjacent to ventricular cells (Fig. 6D). This prompted us to ask whether defects in differentiation of aSHF cells was affected by exogenous RA. A neighboring mesenchymal cell cluster (C6) was strongly over-represented in RA embryos, suggesting that exogenous RA may cause defects in the normal lineage trajectory of aSHF progenitors resulting in a shunt of these cells toward a pharyngeal/mesenchymal cell type, consistent with the multipotent nature of the aSHF (Fig. 7A).

This led us to interrogate transcriptional changes in the aSHF as a result of RA exposure, finding that exogenous RA caused a downregulation in processes involved in ‘myotube development’, ‘ventricular cardiac muscle development’ and ‘cardiac myofibril assembly’ (Fig. 7B). These were accompanied by an increase in processes involved in ‘head morphogenesis’ and ‘head development’, consistent with the subtle head deformities observed previously. We further found increased negative regulation of FGF signaling, which is known to be required for differentiation of aSHF cells towards cardiomyocytes (Pradhan et al., 2017; de Soysa et al., 2019). Consistent with this, aSHF cells from RA-exposed embryos had decreased expression of cardiomyocyte markers such as Myl2 and Nkx2-5, with downregulation of Bmp4, which plays a key role in differentiation of the aSHF progenitors. These changes were accompanied by increases in Pbx1 and Sox9, both genes involved in development of the pharyngeal lineages (Cheung and Briscoe, 2003; Welsh et al., 2018). These data suggest that exogenous RA causes defects in different signaling pathways that subsequently may result in a block in differentiation of aSHF towards the ventricular lineage and a shunt towards other lineages within its multipotent capacity, such as the pharyngeal mesoderm.

In order to better visualize this shunt in differentiation, we subset the aSHF populations as well as the ventricular cardiomyocytes and mesenchymal/pharyngeal mesoderm populations and reclustered them (Fig. 7C,D). We observed diverging fates from the aSHF toward the ventricular cardiomyocytes and the pharyngeal/mesenchymal cell types (Fig. 7E). In order to better visualize these transition states between populations, we constructed trajectories in URD, using the aSHF as a starting population and the mesenchymal/ventricular cardiomyocytes/pharyngeal mesoderm clusters as the terminal population. We observed that cells from RA-treated embryos had impaired differentiation toward a Tnnt2+ ventricular cardiomyocyte identity at later time points, and were over-represented in the differentiation arm towards mesenchymal cell types, suggesting again that RA-treatment had caused impaired differentiation toward cardiomyocyte identity and instead resulted in the formation of mesenchymal/pharyngeal cell types (Fig. 7F-H).

When examining the effect of RA on other relevant populations, we discovered increased expression of processes involved in biosynthetic machinery such as ‘tRNA metabolism’ and ‘amino acid synthesis’ within early cardiomyocyte progenitors (C1) (Fig. 7I). This was coupled with decreased sarcomere organization and differentiation toward cardiomyocytes accompanied by dysregulation in Wnt signaling, potentially suggesting that these cells were arrested during differentiation towards a cardiac lineage. This was further supported by RNA velocity projections, which demonstrated cycling cardiomyocyte progenitors that may be inhibited in their ability to contribute further to differentiated cardiomyocytes (Fig. 7A). The progenitors in RA-exposed embryos demonstrated decreased Nkx2-5 and Tnnt2 expression but interestingly had increased expression of Hopx, perhaps indicative of an arrest of differentiation at the cardiomyoblast stage that could not be overcome to form functional cardiomyocytes (Jain et al., 2015).

Finally, we investigated changes in RA-exposed differentiated ventricular cardiomyocytes, where we uncovered decreases in processes involved in ‘striated muscle development’ and ‘ventricular muscle morphogenesis’, but also found that RA caused defects in metabolic and energetic components such as the respiratory electron transport chain (Fig. 7J). Other genes that were downregulated in the RA condition included sarcomeric genes and cardiac transcription factors such as Nkx2-5 and the Foxa2 lineage-tracing marker EGFP, as well as ventricular markers such as Ckb (Fig. 7J). Collectively, this indicated that not only were ventricular cells present at lower frequency due to the improper differentiation of progenitors, but those that did exist may also have decreased fitness due to perturbations in their maturation and metabolism.

Here, we have generated a single cell transcriptomic atlas to interrogate dynamic changes of differentiation over time, building on existing work defining the heterogeneity of heart field progenitors (Ivanovitch et al., 2021; de Soysa et al., 2019; Tyser et al., 2021; Zhang et al., 2021). We uncovered broadly similar cell populations and investigated in depth their relationship to each other through a variety of computational lineage trajectory tools. We show that FHF cells integrated with differentiating cardiomyocytes when analyzed across real developmental time, adding to a growing understanding of the FHF representing a differentiating early cardiomyocyte population that is altogether different from its SHF progenitor counterpart, which undergoes replication and retains multipotent capacity. We were intrigued to find clusters of SHF progenitors at multiple stages of differentiation, including their ‘stable’ replicative state as well as intermediary states during differentiation. Visualization of gene expression across differentiation time uncovered comprehensive information about the temporal ordering of genes involved in cardiac differentiation and the dynamics of their expression in a cell type-specific manner. This identified several candidates previously described in genetic mutants – for example, the upregulation of Fgf8 and Bmp4 in differentiation of the aSHF toward the OFT/RV (De Zoysa et al., 2020; Hutson et al., 2010; Ilagan et al., 2006; McCulley et al., 2008; Park et al., 2006; Watanabe et al., 2012) – but also uncovered new concepts such as the differences in canonical/non-canonical Wnt signaling across aSHF and LPM/pSHF differentiation. Although the importance of Wnt signaling in anterior/posterior patterning and heart formation has long been recognized (Ai et al., 2007; Gessert and Kühl, 2010; Marvin et al., 2001; Sinha et al., 2015; Yamaguchi, 2001), we have shown here that different components of Wnt signaling are used in distinct patterns during differentiation from the arterial/venous poles of the heart. Our analysis also showed that differentiation of the pSHF and aSHF progenitors involves deployment of separate but sometimes overlapping components of the cardiac GRN such as Tbx5/Tbx20 and Robo/Slit (Bruneau, 2013; Paige et al., 2015; Srivastava, 2006). This raises questions about how different cardiac progenitors might use discrete components of the GRN during differentiation and what upstream mechanisms are driving these differences.

We identified the genes with highest velocity for each cluster and across the course of differentiation from aSHF and pSHF/LPM progenitors to cardiomyocytes. The highly dynamic regulation of these genes suggests that some members of this list may act as putative drivers of differentiation, though follow-up studies would be needed to determine this experimentally. These included a number of candidates for which the role in cardiac development has not previously been shown, such as Prtg, or candidates that have been implicated in cardiac disease but have an unclear role in heart development, such as Cap2 and NebI/Nexn. Importantly, our data revealed not only the ordered cascade of these dynamically regulated genes and other known markers across cardiac differentiation, it also showed how the dynamic regulation of these genes can occur in a cell type-specific context and in intermediate populations during development. This information should be used to inform the design of functional follow-up studies of these candidates, for example in the selection of genes for conditional deletion/knock out studies, or the stages that should be profiled to capture any potential phenotypes.

We further characterized the atrial and ventricular cardiomyocyte populations and identified putative candidates of atrial/ventricular fate and identity. Use of Foxa2 lineage tracing identified genes and processes divergent at the earliest stages of differentiation toward a ventricular lineage, including a number of lipoproteins such as Apom and Apoe and mRNA binding proteins such as Apobec2. These have not been studied in depth during cardiac development but play roles in myoblast differentiation in skeletal muscle and regulation of left-right axis specification in other systems (Ohtsubo et al., 2017; Vonica et al., 2011). The differences in metabolic components between atrial and ventricular cells may be related to expected differences in hemodynamic load between these chambers and the increased energetic demands on the ventricle as it continues to develop, and mirrors changes in metabolism that drive maturation of ventricular cells at later stages of development (Chattergoon, 2019; Murphy et al., 2021; Wickramasinghe et al., 2022). Alternatively, these may represent genes involved in early differentiation, as our previous work has indicated that isolation of Foxa2 lineage-traced cells during migration of cardiac mesoderm can identify early migrating cells that experience different signaling environments than their posterior counterparts (Bardot et al., 2020 preprint).

In-depth single cell data from developmentally relevant stages promise to provide reliable and easy-to-compare-to blueprints for newly emerging in vitro models. When compared with our embryo data, organoid-derived cardiomyocytes mapped across several progenitor subtypes but primarily overlapped with intermediate cell types of differentiating pSHF/aSHF as opposed to chamber-specific cardiomyocytes. These data highlight the utility of organoids for capturing transient differentiation events during early development, which are thought to be particularly affected in CHDs (Bruneau, 2008, 2020; Kloesel et al., 2016). It remains unclear, however, whether the mapping to discrete progenitors also informs the distribution of more differentiated cardiomyocytes that may arise at later stages in culture. This analysis will be important for understanding both the limitations of these models as well as providing strategies for their further development. In addition, this underscores the power of cross-model comparisons for a deeper understanding of these questions, and should be further explored across data from other animal models (Asp et al., 2019; Cao et al., 2020; Cui et al., 2019; Holowiecki et al., 2020; Mantri et al., 2021; Wang et al., 2019).

The use of scRNA-seq data to understand cell type-specific effects of genetic defects has been explored previously (Kathiriya et al., 2021; de Soysa et al., 2019). We extended this to an in utero exposure model of RA, as RA has been well-characterized as a teratogen (Lammer et al., 1985; Piersma et al., 2017). The hypomorphic ventricular phenotype observed in RA-exposed embryos is consistent with classical studies in the chick (Osmond et al., 1991; Yutzey et al., 1994), zebrafish (Stainier and Fishman, 1992) and mouse (Xavier-Neto et al., 1999), demonstrating that exogenous RA results in changes to atrial and ventricular chamber size in a concentration-dependent fashion (Bernheim and Meilhac, 2020; Perl and Waxman, 2019; Waxman and Yelon, 2009). Further studies have shown that depletion of RA receptors results in a low increase in RA concentration and an increase in the number of atrial cardiomyocytes, and combining a Cyp26a1 depletion with it results in an intermediate increase in RA concentration, with reduced ventricular cardiomyocytes and no effect on atrial cells (D'Aniello et al., 2013). The work presented here builds on these studies and aims to determine the transcriptional effects induced directly by RA signaling in a cell type-specific manner (Hochgreb et al., 2003; Lin et al., 2010; Perl and Waxman, 2020). We demonstrate that in utero exposure to RA at E7.5 did not affect the distribution of heart field progenitors between control and RA-exposed conditions, which was interesting as previous work has shown a role for RA in restriction of the cardiac progenitor pool (Duong et al., 2021; Keegan et al., 2004, 2005). This may be because we elected to inject at a time point following migration of the LPM, in an effort to focus on the effect of RA on cardiac mesodermal cell differentiation and expansion. Future work should examine whether earlier injection of RA at similar concentrations might recapitulate previously observed defects in progenitor size, and determine whether this occurs in an atrial/ventricular-specific manner.

Our observation that aSHF differentiation towards a myocardial fate was negatively affected by RA exposure is consistent with previous work demonstrating posterior expansion of Fgf8+ and Isl1+ SHF progenitors and defects in OFT septation in RA-deficient mouse mutants. Antagonism of Fgf8 and RA has been well-characterized in the literature in various contexts, placing RA at the top of an important signaling cascade driving growth and differentiation of the developing heart (del Corral et al., 2003; Cunningham et al., 2013; Pasini et al., 2012). Interestingly, recent work has indicated that RA-deficient zebrafish embryos display defects in differentiation of SHF-derived cardiomyocytes. It is unclear whether this represents differences across species, or perhaps concentration differences between KO studies and our low-dose RA exposure model. The increase in processes involved in head morphogenesis is also consistent with the role of RA as a teratogen, and a cause for craniofacial malformations and microcephaly (Petrelli et al., 2019). Future work should determine whether an RA concentration-dependent mechanism fine tunes control of SHF progenitor expansion/differentiation across various lineages, or whether the increase in pharyngeal mesoderm specification is an otherwise aberrant differentiation event driven by restriction of cardiomyocyte fate.

Foxa2Cre mice (C57BL/6J) mice were shared with us by Dr Heiko Lickert (Uetzmann et al., 2008) and maintained on a mixed background. Rosa26-mTmG mice were obtained from the Jackson Laboratory. All murine embryo scRNA-seq studies were performed on Foxa2Cre;R26R:mTmG embryos derived from crossing homozygous Foxa2Cre males with homozygous Rosa26-mTmG reporter females between 2-3 months of age. For timed matings, the day of plug identification was considered to be E0.5. All animals were housed in the Center for Comparative Medicine and Surgery (CCMS) facilities at Icahn School of Medicine and all experiments were conducted in accordance with the guidelines and approval of the Institutional Animal Care and Use Committee at Icahn School of Medicine at Mount Sinai. For studies on the effect of RA, pregnant mTmG females previously crossed with Foxa2Cre males were injected intraperitoneally at E7.25 with 65 mg/kg, 16.25 mg/kg or 3.25 mg/kg of RA dissolved in corn oil. Single cell analysis on RA-exposed embryos was performed following injection of 3.25 mg/kg RA.

Foxa2Cre;mTmG embryos were dissected at E8.25, E8.75 and E9.25. The cardiac region was sub-dissected at each stage (Fig. S1), including the surrounding endoderm, head folds and pharyngeal structures. This was carried out in order to ensure complete dissection of the cardiac structures and also to include transcriptional information within key tissues that co-develop with the cardiac structures at each stage. Embryos were dissected in cold PBS with 0.1% bovine serum albumin (BSA) and kept on ice until dissociation. Three stage-matched littermates of unknown sex were incubated separately in 200 µl of 0.25% Trypsin-EDTA in a 37°C water bath for 5 min or until the tissue dissociated completely, vortexing once halfway. Trypsin solution was quenched with 1 ml DMEM with 10% fetal bovine serum (FBS), centrifuged at 300 g for 5 min, and resuspended in 50 µl PBS with 0.1% BSA for scRNA-seq preparation.

Quantification of the lineage-traced Foxa2Cre:YFP contributions to the CC was performed using previously published methods (Bardot et al., 2017b). Briefly, whole E8.25 embryos were stained for GFP and Nkx2-5 and mounted for confocal microscopy [Nkx2-5 antibody: Santa Cruz Biotechnology, sc8697, 1:100 (IF); GFP antibody: Abcam, ab13970, 1:500]. Masks were generated for the Nkx2-5 and YFP regions in Imaris Viewer and volumes were calculated for each region. We then quantified the percentage of the Nkx2-5+ volume that was occupied by YFP+ cells and compared these across control and RA-injected embryos.

Individual samples for each time point were hashed using the Chromium platform in order to determine that downstream analysis was not affected by biological variability between embryos. After hashing, libraries were generated using the Chromium platform with the 3′ gene expression (3′ GEX) V3 kit, using an input of ∼10,000 cells. Gel-Bead in Emulsions (GEMs) were generated on the sample chip in the Chromium controller. Barcoded cDNA was extracted from the GEMs using Post-GEM RT-cleanup and amplified for 12 cycles. Amplified cDNA was then fragmented and subjected to end-repair, poly A-tailing, adapter ligation and 10x-specific sample indexing following the manufacturer's protocol. Libraries were quantified using Bioanalyzer and QuBit analysis. Libraries were sequenced in paired end mode on a NovaSeq instrument targeting a depth of 50,000 reads per cell.

Sequencing data were aligned and quantified using the Cell Ranger Single-Cell Software Suite (version 3.0, 10x Genomics) against the provided GRCm38 (mm10) mouse reference genome. Lineage tracing information was determined by further alignment to a custom reference for EGFP using provided maps of constructs used in creation of mTmG mice from Jackson Laboratory. Complementary Tomato expression information could not be obtained owing to increased distance from the polyA tail compared with EGFP.

Downstream differential expression and clustering analysis was performed using the Seurat V.4.0 package, as described in the tutorials (http://satijalab.org/seruat/). cellRanger matrices were imported for each sample, and distributions of nCountRNA, nFeatures and percentage mitochondrial gene expression were examined for each sample in order to filter out doublets or cells of low quality. Cells with greater than 20% of genes coming from mitochondrial genes were selected against, as well as those with fewer than 200 genes.

We then normalized the resulting subset Seurat objects using the scTransform workflow and further scaled and normalized the RNA assay in order to perform downstream differential expression analysis. We performed principal component analysis using the highly variable genes for each sample. The most significant principal components (20-30 depending on sample) were used for graph-based unsupervised clustering (FindClusters and FindNeighbor functions). Uniform manifold approximation and projection (UMAP) was performed using a standard resolution parameter of 0.5 and iteratively modified after performing marker gene expression and examining expression of key markers.

In order to generate a combined Seurat object encompassing all three time points, we visualized expression of cardiac and mesoderm genes (Pdgfra, Nkx2-5), endothelial/endocardial (Pecam1, Nfatc1), endoderm (Epcam, Sox17) and ectoderm (Sox2, Pou3f1). We subset clusters with positive expression of mesoderm/cardiac markers, and more completely annotated clusters that did not express these markers but clustered closely through GO/KEGG analysis and DGE analysis. These represented related lineages such as the neural crest or early/related mesoderm progenitors. These were subset together with the cardiac lineages for each stage.

Seurat objects for the CC, PHT and HT stages were labeled with their original sample ID and merged together using the scTransform-based integration workflow (Hafemeister and Satija, 2019). For analyses involving comparison of control and RA-injected embryos, the above steps were repeated by combining objects from both conditions either for corresponding stages or across all samples. Downstream clustering and dimensionality reduction was performed using similar methods described above. Cell cCycle quantification and annotation was performed on merged samples using the CellCycleScoring vignette available from the Satija lab (https://satijalab.org/). Quantification of time point, condition or sample contribution to clusters was determined by normalizing the frequency of a given metadata classification to the total cell number from each sample in order to control for cell number differences across samples.

In most cases, differential expression was determined using the FindMarkers function on the ‘RNA’ assay, using a negative binomial model with parameters for min.pct=0.2 and a P-value cut off of <0.01 either by comparing a single cluster with all others for annotation, or in pairwise fashion for clusters of interest. Differential expression of Foxa2-lineage-traced positive and negative cells was analyzed by classifying each cell with a binary variable based on expression of EGFP>1.0 and then performing differential expression on subsets of cells across the EGFP+/− metadata information. A similar strategy was employed to query differentially expressed genes from embryos post RA injection on a cluster-by-cluster basis across control and RA condition metadata information.

Differential gene expression was performed in Seurat through the FindMarkers function in Seurat using a negative binomial test using a P-value cutoff of <0.05. Gene set enrichment for GO terms and KEGG pathways were performed using the ClusterProfiler tool (https://guangchuangyu.github.io/software/clusterProfiler/) using a cutoff P-value of <0.05.

RNA velocity was calculated using Velocyto and scVelo according to the developer's manual (Bergen et al., 2020; La Manno et al., 2018). Briefly, .loom files containing exon and intron information of all samples were created using the Velocyto python package ‘velocyto run10x()’ command and merged using ‘loompy.combine()’ in the loompy package. Using the scVelo python package ‘scv.pp.filter_and_normalize()’ and ‘scv.pp.moments()’ functions, the data was normalized and moments computed. Subsequently the ‘scv.tl.velocity(stochastic)’ and ‘scv.tl.velocity_graph()’ functions were used to estimate RNA velocity, which was then projected onto Seurat-derived UMAP coordinates in conjunction with the ‘scv.pl.velocity_embedding_stream()’ function. To analyze the dynamics of individual genes of interest, phase portraits, gene-specific RNA velocity UMAP projections and expression plots were created through ‘scv.pl.velocity()’. The length and coherence of the velocity vectors, which indicate differentiation speed and directional confidence, respectively, were calculated using ‘scv.tl.velocity_confidence()’. To visualize individual cell connections in more detail, high confidence transitions were plotted through ‘scv.pl.velocity_graph()’ with a threshold of 0.5. A dynamical model was applied to analyze transcriptional states and cell-internal latent time by running ‘scv.tl.recover_dynamics()’ and subsequently recomputing RNA velocities [‘scv.tl.velocity(dynamic)’ and ‘scv.tl.velocity_graph()’ functions]. The latent time, based solely on the transcriptional dynamics of a cell, was then determined using ‘scv.tl.latent_time()’ and plotted with ‘scv.pl.scatter()’. Finally, the dynamical model also uncovered putative driver genes, as these show increased dynamic behavior [‘scv.tl.rank_dynamical_genes()’].

To reconstruct branching trajectory trees of the aSHF subset (Fig. 7C) in the control/RA merged dataset we used URD (v1.1.0). Seurat information was imported into URD, including relevant count matrices and metadata information using custom functions. We calculated a diffusion map using the calcDM function from URD with knn=200 and sigma.use=10. As the root, we assigned both aSHF cell clusters from the original dataset (C8, C17). Cells were then ordered in pseudotime by simulating diffusion from the root to calculate the distance of each cell from the root. For this, we used the floodPseudotime function with n=10 (number of simulations) and minimum.cells.flooded=2. In total, 200 simulations were performed. Mesenchymal cells, ventricular cardiomyocytes and pharyngeal mesoderm populations were defined as tips for the RNA full tree, using the RNA velocity information to guide choice of endpoints between related populations. To apply URD, we used pseudotimeWeightTransitionMatrix with parameters optimal.cells.forward=40 and max.cells.back=80 to determine the slope and inflection point of the logistic function used to bias the transition probabilities. We simulated random walks on the cell–cell graph from each tip to the root using processRandomWalks function from URD. In total, 35,000 random walks were performed per tip for the RNA full tree. Finally, trees were built using buildTree function, which starts from each tip and joins trajectories that visited the same cells and compares all predefined tips in a pairwise manner. The following parameters were used to build the tree: bins.per.pseudotime.window=8, cells.per.pseudotime.bin=80, divergence.method=‘preference’, p.thresh=0.0.

Machine learning was used to compare our dataset to others by training a random forest classifier (singleCellNet R package) on fully cell type-annotated external data (Rossi et al., 2021) and assessing it on a withheld subset (Tan and Cahan, 2019). The classifier was then applied to the current dataset and a matching cell type was predicted for each cell. In turn, a classifier was trained on our data and applied to the external datasets.

We thank Bhavana Shewale and Vera Soeder for critical feedback on the manuscript. We thank the Icahn School of Medicine at Mount Sinai (ISMMS) Microscopy Core (Shilpa Dilipkumar, Katarzyna Cialowicz, Nikolaos Tzavaras, Glenn Doherty), the ISMMS Genomics Core (Sanjana Shroff, Mondale Wang) and the ISMMS Center for Comparative Medicine and Surgery for their technical assistance. Evan Bardot provided essential training and helped establish the protocol to isolate cardiac regions from early embryos. Will Zhao provided assistance with aspects of scRNA-seq data analysis in Seurat.

The peer review history is available online at https://journals.biologists.com/dev/article-lookup/doi/10.1242/dev.200557