Abstract Background Single‐cell RNA sequencing is widely used in cancer research and organ development because of its powerful ability to analyze cellular heterogeneity. However, its application in cardiomyocytes is dissatisfactory mainly because the cardiomyocytes are too large and fragile to withstand traditional single‐cell approaches. Methods and Results Through designing the isolation procedure of neonatal mouse cardiac cells, we provide detailed cellular atlases of the heart at single‐cell resolution across 4 different stages after birth. We have obtained 10 000 cardiomyocytes; to our knowledge, this is the most extensive reference framework to date. Moreover, we have discovered unexpected erythrocyte‐like cardiomyocyte–terminal cardiomyocytes, comprising more than a third of all cardiomyocytes. Only a few genes are highly expressed in these cardiomyocytes. They are highly differentiated cardiomyocytes that function as contraction pumps. In addition, we have identified 2 cardiomyocyte‐like conducting cells, lending support to the theory that the sinoatrial node pacemaker cells are specialized cardiomyocytes. Notably, we provide an initial blueprint for comprehensive interactions between cardiomyocytes and other cardiac cells. Conclusions This mouse cardiac cell atlas improves our understanding of cardiomyocyte heterogeneity and provides a valuable reference in response to varying physiological conditions and diseases. Keywords: cardiac cell atlas, cardiomyocyte heterogeneity, single‐cell RNA sequencing, terminal cardiomyocyte Subject Categories: Basic Science Research, Genetics __________________________________________________________________ Nonstandard Abbreviations and Acronyms aCM atrial cardiomyocyte BMP bone morphogenetic protein MIF migration inhibitory factor scRNA‐seq single‐cell RNA sequencing tCM terminal cardiomyocyte UMI unique molecular identifier vCM ventricle cardiomyocyte Research Perspective. What Is New? * We have constructed the single‐cell transcriptome of thousands of cardiomyocytes; to our knowledge, this is the most extensive reference framework to date. * We have discovered unexpected erythrocyte‐like cardiomyocyte–terminal cardiomyocytes, comprising more than a third of all the cardiomyocytes. * We provide a substantial blueprint for comprehensive interactions between cardiomyocytes and other cardiac cells. What Question Should Be Addressed Next? * The terminal cardiomyocytes were newly identified erythrocyte‐like cardiomyocytes, and there are 2 questions to be addressed: “What is their origin?” and “What are the potential clinical applications of these cells?” The heart is the first organ to develop in the embryo that works without rest, ceaselessly beating >2 billion times during the human lifetime.[31] ^1 Given that much of cardiac research has been focused on the cardiomyocyte, knowledge of the entire repertoire of cardiomyocytes is still far from enough because of cellular and functional heterogeneity.[32] ^2 , [33]^3 Single‐cell RNA sequencing (scRNA‐seq), a powerful technology to solve heterogeneity problems, is widely used in cancer research and organ development.[34] ^4 , [35]^5 , [36]^6 , [37]^7 Unfortunately, the cardiomyocytes are too fragile and large.[38] ^2 Therefore, they are unable to withstand traditional single‐cell approaches. As a result, most studies have tended to an alternative strategy: single‐nucleus RNA sequencing.[39] ^8 , [40]^9 , [41]^10 For example, the map of the dynamic transcriptional landscape of 5 distinct cardiomyocyte populations in healthy, injured, and regenerating mouse hearts was constructed by single‐nucleus RNA sequencing.[42] ^10 Likewise, through characterizing the transcriptome of 880 000 nuclei of ventricular cells, Reichart et al have identified 5 subgroups of ventricular cardiomyocytes.[43] ^9 However, single‐nucleus RNA sequencing captured a diversity of types not represented in the scRNA‐seq data set because the nucleus RNA is immature RNA.[44] ^11 Therefore, there is a dire need for a comprehensive analysis of single‐cell expression from high‐quality cardiomyocytes. By designing the isolation procedure of neonatal mouse hearts, we analyzed the whole heart from 4 different developmental stages, days 1, 2, 4, and 6 after mice were born, at single‐cell resolution. Through collecting >6 heart samples in each sample, we have obtained >10 000 cardiomyocytes; to our knowledge, this is the largest reference framework to date. Notably, we have identified a new subgroup of cardiomyocytes: terminal cardiomyocytes (tCMs). We uncovered cell constitution and reconstitution in the heart and provided strategic hints in response to varying physiological conditions and diseases. Methods scRNA‐seq data that support the findings of this study have been deposited in the Gene Expression Omnibus under accession code [45]GSE232466 ([46]https://www.ncbi.nlm.nih.gov/search/all/?term=%20GSE232466). The other data that support the findings of this study are available from the corresponding author on reasonable request. Animals A total of 80 C57BL/6 mice (days 1, 2, 4, and 6 after birth) were purchased from Shanghai Slake Experimental Animal Co, Ltd. All animal experiments in this study have been approved and conducted under the oversight of the Shanghai University of Medicine and Health Sciences. Sex was not determined for neonatal pups in the study. Hematoxylin and Eosin Staining For hematoxylin and eosin staining, the heart tissues were fixed with 4% paraformaldehyde, then embedded in paraffin and sectioned (6 μm). The tissues were stained with hematoxylin and eosin and imaged under high (×400) magnification using a light microscope (Leica, Germany) to determine the extent of tissue damage. Immunofluorescence Staining The paraffin sections were deparaffinized in water, and then antigen retrieval was performed without drying the sections during this process. After natural cooling, the slides were detained. After the slices were dry, we used a histochemical pen to draw circles around the tissue. BSA was added dropwise to seal, and then primary antibodies against titin (ab307447, Abcam), glycophorin A (bs‐10017R, Bioss), ribosomal protein S3 (ab248249, Abcam) prepared with PBS in a certain proportion were added dropwise. Sections were then incubated overnight at 4 °C in the box. After the primary antibodies were washed with PBS, we added the secondary antibody and incubated at room temperature for 50 minutes in the dark. After washing with PBS, the sections were incubated with 4′,6‐diamidino‐2′‐phenylindole staining solution at room temperature in the dark; then autofluorescence quencher was added dropwise for 5 minutes, and washed with running water for 10 minutes. Finally, the sections were mounted, the sections were observed under a fluorescence microscope, and images were collected. Fluorescence In Situ Hybridization After the tissue was taken out and washed, it was immediately placed in fixative for >12 hours. After fixation, the samples were dehydrated with graded alcohol, dipped in wax, and embedded. The paraffin was sliced by a microtome, picked out by a slicer, and baked in a 62° C oven for 2 hours. Paraffin sections were deparaffinized to water, and finally soaked in diethylpyrocarbonate water. The slices were then boiled in the repair solution for 10 to 15 minutes and allowed to cool naturally. According to the different index characteristics of different tissues, proteinase K (20 μg/mL) was added dropwise to digest at 37 °C, and after cleaning, the prehybridization solution was added dropwise and incubated at 37 °C for 1 hour. Then, the probe‐containing hybridization solution was added dropwise, and the slides were allowed to hybridize overnight in an incubator. The hybridization solution was washed off, then 4′,6‐diamidino‐2′‐phenylindole staining solution was added dropwise, and the slides were incubated in the dark for 8 minutes. After washing, antifluorescence quenching mounting medium was added dropwise to mount the slides. Sections were observed under a Nikon upright fluorescence microscope, and images were collected. scRNA‐seq of Neonatal Mouse Hearts For each scRNA‐seq sample, 6 to 9 heart samples of each group from neonatal mice were dissociated into single‐cell suspensions. The tissues were minced with a razor blade, and then digested in a digestion buffer with 0.75% Collagenase/Dispase (Roche, 11097113 001). The cell suspensions were filtered with cell strainers, and centrifuged. Subsequently, the cell pellet was resuspended with PBS and treated with 1× Red Blood Cell Lysis Solution (Thermo Fisher, catalog No. 00‐4333‐57) to remove red blood cells. Trypan blue (Thermo Fisher, catalog No. 14190144) was used to check the viability of cells. Then, 10× Library Preparation and Sequencing Beads with unique molecular identifiers (UMIs) and cell barcodes are loaded to near saturation so that each cell is paired with a bead in a gel bead emulsion. After exposure to cell lysis buffer, polyadenylated RNA molecules hybridize to the beads. Beads were collected in 1 tube for reverse transcription. In cDNA synthesis, each cDNA molecule is labeled at the 5′ end (ie, the 3′ end of the messenger RNA transcript), and the UMI and cell tag indicate its cellular origin. Briefly, 10× beads are followed by second‐strand cDNA synthesis, adapter ligation, and universal amplification. Sequencing libraries were prepared using randomly interrupted whole‐transcriptome amplification products, enriched for 3′ ends of transcripts linked to cellular barcodes and UMIs. All remaining steps, including library construction, were performed according to the standard manufacturer's protocol ([47]CG000206 RevD). Sequencing libraries were quantified using the High Sensitivity DNA Chip (Agilent) and Qubit High Sensitivity DNA Assay (Thermo Fisher Scientific) on the Bioanalyzer 2100. Libraries were sequenced on a NovaSeq6000 (Illumina) using 2×150 chemistry. scRNA‐seq data processing reads were processed using the Cell Ranger 2.1.0 pipeline with default and recommended parameters. FASTQs generated from Illumina sequencing output were aligned to the mouse genome GRCm38 using the STAR algorithm.[48] ^12 Next, a gene‐barcode matrix was generated for each sample by counting UMIs and filtering non–cell‐associated barcodes. Finally, a gene barcode matrix containing barcoded cell and gene expression numbers was generated. This output was imported into the Seurat, version 3.0.2, R toolkit for quality control and downstream analysis of scRNA‐seq data.[49] ^13 All functions are executed with default arguments unless otherwise specified. Cells with <200 or >6000 genes detected were excluded (each gene must align with at least 1 UMI in at least 3 cells). Calculation of mitochondrial gene expression was performed using the percentage feature set function Seurat Bag. To eliminate cells with low viability, cells with >10% mitochondrial gene expression were excluded. Data normalization (Seurat's Normalize Data function package) was used to extract a subset of mutable genes. Variable genes were identified, controlling for a strong relationship between variability and mean expression. We then use Find Integration Anchors and integrate data from the Seurat package to identify anchors between data sets, and then integrate data from different samples.[50] ^13 , [51]^14 We then performed a principal component analysis and, after scaling the data, reduced the data to the top 30 principal component analysis components. The louvain algorithm was used to perform cluster analysis on the normalized data. Visualize the cluster results on a 2‐dimensional map generated using t‐distributed stochastic neighbor embedding or uniform manifold approximation and projection.[52] ^15 , [53]^16 Cell type and subtype identification are obtained by t‐distributed stochastic neighbor embedding. For subclustering, we apply the same steps for scaling, dimensionality reduction, and clustering on a given data set (usually limited to 1 type of cell). For each cluster, the Wilcoxon rank‐sum test was used to find well‐expressed significant genes compared with the rest of the cluster. Single R and known marker genes were used to identify cell types.[54] ^17 Cell‐Cell Communication Analysis Cell‐cell communication analysis was conducted with the scRNA‐seq data by using Cellchat ([55]http://www.cellchat.org/cellchatdb/).[56] ^18 Only receptors and ligands expressed in >10% of cells of any type from the different groups were evaluated, whereas cell‐cell communication was considered nonexistent if the ligand or the receptor was undetectable. Averaged expression of each ligand‐receptor pair was analyzed between different cell types, and only those with P<0.01 were chosen for the prediction of cell‐cell communication between any 2 cell types. Pseudotime Trajectory Analysis The R package monocle2, version 2. 14. 0,[57] ^19 was used to conduct the pseudotime series analysis. Briefly, differential expression analysis was performed to identify the top significantly differentially expressed genes (false discovery rate, <0.05) between healthy and disease conditions to build the disease trajectory, and then every single cell was assigned a numeric pseudotime value and then ordered along the disease trajectory. Transcription factors and cytokines are included in the analyses if their expressions were detected in at least 10% of cells within a cluster and are significantly differentially expressed (false discovery rate, <0.25) along the pseudotemporal relationship among cells. Analysis of Gene Differential Expression and Enrichment Analysis Wilcoxon algorithm was used to analyze all marker genes of each cell type, and group 1 versus rest was used to score marker genes. Gene Ontology pathway enrichment analysis of differentially expressed genes was performed by cluster Profiler, version 3.14.0, package. Results Overview of the Cell Composition of the Neonatal Mouse Hearts Efficient scRNA‐seq of primary cardiomyocytes is hindered by the large size and fragility of these cells. To overcome the size limitation of the 10X Chromium Controller, we have performed 3 things. First, we chose the neonatal mouse heart whose cardiomyocytes are small.[58] ^20 Indeed, most of the cardiomyocytes at day 1 are small than (50 μm), which is small enough to pass through the 10X Chromium Controller. Second, we used the entire heart instead of specific parts, such as the left ventricle, thus obtaining all types of cardiomyocytes. Furthermore, each sample contained >5 whole hearts, resulting in a collection of millions of cardiomyocytes, sufficient for processing through the 10X Chromium Controller (Figure [59]1A). The project comprised 4 groups: days 1, 2, 4, and 6. Day 1 incorporated 9 hearts, and the others contained 6 hearts each. Using this approach, we successfully collected a total of 50 955 single cells (Figure [60]S1A) from the 4 samples (Table [61]S1). A median depth of 56 487 reads and 2776 genes were detected per cell (Figure [62]S1B). Previously, 23 clusters were acquired (Figure [63]S2A and S2B). Subsequently, we annotated each cluster according to the known cell subtypes and their markers. We classified the entire population into 14 major cell clusters: cardiomyocyte, endothelial cell, fibroblast, smooth muscle cell, B lymphocyte, lymphoid, macrophage, myeloid, neutrophil, mesothelial cell, pericyte, neuronal cell, early erythrocyte, and proliferating cell, according to their respective molecular features (Figure [64]1B). Fortunately, a large number of cardiomyocytes were identified (Figure [65]1B). Indeed, cardiomyocyte, endothelial cell, and fibroblast were the major subgroup, counting for almost 83% of the total cells, which is familiar with previous study[66] ^21 (Figure [67]1C). Intriguingly, early erythrocyte, a member of erythroid lineages, was identified (Figure [68]1B and [69]1C). Notably, we found that early erythrocytes were close to a portion of cardiomyocytes, suggesting some similarity between them. The expression of representative biomarkers of each cluster was displayed (Figure [70]1D and Figure [71]S3A). Afterward, we validated the expression of cell markers titin (cardiomyocyte) and glycophorin A (early erythrocyte) by immunofluorescence (Figure [72]1E). Cell number of cardiomyocytes decreased sharply during this period (Figure [73]1F and [74]1G), probably because of the continuing growth of cardiomyocyte size (Figure [75]S3B). Figure 1. A single‐cell atlas of neonatal mouse heart. Figure 1 [76]Open in a new tab A, Schematic of workflows for single‐cell RNA sequencing (scRNA‐seq) of neonatal mouse heart. B, Uniform manifold approximation and projection (UMAP) of 50 955 sequenced cardiac cells split by cohort. C, Pie map shows the proportion of each cell population in the cardiac cells. D, Marker genes for each cell population in the cardiac cells. E, Immunofluorescence (IF) of cell type–specific markers in atrium and ventricles: titin (TTN) (red, cardiomyocytes) and glycophorin A (GYPA) (red, early erythrocyte). Nuclei are counterstained with 4′,6‐diamidino‐2′‐phenylindole (DAPI) (blue). Scale bars, 50 μm. F, UMAP shows the changes of each cell population during different stages after birth. G, The relative variation of the proportion of each cell population during different stages after birth. A New Major Subgroup of Cardiomyocyte As cardiomyocytes play a crucial role in the heart's function, we directed our subsequent research toward them. Applying unsupervised clustering to partition all the cardiomyocytes, we found the cardiomyocytes were roughly divided into 3 subgroups: cardiomyocyte I, cardiomyocyte II, and cardiomyocyte III (Figure [77]2A). Myl7 and Nppa, 2 classic atrial cardiomyocyte markers, were highly expressed in cardiomyocyte III, indicating that cardiomyocyte III is the atrial cardiomyocyte (Figure [78]2B). Similarly, Myh7 and Myl2, 2 typical ventricle cardiomyocyte markers, were highly expressed in cardiomyocyte II, suggesting that cardiomyocyte II is the ventricle cardiomyocyte (Figure [79]2B). Notably, the remaining cluster, cardiomyocyte I, counting more than a third of total cardiomyocytes, decreased sharply during the development, confirming that they were cardiomyocytes because their number decreased sharply during development (Figure [80]2C). In addition, this result suggests that these cardiomyocytes grow faster than the other 2 clusters during this period (Figure [81]2C). More important, they were next to the early erythrocyte (Figure [82]2D), a highly differentiated cell with low RNA expression.[83] ^22 Although these cardiomyocytes are close to erthyocytes, their expression patterns are different. The erthyocytes highly expressed Gypa, Alas2, and Hbb‐bt genes, whereas the tCMs highly expressed Myl2, Tnnt2, and Myl3. Nevertheless, their similarities are obvious. The average number of genes in cardiomyocyte I was small, which was similar to early erythrocytes (Figure [84]2E). Indeed, the average number of genes in cardiomyocyte I was even lower than that in early erythrocytes (Figure [85]2E). Furthermore, the number of genes and UMIs in cardiomyocyte I was much lower than that in cardiomyocyte II and cardiomyocyte III (Figure [86]2E). In addition, several cardiomyocytes with high‐level proteins involved in contractile force generating showed a low mRNA level (Figure [87]2F), indicating that they were cardiomyocyte I. Just like early erythrocytes, cardiomyocyte I might be a highly differentiated terminal cell. Hence, we named it tCM. Figure 2. Overview of cardiomyocytes (CMs) reveals an unexpected erythrocyte‐like cardiomyocyte. Figure 2 [88]Open in a new tab A, Uniform manifold approximation and projection (UMAP) of clusters formed by all the cardiomyocytes. B, Marker genes for the cluster CM II and CM III in the cardiomyocytes. C, UMAP projection of the 3 cardiomyocyte subclusters during different times. D, UMAP projection shows the relative location of CM I and early erythrocyte. E, Left: comparison of CM I with early erythrocyte on the number of unique molecular identifiers (UMIs) and genes. Right: comparison of CM I with CM II and CM III on the number of UMIs and genes. F, Immunofluorescence (IF) and fluorescence in situ hybridization (FISH) of cardiomyocyte marker titin (TTN) in atrium and ventricles during different stages after birth: IF (green, cardiomyocytes) and FISH (red, cardiomyocytes). Nuclei are counterstained with 4′,6‐diamidino‐2′‐phenylindole (DAPI) (blue). Scale bars, 20 μm. tCM is a Highly Differentiated Cell To explore whether the tCMs are highly differentiated, we detected the character of cell proliferation. Cyclin‐dependent kinases 4 and 6 play critical roles in cell proliferation, and help to drive the G[1]/S transition.[89] ^23 , [90]^24 Intriguingly, we found that these 2 genes were expressed in the atrial and ventricle cardiomyocytes in the neonatal heart, whereas they were almost not expressed in tCMs (Figure [91]3A). Pcna and Kras, another 2 genes related to cell proliferation,[92] ^25 , [93]^26 , [94]^27 were also not expressed in tCMs, confirming the worst proliferative capacity of these cells (Figure [95]3A). Then, we analyzed the genes that were highly expressed in the tCMs. Consistent with previous results, most tCMs express <1000 genes (Figure [96]3B). Interestingly, through analyzing the expression of the 300 most highly expressed genes, we observed an entirely different expression pattern of the tCM (Figure [97]3C). We next analyzed their biological function of these highest expressed genes. We illustrated that the top 10 genes were enriched in muscle contraction; genes ranking in 11 to 50 were enriched in heart contraction, nicotinamide adenine dinucleotide metabolic process, and mitochondrial ATP synthesis coupled electron transport; and genes ranking 51 to 300 were involved in cardiac muscle cell development, ATP biosynthetic process, and fatty acid oxidation (Figure [98]3D). This result confirmed that the tCM functions as a contraction apparatus. When aligning atrial cardiomyocytes, ventricle cardiomyocytes, and tCMs in pseudotime, we observed a relatively old‐age distribution of tCMs, suggesting them as terminally differentiated cells (Figure [99]3E). Therefore, these results imply tCMs were highly differentiated cardiomyocytes that function as contraction pump (Figure [100]3F). Because RNA polymerase II plays key roles in RNA synthesis,[101] ^28 we detected the expression of its subunits. As expected, the expression levels Taf1 and Polr2a were low in a portion of cardiomyocytes (Figure [102]3G). And a portion of cardiomyocytes with a high level of titin show a low level of TAF1, indicating that they belong to the tCMs (Figure [103]3H). Figure 3. The characteristics of terminal cardiomyocytes. Figure 3 [104]Open in a new tab A, Plot map shows the expression of genes involved in cell proliferation in the cardiomyocytes. B, Plot map shows the correlation of unique molecular identifiers (UMIs) and genes in each cell of terminal cardiomyocytes (tCMs). C, Comparison of the relative expression of the highest 300 genes in the 3 cardiomyocytes. D, Gene Ontology (GO) analyses of the 300 most highly expressed genes in the terminal cardiomyocytes. E, Monocle analyses showing the ordering of atrium cardiomyocyte (top), ventricle cardiomyocyte (bottom), and terminal cardiomyocyte (middle) along pseudo‐time trajectories. F, BioRender infographic summarizes the predicted relationship of the 3 cardiomyocyte clusters. G, Plot map shows the expression of genes involved in RNA polymerase II in the cardiomyocytes. H, Immunofluorescence (IF) of TTN (titin) and TAF1 (TATA‐box binding protein associated factor 1) in the atrium and ventricles: TTN (red, cardiomyocyte marker) and TAF1 (green, RNA polymerase II marker). Nuclei are counterstained with 4′,6‐diamidino‐2′‐phenylindole (DAPI) (blue). Scale bars, 50 μm. nGene indicates number of genes; nUMI, number of unique molecular identifiers; and RPS3, ribosomal protein S3. Seven Subclusters of Atrial and Ventricle Cardiomyocytes Next, we examined atrial cardiomyocytes (aCMs) and ventricle cardiomyocytes (vCMs). A total of 7 cell clusters were identified following unsupervised clustering (Figure [105]4A). Interestingly, these cells were perfectly divided into 2 groups: aCMs located on the right side, and vCMs located on the left side (Figure [106]4A). Both aCMs and vCMs could be divided into 4 subclusters (Figure [107]4A). aCMs included aCM_Bmp10, aCM_Fabp5, aCM_Hcn4, and CM_Mki67; and vCMs contained CM_Mki67, vCM_Cntn2, vCM_Hmgcs2, vCM_Nppb (Figure [108]4A), and cell subclusters that belong to the vCMs decreased sharply during development (Figure [109]S4A). aCM_Hcn4, a subgroup of aCM, abundantly expressed conduction‐related genes, including Hcn4, Vsnl1, and Shox2 (Figure [110]4B and [111]4C and Table [112]S2), indicating that they belong to the sinoatrial node pacemaker cells.[113] ^29 On the basis of Gene Ontology enrichment analysis, the conserved markers of aCM_Hcn4 were enriched in 3 biological processes: cellular response to cAMP, regulation of heart rate by cardiac conduction, and cardiac conduction system development (Figure [114]4D). Similarly, vCM_Cntn2 also showed robust levels of molecules with conductional functions (Cntn2, Slit2, Sema3c, and Irx2) (Figure [115]4B and [116]4C), suggesting that they belong to the His‐Purkinje system.[117] ^30 Likewise, the conserved markers of vCM_Cntn2 were involved in conduction: cell‐cell signaling involved in cardiac conduction, cation transmembrane transport, and sodium ion export across the plasma membrane (Figure [118]4D). CM_Mki67, which existed in both in atrial and ventricle cardiomyocytes, showed enriched expression of Mki67, Cenpf, and Cdca8 (Figure [119]4B and [120]4C). The number of these proliferative cells gradually decreased during development (Figure [121]4E). Figure 4. Seven subclusters of atrial and ventricle cardiomyocytes. Figure 4 [122]Open in a new tab A, Uniform manifold approximation and projection (UMAP) of subclusters of atrial and ventricle cardiomyocytes. B and C, Marker genes for each subpopulation in the 7 subclusters. D, Gene Ontology (GO) analyses of the 30 most highly expressed markers in each cell subcluster of the 7 cardiomyocytes (the blue bar represents the fold enrichment, which refers to the left y axis; the orange line represents the ‐log [P value], which refers to the right y axis). E, Representative images of immunofluorescence (IF) with cardiomyocyte marker TTN (titin) and MKI67 (proliferation marker antigen identified by monoclonal antibody Ki 67) in atrium and ventricles: TTN (red) and MKI67 (green). Nuclei are counterstained with 4′,6‐diamidino‐2′‐phenylindole (DAPI) (blue). Scale bars, 50 μm. aCM indicates atrial cardiomyocyte; and vCM, ventricle cardiomyocyte. There are 2 major subgroups in the aCMs: aCM_Bmp10 and aCM_Fabp5. aCM_Bmp10 displayed functional enrichment in biological behaviors related to the regulation of heart contraction, cardiac muscle tissue development, and regulation of blood circulation, whereas aCM_Fabp5 exhibited enrichment in cardiac muscle contraction, response to abiotic stimulus, and mitogen‐activated protein kinase cascade. Furthermore, 2 major subgroups (vCM_Hmgcs2 and vCM_Nppb) were identified in the vCMs as well. Gene Ontology analysis revealed that vCM_Hmgcs2 was enriched in functions related to metabolism, including lipid oxidation, fatty acid β‐oxidation, and long‐chain fatty acid import into cell; on the contrary, vCM_Nppb was associated with oxidative phosphorylation, regulation of heart rate, and muscle hypertrophy. We next examined differentially expressed genes between them (Figure [123]S4B). Bmp10, a conserved marker of aCM_Bmp10, is predominantly distributed in the heart and preserves cardiac function.[124] ^31 In contrast, the transcription factor Pitx2, which is expressed in aCM_Fabp5, is implicated in atrial fibrillation predisposition.[125] ^32 Besides Nppb, Ankyrin repeat domain 1 was found predominantly within vCM_Nppb, playing crucial roles in cardiogenesis, mechanosensing, and intracellular signaling.[126] ^33 Cell‐Cell Interactions Then, we validated cell‐cell interactions with cellchat to identify receptor‐ligand circuits among all clusters of cells.[127] ^18 We identified that fibroblast cell populations were the most prominent sources of cell‐cell interaction (Figure [128]5A). On the contrary, early erythrocyte showed the fewest interaction with other cells, followed by tCMs (Figure [129]5A). Subsequently, we detected detailed communications for individual pathways. Five patterns for incoming signaling and 5 patterns for outgoing signaling were uncovered (Figure [130]5B). The incoming signaling, for example, reveals that a large portion of incoming aCM and vCM signaling was characterized by pattern 4, which represents multiple pathways, including but not limited to IGF (insulin‐like growth factor), EGF (epidermal growth factor), and GRN (granulin) (Figure [131]5B). On the other hand, the communication patterns of secreting cells shows that outgoing aCM and vCM signaling were dominated by pattern 1, which include signaling pathways such as BMP (bone morphogenetic protein) 10, WNT (Wingless/Integrated), as well as vascular endothelial growth factor (Figure [132]5B). This result reveals that signaling network in the heart is complex and highly redundant with multiple ligand sources and a large portion of target cells. Figure 5. Cell‐cell interaction. Figure 5 [133]Open in a new tab A, Predicted cell‐cell interactions with cellchat identify receptor‐ligand circuits among all clusters of cardiac cells (left: interaction weight/strength; right: number of interactions). The edge width is proportional to the indicated ligand‐receptor (L‐R) pairs. B, The inferred incoming communication patterns of the target cell (left) and the inferred outgoing communication patterns (right) of secreting cells, which show the correspondence between the inferred latent patterns and cell groups, as well as signaling pathways. The thickness of the flow indicates the contribution of the cell group or signaling pathway to each latent pattern. C, Violin plot shows the expression of genes involved in the inferred MIF (macrophage migration inhibitory factor) signaling network in each subgroup based on single‐cell RNA sequencing (scRNA‐seq) data. D, Heat map shows the relative importance of each cell group based on the computed 4 network centrality measures of the MIF signaling network. E, Relative contribution of each MIF ligand‐receptor pair. F, Hierarchical plot shows the inferred intercellular communication network for MIF signaling. It consists of 2 parts: the left portion highlights the MIF signaling to lymphoid and early erythrocyte clusters, and the right portion highlights MIF signaling to other cell clusters. Solid and open circles represent source and target, respectively. Circle sizes are proportional to the number of cells in each cell group, and edge width represents the communication probability. Edge colors are consistent with the signaling source. G, BioRender infographic shows predicted MIF signaling between cardiomyocytes and other cells. MIF is secreted by all the cells, and endothelial cells and immune cells, including macrophage, myeloid, and B lymphocyte, have the receptors. H, Predicted cell‐cell interactions with cellchat identify receptor‐ligand circuits among 7 clusters of cardiomyocytes (left: interaction weight/strength; right: number of interactions). I, Violin plot shows the expression of genes involved in the inferred BMP (bone morphogenetic protein) signaling network in each subgroup based on scRNA‐seq data. J, Hierarchical plot shows the inferred intercellular communication network for MIF signaling. K, BioRender infographic shows predicted BMP signaling between cardiomyocytes. BMP is accepted by all the cells, and aCM_Hcn4 and vCM_Nppb are the secretors. aCM indicates atrial cardiomyocyte; tCM, terminal cardiomyocyte; and vCM, ventricle cardiomyocyte. We next analyzed the MIF (macrophage migration inhibitory factor) signaling network. MIF signaling was at the top of incoming signaling (Figure [134]5B). Interestingly, approximately all of the 14 cardiac cells showed high MIF expression except early erythrocytes (Figure [135]5C). B lymphocyte and myeloid expressed all the receptors except Ackr3, whereas macrophage expressed all the receptors except Cxcr4 (Figure [136]5C). As expected, immune cells, including B lymphocyte, macrophage, and myeloid, are the most prominent mediators, suggesting their roles as gatekeepers of cell‐cell communication (Figure [137]5D). Subsequently, we detected the relative contribution of each ligand‐receptor pair to the overall communication network of MIF signaling pathway. MIF‐(Cd74+Cd44) was the majority of MIF interactions (Figure [138]5E). Furthermore, network centrality analysis of the MIF signaling network demonstrated that almost all the clusters were senders, whereas fibroblast, endothelial cell, and immune cells, including B lymphocyte, macrophage, and myeloid, were important target cells (Figure [139]5F). Collectively, MIF, secreted by all the cardiomyocytes and other cardiac cells, targets fibroblast, endothelial cell, and immune cells (Figure [140]5G). Then we analyzed the cell‐cell communication among the 7 cell clusters of cardiomyocytes. aCM_Bmp10 was the most prominent source of cell‐cell interaction (Figure [141]5H). Five patterns for incoming signaling and 5 patterns for outgoing signaling were identified (Figure [142]S5A and S5B). For example, we revealed that a large portion of incoming aCM_Bmp10 signaling was characterized by pattern 3, which represents EGF and GRN pathways (Figure [143]S5A). On the other hand, the communication patterns of secreting cells show that outgoing aCM_Hcn4 signaling was dominated by pattern 3, which include signaling pathways such as BMP and AGT (angiotensinogen) (Figure [144]S5B). Then, we focused on pivotal signaling on cardiovascular system–BMP signaling. In contrast to MIF signaling that ligand MIF was expressed in almost all clusters, and 2 receptors of BMP signaling, Acvr1 and Bmpr2, were highly expressed in all 7 clusters (Figure [145]5I). aCM_Hcn4 and vCM_Nppb were the prominent crossroads of this network (Figure [146]S5C). Bmp2‐(Bmpr1a+Bmpr2) was the majority of BMP interactions (Figure [147]S5D). Network centrality analysis of the BMP signaling confirmed that aCM_Hcn4 and vCM_Nppb played crucial roles as the signal sender, and all the clusters contained signal receptors, as shown in Figure [148]5J and Figure [149]5K. Discussion The study has constructed a comprehensive map of the transcriptional landscape in neonatal mouse hearts composed of scRNA‐seq for >50 000 cells. This work provides at least 2 novel advances to enhance our understanding of cardiovascular biology. First, we have generated the largest collection of single‐cell transcriptomes of cardiomyocytes to date. Furthermore, we identified an unexpected major cluster: tCMs. Second, with this robust data set, we defined 2 subclusters of cardiomyocytes, aCM_Hcn4 and vCM_Cntn2, which are hard to identify in the single‐nucleus RNA sequencing data set.[150] ^10 , [151]^34 Together, this cardiac cell atlas generated an analytic and statistical framework for our understanding of the heart and provided a valuable reference to the scientific community. Although heterogeneous cardiomyocytes have been extensively studied, studies have not been of sufficient scale to characterize them because of limited technical difficulties in isolating high‐quality primary cardiomyocytes.[152] ^35 After digesting and sequencing several neonatal mouse hearts together, we have successfully obtained tens of thousands of single‐cell transcriptomes of cardiomyocytes. Notably, the unexpected major subcluster, tCM, was identified. They are similar to the early erythrocytes, suggesting that they were highly differentiated cardiomyocytes. Moreover, tCMs have only one‐tenth of UMIs of other cardiomyocytes, indicating that they lost most of their cellular function, except as contraction apparatus, including contractile force generating, mitochondrial ATP synthesis, and fatty acid oxidation. Interestingly, Tucker et al have uncovered 2 particular cardiomyocyte clusters, which they removed from subsequent analyses because nonnuclear regions were the source of these cells.[153] ^1 We speculated that these 2 clusters might be similar to tCMs because of a lack of nuclear mRNA. Therefore, in future studies, it will be necessary to systematically determine the existence and function of tCMs in the human heart. Each heartbeat is triggered by an electrical impulse spontaneously released by the sinoatrial node and then conducted by the His‐Purkinje system.[154] ^36 However, the nature of conduction system cells (cardiomyocyte‐ or neuron‐like cells) remains incompletely known.[155] ^37 Besides neuronal‐like cells that highly expressed Gfra3, Plp1, and Kcna1, we have also discovered hundreds of specialized cardiomyocytes that encode contractile force‐generating sarcomere genes, including Ttn and Mybpc3, but highly expressed genes involved in conduction, such as Hch4 and Slit2. For example, aCM_Hcn4s were sinuatrial node pacemaker cell‐like cells that showed robust levels of Hcn4 and Vsnl1, whereas vCM_Cntn2s were Purkinje cell–like cells that showed enriched expression of Slit2 and Sema3c. These results may, at least in part, demonstrate that a portion of cardiomyocytes plays critical roles in the cardiac conduction system. In addition, Chen et al revealed that sinuatrial node pacemaker cells are glutamatergic neuron‐like cells in adult mice by single‐cell transcriptome analysis.[156] ^38 The probable reason for these differences may be that sinuatrial node pacemaker cells comprise both cardiomyocyte‐ and neuron‐like cells. Further studies on other species, such as monkeys and humans, may help to solve this puzzle. We recognize that this study was subject to several limitations. First, we have only investigated 4 different stages when most of the cardiomyocytes were mature, lacking the data of the development of cardiomyocytes and the origin of the tCMs. Second, the number of cardiomyocytes decreased sharply in the elder samples, and the characters of these missing cardiomyocytes is still unknown. Nevertheless, we expect that our results will inform studies of tCMs in other species (including rat and human), as well as studies with large cohorts to determine the origin and differentiation of tCMs, and afford pivotal insights to understand the cellular mechanism of heart diseases. Sources of Funding This study was supported by the Natural Science Foundation of Shanghai (21ZR1428400) and the Shanghai Medical Science and Technology Support Project (21S11901700). Disclosures None. Supporting information Tables S1–S2 Figures S1–S5 [157]Click here for additional data file.^ (782.2KB, pdf) Acknowledgments