We acquired the available snRNA-seq data derived from seven normal human heart samples from the Broad Institute's Single Cell Portal¹ under the Study ID SCP498 online. In which, “healthy_human_4chamber_map_unnormalized_V3.h5ad” was used in the downstream analyses. We used the gene expression matrix and cell type annotation provided by the original article (26).

Since the low-quality cells were excluded from the original publication, we used all the data available. The downloaded unique molecular identifier count matrix was converted to Seurat object using the R package Seurat v.3.2.2 (27). We then normalized the raw gene expression matrix using the NormalizeData function and visualize the expression level using the Violin plot function in Seurat. Uniform manifold approximation and projection (UMAP)² was used to visualize the results.

Downstream analyses, including highly variable gene detection (FindVariableFeatures, method= “vst,” nfeatures = 2,000), data feature scaling (ScaleData), PCA (Principal Component Analysis, RunPCA, from highly variable genes), an integration procedure (RunHarmony, npcs = 23), neighbor construction (FindNeighbors, using the output of RunHarmony), and clustering (FindClusters, resolution = 0.4), were performed using the R package Seurat (v.3.2.2) and harmony (v1.0). Cell type labels were assigned manually using the annotation provided in the original article (26). We further investigated clusters for which the genes indicated additional diversity and specificity.

Unless otherwise stated, tests for differentially expressed genes (DEGs) were performed using the FindAllMarker or FindMarker function in the Seurat package (with the default test method: Wilcoxon rank-sum test). To compute for the DEGs, all genes were included and expressed in at least 25% of cells in either of the two populations compared.

Cell–cell interaction weights for each ligand-receptor pair were computed as the product of the fold change of ligands in sender-cell types and the fold change of the corresponding receptors in receiver cell types. We used the “layout_nicely” function in the R package “igraph” to visualize the cell-cell interaction.

The R and Python scripts used to analyze the data in this study are available from the corresponding authors upon reasonable request (Table 1).

All statistical tests were performed using the R statistical programming language (V.3.6.2, R Foundation for Statistical Computing, Vienna, Austria). For analysis of snRNA-seq data, we used Wilcoxon rank-sum tests (between two groups) or Kruskal–Wallis test (for comparisons of multiple groups) to detect differential genes, with P-values adjusted based on Bonferroni correction. Statistical significance was set at P < 0.05.

We obtained single-nucleus transcriptome data of seven human donor hearts containing 287,269 nuclei sequenced by 10X Genomics from the previous study (26). The basic characteristics of the human donors are shown in Supplementary Table 1. As shown in Figure 1A, unsupervised graph-based clustering revealed 11 cell types in the human heart, of which multiple cell types expressed the ACE2 gene (Figure 1B). Clear differences in the proportion of ACE2-positive cells were found among the different cell types (Figure 1C). Furthermore, high percentages of ACE2-positive cells were observed in the pericytes (8.3%), fibroblasts (5.1%), and adipocytes (4.4%) in the heart (Table 2). These results are slightly different from those reported in the literature (17). In addition to pericytes, we predicted that fibroblasts and adipocytes could be vulnerable targets for SARS-CoV-2 infection. Similar to previous reports, the pericyte showed a significant intercellular interaction with the endothelium (Figure 1D). In addition, fibroblast demonstrated significant reciprocal interactions with vascular smooth muscle cells (Figure 1D). These results suggest that vascular injury plays an important role in the occurrence and development of cardiac injury in COVID-19. After binding to the cell surface, SARS-CoV-2 requires S protein cleavage proteases, such as TMPRSS2, FURIN, CTSL, and CTSB, to facilitate cell entry by inducing the fusion of cellular and viral membranes. TMPRSS2/4/11A, the well-known co-receptor for SARS-CoV-2 in other organs, demonstrated almost no co-expression (<1%) with ACE2 in various cell types of the human heart (Figure 1E; Supplementary Table 2). However, FURIN and CTSL/B seemed to play significant roles as co-receptors for SARS-CoV-2 entry into the human heart (Figure 1F). Further studies are required to elucidate the exact mechanism.

It has been proven that the male sex is a risk factor for disease severity in COVID-19 (28). The heart is composed of different chambers, which have distinct structural and functional characteristics. Hence, gaining a deeper understanding of the expression characteristics of SARS-CoV-2 cellular receptors in the human heart at spatial and sex-specific levels is of great importance. No significant differences in the overall distribution of cardiac cells between the sexes were observed (Figure 2A). The proportion of ACE2-positive cells was further compared between different genders in each cell type. As shown in Figure 2B, the ratio of ACE2-positive cells had a gender-differentiated effect in pericytes, fibroblasts, and the endothelium. However, it is interesting to note that even with a lower ACE2 proportion, higher expression levels of ACE2 were observed in pericytes, fibroblasts, and endothelium in male hearts (Figure 2C). This might be a reason why there is no significant difference observed in adverse cardiac events between patients with COVID-19 of different genders.

However, there were remarkable differences in the gene expression between the atria and ventricles of the human heart (Figure 2D). The frequency of ACE2-positive cells in each cell type from the ventricles was significantly higher than that in the atria (Figure 2E). For instance, the proportion of ACE2-positive cells in ventricular cardiomyocytes (2.1%) were significantly higher than that in the atrial cardiomyocytes (0.5%). This suggests that human ventricular cells are more susceptible to SARS-CoV-2 infection.

Previous studies have reported potential receptors involved in SARS-CoV-2 binding and infection other than the ACE2 (20–25). Each potential receptor exhibited a distinct expression pattern (Figure 3A). The human cardiac cells showed nearly no expression of ASGR1, while the other receptors (BSG, HSPA5, KREMEN1, NRP1, ANPEP, and AXL) showed scattered expression in each cell type (Figure 3A; Supplementary Tables 3–5). The expression distributions and correlations of the potential receptors and ACE2 were analyzed using feature plots. There was no significant association between the distribution of ACE2 and the potential receptors (Figures 3B,C). Dissimilarities among the expression characteristics of the potential receptors, compared with ACE2, might be the reason for the different clinical features of COVID-19 and SARS. Significantly different from the distribution pattern of ACE2, positive expression of BSG, KREMEN1, and NRP1 were highest in the ventricular cardiomyocytes (Figures 3D,E). These results may shed light on the susceptibility of human ventricular cardiomyocytes to SARS-CoV-2.

First, our results and conclusions are data-driven, which needed further laboratory verifications. Additionally, the original sample size was relatively small, which restricted its further implications. Finally, the snRNA-seq data we included was from normal human hearts, which could, therefore, limit further analysis of the correlation between gene expression features and clinical characteristics in patients with COVID-19.