Introduction

Heart development requires coordinated specification of multiple lineages, each characterized by serial cell fate determination events. Identifying the developmental trajectories is the key to understanding heart formation (Buckingham et al. 2005; Kelly 2023). In past decades, a stepwise determination of early cardiac lineage hierarchy has been primarily established (Parameswaran and Tam 1995; Abu-Issa and Kirby 2007). The current model designates cardiac progenitor cells into discrete pools, including the first and second heart field (FHF and SHF) (Buckingham et al. 2005), and a newly classified juxta-cardiac field (JCF) (Tyser et al. 2021). FHF mainly contributes to the left ventricle (LV) and the atria. SHF, located in a dorsal-medial region to FHF, progressively develops into cells in the right ventricle (RV), the outflow tract (OFT) and the atria (Cai et al. 2003). Interestingly, JCF contributes to not only epicardium and pericardium, but also cardiomyocytes (CMs) of LV and the atria (Tyser et al. 2021; Zhang et al. 2021).

Cardiac progenitors are mostly generated from the Mesp1-expressing (Mesp1⁺) nascent mesoderm (NM) cells during gastrulation (Saga et al. 1999; Devine et al. 2014; Lescroart et al. 2014). Pioneering work has revealed that the developmental capacities of each cardiac progenitor pool are highly related to the spatial-temporal constriction during the specification of NM cells (Lescroart et al. 2018; Ivanovitch et al. 2021; Zhang et al. 2021). Temporally inducible lineage tracing indicates that E6.5 Mesp1⁺ cells mostly contribute to LV, whereas E7.25 Mesp1⁺ cells give rise to RV, atria, OFT, and inflow tracts (IFT) (Lescroart et al. 2018). It seems that early-streak stage NM cells differentiate into FHF pools, while late-streak stage NM cells relate to SHF progenitors. Interestingly, JCF population is also derived from the Mesp1⁺ NM cells in the gastrula (Zhang et al. 2021). Recent studies on single-cell transcriptomic data of the late headfold stage embryos have revealed that JCF shares a number of molecular markers with FHF, but lacks Nkx2.5 expression and exhibits specific Mab21l2 expression (Tyser et al. 2021). However, unlike FHF cells, JCF cells are largely located at the embryonic-extraembryonic mesodermal interface, as revealed by the Mab21l2 expression, rostrally to the Nkx2-5 positive cardiac crescent region (Tyser et al. 2021). It remains unclear the molecular signaling underlying the early specification of NM cells into JCF and FHF population.

In this study, by bridging the transcriptional landscapes between the gastrula and the headfold stage in early mouse embryos, we systematically depict the transcriptional trajectories leading to CMs during early mouse development, and characterize the epigenetic landscapes that underlie early mesodermal lineage specification. The analyses reveal two dinstinct developmental trajectories towards CMs, namely the MJH trajectory (the Hand1-expressing early extraembryonic mesoderm – JCF and FHF – CM) and the PSH trajectory (the pharyngeal mesoderm cells – SHF – CM). Through single-cell multi-omics analysis, the core regulatory network (CRN) in MJH is identified, consisting of transcription factors (TFs) GATA4, TEAD4, HAND1 and FOXF1. Further functional analysis indicates that HAND1 and FOXF1 are activated sequentially, and both required for mesodermal specification and the expression of the MJH specific genes. Taken together, our study unveils the transcriptional and epigenetic dynamics during early cardiac specification, demonstrates the crucial roles of HAND1 and FOXF1 in driving early cardiac specification, and provides insights for the investigation of congenital heart defects.

Results

Transcriptional dynamics during the specification of early cardiac progenitors

To delineate the origin of cardiac progenitors, we constructed the E6.5-8.5 developmental trajectories using the published mouse gastrulation cell atlas by performing the Waddington-Optimal-Transport (WOT) analysis to infer ancestor-descendant fates of cells (Schiebinger et al. 2019), (Pijuan-Sala et al. 2019). CMs were used as trajectory endpoints and traced back to the E6.5 epiblast. Clear trajectory separation was observed within E7.5-7.75 (Figure 1A). Besides the pharyngeal mesoderm (PM) cell cluster at E7.75, a subset of E7.5 early extraembryonic mesoderm (EEM) cells (Zhang et al. 2021) was specifically identified in a distinct developmental trajectory. PM population was marked by the expression of Isl1, Sfrp1, Tcf21, Tbx1 and Irx3, suggesting its relationship with the SHF progenitors; the EEM cells exhibited highly expressing Hand1, Pmp22, Foxf1 and Spin2c (Figure S1A). In addition, cells along the EEM trajectory also expressed higher levels of the FHF signature genes, including Tbx5 and Hcn4, suggesting their contribution to the FHF (Zhang et al. 2021). Compared with the PM trajectory, the EEM trajectory was mainly composed of cells in the later stages of cardiac development (Figure 1B and S1B). By the heart looping stage (E8.5), the two trajectories indeed exhibited distinct contributions to cardiac structures LV and OFT, respectively (Figure 1C).

Figure 1:

Thus, we here refer to these two transcriptional trajectories as EEM-JCF/FHF (MJH) and PM-SHF (PSH) (Figure 1A). The spatial correlation of these two trajectories with heart fields was confirmed using the dataset from manually micro-dissected mesodermal cells in cardiac regions of E7.75-8.25 mouse embryos (Figure S1C) (Tyser et al. 2021). The MJH and PSH trajectories contained overlapping but temporally distinct cell types (Figure 1B). The separation of the two trajectories can be observed as early as E7.0 (Figure 1B). NM cells, marked by Mesp1 and Lefty2, at E7.0 were more likely to be the multipotent progenitors of the PSH trajectory; whereas mixed mesoderm (MM) cells, marked by Hand1 and Msx2, in later developmental stage and EEM cells tended to contribute to the MJH trajectory (Figures 1D and 1E).

We then systematically characterized the stage specific marker genes in the MJH and PSH trajectories. These two trajectories exhibited discrete and dynamic gene expression profiles during development (Figure 1F). We also observed relatively similar inhibitory Wnt and Nodal, as well as active Fgf and Notch, signaling activities along the two trajectories (Figures 1G and S2A). Interestingly, the early stage of the MJH trajectory seems show higher Bmp and Yap signaling activities (Figure 1G). Temporal expression profiles of the Bmp genes indicated that Bmp4/5/7 were dynamically expressed during cardiac specification, with Bmp4 demonstrating higher MJH specificity and at least 0.5 days earlier activation (Figure S2B). Geo-seq data analysis indicated that Bmp4 was highly specific to mesoderm, and enriched at the proximal mesodermal ends (layer 11 at E7.0, layer 9-10 at E7.5) with distinct anterior-posterior preference at E7.5 (Figure S2C). For the target genes of Bmp signaling, several genes (Hand1, Car4, Arl4c and Pmp22) showed MJH specific activation-to-repression dynamics, similar to Bmp4 (Figures 1F and S2D).

Epigenetic signatures of the early MJH and PSH progenitors

To investigate the epigenetic regulation during cardiac cell fate decisions, we performed multi-omic analysis by combining single-nucleus RNA-sequencing (snRNA-seq) and scATAC-seq to generate paired, cell type specific transcriptome and chromatin accessibility profiles of 13,226 cells in E7.0 mouse embryos. The single-cell transcriptomic data were integrated with the published E7.0 mouse embryo cell atlas data, followed by label transfer and gene expression-based cell type identification (Figures S3A and S3B). For the scATAC-seq data, we scored the genome-wide ATAC activities with bin sizes of 10 kb prior to UMAP analysis, which yielded cell clusters similar to transcriptome-based analysis (Figure S3A).

9 clusters of cells were identified through clustering analysis at both transcriptional and epigenetic levels, which are NM cells (Clusters 0, 1 and 2), primordia germ cells (PGC, Cluster 4), hematoendothelial progenitors (Haem, Clusters 5 and 6), and EEM cells (Clusters 3, 7 and 8) (Figures 2A, 2B and S3C). RNA velocity also supported the four possible trajectories mentioned above for the earliest mesodermal cell specification (Figure 2C). WOT analysis revealed that Clusters 3, 7, and 8 showed intermediate to high probabilities of belonging to the MJH trajectory; pseudotime analysis indicated that Cluster 8 represented the late differentiated EEM populations (Figure S4A). Although Cluster 2 represented the relatively late stage of NM cell population by pseudotime analysis, Cluster 1, 0 and 2 demonstrated similar probabilities of belonging to the PSH trajectory (Figure S4A). Thus, the analyses further indicated that EEM cells and hematoendothelial progenitors were clearly separated from the mesoderm, while the PSH trajectory related cells still remained at the early NM stage at E7.0.

Figure 2:

We further analyzed the chromatin accessibilities of these cell clusters. Total 90,661 chromatin accessible elements (CAEs) were detected, 7,206 of which were differentially accessible elements (DAEs) across the 9 clusters (Figure 2d). The DAEs were annotated to their target genes by enhancer-promoter (EP) pairing analysis. Consistent with the clustering analysis based on gene expression (Figure S3C), DAVID functional term analysis revealed that the DAE target genes in Clusters 3, 7 and 8, such as Hand1, Foxf1, Bmp4 and Msx1, were mainly associated with heart morphogenesis, that those in Clusters 5 and 6, like Tal1, Lmo2 and Fli1, were related to angiogenesis and vasculature development, and that those in Clusters 0, 1, 2, for examples, T, Zic2, Lhx1 and Gata6, were enriched in gastrulation and mesoderm development (Figures 2D and 2E).

To characterize the spatiotemporal chromatin dynamics of the DAEs in MJH and PSH, we quantified the occupancies of the enhancer marks H3K4me1 and H3K27ac, as well as the promoter mark H3K4me3 at these DAEs across the E6.5-7.5 developmental stages(Yang et al. 2019). MJH/Cluster 8 and PSH/Cluster 2 specific DAEs could potentially function as enhancers and become activated at anterior regions of E7.5 embryos, as the DAEs were generally marked by H3K27ac and H3K4me1, but not H3K4me3 (Figure 2F). The enrichment of H3K4me1 at E7.0 and even earlier at E6.5 stage, along with the higher levels of H3K27ac at these DAEs at E7.5 stage, suggested that many of these DAEs could be dormant or inactive enhancers during earlier stages like E6.5, but primed for later activation during lineage specification (Figure 2F). Indeed, the integrated analysis on chromantin accessibility of the DAEs, shown by ATAC-seq, and their target gene expression levels supported that a large portion of the MJH/Cluster 8 DAEs were primed before the full activation of their target genes, like Bmp4, Hand1 and Foxf1 (Figure 2G). For example, seven DAEs associated with the Bmp4 gene were identifited, three of which were primed at E7.0 as marked by low levels of H3K27ac but high levels of H3K4me1 (Figure 2H). Taken together, the combined transcriptome and chromatin accessability analysis further supported the early lineage segregation of MJH and the epigenetic priming at gastrulation stage for early cardiac genes.

Identification of lineage specific key TFs

An integrated analysis of motif enrichment at the DAEs and TF expression data allowed us to identify potential lineage specific key TFs. The PSH/Cluster 2 specific DAEs showed motif enrichment similar to the recognition sequences of known NM specific TFs, including GATA4, ZIC3, EOMES, OTX2, and LHX1 (Figure S4B). Binding sites for hematoendothelium-related TFs, such as GATA2, FLI1, JUNB, and SOX7, were enriched in the hematoendothelial progenitors/Cluster 6 DAEs. In contrast, the binding motifs of GATA4, HAND1, FOXF1 and TEAD4 were highly over-represented in the MJH/Cluster 8 specific DAEs (Figure 3A). Among those MJH-related TFs, GATA4 and TEAD4 showed similar expression and motif activities at Clusters 2, 7, 8 of both MJH and PSH lineages. HAND1 and FOXF1 demonstrated both strong motif activities and specific expression at Cluster 7 and Cluster 8. Interestingly, the expression of HAND1 and FOXF1 seemed relatively transient at NM and EEM cells, and then became downregulated at CM of E7.75 (Figure S4C).

Figure 3:

Based on E7.0 single-cell multi-omics data analysis, we predicted a core regulatory network (CRN) centering on the four TFs (Figure 3b). Functional enrichment analyses indicated that this CRN could control key aspects of MJH specification, including Wnt signaling, epithelium cell migration, cell number maintenance, mesenchyme development, and cardiac development. In the CRN, GATA4 and TEAD4 controlled larger gene sets related to the transition from epiblast to mesodermal status, necessary for both MJH and PSH. HAND1 and FOXF1 co-regulated functionally more specific gene sets critical for differentiation to EEM status (Figure 3b). Consistently, most HAND1 and FOXF1 target genes were specifically expressed in EEM, in contrast to GATA4 and TEAD4 target genes (Figure 3b).

We also performed Chromatin immunoprecipitation sequencing (ChIP-Seq) analysis to profile the chromatin occupancies of HAND1 and FOXF1 in mesoderm (MES) and cardiac progenitor (CP) cells derived from the step-wise directed cardiomyocyte differentiation of mouse embryonic stem cells(Wamstad et al. 2012). We also collected published GATA4 ChIP-seq data(He et al. 2014) of E12.5 mouse embryonic heart, FLI1, ZIC2, ZIC3 and MESP1 ChIP-seq data(Lin et al. 2022) of 2.5 days EB differentiation. Direct comparison of ChIP-seq occupancy profiles with DAEs confirmed the specific enrichment of GATA4, HAND1, and FOXF1 at clusters 7 and 8 of MJH lineages, Fli1 at Clusters 5 and 6, while MESP1 is specifically enrichment at mesoderm cell clusters (Figure 3C). We also noticed the enrichment of HAND1 at early Cluster 3-specific DAEs, while FOXF1 tends to show more speicifc enrichment at cardiac specific enhancers at later CP cells.

HAND1 and FOXF1 regulate the MJH specific genes

In order to further investigate the molecular roles of HAND1 and FOXF1 in MJH specification, we generated Hand1 and Foxf1 KO mouse embryonic stem cell (mESC) lines (Figures S5A-D), followed by in vitro mesoderm differentiation and RNA-seq analyses. 2,331 down-regulated and 1,714 up-regulated genes in Hand1 KO mesoderm (MES) cells, and 870 down-regulated and 970 up-regulated genes in Foxf1 KO MES cells with fold-change (FC) > 1.5 and P-value < 1e-5 were identified (Supplementary Table 2). To explore whether HAND1 and FOXF1 are required for the proper expression of the MJH related genes, we examined the expression of the signature genes of the 9 cell clusters in control, Hand1 KO and Foxf1 KO MES cells. First of all, over 90% of the cluster specific genes were detected in the control MES transcriptome, indicating that the in vitro differentiation model could be a reliable tool to study the regulation of these cluster specific signature genes (Figure 4A). Indeed, whole transcriptome cosine similarity analysis revealed that the in vitro differentiated MES cells were more close to MM and EEM cell state transcriptome-wide (Figure 4B). Hand1 KO and Foxf1 KO led to down-regulated expression of a large portion of the Clusters 3, 7 and 8 MJH marker genes, but up-regulated expression of many of the Clusters 0 and 1 PSH marker genes in MES cells (Figure 4A). Consistently, the cosine similarity suggested that Hand1 and Foxf1 depletion could lead to EEM differentiation defects (Figure 4B). In addition, many of the down-regulated MJH specific genes and up-regulated PSH specific genes were directly bound by HAND1 and FOXF1 (Figures 4A and S5E). The analysis indicated that HAND1 and FOXF1 could dually regulate MJH specification through directly activating the MJH specific genes and inhibiting the PSH specific genes.

Figure 4:

Mutual regulation between HAND1 and FOXF1 in driving cardiac specific gene expression

We found that around 50% of the dysregulated genes in Foxf1 KO MES cells were also dysregulated in Hand1 KO MES cells, suggesting their synergistic function in transcriptional regulation (Figure 4C). For example, the key MJH specific genes, including Hand1, Foxf1, Bmp4, Tbx20 and Pmp22 were significantly down-regulated, while the epiblast genes Dnmt3b, Sema6a and Fst, and the NM specific genes Otx2 and Zic2 were substantially up-regulated in both KO cells (Figures 4C and S5E). Depletion of Hand1 blocked the activation of Foxf1 during cardic progenitor differentiation, and vice versa, while Hand1 overexpression was able to activate Foxf1 and the other MJH specific genes, and vice verse (Figures S6A-D).

The functional relevance between HAND1 and FOXF1 in target gene regulation could be attribuited to the enrichment of HAND1 at the Foxf1 enhancers, and vice versa (Figure 4D). We then used CRISPR-Cas9 technology to delete the putative enhancers of Hand1 and Foxf1 (Figures S6E-F). Indeed, deletion of the HAND1-bound Foxf1 enhancer (Foxf1-eH) abolished the activation of Foxf1; deletion of the FOXF1-bound Hand1 enhancer (Hand1-eF) also significantly reduced the expression of Hand1 during cardiac differentiation (Figures 4E and 4F). Importantly, the induction of the MJH specific genes Bmp4, Pmp22 and Spin2c was severely impaired after deletion of either Foxf1-eH or Hand1-eF. Overexpression of HAND1 in the Hand1-eF KO cells and FOXF1 in the Foxf1-eH KO cells were able to rescue the levels of these MJH specific genes during cardiac differentiation. To further investigate the role of Foxf1-eH in Foxf1 expression, we performed CRISPR activation (CRISPRa) assay to activate Foxf1-eH (Figures 4G). CRISPRa of Foxf1-eH led to a specific increase the expression of Foxf1 and its downstream target genes, but not its neighboring genes (Figure 4H). Together, our data indicated that mutual regulation between HAND1 and FOXF1 could play a key role in MJH cardiac progenitor specification.

Hand1 KO leads to MES overproliferation but cell death after exiting from the MES status

Meta analysis demonstrated that HAND1 occupied the early Cluster 3 specific DAEs, while FOXF1 tended to show more specific enrichment at the late cardiac specific Clusters 7 and 8 enhancers (Figures 3C and 5A). We also noticed that the Hand1 KO MES colonies were evidently much larger than those of WT control and Foxf1 KO, though the same number of EB cells were seeded at equal density for MES differentiation (Figure 5B). While Foxf1 KO barely affected cell proliferation rate, the count of the Hand1 KO cells was substantially increased when the cells were differentiated from EB towards MES state (Figure 5C). The Hand1 KO cells gradually lost viability upon in vitro cardiac lineage induction. Consistent with the increased proliferation rate observed, gene ontology (GO) analysis revealed that the genes involved in negative regulation of proliferation and positive regulation of cell migration were specifically down-regulated in Hand1 KO, but not Foxf1 KO, MES cells (Figures 5D and S6G).

Figure 5:

The Foxf1 KO MES colonies derived from MES differentiation appeared phenotypically normal, but were not able to further differentiate into beating CMs (Figure 5B and Supplementary Video 1). The expression of the mature CM markers Tnnt2, Myh6 and Myh7 were also substantially lower after Foxf1 KO (Figure S6H). To further examine the function of FOXF1 in cardiac progenitor specification, we performed RNA-seq analysis in control and the CP cells derived Foxf1 KO mESCs. Gene Set Enrichment Analysis (GSEA) revealed the expression levels of the NM specific genes (like Zic3, Pou5f1 and Cited1) and the PSH specific genes (like Isl1, Ifitm1 and Foxc2) were remarkably up-regulated, while the expression levels of the MJH specific genes (like Hand1, Spin2c and Tdo2) and the early CM specific genes (like Acta2, Tnnc1 and Tnnt2) were significantly reduced (Figure 5E). Thus, our data further supported the specific and synergistic roles of HAND1 and FOXF1 in MJH cardiac progenitor specification.

Genetic loss of Hand1 blocks the specification of mesoderm along MJH

To further invesigate the roles of HAND1 in MJH specification in vivo, we generated floxed allele of Hand1 by inserting loxP sites flanking the exon 1 of the Hand1 gene. Genetic crosses of the Hand1^fl/fl mice with Mesp1-Cre mice(Saga et al. 1999) allowed specific deletion of Hand1 in mesodermal cells (Figure S7A). Consistently, Foxf1 is also co-expressed in HAND1-positive EEM cells and its expression was also drastically down-regulated in MESP1-CRE driven Hand1 conditional KO (Hand1 CKO) embryos (Figure 6A). Mesodermal deletion of Hand1 generally led to smaller embryos at E7.0 and also later stages (Figure S7B). The Hand1 CKO embryos appeared to die by E9.5 due to embryonic turning failure and heart looping abnormality, mirroring the previously reported phenotype of Hand1 KO mice (Firulli et al. 1998; Riley et al. 1998) (Figure 6B).

Figure 6:

To assess how Hand1 loss affects early mesoderm development in vivo, we analyzed the single-cell transcriptomics of control and Hand1 CKO embryos at E7.0. Integrated analysis indicated the loss of EEM cells, but the abnormal accumulation of primitive streak (PS), NM and MM cells in Hand1 CKO embryos (Figures 6C and 6D), which suggested that specific deletion of Hand1 in mesodermal cells strongly affected early mesoderm differentiation. Detailed analysis regarding the percentage of each cell type revealed the specific reduction of cell numbers from EEM, hematoendothelial progenitors, and ExE ectoderm cell clusters in Hand1 CKO embryos (Figure 6D). We then compared the gene expression pattern for each cell cluster and calculated the significantly affected genes between Control and Hand1 CKO embryos (Figure 6E). Consistently, EEM cell cluster was the most affected clusters upon Hand1 loss (Figures 6E and S7C). Although the percentages of hematoendothelial progenitors and ExE ectoderm cells were reduced, the expression of both cell type specific marker genes did not seem affected drastically. To further illustrate the developmental progression of the mesodermal lineage, we performed paired URD lineage inference analysis(Farrell et al. 2018), further confirming the specific development block of NM specification towards EEM in MJH (Figure 6F). We also performed label transfer analysis to identify MJH and PSH cells from the E7.0 scRNA-seq data (Figure S7D). Our analysis showed that the fraction of MJH cells increased by over 2-fold, whereas the fraction of PSH cells remained unchanged (Figures S7E and S7F), supporting that Hand1 deletion driven by MESP1-CRE blocked the MJH direction of mesoderm differentiation.

The abnormal accumulation of PS, NM and MM cells in Hand1 CKO embryos appeared to be consistent with the phenotypes observed in in vitro mesoderm differentiation of Hand1 KO mESCs (Figures 4B and 5B). Notably, the genes involved in negative regulation of proliferation and down-regulated in Hand1 KO MESs were also reduced in the MJH, but not PSH, trjectory of E7.0 Hand1 CKO embryos (Figure 7A). In addition, cell migration related genes were also affected in Hand1-depleted MES and embryos. To further validate this phenotype, we performed sequential DAPI staining on cryo-sectioned E7.0 control and Hand1 CKO embryos, followed by cellular segmentation and cell density measurement (Figures 7B and 7C). The analysis revealed that the mesoderm cells near the extraembryonic region in Hand1 CKO embryos were more compacted (Figures 7B’’, 7C’’, and 7D), while the distal region of the Hand1 CKO embryos showed no obvious difference from the control embryos (Figures 7B’’’ and 7C’’’). In addition, reduced exocoelomic cavity (EC) size and increased number of mesodermal cells in extraembryonic region were also observed in Hand1 CKO embryos (Figures 7B’ and 7C’). Hematoxylin and eosin (H&E) staining of E7.5 embryos also supported reduced embryo size and accumulated cells at mesodermal regions upon loss of Hand1 (Figure 7E). These data together establish HAND1 as a factor in promoting the specification of mesodermal cells toward MJH.

Figure 7:

Discussion

In this study, we described the transcriptional trajectories and epigenetic landscapes of early cardiac specification event in mouse embryos, and identified the CRNs underlying the early mesodermal lineage specification. Mechanistically, we demonstrated that this earliest cardiac MJH specification event was tightly regulated by HAND1 and FOXF1. HAND1 and FOXF1 were mutually regulated, but also performed their respective functions in MJH cell fate determination. In vitro differentiation and in vivo mouse model analyses indicated that HAND1 was essential for exiting from NM program during cardiac specification. Depletion of FOXF1 impaired the capacity of mesodermal cells in cardiac specification (Figure 7F). In sum, our findings provided new insights into the transcriptional determinants that specify the early cardiac progenitors, paving the way for the identification of potential therapeutic targets for treating congenital heart defects.

Recent studies using scRNA-seq have reported the roadmaps of mammalian embryonic lineage development(Pijuan-Sala et al. 2019; Mittnenzweig et al. 2021; Qiu et al. 2022). Together with studies which focus on cardiac progenitors(de Soysa et al. 2019; Tyser et al. 2021; Zhang et al. 2021), several models have been proposed to explain the multi-lineage process of heart formation. In this study, we found significant distinction between the MJH and PSH trajectories at E7.5-7.75 and performed forward-backward tracing to achieve two non-overlapping trajectories. The two trajectories are consistent with the previous clonal analysis and also the HAND1+ cells lineage tracing by Zhang et al. (Meilhac et al. 2004). Importantly, we here identified specific trajectory-related gene cohorts throughout the whole process, shedding light on the molecular mechanisms by which HAND1 is crucial for promoting the exit of NM program during the early MJH specification process.

MESP1 serves as a master regulator in the establishment of the cardiac lineage. An analysis of previously published MESP1-regulated genes revealed that GATA4, HAND1 and FOXF1 were directly controlled by MESP1. GATA4 was rapidly induced by MESP1 within 12 hours, while HAND1 and FOXF1 were activated after 24 hours of MESP1 induction. On a pseudotime scale, we noticed a sequential-temporal expression pattern of GATA4, HAND1 and FOXF1 in MJH. Our transcriptional and epigenomic analyses suggested a sequential activation model of MESP1-GATA4-HAND1-FOXF1. However, it should be noted that MESP1 and GATA4 were activated in both MJH and PSH. The process of cardiac lineage segregation is a complex one that may involve TF regulatory networks and signaling pathways. It has been reported that BMP signaling activated the expression of Hand1 during heart formation (Zheng et al. 2021). Since Bmp4 displayed higher specificity for MJH, it might explain the activation of Hand1, and subsequent Foxf1 activation in MJH, but not the PSH. Thus, HAND1 and FOXF1 might be subject to feed-forward activation by MESP1 and GATA4 in concert with BMP signalling, thereby promoting MJH cardiac specification. Future studies are worthwhile to further investigate how TFs work together with signaling pathways to promote cardiac lineage segregation.

We further demonstrated the synergistic roles of HAND1 and FOXF1 in early cardiac specification. Foxf1-/- mice were not able to survive beyond E10.0, with incomplete separation of splanchnic and somatic mesoderm, leading to abnormal coelom formation. Interestingly, FOXF1 marked the mesothelium lining (a monolayer epithelial lining of the exocoelomic cavity), the anterior distal part of which extends to the presumptive cardiac mesoderm (Fleury et al. 2015). Such an expression pattern is coincident with that of HAND1 in JCF (Tyser et al. 2021). Our studies showed that FOXF1 was indispensable for the specification of mesoderm toward early cardiac progenitors, through binding the MJH progenitor specific enhancers for their later activation. Together, the combined analyses on the developmental route of CMs and their transcriptional determinants will further enhance our understanding of the etiology behind congenital heart defects, ultimately providing insights into potential regenerative strategies for heart disease treatment.

Methods

Mice

All mouse experiments were approved by the Animal Care and Use Committee at Southeast University, and performed in accordance with institutional guidelines. Mice were housed in cages under SPF conditions and had free access to water and food. Hand1^fl/fl mice were produced by Cyagen with loxP sites inserted into the regions surrounding exon 1 in the Hand1 locus. Hand1^fl/+; Mesp1-Cre mice were generated by crossing Hand1^fl/fl mice with Mesp1-Cre mice (Nie et al. 2015).

The Hand1^fl/fl or Hand1^fl/+ female mice caged with male Hand1^fl/+; Mesp1-Cre mice. Females were screened for vaginal plugs following morning (E0.5). To obtain post-implantation embryos, female mice at 7.0 or 9.5 days post-coitum (d.p.c.) were executed and the uteri were dissected and transferred to a petri dish with PBS. For 7.0 d.p.c embryos, each decidua was carefully freed from the uterine muscle layers using properly sharpened forceps. Then the embryos were carefully separated from decidua. Reichart’s membrane and the ectoplacental cone were also removed from the 7.0 d.p.c embryos. The uteri surrounding the 9.5 d.p.c embryos were cut open with a small incision, the embryos were genetly squeezed out, and then the amniotic membranes were removed. Images of the embryo were acquired on a stereo microscope (Mshot, MZ62; Olympus, SC180).

Genomic DNA for mouse genotyping was obtained from mouse tail biopsies. Genotyping of mouse embryos was performed with genomic DNA, which was obtained by collecting extra-embryonic regions subsequently digestion using Mouse Direct PCR Kit (Bimake, B40013). PCR reactions were used to detect the Cre transgene and the Hand1 loxP site. Thermal cycle reactions were as follows: 3 min at 94 °C, 35 cycles of 30 s at 94 °C, 35 s at 60 °C, 45 s at 72 °C and a final 5 min extension at 72 °C. Primer sequences used in this study are listed in Supplementary Table 3.

Antibodies

Antibody against HAND1 (sc-390376) (WB: 1:1000, IF for embryos: 1:100, IF for cells: 1:200) was purchased from Santa Cruz. Antibody against FOXF1 (Abclonal, A13017) (WB: 1:2000, IF for embryos: 1:100, IF for cells: 1:500) was purchased from Abclonal. Antibody against HA (ChIP: 3 μg) was generated in house. Antibody against TUBULIN (66031-1-Ig) (WB: 1:50, 000) was purchased from Proteintech. Goat anti-mouse IgG Alexa Fluor 488(A-11001) and goat anti-rabbit IgG Alexa Fluor 546 (A-11035) were purchased from Thermo Fisher Scientific.

Mouse ESC culture

Mouse E14 ESCs were maintained in DMEM (Hyclone, SH30243.01) supplemented with 15% FBS, 1 × nonessential amino acids (STEMCELL, 07600), 1 × GlutaMAX (Gibco, 35050-061), 1 × penicillin streptomycin solution (Sangon Biotech, E607011-0500), b-mercaptoethanol (ALDRICH), 0.1 μg/ml leukemia inhibitory factor (Novoprotein, C690), 3 μM CHIR99021 and 1 μM PD0325901 in gelatin-coated plates at 37 °C under 5% CO₂.

CRISPR-Cas9 guided KO

SgRNAs targeting Hand1, Foxf1 and their enhancer sites were cloned into lentiCRISPR v2. These constructs were transfected into 70% confluent 293T cells together with 6 μg of psPAX2 packaging plasmids and 2 μg of pMD2.G envelope plasmids using Highgene (Abclonal, RM09014). The media was half-replaced with fresh DMEM supplemented with 10% FBS 6 h after transfection. The lentiviral supernatants were harvested at 24 h, 48 h, and 72 h post-transfection, filtered through 0.45 μm filters, and concentrated at 50,000 × g for 0.5 h. Mouse E14 ESCs were infected with concentrated lentiviral particles with polybrene (Sigma) at the concentration of 8 μg/ml. 2 μg/ml puromycin was used for selection for 48 h and individual colonies were picked and expanded in 48-well plates. The clones were screened with genomic PCR, and confirmed by TA cloning and Sanger sequencing. Oligo sequences used in this study are listed in Supplementary Table 3.

In vitro cardiac differentiation

Mouse E14 ESCs were differentiated into embryoid bodies (EBs) at a density of 75,000 cells/ml in 6 cm dishes for 48 h in serum-free media (DMEM (Hyclone, SH30243.01), DMEM/F12 (1:1) (Gibco, 11320-033), 0.05% BSA (Sigma, A1933), 1 × GlutaMax (Gibco, 35050-061), B27 supplement (Gibco, 12587010), N2 supplement (Gibco, 17502048), supplemented with 50 mg/ml ascorbic acid (Alfa Aesar, 50-81-7) and 4.5×10^-4 M monothioglycerol). EBs were dissociated into single cells and re-aggregated as MES cells for 40 h at 50,000 cells/ml in the presence of 5 ng/mL human VEGF (Novoprotein, C083) and 10 ng/ml human Activin A (Peprotech, 96-120-14-10) and 0.3 ng/ml human BMP4 (R&D, 5020-BP-010). MES cells were dissociated and plated as a monolayer in gelatin-coated 12-well plate at 50,000 cells/ml in StemPro-34 (Gibco, 10639011) supplemented with 5 ng/mL VEGF (Novoprotein, C083), 10 ng/mL human basic FGF (Gibco, PMG0035) and 25 ng/mL FGF10 (R&D, 345-FG-250). CP cells were harvested after differentiation for 32 h. Contracting CMs can be observed after another 5 days. Differentiation stages were confirmed by the expression of marker genes.

Immunofluorescence

Mouse Embryos: Embryos were fixed in 4% PFA, then rinsed in PBS and incubated in 30% sucrose at 4 °C until sinking. Embryos were then embedded in OCT medium and snap-frozen in liquid nitrogen and stored at −80 °C. Sections were cut at a thickness of 10 μm. To perform IF, embryo sections were washed three times for 5 min each, then antigen repair was performed using citrate. The samples were permeated for 40 min with 0.3% TritonX-100 followed by washing three times with PBS for 5 min each. The samples were then blocked with ReadyProbes 2.5% Normal Goat Serum (Thermo) for 1 h at room temperature. The samples were incubated overnight at 4 °C with primary antibodies followed by washing three times with PBS for 10 min each. The samples were incubated with secondary antibodies for 1 h at room temperature, then washed three times with PBS for 10 min each. The samples were then incubated with DAPI at room temperature for 10 min followed by washing one time with PBS for 5 min. After the final wash, the samples were mounted with mounting buffer. Images were captured using a Zeiss 700 laser confocal microscope.

MES and CP cells were dissociated into single cells and fixed on 12-well plate cell slide with 4% PFA at room temperature for 15 min. The fixed cells were washed 3 times with PBS and permeabilized with 0.2% Triton X-100 for 5 min at room temperature. The cells were then blocked with PBS containing 2% BSA and 0.3% TritonX-100 for 1 h at room temperature. Appropriate dilution of primary antibody was added and cells were incubated overnight at 4 °C. After washing with PBS for 3 times, the cells were stained with secondary antibodies and DAPI for 1 h at room temperature followed by mounting on slides. Images were captured using a Zeiss 700 laser confocal microscope.

H&E staining

Whole E7.5 embryos were fixed overnight with 4% PFA, embedded vertically in clean paraffin, then sliced to obtain 7 µm paraffin sections. Standard hematoxylin and eosin (H&E) staining methods were used to stain the sections. The images of H&E staining were obtained on a microscope (Olympus, IX73).

Western blotting

Cells were washed twice in ice-cold PBS and incubated in 1 × SDS lysis buffer for 15 min at 95 ℃. Lysates were pre-cleared by maximum speed centrifugation for 3 min, then separated by SDS-PAGE and transferred to polyvinylidene fluoride (PVDF) membrane. 5% non-fat dry milk was used for blocking and primary antibody was added. The membrane was incubated overnight at 4 °C, then washed, incubated with secondary antibody for 1 h at room temperature. ECL substrate was used for imaging by autoradiography.

Quantitative RT-PCR and bulk RNA-seq

Total RNA was extracted from cells using RNA isolater Total RNA Extraction Reagent (Vazyme). 500 ng of total RNA was used to synthesize cDNA using ABScript II RT Mix (Abclonal, RK20403). Resultant cDNA was diluted in water and 12.5 ng cDNA was used in each qRT-PCR reaction. Reactions were run on CFX96 (Bio-Rad) using 2 × Universal SYBR Green Fast qPCR Mix (Abclonal, RK21203). The relative expression levels of genes of interest were normalized to the housekeeping gene Actin. Each experiment contains at least three biological replicates. Primer sequences used in this study are listed in Supplementary Table 3. For RNA-seq, 1 × 10⁶ MES or CP cells were harvested and RNA extraction was performed using Rneasy mini plus kit (Qiagen). 1 μg of total RNA was used for the construction of sequencing libraries and sequencing.

snRNA-seq

To prepare single cells for snRNA-seq, embryos at E7.0 were dissected in cold sterile 1 × PBS without Ca²⁺, Mg²⁺ under a stereo microscope. Embryos were staged based on their morphology. The Reichert’s membrane and ectoplacental cone were removed, and genotypic identification of embryos was carried out. Embryos were placed into a 1.5 ml microfuge tube and digested into single cells using TrypLETM Express (Thermo) at 37 °C. Single-cell RNA sequencing was performed on single-cell suspensions using 10 × Genomics.

10 × Multiome library preparation and high throughput sequencing

Embryos staged at E7.0 were collected and washed with cooled 0.5% BSA/PBS for twice. Sufficient embryos were pooled and then subjected to 200 μL TrypLE for cell dissociation for 10 min at 37 °C with frequent gentle mixture. Single-cell suspension of embryos were then quenched and washed with 0.5% BSA/PBS, and finally filtered using 40 μm Flowmi cell strainer. The acquired single cell suspension were then subjected to nuclei isolation, library preparation by following the manufacturer’s instruction. The library was sequenced on Novaseq 6000 platform with recommended sequencing depths and read lengths.

ChIP, ChIP-seq library preparation

5 × e⁶ cells were crosslinked with 1% formaldehyde for 10 min and quenched with 0.125 M glycine for 5 min at room temperature. Fixed cells were incubated in 0.5 ml lysis buffer (50 mM Tris-HCl [pH 8.0], 1 % SDS, 5 mM EDTA) for 10 mins on ice. After adding 1 ml dilution buffer (20 mM Tris-HCl [pH 8.0], 150 mM NaCl, 2 mM EDTA, 1% Triton X-100), chromatins were sonicated into 200-800 bp fragments using Bioruptor (Diagenode) and immunoprecipitated with protein A agarose beads and specific antibody at 4 ℃ for 12 h. Immunoprecipitates were washed with Wash Buffer I (20 mM Tris-HCl [pH 8.0], 150 mM NaCl, 2 mM EDTA, 1% Triton X-100, 0.1% SDS), Wash Buffer II (20 mM Tris-HCl [pH 8.0], 500 mM NaCl, 2 mM EDTA, 1% Triton X-100, 0.1% SDS), Wash Buffer III (10 mM Tris-HCl [pH 8.0], 0.25 M LiCl, 1 mM EDTA, 1% deoxycholate, 1% NP-40) and TE, respectively. After the final wash, DNA was eluted and reverse-crosslinked at 65 °C for at least 6 h. DNA was then purified and used for PCR amplification or ChIP-seq library preparation. ChIP-seq libraries were prepared with VAHTS Universal DNA Library Prep Kit for sequencing.

Data processing and quality control

Single cell multi-omics data: sequence alignment

We used Cell Ranger ARC software suite for the analysis of single cell multi-omics (ATAC & Gene Expression) data (https://support.10xgenomics.com/single-cell-multiome-atac-gex/software/pipelines/latest/what-is-cell-ranger-arc). We applied the “cellranger-arc count” function for barcode counting, adapter/primer removal and sequence alignment. The processed reads were aligned to mm10 using the BWA-MEM algorithm. For ATAC sequencing data, BAM files were generated for downstream analysis. For Gene Expression (GEX) data, “cellranger-arc count” generated gene-barcode matrices were used for further analysis. Version of Cell Ranger ARC genomic sequence and gene annotation: refdata-cellranger-arc-mm10-2020-A-2.0.0.

snRNA-seq data

snRNA-seq datasets of the current study include: (1) E6.5-8.5 whole mouse embryos; (2) micro-dissected anterior cardiac regions of mouse embryos at early crescent to linear heart tube (∼E7.75-8.25) stages. (3) GEX part of the single cell multi-omics data of E7.0 mouse embryos. (4) Control and Hand1 CKO E7.0 mouse embryos.

We used R package ‘Seurat’ for QC and normalization purposes (Hao et al. 2021). Cells with abnormal sequencing depth (nFeature_RNA < 2000 or nCount_RNA > 1e5) or with high mitochondrial ratio (percent.mt > 5) were excluded. We used the ‘SCTransform’ to perform normalization, variance stabilization, and regression of cell cycle scores (using the ‘CellCycleScoring’ function). For dataset (1), we selected E7.0 samples, removed cells labelled as ‘doublet’ or ‘stripped’ and regressed out sequencing batches.

snATAC-seq data

Bam files of snATAC-seq from multiome data were processed using Snaptools (Fang et al. 2021) function ‘snap-pre’ to remove low-quality fragments (MAPQ < 30), over-sized fragments (length > 1,000 bp), secondary alignments and PCR duplicates. This function also generates cell-by-bin matrices of variable resolution (bin sizes: 1kb, 5kb, 10kb) for downstream analysis. The range of fragment coverage for bin selection was set between 500 and 20,000.

ChIP-seq data analysis

Adapter sequences were removed from Fastq files using trim_galore (v0.6.7) with default parameters. After trimming, ChIP-seq reads were aligned to the mm10 mouse genome using Bowtie2 (v2.3.5.1) (Langmead and Salzberg 2012) with ‘--no-mixed’ and ‘--no-discordant’ parameters. The aligned files were sorted and converted to BAM format using SAMtools (v1.9). BigWig files were subsequently generated using deepTools bamCoverage (v3.5.0) (Ramirez et al. 2016), employing CPM normalization and ignoring duplicates. Peak calling was performed with MACS3 (3.0.0a6) (Zhang et al. 2008) for HAND1 (Q-value < 1e-5) and FOXF1 (Q-value < 1e-10) ChIP-seq data.

Tracing the CM, MJH and PSH trajectories

The trajectories of each in Figure 1 were inferred using the Waddington-OT Python package (v1.0.8) (Schiebinger et al. 2019) with a predefined starting cell set (E8.5 CM cells for tracing the CM trajectory, E7.5 EEM cells in the CM trajectory for tracing MJH trajectory, and E7.75 pharyngeal mesoderm cells in the CM trajectory for tracing the PSH trajectory), and cells of each trajectory are selected if WOT score > 0.0001. To specifically distinguish between MJH and PSH, the difference in the WOT score of the two trajectories was calculated, and if the difference is greater than 0, it belongs to the MJH trajectory and vice versa to the PSH trajectory. The WOT package was used with default parameters as in the Waddington-OT online tutorial (https://broadinstitute.github.io/wot/tutorial/).

Signaling pathway enrichment analysis

Potential signaling-activated/inhibited genes of Bmp, Yap, Wnt, Nodal, Notch and Fgf were collected from Peng et al. 2019 (Peng et al. 2019). Enrichment of signaling target gene sets was determined by the average expression levels of activated/inhibited genes.

Spatial mapping of single cells

Locations of mouse E7.0/7.5 mesodermal cells were inferred by comparison with the GEO-seq data, where each sample represents 5–40 cells with defined spatial locations. Transfer component analysis (TCA) (Pan et al. 2011) was performed to achieve shared representation of scRNA-seq and GEO-seq samples. Single cells were mapped to GEO-seq locations with the highest correlation coefficients.

Identification of pseudotime-dependent genes

Single-cell pseudotime trajectory analysis was performed using R package Monocle 2 (v2.22.0) (Trapnell et al. 2014) according to the online tutorials (http://cole-trapnell-lab.github.io/monocle-release/). Monocle object was directly constructed using Monocle implemented new Cell Data Set function from Seurat (v4.1.1) object, and Monocle implemented differentialGeneTest function was used to find highly variable genes for ordering. Based on these, we selected key and specific genes in MJH and PSH trajectory, respectively, and further visualized them in the heatmap using cell bin-by-gene matrix.

snRNA-seq data integration and label transfer

This research included data integration for the following datasets (dataset numbers in ‘Data processing and quality control: snRNA-seq data’): (1) and (2) for comparison of the predicted trajectories and MJH cells; (1) and (3) for cell-type annotation of multiome single cells; (1) and (4) for cell-type annotation of Ctrl/Hand1 CKO single cells. For multiome dataset, integration was performed by two steps: 1) whole embryo integration; 2) according to the predicted cell types of (1), mesoderm specific integration was done by selecting relevant cell types (NM, MM, EEM, Haem, PGC).

The integration and label transfer process follows the pipeline provided by Seurat (Hao et al. 2021): https://satijalab.org/seurat/articles/integration_introduction.html. 4,000 genes were selected for integration based on variable gene sets of individual datasets. We selected “SCT” as the normalization method for identification of integration or transfer anchors.

Clustering analysis of snRNA-seq data

We performed clustering analysis for cells of the integrated mesodermal lineage. Seurat (Hao et al. 2021) function “FindNeighbors” and “FindClusters” were applied to the Top 20 principle compoments with the following parameter setting: annoy.metric = “cosine”, resolution = 0.95. This analysis resulted in nine clusters (C0-8), which is also used for snATAC-seq analysis.

FHF and SHF gene signature scores

For each gene set, we defined the average z-score normalized expression levels as the signature score per cell. The FHF and SHF gene sets were collected from Soysa et al. 2019 (de Soysa et al. 2019) Supplementary Table 1 (FHF: “FHF” genes of tab “Figure1de_n = 21,366 cells”; SHF: “MP” genes of tab “Fig 2ab_n = 2,103 cells”).

snATAC-seq data analysis

snATAC-seq bam files were combined for each cluster, C0-8, for identification of accessible elements (AE). AEs were then collected as a basic set of peak annotation. We used Snaptools (Fang et al. 2021) function ‘add_pmat’ to generate a cell-by-peak matrix (pmat), where the value of each element represent the accessibility per cell per peak. Diffusion map followed by UMAP analysis was performed to generate the snATAC-seq version of 2D data visualization.

Cluster specific AEs were identified using SnapATAC function ‘findDAR’. For analysis of C2-8 specific AEs, we used C1, which represents the least differentiated cell group, as the background cluster. For C0-1, we used cells from all other clusters as background. Cluster specific AEs were provided to HOMER (Heinz et al. 2010) function ‘findMotifsGenome.pl’ for motif analysis, and to deepTools (Ramirez et al. 2016) for plotting heatmaps of ChIP-seq data.

Definition of enhancer regulated target genes

Enhancer-promoter (EP) pairs were predicted by SnapATAC function ‘predictGenePeakPair’ (Fang et al. 2021). This method performs logistic regression using peak accessibility (snATAC-seq) and expression (snRNA-seq) of neighboring genes to identify EP links. Candicate EP pairs need to be less than 50kb apart. Cluster specific target genes were defined as EP-associated genes of cluster specific AEs. Functional enrichment analysis was performed using R package clusterProfiler (Wu et al. 2021) (‘compareCluster’ function, ontology was selected as ‘BP’).

Analysis of TF activity

We analyzed TF activity, defined as the relative accessibility of AEs containing the motif sequence(Schep et al. 2017) of a TF, using R package chromVAR (https://greenleaflab.github.io/chromVAR/articles/Introduction.html). SnapATAC-generated cell-by-peak matrices were converted to ‘SummarizedExperiment’ format as the input of chromVAR. The TF activity is calculated by the ‘computeDeviations’ function by comparing the accessibility of motif-containing AEs with background peaks with similar GC content. Genomic coordinates of motifs were acquired by the ‘getMatrixSet’ function and R package ‘JASPAR2020’. Synergy scores between TF pairs were computed by the ‘getAnnotationSynergy’ function.

Classification of activation modes of gene-enhancer pairs

We performed pseudotime trajectory analysis for E7.0 MJH trajectory (C0, C3, C7 and C8) in multiome using the abovementioned method. To further explore the underlying mesoderm developmental mechanism of enhancer clusters in regulating gene expression, we first built gene-enhancer pairs for preselected clusters based on snATAC-seq and snRNA-seq datasets. Secondly, we partitioned cells from E7.0 mesoderm lineage into pseudotime bins and calculated the cell bin-by-gene matrix of snRNA-seq and the cell bin-by-peak matrix of snATAC-seq, respectively. Finally, we divided activation modes of gene-enhancer pairs into three groups, fast, sync., and slow, according to the order of bins corresponding to their activation time.

Inference for URD lineage tree

To reconstruct branching developmental trajectory trees for E7.0 mesoendodermal cells (Ctrl and Hand1 CKO), we used URD (v1.1.1) (Farrell et al. 2018). We used all cells assigned as epiblast as the root of the tree, and cells of clusters that contained the most differentiated cell-types at the latest pseudotime as the tips, where the CM branch was divided into MJH and PSH trajectories. By the movement of coordinates, the cells of Ctrl and Hand1 KO are located on either side of the branch, respectively.

CRN construction and in silico KO

The CRNs based on multi-omics (scRNA + snATAC) data was constructed by SCENIC+ (v1.0.1) (ref) (details available via https://scenicplus.readthedocs.io/en/latest/pbmc_multiome_tutorial.html). CRN-based in-silico KO was conducted by SCENIC + module “scenicplus.simulation”. In this step the expression level of a TF was set to zero and SCENIC + propagates the effect through CRN in an iterative manner to obtain the perturbed expression matrix. Cell-state transition vectors were generated and projected to the tSNE map using “plot_perturbation_effect_in_embedding” function.

Cell nucleus segmentation

We employed the cellpose (V2.1.1) (Stringer et al. 2021)algorithm for segmentation of DAPI staining for embyo transverse sections. The pre-trained CP model (--pretrained_model CP) was used to obtain the masks of images. We used option “--diameter 8” to resize the image to conform to the input parameters of the pre-trained model. We utilized the MorphoLibJ plugin of Fiji to obtain morphological metrics of the segmented nuclei. Finally, we utilized the VAA3D (2.938) (Peng et al. 2014) software to perform partitioning of the embryo, to identify cells belonging to the endoderm, mesoderm or epiblast.

Statistics and reproducibility

Numbers of biological replicates, statistical tests and P-values are reported in the figure legends. If not mentioned otherwise in the figure legend, statistical significance was determined using two-tailed Student’s t-test, provided by GraphPad Prism9 statistical software.

Data availability

Sequencing (sc/snRNA-seq, snATAC-seq, ChIP-seq, and RNA-seq) data that support the findings of this study have been deposited in the Gene Expression Omnibus under accession numbers GSE245713. Previously published ChIP-seq data that were re-analyzed here are available under accession codes GSE165107 (MESP1, ZIC2 and ZIC3, 2.5-day EB), GSE47085 (HAND1, FLK+MES cells), GSE52123 (GATA4, E12.5 mouse heart) and GSE69099 (FLI1, ES derived hemogenic endothelium). Previously published scRNA–seq data that were re-analyzed here are available under accession codes E-MTAB-6967 (ArrayExpress, E6.5-8.5 mouse embryos) and E-MTAB-7403 (ArrayExpress, micro-dissected heart-related samples of E7.75-8.5 mouse embryos). Previously published GEO-seq data that were re-analyzed here are available under accession codes GSE171588. All other data supporting the findings of this study are available from the corresponding author on reasonable request.

Acknowledgements

The authors are grateful to the Lin and Luo lab members for helpful discussion of this study. Studies in this manuscript were supported by funds provided by National Key R&D Program of China (2018YFA0800100 and 2018YFA0800101 to C.L.; 2018YFA0800103 to Z.L.), the National Natural Science Foundation of China (32030017 and 31970617 to C.L.; 31970626 to Z.L.; 32100529 to P.X.), Shenzhen Science and Technology Program (JCYJ20210324133602008 and JCYJ20220530160417038 to C.L.; JCYJ20210324133601005 and JCYJ20220530160416037 to Z.L.)

Declaration of interests

The authors declare that they have no conflicts of interest.

Supplemental Tables

Table S1.

Table S2.

Table S3: