The single-cell RNA-seq data of two normal left ventricle samples were collected from Gene Expression Omnibus (GEO) with accession number GSE134355 (Han et al., 2020). To identify cell types and key genes related to heart failure, we downloaded the single-cell RNA-seq data of two normal, four dilated cardiomyopathy (DCM), and two ischemic cardiomyopathy (ICM) hearts samples (accession number: GSE121893 (Wang et al., 2020)), one scRNA-seq data of one normal, two DCM and two ICM hearts (accession number: GSE145154 (Rao et al., 2021)) for validation, and bulk RNA-seq data of 14 non-failure (NF), 37 DCM, and 13 ICM samples from GEO database (accession number: GSE116250 (Sweet et al., 2018)). The RNA-seq data of fibroblasts induced by TGFβ1 and control samples, and the microarray-based gene expression data for validation were downloaded from GEO with accession numbers GSE97358 (Schafer et al., 2017) and GSE5406 (Hannenhalli et al., 2006), respectively.

The unique molecular identifiers (UMIs) count-based scRNA-seq data of the two normal left ventricle samples were used for the cell clustering analysis, which was implemented in R Seurat v3.2.3 package. Cells with less than 500 UMIs were eliminated and features detected in less than 3 cells were filtered. The two hearts were integrated using the anchors by Reciprocal PCA. The expression data was normalized by LogNormalize method with scale factor = 1000,000, and top 2000 highly variable features were selected by FindVariableFeatures with dispersion method. The clusters were found at a resolution of four by FindClusters, and T-distributed Stochastic Neighbor Embedding (t-SNE) was applied to reduce the dimensionality. The cell-type marker genes were detected by FindAllMarkers function at adjusted p-value < 0.05, minimal percentage >0.25, and log2 fold change >0.25. All the marker genes of the cell clusters were collected from the earlier study (Han et al., 2020). This analysis was implemented by R Seurat v3.2.3 package (Stuart et al., 2019).

The bulk RNA-seq data was downloaded from GEO database (GEO accession number: GSE116250 (Sweet et al., 2018)). The FPKM-based gene expression data were used for PCA analysis. Specifically, gene expressions higher than 1 FPKM in more than five samples were transformed to log2 (FPKM + 1), and the principal components were calculated by R FactoMineR package (Le, Josse, and Husson, 2008) and visualized by R factoextra package.

The pre-normalized microarray data and the RNA-seq data normalized to log2 (FPKM or RPKM +1) were tested by student t test and fold change. The count-based RNA-seq data was processed in R/Bioconductor DESeq2 package (Love, Huber, and Anders, 2014). All p-values were adjusted using the Benjamini and Hochberg approach. Genes with an adjusted p-value less than 0.05 and a fold change more than two were deemed as differentially expressed genes. Those genes could be ranked by the student t statistic to measure the differential expression levels.

The upregulated or downregulated genes in DCM/ICM samples were used for the identification of cell types significantly altered in HF. The gene set overrepresentation enrichment analysis (Fisher’s exact test) was employed to evaluate the significance of the differentially expressed genes (DEGs) against the cell type specific marker genes, which was implemented in R clusterProfiler (Yu, Wang, Han, and He, 2012) package.

The gene set enrichment analysis (GSEA) was used to calculate the enrichment degree of those upregulated genes involved in HF or cardiac fibrosis in endothelial cells. Specifically, all the genes were pre-ranked by the t statistic, which represented the differential expression levels. The GSEA analysis was implemented in R clusterProfiler (Yu et al., 2012), and the genes identified as core enrichment in this analysis were considered as key components.

The gene set overrepresentation enrichment analysis (ORA) was employed to identify the Reactome pathways enriched by previously detected endothelial cell specific marker genes and cardiac fibrosis-related genes in HF. This analysis was implemented in R ReactomePA package and visualized by R clusterProfiler (Yu et al., 2012) package.

The cell activity was estimated using single-sample Gene Set Variation Analysis (Hanzelmann, Castelo, and Guinney, 2013) (GSVA). Specifically, gene expression profiles and cell type specific marker genes were used as the input for GSVA to estimate the relative activities for each cell type and each sample.

The two-sample comparison was conducted by student t test, and the multiple-sample comparison was implemented by analysis of variance (ANOVA). The p-values for multiple-sample comparisons were adjusted to q-values by the Benjamini and Hochberg method. Any p-values or q-values less than 0.05 were considered as statistically significant.

To identify and characterize the cell types in the human left ventricle (LV), we collected two single-cell RNA sequencing datasets (scRNA-seq) of left ventricle provided by earlier study (Han et al., 2020). Subsequently, we eliminated the cells with low quality and retained 1,324 and 1,480 cells for further analysis (Materials and methods). As shown in Figure 1A, the cells from the two hearts were clustered into 18 clusters by the T-distributed Stochastic Neighbor Embedding (t-SNE) analysis, respectively. Using scHCL method, we successfully annotated 11 cell types for the two hearts (Figure 1A). Notably, the marker genes were specifically expressed in the cell types (Figure 1B). These results indicated that the cell types in the human left ventricle tissues could be identified and well-characterized by the scRNA-seq data.

With the cell types and marker genes in the left ventricles, we aimed to identify the cell types altered in the left ventricles of heart failure. We analyzed the gene expression profiles of 14 NF, 37 DCM, and 13 ICM samples from previous study (Sweet et al., 2018). The PCA and differential expression analysis revealed that the samples from the three groups exhibited significantly different expression patterns (Figures 2A,B). Furthermore, we also conducted GSEA on the marker genes of cell types to test whether those marker genes were clustered within the upregulated or downregulated genes of ICM or DCM. Specifically, the marker genes of fibroblast and endothelial cell were significantly enriched within the upregulated genes in both DCM and ICM (Figure 2C, adjusted p-value < 0.05), suggesting that the dysfunction of the two cell types might be associated with both DCM and ICM. Moreover, marker genes of dendritic cell, M1/2 macrophage, neutrophil, and smooth muscle cell were more specifically enriched within the upregulated genes in ICM (Figure 2C, adjusted p-value < 0.05). These results indicated that DCM and ICM had both similarity and specificity in the pathogenesis of heart failure based on these disease-related cell types.

As the endothelial cell and fibroblast could be activated in response to HF (Colombo et al., 2005), we then investigated the key regulators in the ECs and fibroblasts of HF, and collected scRNA-seq data of 1,082 endothelial cells from the left ventricles of NF, DCM, and ICM samples (Wang et al., 2020). The comparison of DCM and ICM samples with NF samples revealed that the endothelial cell specific marker genes were highly enriched in the upregulated genes of HF endothelial cells (Figure 3A, FDR <0.05). Specifically, a total of 24 EC marker genes were found to be upregulated in both HF tissues (bulk RNA-seq) and the endothelial cells of HF samples (scRNA-seq) (Figure 3B, p-value < 0.05). The pathway enrichment analysis identified inflammation-related cell adhesion molecules (CAMs) as key regulators, including CD74, HLA-B, HLA-E, HLA-DRB1, HLA-DQA1, HES1 and CLDN5, involved in the pathogenesis of HF (Figure 3C, FDR <0.05).

Furthermore, as transforming growth factor β1 (TGFβ1) is the principal pro-fibrotic factor in fibroblast activation (Akhurst & Hata, 2012), (Davis & Molkentin, 2014), which played vital roles in cardiac fibrosis (Ma, Iyer, Jung, Czubryt, & Lindsey, 2017), we examined whether the upregulated fibroblast marker genes in HF were involved in cardiac fibrosis. Consistently, we identified a large proportion of fibroblast marker genes upregulated in TGFβ1 induced cardiac fibroblast by differential expression analysis and GSEA (Figure 4A, FDR <0.05). Among these fibroblast marker genes, 29 were also upregulated in both HF tissues and fibroblast with TGFβ1 treatment (Figure 4B, FDR <0.05). The functional characterization of these genes revealed that LTBP2, LTBP1, COL3A1, MFAP4, COL12A1, COL1A1, COL1A2, MMP2, TIMP2, and PCOLCE2 were primarily involved in extracellular matrix (ECM) organization and collagen biogenesis/formation/degradation (Figure 4C, FDR <0.05). Collectively, these results indicated that inflammation-related CAMs and ECM proteins such as collagens were specifically secreted by endothelial cell and fibroblast, respectively, and might induce cardiac inflammation and fibrosis during heart failure.

As ICM had more specific immune cell types, such as macrophage and dendritic cell (DC), than DCM, we then estimated the activities of immune cells including macrophage, DC, and neutrophil. Neutrophil and macrophage appeared to have higher activities in ICM than DCM and NF (Figure 5A, p-value < 0.05). Consistently, the marker genes of neutrophil and macrophage were also observed to be specifically upregulated in ICM (Figure 5B, p-value < 0.05). The cell-cell communication analysis revealed that the autocrine ligand-receptor interaction induced chemokine signaling activation in neutrophil and macrophage might be responsible for the immune response in ICM (Figure 5C). Particularly, the ligands, CCL3, and CCL4, and the receptor CCR5 were specifically upregulated in ICM as compared with DCM and normal controls (Figure 5D). These results indicated that higher inflammation in ICM might be associated with autocrine CCL3/CCL4–CCR5 interaction induced chemokine signaling activation.

We collected an independent gene expression dataset from previous study (Hannenhalli et al., 2006) for validation. The inflammation-related CAMs such as HLA-E, HLA−DQA1, HLA−DRB1, and CD74, and all the ECM genes were upregulated in the HF samples of bulk RNA-seq dataset (GSE121893, Figure 6A, p-value < 0.05). Notably, the ECM genes were also upregulated in the fibroblasts of HF from an independent scRNA-seq dataset (Figure 6B). Furthermore, neutrophil and macrophage activities also appeared to be higher in ICM compared with NF and DCM, and the upregulation of autocrine ligand-receptor pairs in ICM, CCL3/CCL4 –CCR5, was also observed in the validation dataset (Figures 6C,D, p-value < 0.05). Consistently, the CCL3 and CCL4 were expressed higher in the macrophages of ICM than the DCM and normal hearts (Figure 6E). These results further indicated that inflammation-related CAMs and ECM proteins, which were specifically secreted by endothelial cell and fibroblast, respectively, and chemokine signaling activation in neutrophil and macrophage might induce cardiac inflammation and fibrosis during heart failure.