The GSE106118 dataset in the GEO database (https://www.ncbi.nlm.nih.gov/geo/) contains human embryonic heart development data using STRT-seq. A total of 4,948 cells from 20 human fetal hearts were collected from 5 W to 25 W of gestation (Cui et al., 2019). Mouse dataset1 (GSE76118) was obtained using Smart-seq2 and contained the data of three time points: e8.5, e9.5, and e10.5, with a total of 2,967 cells (Li et al., 2016); two time points of mouse embryonic heart development: e13.5 and e14.5, a total of 11,532 cell data were obtained, which was named mouse dataset2 (GSE100861) (Xiao et al., 2018). The above data were imported into the Seurat v4.0.6 R package (https://satijalab.org/seurat/) for stringent quality control.

To avoid the impact of low-quality data on subsequent analyses and to maximize the retention of important and rare cell subpopulations or transcripts, we observed the distribution of samples, eliminated outliers, and retained as many transcripts or cells as possible. For cell filtering, the number of captured transcripts per gene of human embryonic hearts was inferred based on UMIs. Only cells with a minimally detected gene number of 1,000 were kept for downstream analysis. The number of transcripts detected less than 5,000 were also filtered. For mouse dataset1, only cells with the number of genes detected greater than 1,800 and the number of transcripts detected at least 1E + 06 were kept. Mouse dataset2 was processed with the same strategy. In addition, the proportion of mitochondrial genes to all genetic material may indicate whether a cell is in homeostasis or not. The cells with high proportions of mitochondrially-derived genes might indicate damaged or dead since the cell could lead to leakage of cytoplasmic RNA molecules, while mRNA in the mitochondria might be retained. Therefore, cells with high mitochondrial gene expression need to be filtered (Luecken and Theis, 2019; Andrews et al., 2021). In the heart, about 30% of the total mRNA comes from mitochondria due to the high energy demand of cardiomyocytes (Mercer et al., 2011; Galow et al., 2021). Thus, the cells with a mitochondrial gene percentage of >30% were moved out in both species. DoubletFinder was used to identify and remove doublets (McGinnis et al., 2019). Since the two mouse datasets were generated using different sequencing technologies (Drop-seq and Smart-seq2), the different sequencing depths and normalization means prevented the two datasets from being directly combined for analysis. A geometry-based sampling of mouse dataset2 from the Drop-seq platform was performed using Python’s geosketch to ensure equilibrium between the two mouse datasets. This geometry-based sampling approach was able to retain rare cell types and accelerated the integration of large scRNA-seq datasets, thereby improving and accelerating downstream analysis. 2,402 cells were finally sampled from mouse dataset2. For each dataset, the counts of each cell were normalized using the “NormalizeData” function in Seurat with a scale factor of 10,000. ScaleData function was used then for standardization. The top 2000 highly variable genes were selected using the FindVariableFeatures function. To correct potential batch effects in mouse datasets from different platform sources, we integrated them using the R package Harmony (https://github.com/immunogenomics/harmony). In the analysis, following the examples provided in the Harmony package in the R language environment, we performed Harmony within the Seurat 2 workflow with the maximal number of clusters (50) and the maximal number of iterations (100). The top 20 normalized Harmony vectors in PCA space were used as input to the assessment methods.

The pre-processed data were unsupervised and clustered using the FindClusters function in Seurat to identify different cell types based on the abundance of cell subpopulations and the representation of cell subpopulations in cardiac function based on known marker genes. Cells were classified into 20 clusters for humans and 17 clusters for mice. Next, we used the FindAllMarkers function to find differentially expressed genes (DEGs) of each type of cell.

The cross-species comparison was based on one-to-one orthologues annotated by the Ensembl genome annotation system (http://www.ensembl.org/index.html). The average expression levels in major cell types were calculated. Pearson correlations of major cell types in humans and mice were used to determine cross-species similarities.

Reconstructing pseudotime trajectories of cells were constructed using the Monocle2 R package (version 2.18.0, http://cole-trapnell-lab.github.io/monocle-release/) (Trapnell et al., 2014). First, the UMI of CMs, FBs, and ECs/FBs were respectively imported into Monocle2. Genes that were expressed in at least 10 cells were used, and only genes expressed in 5% of cells were kept. We used the differentially expressed genes between time points as the set of ordering genes and performed the dimension reduction and the trajectory analysis. Once we established a trajectory, we used the differentialGeneTest function to find genes that had an expression pattern that varied according to pseudotime.

GO characteristics of gene clusters were determined using the clusterProfiler package (v3.10.0.94) (Yu et al., 2012) for all DEGs with an average logFC value above zero and an adjusted p-value below 0.01.

According to a previous work (DeLaughter et al., 2016), PCA profiles for each time point of the cell types in human hearts were produced first and a density map based on the principal component PC1 that could better reflect the developmental stage was generated. Then, the loading value of PC1 was calculated. After that, the loading value and the expression of homologous genes in mouse corresponding cell types were used to calculate the PC1 value of the mouse. Finally, the PC1 values of the mouse were mapped to the human density plot.

Stringent quality control was performed using Seurat on the abovementioned three datasets (Supplementary Figures S1A–F). 3,757, 2,402, and 9,757 high-quality cells were further analyzed (Satija et al., 2015). To ensure the integrity of the mouse dataset at time points and comparability with the human dataset, the two mouse datasets were then sampled and integrated. We first used a geometry-based sampling approach for mouse dataset2, from which we took a comparable number of cells to dataset1 for subsequent integration (Hie et al., 2019) (Supplementary Figure S1G). After using Harmony integration (Korsunsky et al., 2019), cells from both datasets were overlapped well together (Supplementary Figure S1H).

Here, we used the same software and strategies to analyze the human and mouse datasets. After dimensional reduction, doublets removal, and unbiased clustering, 20 cell clusters were identified in human hearts. According to the expression status of marker genes, we further classified cells into major cell types, including 5 cardiomyocyte clusters (CMs) with high expression of cardiomyocyte-specific marker genes, such as TTN, MYH6, and TNNT2; 5 fibroblast clusters (FBs) with high expression of fibroblast-specific marker genes, such as COL3A1, COL1A2, and FN1; 4 endothelial cell/fibroblast clusters (ECs/FBs) with co-expression of endothelial cell-specific marker genes (NKX2-5, CDH5, and EMCN) and fibroblast-specific marker genes; 2 epicardial cell clusters (EPs) with high expression of epicardial cell-specific marker genes (UPK3B, MSLN, and WT1); macrophage cluster (MACs) with high expression of MS4A4A and SEPP1; macrophage/monocyte cluster (MACs/MONOs) with high expression of both macrophage marker genes and monocyte marker genes (CD68, LYZ, S100A8, and S100A6); natural killer cells cluster (NKTs) with high expression of NKG7, GNLY, and a cluster of blood cells (Figures 1A and B). A total of 17 clusters were identified in mouse hearts, including 5 CM clusters with high expression of Tnnt2, Tnnc1, and Actn2; 5 FB clusters with high expression of Col3a1, Col1a2, and Fn1; 4 EC/FB clusters co-expressing endothelial cell maker and fibroblast marker, such as Pecam1, Cd93, Tek, and Cdh5; 1 epicardial/fibroblast cluster (EPs/Fibroblasts) and a cluster of MACs/MONOs and a blood cells cluster (Figures 1A and B). Major cell types, including CMs, FBs, ECs, and immune cells (MONOs/MACs), were identified in both human and mouse hearts, indicating overall relative conservation of cell types between the two species.

Based on the defined cell types, we further identified cell-type–specific genes by performing differential expression analysis between each cell type and the others. Then GO enrichment analysis (biological process, BP) was performed on the DEGs (Supplementary Table S1). As shown in Figure 1C, in the two species, the cell-type-specific genes of CMs were mainly enriched with muscle tissue development and cardiac morphogenesis. Cell-type-specific genes of FBs were mainly associated with the organization of extracellular matrix, actin filament, and collagen fibril. In the EC/FB populations, in addition to functions such as extracellular matrix organization and cell-substrate adhesion, some functions such as angiogenesis regulation and positive regulation for vasculature development were also enriched. Since FBs exert an important role in cardiac development through the deposition of collagen and other extracellular matrix (ECM) components, and ECs involve maintaining vascular homeostasis and cardiac function, the functional enrichment analysis shows a well-characterized biological profile for each cell type in this work.

Based on the defined cell types in human and mouse (Figure 2A), we further explored the expression of marker genes in each cell type in both species. We found that some previously-unreported genes conserved expression in both human and mouse, such as CSRP3, ENO3, and COX6A2 in CMs; PTN, LTBP4, and MFAP4 in FBs; KDR and ESAM in ECs/FBs. Meanwhile, some species-specific marker genes, such as CKMT2 and MYZAP, which were only expressed in human CMs, and ATP1B1 was only expressed in mouse CMs. As for FBs, MATN2, and COL12A1 genes were highly expressed only in humans, while PAPSS2 expression showed a specificity with mouse (Figure 2B).

To understand whether the human concordant cell types are similar to mouse ones, an analysis based on homologous genes was performed. We found a high correlation between cardiomyocytes population, indicating that identical cell types in two species are more relevant than non-identical ones (Figure 2C).

Next, we observed the composition of the major cell types at different developmental stages in human and mouse hearts. The result showed that the proportion of CMs was greatly reduced, but the proportion of non-cardiomyocytes was increased in late developmental stages compared to early developmental stages (Figure 2D).

In order to further understand the conservation and differences of cellular composition and gene expression pattern at various developmental stages of human and mouse hearts, a comparative analysis of gene expression was applied at different time points for each cell type in both species. Although the process of heart development is synergized by multiple cell lineages, the function of cardiomyocytes is more crucial to maintain heart development (Kelly et al., 2014). Therefore, we first further subpopulated the cardiomyocyte population. An unbiased clustering was used to classify cardiomyocytes into five subpopulations: CM0, CM1, CM2, CM3, and CM4 (Figure 3A, Supplementary Figure S2A). To figure out the progressive molecular changes of CMs during the maturation process, we performed a pseudotime analysis on these cardiomyocyte subpopulations to reconstruct the trajectory of cardiomyocyte development (Figure 3B, Supplementary Figure S2B). In this process, only the genes expressed by at least 5% of the CMs were retained. We used the DEGs between time points as sorted gene sets for dimension reduction and trajectory analysis. For both human and mouse CMs, a pseudotime trajectory was constructed starting from 5 to 9 W (human) and e8.5-e9.5 (mouse) cells and gradually progressing to 23–25 W (human) and e13.5-e14.5 (mouse) cells, respectively. Since selected the dataset with developmental stage information for cell labeling, we found that the pseudotime trajectories colored by developmental time points could reasonably reflect the process of CMs development (Figure 3B, Supplementary Figure S2B). In order to investigate the potential role of CMs subpopulations in their development and in their functional stage, we next observed the expression profile of different CMs subpopulations in the pseudotime trajectory. The density curves of cell subpopulation expression generated from the trajectory showed that in human CMs, hCM4 subpopulation was aligned at the beginning of the trajectory, whereas hCM1 and hCM2 ones mainly appeared at the end of the trajectory. While both mCM1 and mCM3 subpopulations were aligned at the start point of the trajectory and mCM0 one was presented around the trajectory end (Figure 3C), suggesting that hCM4, mCM1, and mCM3 may exert their function in early cardiac development, while hCM1, hCM2, and mCM0 may play a role in later cardiac development. Along with heart development, a series of genes with conserved expression patterns were detected in cardiomyocyte subpopulations (Figure 3D), for example, TNNI3 and MYL2 remain highly expressed, LGALS1 and COX7A1 were significantly and progressively upregulated, and both ANXA3 and SCD were found to be downregulated along the trajectory. Consistent with this, the subpopulations of hCM4, mCM1, and mCM3 that were arranged at the beginning of the developmental trajectory also show relatively high expression of ANXA3 and SCD; the subpopulations of hCM1, hCM2, and mCM0 that were arranged at the end of the developmental trajectory also had relatively high expression of LGALS1 (Supplementary Figure S2C). These significantly expressed genes are really crucial in maintaining heart development. For example, LGALS1, as a member of the highly conserved β-galactoside-binding lectin family, has been identified to play an essential role in the control of inflammation and neovascularization (Camby et al., 2006; Liu and Rabinovich, 2010). Gal-1 encoded a major component of the contractile apparatus of cardiomyocytes involved in the regulation of cardiac function (Dias-Baruffi et al., 2010). As a mitochondrial function-related gene, COX7A1 is necessary for both systolic and diastolic functions (Huttemann et al., 2012). In addition, ANXA3 can protect the myocardium from ischemic injury and reduce infarct size (Meng et al., 2019). SCD encodes a central enzyme that synthesizes monounsaturated fatty acids in lipid metabolism. Recent evidence shows that SCD can reprogram cardiac metabolism, thereby regulating cardiac function (Dobrzyn et al., 2015). In the heart, SCD1 deficiency enhances glucose transport and metabolism at the expense of fatty acid uptake and oxidation (Izquierdo-Garcia et al., 2018). The metabolic changes associated with SCD1 deficiency really protect cardiac myocytes against both necrotic and apoptotic cell death and improve heart function.

Based on the developmental trajectory of CMs, we divided the affected-trajectory genes into three gene clusters that are highly expressed at the initial, transitional, and terminal stages of the trajectory. The functional enrichment of different gene clusters closely reflects the molecular signature of CMs from the immature to the mature state (Figure 3E, Supplementary Table S2). GO analysis revealed that in both human and mouse hearts, high expression at the beginning of the trajectory gene (HEBTG) was enriched in functions such as regulation of mitotic nuclear division and cell cycle (Figure 3F, Supplementary Table S2), while high expression at the end of the trajectory gene (HEETG) was enriched in the functions, including muscle system development and muscle contraction (Figure 3F) (Jouandin et al., 2014).

A temporal evolution of CM gene expression has been analyzed in both ventricular and atrial cardiomyocytes in mouse from e9.5 to P21 (DeLaughter et al., 2016). In a similar way, we also observed a temporal evolution of gene expression in hCMs through PCA analysis at all time points based on the genes expressed in at least 75% of hCMs (Figure 3G). PCA component PC1 separated cells in overlapping, stepwise clusters ordered by age. To identify the most similar time points in human and mouse CMs, a density map was generated according to a reported method (DeLaughter et al., 2016) from the most informative components based on PCA (Figure 3H). MCMs gene expression in different developmental stages was then mapped to the human density plot using genes expressed in at least 75% of mCMs. We found that the expression data from e8.5 mCMs corresponded to 5 W hCMs, while data from e10.5 mCMs corresponded to 9 W hCMs (Figure 3H).

Cardiac FBs have central effects on normal cardiac physiology and cardiovascular disease. The cells also play an important role in development via depositing collagen and other extracellular matrix (ECM) components (Tallquist, 2020). Therefore, we further analyzed the dynamics of the cell transcriptome during cardiac development. Unbiased clustering was utilized to classify human and mouse fibroblasts into five subpopulations: FB0, FB1, FB2, FB3, and FB4 (Figure 4A, Supplementary Figure S3A). The trajectory of FBs development was reconstructed to characterize the progressive molecular changes in the cells with the process (Figure 4B, Supplementary Figure S3B). Similarly, for human and mouse FBs, a temporal trajectory was constructed for each cell from early developmental stages to the later stages. To analyze the potential effects of different FBs subpopulations at the stages of fibroblast development, we next observed the expression status of different FBs subpopulations in the pseudotime trajectory. In hFBs, subpopulation hFB2 was aligned at the beginning of the trajectory, while subpopulation hFB1 appeared mainly at the end of the trajectory. For mFBs, mFB2 subpopulation was aligned at the beginning of the trajectory, and subpopulations of mFB3 and mFB4 appeared mainly at the end of the trajectory (Figure 4C). We also found that some genes with conserved expression patterns in FB subpopulations, such as BTG2, RGS5, and TIMP3, were significantly upregulated along the developmental process (Figure 4D). Consistently, subpopulations hFB1, mFB3, and mFB4 arranged at the end of the developmental trajectory also showed relatively high expression of RGS5 and BTG2 (Supplementary Figure S3C). BTG2 encode an anti-proliferative protein of the BTG/Tob family that has been proved to inhibit proliferation and promote terminal differentiation of neurons (Farioli-Vecchioli et al., 2009) and skeletal muscle progenitors (Feng et al., 2007). RGS5 can promote arterial growth during arteriogenesis (Arnold et al., 2014), and TIMP3 may stimulate fibroblast proliferation and phenotypic differentiation into myofibroblasts at cardiac tissue injury sites (Lovelock et al., 2005). Also, a series of genes, such as LUM, showed a trend of being first being upregulated and then downregulated during both human and mouse heart development. (Figure 4D). Lumican (LUM) localized to the ECM as a keratan sulfate small leucine-rich proteoglycan that directly binds to fibrillar collagen and regulates collagen fibrillogenesis production in some tissues such as cornea, tendon, and skin, also is a central pro-fibrotic molecule in the heart (Mohammadzadeh et al., 2020). Moreover, we next divided the genes affecting the FBs development into three gene clusters that are highly expressed at the initial, transitional, and terminal stages of the trajectory (Figure 4E, Supplementary Table S3). The GO analysis of DEGs of FBs along the trajectory showed that HEETG was enriched in the function of extracellular matrix organization in both human and mouse hearts (Figure 4F, Supplementary Table S3), while HEBTG was enriched with cell cycle function. It should be noted that, HEBTG was also enriched in Notch and p53 signaling pathways, while HEETG was enriched in the BMP signaling pathway. These results suggest that the functional enrichment of DEGs along with trajectories closely reflects the molecular features of progressive maturation for FBs.

We next extended our inspection to FBs to observe the corresponding developmental stages of human and mouse hearts using the same strategy as CMs. We found that the expression data from e10.5 mFBs corresponded to 6 W hFBs, the data of e13.5 mFBs corresponded to 7 W hFBs, and the e14.5 -mFBs data corresponded to 9W hFBs, respectively (Figure 4G).

ECs are crucial players in maintaining vascular homeostasis and cardiac function. We identified EC/FB clusters that express both endothelial and fibroblast markers in both human and mouse. Unbiased clustering was applied to divide the human and mouse EC/FB populations into four subpopulations: EC/FB0, EC/FB1, EC/FB2, and EC/FB3 (Figure 5A; Supplementary Figure S4A). Similarly, we also constructed a temporal trajectory for human and mouse ECs/FBs each with a progression from early to later developmental stages (Figure 5B; Supplementary Figure S4B). In hECs/FBs and mECs/FBs, subpopulations of hEC/FB0 and mEC/FB1 were alighted mainly at the beginning of the trajectory, while subpopulations of hEC/FB2, hEC/FB3, and mEC/FB2 were arranged at the end of the trajectory, suggesting that hEC/FB0 and mEC/FB1 may play a role at early stage of heart development; whereas hEC/FB2, hEC/FB3, and mEC/FB2 may have a function at a slightly later stage during the development (Figure 5C). A series of genes showed conserved expression patterns in the EC/FB subpopulation during the development (Figure 5D). For example, the expression of CD36, FABP4, and ID1 was gradually upregulated during cardiac development, while genes such as FABP5 and HAPLN1 remained highly expressed throughout the development. Consistently, subpopulations hEC/FB2 and mEC/FB2 arranged at the end of the developmental trajectory also had relatively high expression of CD36 and FABP4 (Supplementary Figure S4C). These significant-expression genes are involved in many physiological functions, for example, CD36 is closely related to fatty acid transport (Brinkmann et al., 2002), FABP4 plays a key role in fatty acid uptake in the heart and skeletal muscle (Iso et al., 2013), and the ID family (including ID1) is essential for early heart formation. These results indicate that the representative subgroups hEC/FB2, hEC/FB3, and mEC/FB2 in both human and mouse really exert important effects in late cardiac development, especially in energy metabolism. Moreover, HAPLN1 is an important regulator of ECM interactions during development (DeLaughter et al., 2013), and its expression pattern during cardiac development is consistent with extracellular matrix remodeling in the formation of ventricular trabeculae during early cardiac development.

The functional enrichment of DEGs of ECs/FBs along with the trajectory also revealed that (Figure 5E) HEETG were enriched in the extracellular matrix organization, regulation of endothelial cell proliferation, and endothelial cell proliferation (Figure 5F, Supplementary Table S4), and HEBTG was enriched in fibroblast proliferation, positive regulation of ERK1 and ERK2 cascade, negative regulation of Wnt signaling pathway, stem cell differentiation, and other functions. It should be noted that ERK1 and ERK2 are central molecules to regulate cell division, growth, and survival (Sugden et al., 2011). Similar to CMs and FBs, functional enrichment of DEGs along temporal trajectories also closely reflects the transition of EC/FB populations from immature to mature states. The expression data from e9.5 mECs/FBs corresponded to 10 W hECs/FBs is shown in Figure 5G.