The workflow of data processing and analysis is depicted in Figure 1. We first retrieved the bulk RNA-seq data of the heart tissue of a cardiac disease mouse model from NCBI GEO database (http://ncbi.nlm.nih.gov/geo/). By using the keyword combination [cardiac disease (All Fields) OR heart disease (All Fields)] AND [“Mus musculus” (Organism) AND “Expression profiling by high throughput sequencing” (Filter)], we have collected 17 bulk RNA-seq datasets that contain 184 bulk RNA-seq samples in total. It is known that the quantification, gene coverage, and file format of the author-provided processed gene expression matrixes varied from dataset to dataset. Therefore, to reduce laboratory bias and systematic errors, we calculated and corrected the expression values of each sample by starting from its raw sequencing data. First, the raw sequencing data were downloaded from the SRA database and transformed into FASTQ format by SRA-Toolkit (https://trace.ncbi.nlm.nih.gov/ Traces/sra/sra.cgi?view = software). Adapters and low-quality reads were automatically removed by fastp software (Chen et al., 2018) to control the quality of reads. After quality control, the clean reads were aligned to mouse genome (GRCm38) using STAR (Dobin et al., 2013) (version 2.7.3a with default parameter and GENCODE version M24 GTF gene annotation), and the gene expression quantification of each sample was performed by using the RSEM software (version 1.3.1, with default parameter) (Li and Dewey, 2011). The gene expression read counts from each sample were merged into the gene expression matrix by the Perl script modified from RSEM software.

An unexpected, prominent variation between samples from the same disease condition group but different datasets was observed (Supplementary Figure S1), derived by principal component analysis (PCA), where the samples were more likely to be clustered according to the datasets rather than to the disease conditions. Therefore, the batch effect emerged as a serious problem before any further analysis based on the gene expression matrix. To correct these batch differences, the whole gene expression matrix which contained all samples was first processed by quantile normalization, and then the normalized expression matrix was submitted to batch correction through the popular ComBat tool (Johnson et al., 2007). To minimize the batch effects and retain biological variations between samples, the batch covariate and type covariate were both considered when running the ComBat tool. Nevertheless, there were still some outlier datasets that did not cluster according to disease condition after correction. After removing these outlier samples, the corrected bulk RNA-seq expression matrix covering 110 samples from 12 datasets (Supplementary Table S1) was obtained, which included disease model samples in MI, DCM, HCM, and ARVC and their corresponding healthy samples.

As for the reference scRNA-seq data, we first downloaded the single-cell expression data of mouse hearts from Tabula Muris Consortium (Tabula Muris, 2018). To accurately annotate the cell types, we used the Seurat R package (v3.2.0) (Satija et al., 2015) to perform unsupervised clustering of the single-cell read count profiles. The read counts from each cell were divided by the total read count of that cell and then were multiplied by a scaling factor of 10,000 and log-transformed. To improve the clustering stability, genes detected in less than 10 cells were excluded from the downstream analysis. We performed PCA on the normalized expression matrix using highly variable genes identified by the FindVariableGenes function of Seurat. Following the result of the PCA, the top significant principal components were selected for clustering, with a resolution parameter equal to 0.5. The cell clusters were visualized using t-distributed stochastic neighbor embedding method, and their corresponding markers were listed using the FindAllMarkers function of Seurat. Finally, the cell type or subtype of each cluster was annotated based on the known cell type markers from the CellMarker database (Zhang X. et al., 2019). For healthy heart samples from Tabula Muris Consortium, six cell types were annotated by this procedure, and four of them were used to further build a signature matrix, including cardiomyocyte, endothelial cell, fibroblast, and smooth muscle cell.

To better understand the disease-related cell abundance alteration, cell-type-specific pathway alteration, and/or cell subpopulation alteration, we also introduced disease-related scRNA-seq datasets from the GEO database as the alternative scRNA-seq references. These alternative disease-related GEO references include the following:

(1) MI-related scRNA-seq dataset (GSE120064) of Ren et al.: This reference datasets covered cardiomyocytes, endothelial cells, fibroblasts, smooth muscle cells, and macrophages in the MI mouse model’s heart tissues (Ren et al., 2020).

CIBERSORTx is a new machine-learning-based method for estimating the relative abundance of cell types and cell-type-specific expression from the bulk RNA-seq data. To achieve this, the initial step is to construct a signature matrix featuring the gene expression characteristics of each cell type. We here considered two sources of cell type annotations. First, the aforementioned scRNA-seq reference data processed from Tabula Muris reference or the alternative disease-related GEO references were used for the four major cell types of the heart, which was prepared and formatted as a standard signature matrix by the Create Signature Matrix analysis module of CIBERSORTx with suggested parameters. Second, to further investigate the immune microenvironment of heart tissues of different diseased states, the built-in LM22 immune cell signature matrix was used to estimate the immune cell infiltration and expression. To observe consistency with other mouse expression matrixes, we converted the human gene symbols of LM22 into the corresponding mouse homologous gene symbols before calculation. Then, the two prepared signature matrix files along with the batch-corrected bulk RNA-seq gene expression matrix were submitted to the Impute Cell Fraction and Impute Cell Expression analysis modules of CIBERSORTx to estimate the relative abundance and expression for each of the 26 cell types in the signature matrixes, respectively. When imputing cell fraction, the S-mode was applied to remove batch effect between signature matrix and mixture samples. When imputing gene expression, high resolution mode was used to further enable the estimation of cell-type-specific gene expression profiles for each of the 110 bulk RNA-seq samples. A gene subset file which covers a list of 18,584 genes from the bulk expression file was used as input for analysis of as many genes covered by the bulk transcriptomes as possible.

The reconstructed cell-type-specific gene expression matrixes allow the screening and analysis of genes differentially expressed for each cell type between disease (MI, DCM, HCM, and ARVC) and healthy (or control) samples. Taking the cell-type-specific gene expression matrix as the input, cell-type-specific DEGs were calculated by DESeq2 (Love et al., 2014), and those with adjusted p-value <0.05 were considered as the significant DEGs. The g:Profiler (Raudvere et al., 2019) is a web server for simultaneously performing multiple functional enrichment analyses with timely updated reference gene function annotations. To explore the potential biological functions of these cell-type-specific DEGs, the DEG lists were submitted to g:Profiler to investigate their associated significantly overrepresented (adjusted p-value <0.05) Gene Ontology functional terms and WikiPathways biological pathways. As for the disease gene enrichment analysis, the disease gene annotations were firstly retrieved from DisGeNET database (Piñero et al., 2017), and only “heart disease” and its child diseases terms (e.g., myocardial infarction) were retained. We adopted Fisher exact test to test the cell-type-specific heart disease gene associations. More specifically, genes related to each heart disease were treated as one gene set, and the enrichment of these disease genes in each cell-type-specific DEG set was compared to the background (i.e., all genes in the transcriptome) by Fisher exact test. For each set of cell-type-specific DEGs from MI, HCM, DCM, and ARVC, we separately calculated their significance and retained the disease gene terms satisfying the significance threshold of false discovery rate-adjusted p-value <0.05.

Protein–protein interaction network analysis was performed using GeneMANIA (Warde-Farley et al., 2010). The input genes contained the differentially expressed secretory protein genes in monocytes and the overrepresented pathway-related genes in fibroblasts. The prioritized interacting gene pairs related to the input gene list were used to build a regulatory subnetwork, which was further visualized by Cytoscape software (version 3.8.0). CellChat (Jin et al., 2021), an R-based computational analysis tool, was implemented for the analysis of cell–cell communication. The expression profiles of 26 cell types were used as the input of CellChat to construct the cell–cell communication network. Briefly, we first calculated the significant ligand–receptor pairs (LRs) between any two cell types. CellChat covered several scenarios of LRs, and only the typical Secreted Signaling LR sets in CellChatDB were considered. Afterwards, the LR count matrix between any two cells was constructed, and the healthy LR count matrix subtracted from the diseased LR count matrix was used to reflect changes in cell–cell communication between healthy and disease samples. The changing of cell–cell communication network was also visualized by Cytoscape.

Ten-week-old C57BL6/J mice were used in our experiments. The mice were first anesthetized by isoflurane inhalation and intubated with an endotracheal cannula. As for the myocardial infarction group, intercostal thoracotomy was performed to expose the heart, followed by ligation of the left anterior descending coronary artery. As for the sham-operated mice, they underwent the same surgical incision without ligation. All mice were then assessed by ultrasound imaging to ensure the success of disease modeling. All animal experiment protocols complied with the Animal Management Rules of the Ministry of Health of the People’s Republic of China and the Guide for the Care and Use of the Laboratory Animals of Peking University and were approved by the Laboratory Animal Ethics Office of Peking University Biomedical Ethics Committee (LA2020337).

At 1 or 4 weeks after ligation or sham surgery, hearts were harvested from the euthanized mice and perfused with 20 ml of cold phosphate-buffered saline (PBS). The hearts were then minced and digested with shaking for 30 min at 37°C in RPMI 1640 medium containing collagenase I (450 U/ml, Sigma, St. Louis, USA), collagenase XI (125 U/ml, Sigma, St. Louis, USA), DNase I (60 U/ml, Roche, Basel, Switzerland), and hyaluronidase (60 U/ml, Harveybio, Beijing, China) enzymes. The digested tissues were passed through a 40-μm filter, centrifuged at 1,600 rpm for 6 min, and resuspended in PBS containing 2% fetal bovine serum, followed with antibody staining. The antibodies used in the experiment are presented in Supplementary Table S2.

Single-cell suspensions were labeled with antibodies for 30 min at 4°C and washed in FACS buffer. Flow cytometry analysis was performed on LSRFortessa (BD Biosciences), and the results were analyzed with FlowJo software. Analysis of significant differences on cell counts in flow cytometry between healthy and MI samples was performed using Student’s t-test. The differences on calculative composition were analyzed using Mann–Whitney test.

Multiple bulk RNA-Seq data could be retrieved by searching the GEO database (Barrett et al., 2013) with cardiac disease-related keywords. However, the author-provided gene expression matrixes which are readily downloadable on GEO are not suitable for our integrative analysis purpose because of several pronounced problems, like inappropriate sample source (e.g., blood samples), discrepancy of expression quantification methods, and prominent batch effects. Therefore, we instead started from the raw sequencing data and applied an extensive and carefully standardized pipeline for data collection and processing (see “Materials and Methods”). As a result, we assembled a batch-corrected bulk RNA-seq expression matrix covering 110 samples from 12 independent datasets. Four types of cardiac disease samples, including MI, DCM, HCM, and ARVC, and their corresponding healthy controls were included in this matrix. We noted that a reasonable PCA clustering result could be observed based on this corrected expression matrix, where samples that belong to the same disease group were clustered together well even with the variations in source datasets and disease modeling approaches (Figure 2A). This PCA result indicates a successful reduction of technical bias and retention of biological signatures in the bulk expression matrix produced from our standardized data processing pipeline.

To reconstruct cell-type-specific gene expression profiles, the gene expression characteristics should be first extracted from a comprehensive scRNA-seq reference data. To this end, the scRNA-seq reference data of mouse heart was obtained from Tabula Muris Consortium (Tabula Muris, 2018). By cell clustering and marker gene annotation, we identified four major cell types (Supplementary Figure S2) of heart from this data, including cardiomyocyte (with Actc1 and Actn2 markers), fibroblast (Col1a2 and Col3a1markers), endothelial cell (Egfl7 and Emcn markers), and smooth muscle cell (Rgs5 and Acta2 markers). Previous studies have also confirmed these four types as the main cell types in heart tissue (Litviňuková et al., 2020; Wang et al., 2020), while the macrophages and erythrocytes were not consistently detected in other heart scRNA-seq datasets and arouse concern of blood contamination and thus were not considered in the subsequent analyses. Therefore, the scRNA-seq expression matrix for the four major cell types was used as the representative reference for further analysis. Based on this reference, the input bulk RNA-seq expression matrix can be dissected in a cell-type-specific manner via the recently established CIBERSORTx computational approach (Steen et al., 2020). More specifically, for each bulk RNA-seq sample, the relative fraction and the gene expression profile of each cell type was calculated by CIBERSORTx. To preliminarily assess the reliability of the reconstructed cell-type-specific expression matrix, we performed PCA based on the expression matrix of cardiomyocytes (Figure 2B). The clustering result revealed that the cardiomyocyte expression profiles were clearly clustered according to the disease condition, demonstrating the recapitulation of bona fide biological characteristics of cardiomyocytes in cardiac diseases by the reconstructed expression matrix.

We first investigate the alteration of cell composition in cardiac disease condition based on the estimated cell abundance of CIBERSORTx. The average abundances of the four major cell types in different sample groups are shown in Figure 3A. According to the results, cardiomyocytes which constitute over 80% of these four main cell types exhibit the highest abundance in all sample groups. Following cardiomyocytes, the next abundant cell types are endothelial cells and fibroblasts in healthy samples. Interestingly, a recent study (Pinto et al., 2016) based on new genetic tracers and enhanced flow cytometry techniques has confirmed endothelial cells as the most abundant cell population other than myocytes, while other studies considered that fibroblasts are more abundant than endothelial cells (Banerjee et al., 2007; Bergmann et al., 2015). On the one hand, our results also support the considerable proportions of endothelial cells and fibroblasts in heart tissues, where their abundances are largely comparable in healthy heart tissues. On the other hand and more importantly, by comparing the estimated abundances between the disease and healthy sample groups, a substantial alteration of the relative abundances of these two cell types can be observed.

As noted above, the comparison also indicates alteration of cell composition in pathological conditions. In MI samples, the abundances of cardiomyocytes are significantly decreased compared to the healthy samples (p < 0.01, Mann–Whitney test). One straightforward explanation of this observation would be the hypoxia-induced cell death of a vast number of cardiomyocytes in proximity to the damaged area in MI (Frangogiannis, 2015; Martín-Fernández and Gredilla, 2016). On the contrary, a significant elevation of fibroblast abundance can be observed not only in MI samples (p < 0.01, Mann–Whitney test) but also in DCM samples (p < 0.01, Mann–Whitney test), indicating the cardiac fibrosis progression in these disease conditions. To verify these results, we adopted two other scRNA-seq-based transcriptome reconstruction methods, namely, MuSiC (Wang et al., 2019) and Bisque (Jew et al., 2020), to estimate the cell type abundance using the same single-cell reference. It is noteworthy that, unlike CIBERSORTx, MuSiC and Bisque can only calculate the cell type relative abundance (i.e., cell proportions that sum up to 100%). The estimated cell type abundance scores from MuSiC and Bisque are summarized in Supplementary Figure S3. Despite the distinguishing cellular abundance scores between the different methods, the variations of cell type abundance scores across different disease conditions are highly similar to those from CIBERSORTx—for example, the relative score of cardiomyocytes is decreasing in MI and of fibroblast is increasing in MI and DCM compared with the healthy samples (p < 0.01, Mann–Whitney test). Such trend is the same as that observed in the CIBERSORTx results.

To further verify this observation, we applied flow cytometric assay to compare the fibroblast abundances in healthy and MI mouse hearts. Indeed a significantly increased fibroblast fraction is observed in MI condition in comparison with the sham condition, which also confirmed our result (Figure 3B). Cardiac fibrosis is a prevalent pathophysiological process in many myocardial diseases and considered to be the end state of heart injury along with increasing fibroblasts (Kong et al., 2014). Our analysis implies that cell proliferation may not be the sole source of the increased fibroblasts—for example, in DCM samples, the average abundance of endothelial cells is decreased (p < 0.01, Mann–Whitney test) along with the increased fibroblasts, indicating that endothelial-to-mesenchymal transition (EndMT) would be a common event during heart fibrosis in DCM. EndMT is an intricate cellular differentiation process in which endothelial cells lose their properties and acquire mesenchymal features, and this process would give rise to the transition from endothelial cells to fibroblasts (Cheng et al., 2021). Indeed Xie et al. have demonstrated that EndMT might contribute to myofibroblast recruitment in human DCM, characterized by decreased endothelial markers and increased mesenchymal markers in the immunofluorescence co-localization analysis (Xie et al., 2018). They also found that several cells expressed both mesenchymal and endothelial markers, which occurred exclusively in the cardiac sections of DCM patients but not in normal cardiac samples. Moreover, Zhang et al. has observed a significantly elevated EndMT in DCM rats, and this process can be suppressed through inhibition of TGF-β/ERK signaling (Zhang C. et al., 2019). In a mouse model of cardiac fibrosis, transforming growth factor-beta 1 (TGF-beta1) can induce the endothelial cells to undergo EndMT, whereas bone morphogenetic protein 7 can preserve the endothelial phenotype (Hong et al., 2018; Piera-Velazquez and Jimenez, 2019).

Varied gene expression patterns of fibroblasts can be observed between different disease conditions (Figure 3C), suggesting the plausible functional heterogeneity of fibroblasts in different cardiac diseases. To further explore the related biological functions in each individual cell type, the cell-type-specific DEGs were determined following the method described above (“Materials and Methods”). Among the four major cell types, fibroblasts indeed show the largest amount of DEGs, followed by cardiomyocytes and endothelial cells, while DEGs can hardly be detected in smooth muscle cells (Supplementary Figure S4 and Supplementary Table S3). Based on these identified DEGs, WikiPathways and Gene Ontology (GO) enrichment analyses were performed for each cell type in each specific disease condition. The results of the WikiPathways enrichment analysis are summarized in Figure 4 and Supplementary Table S4. Few significant pathways are obtained from ARVC DEGs for any cell types and smooth muscle cell DEGs for any disease conditions (which is also resulting from the limited DEGs identified in these cases, as shown in Supplementary Table S3), so the functional enrichment analysis results for these DEGs are not presented. In summary, the pathways enriched in cardiomyocytes are largely shared across disease conditions (Figure 4), where the pathways related to mitochondrial energy metabolism are most prominent, such as electron transport chain, oxidative phosphorylation, fatty acid beta oxidation, and tricarboxylic acid cycle. Moreover, the GO enrichment term analysis also provides signatures associated with mitochondrial dysfunction, like mitochondrion organization and electron transfer activity (Supplementary Figure S5). The mammalian heart must contract incessantly to sustain life, which is heavily dependent on adenosine triphosphate (ATP) production and consumption (Doenst et al., 2013). The cardiomyocytes had consumed more than 90% ATP generated by the mitochondria, and thus the cardiomyocyte function requires proper mitochondrial energy metabolism. However, our analysis suggests mitochondrial dysfunction in cardiomyocytes accompanied with the development of cardiac diseases like MI, HCM, and DCM, which could be also confirmed by previous experimental researches (Peoples et al., 2019). In addition to cardiomyocytes, endothelial cells also have many enriched pathways related to energy production, suggesting its similarity with cardiomyocytes in cellular response to disease state.

Unlike cardiomyocytes and endothelial cells, the associations with energy metabolism are much less prevalent for fibroblasts. Instead we find that Wnt signaling pathway is specifically enriched in fibroblasts under DCM condition. Wnt signaling is a complex collection of signal transduction pathways, and the aberrant regulation of the Wnt signaling pathway is associated with many diseases, such as cancer, degenerative diseases, and also cardiovascular diseases (Nusse and Clevers, 2017; Foulquier et al., 2018; Lietman et al., 2018). The significant enrichment of the Wnt pathway in fibroblasts suggests its important role in the occurrence and development of DCM, and indeed one recent study has described that canonical WNT/β-catenin signaling activity is impaired in DCM, and reactivation of its activity can improve cardiac contractility and ameliorate intraventricular conduction defects (Lu et al., 2016). The MAPK signaling pathway is another specifically enriched pathway in fibroblasts in MI, which is also recently proven to modulate the apoptosis of myocardial cells in acute MI heart (Zhang et al., 2020). The GO term enrichment analysis in fibroblasts is also related with MAPK cascade, actin cytoskeleton organization, and cell adhesion molecular binding (Supplementary Figure S5), which also supports the pathway enrichment result. Notably, several pathways mapped to fibroblasts have not been reported, including the EGFR1 signaling pathway, the insulin signaling pathway, and the alpha6-beta4 integrin signaling pathway. The involvement and mechanisms of these pathways in cardiac fibrosis and cardiac pathogenic process would be an interesting topic for further investigation according to our cell-type-specific pathway analysis.

Furthermore, to establish a relationship between cell types and the known disease genes of cardiac diseases, we performed a disease gene enrichment analysis for each cell-type-specific DEG set (see “Materials and Methods”). The results (Supplementary Figure S6) demonstrate the diversity of cell-type-specific associations with disease terms between different disease conditions—for example, in MI, the fibroblasts and endothelial cells are associated with the myocardial infarction term. As for HCM and DCM, the HCM-associated cell types include cardiomyocyte and endothelial cells, while the DCM-associated cell types include cardiomyocytes, fibroblasts, and endothelial cells. This observation indicates their essential role in the occurrence of corresponding diseases and would be helpful for further cell-type-specific disease mechanism investigations.

The above-mentioned analyses were based on the single-cell references derived from healthy samples. However, scRNA-seq data derived from the corresponding heart diseases could provide additional valuable information. To this end, we introduced two other public scRNA-seq datasets derived from MI (Ren et al., 2020) and HCM (Zhang Y. et al., 2019) samples as the alternative, disease-related scRNA-seq references. By reconstruction of bulk transcriptome profiles using these two reference datasets, we have re-estimated the cell type’s abundances and cell-type-specific pathways across different disease conditions (Figures 5A,B). Generally, the overall cell abundances and the overrepresented pathways in DEGs are similar between the results based on healthy single-cell reference and those based on the alternative disease-related single-cell references—for example, the abundance of cardiomyocytes decreases in MI and that of fibroblasts increases in MI and DCM along with a decrease of endothelial cells (p < 0.01, Mann–Whitney test). Furthermore, the disease-reference-based results from a cell-type-specific pathway analysis also recapitulate energy metabolism-related pathways of cardiomyocytes and several overrepresented pathways (EGFR1 signaling pathway, MAPK signaling pathway, focal adhesion, alpha6-beta4 integrin signaling pathway, and regulation of actin cytoskeleton) in fibroblasts in MI (Figures 5C,D). On the other hand, we also found several pathways which are not represented in our previous healthy reference-based analysis, like IL-6 signaling pathway and oxidative damage in fibroblasts. It is well accepted that IL-6 is mainly secreted by fibroblasts in the heart and can promote CF proliferation by promoting EndMT in ECs (Zhang et al., 2021). This observation indicates that combining those different datasets is helpful for a more comprehensive knowledge associated with cardiac diseases.

As shown above (Figure 3A and Figures 4A,B), fibroblasts were suggested to be an important participant of the transcriptome alterations of heart tissue in multiple disease conditions like MI and DCM. It is interesting to investigate if particular subpopulation(s) of fibroblasts contribute more to such alterations. To further identify which subpopulation of fibroblasts is activated in the disease conditions, we further introduced a fibroblast subpopulation-annotated scRNA-seq data (Ruiz-Villalba et al., 2020) as the alternative reference. Six major fibroblast subpopulations (i.e., FB_sp1 to FB_sp6) were annotated in this MI-related, fibroblast subpopulation-annotated scRNA-seq reference. As shown in Figure 5E, FB_sp2 that is featured in the high expression of Cthrc1 is observed to be increased in MI, DCM, and HCM (p < 0.01, Mann–Whitney test). Remarkably, FB_sp2 has been experimentally verified to be increased in the infarcted myocardium (Ruiz-Villalba et al., 2020). FB_sp3 (Efhd1 and Cd248 highly expressed) is another fibroblast subpopulation showing an elevated abundance in MI, DCM, and HCM (p < 0.01, Mann–Whitney test). Other subpopulations, like FB_sp1 (Hsd11b1 highly expressed), is the main subtype in ARVC where the amount of FB_sp2 is very few. Taken together, our results suggest that FB_sp2 may serve a crucial role in other cardiac diseases besides MI, and other subpopulations like FB_sp1 and FB_sp3 may also be associated with some specific cardiac diseases.

Immune cells play crucial roles in the occurrence and development of cardiovascular diseases, and many studies have demonstrated that inflammatory response is prevalent under cardiac pathophysiological conditions (Bartekova et al., 2018). Though some inflammation processes are required for tissue repair during injury and are therefore protective (Liao et al., 2018), prolonged inflammation leads to myocardial remodeling and apoptosis of cardiomyocytes (Frangogiannis, 2014). To further compare the differences in immune cell infiltration among different disease types, we used the built-in LM22 immune cell reference of CIBERSORTx to calculate the abundance and cell-type-specific expression profile of various immune cell subtypes in mouse hearts. We find that alteration of immune cell subtype abundance is a common phenomenon across different cardiac diseases (Figure 6A)—for example, the number of macrophages increases in DCM along with reduction of monocytes (p < 0.01, Mann–Whitney test). It is widely accepted that monocytes are the precursor cells of macrophages and are recruited to the inflammatory site by inflammatory factors (Peet et al., 2020), and such increased number of macrophages and decreased number of monocytes would suggest that the monocytes are more likely to differentiate into macrophages in DCM.

When using the LM22 built-in immune cell reference of CIBERSORTx, the result suggests that the abundance of monocytes increases, but the abundance of macrophages remains almost unchanged. However, several previous studies have reported upregulated macrophages in MI (King et al., 2017; Farbehi et al., 2019). In order to explore the reason of this discrepancy, we used three other scRNA-seq data (GSE129175 from HCM, GSE120064 from MI, and GSE109744 from healthy samples) as the alternative references to estimate macrophage abundance in MI (Figure 6B). We found that two out of the three results indicate an elevated macrophage abundance in MI. Furthermore, our result of the flow cytometry assay confirms the significantly elevated macrophage abundance in mouse hearts in the MI group compared with the sham group (Figures 6C,D). This result indicates that the external scRNA-seq data is a better choice for analyzing macrophage transcriptomes than the built-in LM22 references. However, most other immune cell (sub) types are not covered by external scRNA-seq data, and an analysis based on the LM22 reference is able to provide useful information, as exemplified below. A previous study has speculated that monocytes can serve as critical mediators of the inflammatory response and also mediate cardiac repair in the development of MI (Heidt et al., 2014). Therefore, we assume that not just differentiated macrophages but also monocytes can serve as primary myeloid cells that contribute to MI. Although the GEO scRNA-seq datasets did not cover reference single-cell expression profiles for monocytes, the LM22-based analysis indicates the increased monocyte fraction in MI (Figure 6A). A further flow cytometry assay also confirmed the significantly elevated monocyte abundance in mouse hearts in the MI group compared with the sham group using the gating strategy that could distinguish monocytes from macrophages (Figures 6C,E). Since monocytes usually play a biological role through secretory proteins (Liu et al., 2018; McKiernan et al., 2020), we further explored the expression pattern of secretory proteins in monocytes (Figure 6F). Indeed the heat map shows that the expression characteristics of the secreted proteins in MI are different from those of the healthy and other diseased states.

Certainly the above-mentioned analysis is not merely limited to monocytes but covers 22 representative immune cell subtypes. Results about other immune cell subtypes would also provide helpful clues to investigate the immune involvement in cardiac diseases. An example is ARVC, where the abundance of T cells and eosinophils increases while the abundance of macrophages decreases (p < 0.01, Mann–Whitney test). A clinical investigation has shown that inflammatory infiltration in the ventricular myocardium of ARVC samples is associated with severe structural heart changes, and T-lymphocytes are the main infiltrated immune cell type (Campuzano et al., 2012), indicating that T cells may exacerbate the progress of ARVC.

As monocytes are elevated in myocardial infarction along with proliferation of fibroblasts, we speculate that monocytes may regulate cardiac fibrosis through secretory-protein-mediated cell–cell communication. To explore this possibility, we first analyzed the possible protein–protein interaction (PPI) between the differentially expressed secreted proteins from monocytes and the overrepresented pathway-related proteins. The GeneMANIA PPI network analysis confirms the wide interactions between these two groups of proteins (Supplementary Figure S7), supporting the hypothesis that several proteins can be secreted by monocytes to regulate the function pathways in fibroblasts during the pathological cardiac fibrosis process.

We further extended the cell-to-cell communication analysis to all of the 26 cell types (4 major cell types + 22 immune cell subtypes) covered in this study. We used CellChat to calculate changes in cell-to-cell communication patterns for each disease type (Figure 7). In line with the cell abundance and functional pathway analyses, for most cases, fibroblasts are at the center of the cell communication network, especially in MI and ARVC, suggesting its central role in the development of cardiac diseases. In addition, immune cells, such as activated NK cells, show strong connections with many other cell types in HCM, DCM, and MI. A recent study found that NK cells play a significant role in repairing injured tissue and maintaining tissue homeostasis (Ong et al., 2017). NK cells may function in a cardiac immune environment directly through receptor–ligand interactions or indirectly through cytokine secretion. One recent study has demonstrated that NK cells can not only protect against the development of cardiac fibrosis by limiting collagen formation in cardiac fibroblasts and by preventing the accumulation of specific inflammatory populations in the heart but also prevent monocrotaline-induced endothelial damage (Notas et al., 2009). In our result, activated NK cells enhance the interaction with macrophage in DCM and ARVC and show strong interactions with dendritic cells in MI and HCM. Indeed NK cells have been reported to accelerate the maturation of macrophages and dendritic cells, indicating its role in altering the cardiac immune microenvironment, which again suggests that the cell–cell communication analysis would be a reasonable proxy for future research of cardiac disease-associated immune microenvironment shifts.