ECM-producing fibroblasts have similar phenotypes across mouse and human single-cell transcriptome atlas

To comprehensively characterize the subpopulations of fibroblasts in MI, we integrated two previously published scRNA-seq datasets [16,19]. After quality control and filtering, we ended up with a dataset containing 12,0121 cells categorized into 11 different cell types (Fig. 1A). These cells were extracted from mouse hearts 1, 3, 5, 7, 14 and 28 days after MI (Fig. 1B). We extracted fibroblasts from the integrated scRNA-seq dataset and re-clustered them to get higher resolution for heterogeneity and ultimately identified 12 subclusters (Fig. 1C). Analyzing the composition of fibroblasts at different time points we identified three subsets of fibroblasts (cluster 5, cluster 9 and cluster 10) that were absent in controls but appeared during the acute phase of the MI period (Fig. 1D). Next, we analyzed all the fibroblast subsets to identify genes characteristic of different subsets (Table S5). Characterization genes for cluster 5 included Cthrc1, Acta2, Fn1, Col1a1 and Postn (Figure S1A; Table S5). Characterization genes for cluster 9 included Mt2, Ccl2, Timp1, Mt1 and Prg4. Cluster 9 was present in large numbers on day 1 after MI, decreases rapidly on day 5, and exhibited high expression of genes related to immunity. Characterization genes for cluster 10 included Stmn1, H2afz, Cks2, Cenpa and Acta2 (Figure S1A; Table S5). We also found that MKi67 was almost exclusively expressed in cluster 10 (Figure S1B), suggesting that cluster 10 represents a subset of fibroblasts with high proliferative capacity. Of note, we identified cluster 7, which is abundant in the late phase of MI but rare in the control and acute phases of MI. Characterization genes for cluster 7 included Comp, Sfrp2, Cilp, Eln, Wisp2 and Ctgf (Figure S1A; Table S5). This subcluster exhibits the characteristic of matrifibrocytes, which appeared late in the infarction phase as mentioned in the previous study, persisted in the scarred area, expressed the Comp gene, and lost αSMA expression (Figure S1A) [39,40]. The biological processes of cluster 5 were enriched in the collagen fibril organization and ECM organization (Fig. 1F), suggesting that this cluster was closely associated with fibrotic scar formation and may represent the reparative cardiac fibroblasts (RCFs) mentioned in a previous study [23]. To test our hypothesis, we reanalyzed the Ruiz-Villalba et al. generated scRNA-seq, and clustering analysis found Cthrc1 is highly expressed in cluster 3, suggesting that cluster 3 represents the RCFs mentioned by Ruiz-Villalba et al. (Figure S1C and S1D). We then extracted the top 50 DEGs in this subpopulation and mapped them to our fibroblast dataset, and finally found this signature was highly expressed in cluster 5, in concordance with our hypothesis (Figure S1E). Of note, we found that the top 50 DEGs of matrifibrocytes, which appeared in the late phase of MI, were also predominantly enriched with ECM organization (Fig. 1G). DEGs analysis further confirmed that Postn, Cthrc1 and Acta2 were expressed in RCFs, while Comp and Gsn were enriched in matrifibrocytes (Fig. 1H). We next used pySCENIC [38] to explore whether the gene regulatory structure (regulons) in RCFs was deranged compared with the matrifibrocytes. The results showed significant differences in the TFs of these two types of fibroblasts, with RCFs exhibiting much higher expression of transcription factors such as Bcl11a and Sirt6, as well as high expression of SOX4, SOX6 and SOX8 (Fig. 1I). Among them, SOX6 is considered an important target for cardiovascular disease [41]. Matrifibrocytes, in contrast, had high expression of E2f6 and Atf3 (Fig. 1I). GO enrichment analysis of Mki67+ fibroblasts showed that the DEGs were mainly enriched in mitotic spindle organization, microtubule cytoskeleton organization involved in mitosis and mitotic cell cycle phase transition (Figure S1F). DEGs analysis showed that Col1a1, Col3a1 and Sparc were expressed in RCFs, while Stmn1 and Cenpa were enriched in Mki67+ fibroblasts (Figure S1G).

Fig. 1

Activated fibroblasts have similar phenotypes across mouse and human single-cell data. A Overview of the mouse scRNA-seq datasets. The project distribution, time distribution, clusters distribution and cell types distribution are depicted as UMAP plots. B UMAP plots of cells from 0, 1, 3, 5, 7, 14, and 28 days after MI in the integrated mouse scRNA-seq dataset colored by cell type. C UMAP scRNA-seq plot of fibroblasts from the integrated mouse scRNA-seq dataset. A total of 12 clusters were identified. D UMAP plots of fibroblasts from 0, 1, 3, 7, 14, and 28 days after MI colored by cluster. E Gene expression of Postn, Fn1, Col1a1 and Runx1. F GO enrichment analysis of the top 50 genes in cluster 5 (RCFs, mouse). G GO enrichment analysis of the top 50 genes in cluster 7 (matrifibrocyte, mouse). H The volcano plots showing the DEGs between cluster 5 vs. cluster 7. I Heatmap showing TFs activity for cluster 5, cluster 7 and cluster 10 fibroblasts subsets (mouse). J UMAP plot of cells from the integrated human scRNA-seq dataset colored by cell type. K UMAP scRNA-seq plot of fibroblasts from the integrated human scRNA-seq dataset. A total of 12 clusters were identified. L UMAP plots showing mouse RCFs scRNA-seq signatures using the top 10 DEGs from mouse RCF mapped into fibroblasts from the integrated human scRNA-seq dataset. M UMAP plots showing mouse cluster 7 scRNA-seq signatures using the top 10 DEGs from mouse cluster 7 mapped into fibroblasts from the integrated human scRNA-seq dataset. N Gene expression of COMP, SFRP2, CILP and WISP2. O Gene expression of POSTN, FN1, COL1A1 and RUNX1

To compare fibroblast subsets in human heart tissues, we downloaded snRNA-seq datasets of MI and ICM from the Human Cell Atlas and the GEO database and performed quality control and integration and a dataset of 72,515 cells was finally obtained. Clusters were annotated with marker genes from the literature [12,15] and nine major cardiac cell types were identified (Fig. 1J). we extracted fibroblasts from the human scRNA-seq dataset and reclustered them for analysis. Unbiased clustering eventually grouped cardiac fibroblasts into 12 different subsets (Fig. 1K). We first mapped the gene expression signature of mouse RCFs using an AUC score. This signature was most highly expressed in human cluster 3 and cluster 6 (Fig. 1L). At the same time, the marker gene (POSTN, FN1, COL1A1 and RUNX1), which was highly expressed in mouse RCFs, was also highly expressed in human fibroblasts cluster 3 (Fig. 1O). This result suggests that human fibroblasts cluster 3 and mouse RCFs have similar phenotypes. Next, we mapped the gene expression signature of mouse matrifibrocytes using an AUC score. This signature was highly expressed in human cluster 6, cluster 4 and some cells in cluster 1 and cluster 5 (Fig. 1M). SFRP2, CILP and WISP2 were highly expressed in cluster 6, whereas COMP was mainly expressed in cluster 1 and cluster 5 (Fig. 1N). We simultaneously analyzed human cluster 3 (human RCFs) for DEGs analysis with human cluster 6 and human cluster 0, respectively (Figure S1H and S1I). GO enrichment analysis of human cluster 3 showed that the DEGs were mainly enriched in collagen fibril organization and extracellular matrix organization, while the DEGs in human cluster 6 were mainly enriched in cytoplasmic translation and peptide biosynthetic process (Figure S1J and S1K). In summary, by dynamically observing the changes in fibroblast subsets after MI, we observed four subsets of fibroblasts appearing after infarction, including Mki67+ fibroblasts, RCFs, Ccl2+ fibroblasts and matrifibrocytes. The signature of ECM-producing myofibroblasts (RCFs and matrifibrocytes) were also present in human cardiac single-cell transcriptome atlas.

Mki67+ fibroblasts contribute to RCFs during the acute phase of MI

To explore the differentiation trajectories of different subsets of fibroblasts during MI, we utilized the E-MTAB-7895 dataset for RNA velocity analysis since this method requires raw sequencing data. We divided these fibroblasts into three groups based on different periods: acute phase of MI (1, 3, 5, 7 days post-MI), subacute phase of MI (14, 28 days post-MI), and control group (non-surgical) (Fig. 2A and 2B). Here, we used two different models “stochastic” and “dynamical” to estimate RNA velocity, and we also used CellRank to compute a transition matrix based on RNA velocity. We calculated the rate of cell differentiation, which can be measured by velocity length. We first analyzed the fibroblast populations in the subacute and control phases. We found that cluster 0 in the control group may contribute to other fibroblast subsets, with no clear evolutionary sequence between the various subsets (Fig. 2C). This is consistent with the higher rate of differentiation of cluster 0 (Fig. 2D). Similarly, there was no significant differentiation trajectory at 14 and 28 days after MI (Fig. 2E). Of note, cluster 7 (matrifibrocytes), a subset of fibroblasts appearing mainly in the subacute group exhibited a higher rate of differentiation (Fig. 2F). In sum, our RNA velocity analysis results showed that the transitions between fibroblast subsets in the later stages of MI resembling the dynamics between fibroblast subsets observed in healthy myocardium.

Fig. 2

Mki67+ fibroblasts contribute to RCFs during the acute phase of MI. A UMAP plots of fibroblast clusters showing time distribution. B Proportion of fibroblast clusters showing in bar plots from control, stable and acute groups. C RNA velocity analysis showing the transition potential among fibroblast clusters at the control group, transcription dynamics based on the dynamical model (left), transcription dynamics based on the stochastic model (middle). Using CellRank’s VelocityKernel to compute transition matrix based on RNA velocity (right). D The velocity length (left) and velocity confidence (right) of the fibroblasts in control group. E RNA velocity analysis showing the transition potential among fibroblast clusters at subacute group, transcription dynamics based on dynamical model (left), transcription dynamics based on stochastic model (middle). Using CellRank's VelocityKernel to compute transition matrix based on RNA velocity (right). F The velocity length (left) and velocity confidence (right) of the fibroblasts in subacute group. G RNA velocity analysis showing the transition potential among fibroblast clusters at the acute group, transcription dynamics based on the dynamical model (left), transcription dynamics based on the stochastic model (middle). Using CellRank's VelocityKernel to compute transition matrix based on RNA velocity (right). H The velocity length (left) and velocity confidence (right) of the fibroblasts in acute group. I UMAP analysis of all fibroblast clusters at acute group embedded with PAGA connectivities for trajectory inference. J UMAP plots colored by the recovered latent time in the scRNA-seq datasets at acute group (left) and expression heatmaps for top 100 genes with fibroblasts ordered by latent time values. K UMAP plots of differentiation levels and the distribution of fibroblasts clusters (left). Boxplot showing the comparison of CytoTRACE score between different fibroblasts clusters (right)

We next focused on RNA velocity analysis of fibroblast subsets in the acute phase of MI. We found that the results of the RNA velocity analysis of the “stochastic” model were similar to those of the “dynamical” model and the transition matrix based on RNA velocity. In the period from 1 to 7 after MI, we found a transition from cluster 10 (Mki67 + fibroblasts) to cluster 5 (RCFs), immediately followed by a transition from RCFs to matrifibrocytes, cluster 4 (Fig. 2G). This is consistent with a previous study that matrifibrocytes are transformed from myofibroblasts [39]. We also observed a transition from cluster 9 to cluster 0 and cluster 2. Cluster 0, cluster 1 and cluster 3 are supposed to be the terminal subgroups of differentiation. Whereas there were no significant differences in differentiation rates between subsets of fibroblasts in the acute group (Fig. 2H). PAGA of RNA velocity provides a graph-like map of the data topology with weighted edges corresponding to the connectivity between two clusters [42]. The PAGA velocity graph also showed a direction from Mki67+ fibroblasts to RCFs (Fig. 2I). The above results suggested that Mki67+ fibroblasts, which were rapidly activated after MI injury, were converted to RCFs and participated in subsequent cardiac wound healing. However, latent time analysis showed no significant difference between RCFs and Mki67+ fibroblasts (Fig. 2J). We also used CytoTRACE to assess the differentiation potential of different fibroblast subsets and obtained similar results to the RNA rate analysis, with cluster 10 (Mki67+ fibroblasts), cluster 5 (RCFs), cluster 9 (Ccl2+ fibroblasts) and cluster 7 (matrifibrocytes) being the subsets with the highest differentiation capacity, representing an intermediate differentiation state (Fig. 2K). Finally, we used Slingshot to track the differentiation trajectory of fibroblast subsets. We set Mki67+ fibroblasts as our starting point because Mki67+ fibroblasts have a higher differentiation capacity compared to other subsets. Consistent with the results of the above analysis, RCFs were derived from Mki67+ fibroblasts and served as an intermediate for the eventual emergence of other fibroblast subsets (Figure S2A through S2E).

CTHRC1+ fibroblasts represent an activated fibroblast population in chronic disease state

Our results indicate that Cthrc1 was the top marker gene in RCFs, but RCFs were predominantly seen in the acute phase after MI injury. Previous studies have shown that CTHRC1 also appears in the late infarct stage and is expressed in human cardiac tissues from ICM and DCM patients [23]. To further explore the gene expression level of CTHRC1 in cardiac fibrosis and its correlation with the chronic disease state, we obtained the scRNA-seq datasets GSE155882 and GSE185265 related to cardiac fibrosis in mice from the GEO database. The GSE155882 dataset includes mice treated in four different subgroups, while GSE185265 contains three subgroups (Figure S3A). We first integrated the data from different treatment groups in GSE155882 and then extracted fibroblast subpopulations based on the marker gene of fibroblasts and performed dimensionality reduction clustering analysis (Fig. 3A). Fibroblast populations with high expression of Cthrc1 also had high expression of Postn, a classical marker of activated fibroblasts (Fig. 3B). DEGs analysis revealed a significant upregulation of Postn and Cthrc1 in cluster 6 (Fig. 3C). The expression of Postn was reduced in fibroblasts after treatment with JQ1, a drug that effectively reduces fibrosis and improves cardiac function [22], and this effect was attenuated after the withdrawal of JQ1 midway. Similar to the results for Postn, Cthrc1 expression in fibroblasts was reduced after JQ1 treatment (Fig. 3D). We used the same approach to integrate and process another scRNA-seq dataset GSE185265 (Figure S3B). Similar to the results of the TAC model, Cthrc1+ fibroblasts highly express Postn (Figure S3C). The expression level of Cthrc1 was significantly reduced in the late stage of MI, but it was still evident that Postn and Cthrc1 expression levels were reduced in TTg mice [21] (better heart function and less cardiac fibrosis) compared to control mice (Figure S3D and S3E). Our immunofluorescence results showed that some Postn-positive fibroblasts also expressed Cthrc1 at 7 days after MI and in TAC model hearts, but such fibroblasts were rarely present in cardiac tissues in the control group and at 28 days after MI (Fig. 3E). Taken together, the analysis of the two mouse scRNA-seq datasets suggested that in chronic disease states, fibroblast populations with high Cthrc1 expression (Cthrc1+ fibroblasts), a subset of Postn+ fibroblasts, are associated with worsened cardiac function and exacerbated fibrosis.

Fig. 3

Characterization of CTHRC1+ fibroblasts in a mouse model of heart failure. A UMAP plot showing the fibroblast and their subclusters. A total of 8 clusters were identified. All other types of cells were colored in gray. B Gene expression of Postn and Cthrc1. C Expression by cluster of known activated fibroblast related genes shown as violin plots in fibroblasts. D Expression by sample of Postn and Cthrc1 shown as feature plots in fibroblasts. E Spatial location of POSTN, CTHRC1 in the control group, the infarct zone at 7, 28 days post MI and 28 days post TAC. POSTN (green), CTHRC1 (red), Nuclei (DAPI, blue). Co-localizations are in yellow (arrows). Scale bars: 100 µm

To explore whether there are similar results in human heart tissue, we collected snRNA-seq datasets of human DCM and HCM (Figure S4A). We first integrated multiple DCM samples and then extracted fibroblasts which were finally divided into 5 subsets (Figure S4B). Compared with other subsets, cluster 5 and cluster 6 highly expressed POSTN, and were mainly present in the DCM group (Figure S4C and S4D). In addition, we found that CTHRC1 was mainly expressed in cluster 5 and highly overlapped with POSTN+ fibroblasts. For HCM, unbiased clustering grouped the fibroblasts into six clusters. Cluster 1 and cluster 5 increased considerably in HCM vs. CTRL (Figure S4F and S4G). Consistent with this, cluster 1 and cluster 5 were highly expressive of genes associated with activated fibroblasts such as POSTN and MEOX1 [22], meanwhile CTHRC1 was also predominantly expressed in relation to these two subsets (Figure S4H). Compared to controls, fibroblasts from the HCM group exhibited higher expression of CTHRC1 (Figure S4I). The results of GO enrichment analysis for cluster 1 and cluster 5 also showed that the DEGs were mainly enriched in the ECM organization, extracellular structure organization and collagen fibril organization (Figure S4J). Similar to the results obtained from the mouse scRNA-seq dataset, CTHRC1+ fibroblast in DCM and HCM often also have high expression of POSTN and represents activated fibroblasts.

Activator protein-1 is a common transcription factor in CTHRC1+ fibroblasts in chronic disease state

To better understand the characteristics of CTHRC1+ fibroblasts, we subjected CTHRC1+ fibroblasts to DEG analysis with all other fibroblast subpopulations. We took the intersection of the genes upregulated in the dataset of TAC mice with the genes whose expression was upregulated in the dataset of MI 3 month mice, and finally obtained 16 genes (Fig. 4A and 4B). These intersecting genes included Comp, Fn1, Cilp, Spp1 and Postn, which were associated with fibrosis. Similarly, the common genes that were highly expressed in CTHRC1+ fibroblasts in DCM and HCM compared to other fibroblast subpopulations also included profibrotic genes such as POSTN, COL1A1 and COL1A2 (Fig. 4C and 4D). We transfected the constructed adenovirus into cardiac fibroblasts to knock down the Cthrc1 gene in the cells, and the results of qRT-PCR showed that the Cthrc1 gene was successfully knocked down after infection (Fig. 4E and 4F). We found that expression of collagen-associated mRNA levels was reduced in Cthrc1 knockdown cells under Ang II stimulation compared to shCtrl group (Fig. 4G). Similarly, immunofluorescence results showed that the fluorescence intensity level of Col1a1 was lower in the shCthrc1 group compared to the shCtrl group, but the fluorescence intensity of Col3a1 was not significantly different between the two groups (Fig. 4H and I). We next analyzed the TFs of CTHRC1+ fibroblasts using pySCENIC in an attempt to discover their shared transcription factor profile. Common TFs for Cthrc1+ fibroblasts in the TAC model and 3 months after MI included Klf4, Klf2, Npdc1, Xbp1, Jund, Maff, Junb and Tfdp1 (Fig. 4J). TFs shared by CTHRC1+ fibroblasts in DCM and HCM were PKNOX1, FOXO1, CREB5, BBX, RFX2, BRF2, RXRG and FOS (Fig. 4K). Combining the results of these analyses we found that activator protein-1 (such as JUND, JUNB and FOS) may be a common transcription factor in CTHRC1+ fibroblasts. Altogether, the above results suggested that, CTHRC1+ fibroblasts in the chronic disease state tend to represent persistently activated fibroblasts and activator protein-1 may be involved in the regulation of CTHRC1+ fibroblasts.

Fig. 4

Gene expression profile of CTHRC1+ fibroblasts. A Volcano plot of DEGs in Cthrc1+ fibroblasts versus all other types of fibroblasts. B Shared differentially expressed genes and their PPI network. C Volcano plot of DEGs in Cthrc1+ fibroblasts versus all other types of fibroblasts. D Shared differentially expressed genes and their PPI network. E–F GFP expression level and mRNA expression level of CTHRC1 in the shCtrl and shCthrc1 groups. Scale bars: 100 µm G The mRNA expression level of Col1a1, Col3a1, Ctgf, Postn, Fn1 in the shCtrl and shCthrc1 groups (n = 3). H Col1a1 expression detected by cellular immunofluorescence staining (n = 4). Scale bars: 50 µm. I Col3a1 expression detected by cellular immunofluorescence staining (n = 4). Scale bars: 50 µm. J Heatmap showing top 50 TFs activity for TAC-fibroblast cluster 6 (Cthrc1+ fibroblasts) and MI-fibro cluster 8 (Cthrc1+ fibroblasts). Venn diagram showing TFs common to both disease models. K Heatmap showing top 50 TFs activity for DCM-fibroblast cluster 5 (CTHRC1+ fibroblasts) and HCM-fibro cluster 1 (CTHRC1+ fibroblasts). Venn diagram showing TFs common to DCM and HCM. *P < 0.05, **P < 0.01

Gpnmb+Fabp5+ macrophages possess characteristics of SAMs and SPP1hi macrophages present in the early stages of MI

Immune-fibroblast interaction is an important mechanism of MI. To study macrophage alterations after MI, we extracted macrophages from the integrated mouse scRNA-seq dataset, and unbiased clustering grouped the macrophages into 8 subclusters (Fig. 5A and 5B). Cluster 3 high expression of Lyve1 may represent a resident Lyve1+ macrophage (Fig. 5C). It has recently been suggested that SAMs express TREM2 and/or CD9 in liver fibrosis and lung fibrosis [7]. We therefore validated the expression levels of these genes in different macrophage subsets. We found that Cd9 and Trem2 were expressed in several subsets including cluster 0, cluster 1, cluster 2, and cluster 4 (Fig. 5C). Spp1+ macrophages, which have recently been suggested to be strongly associated with tissue fibrosis [43,44]. We found that Spp1 was highly expressed in cluster 0, cluster 1 and cluster 4 (Fig. 5C). Cluster 0 and cluster 4 began to expand mainly at 3 days after MI, peaking at day 7, but decreased dramatically on day 14, whereas cluster 1 expanded rapidly on day 1 after MI but decreased substantially on day 5 (Fig. 5D). DEGs analysis of the different subsets revealed that genes highly expressed in cluster 0 included Gpnmb, Fabp5, Ctsd, Syngr1, Trem2, Cd63 and Spp1 (Fig. 5E; Table S6). In concordance with this, compared with the cluster 1 and cluster 3, cluster 0 highly expresses Apoe, Gpnmb, Lgals3 and Spp1 (Fig. 5F and 5G). The temporal sequentiality of the appearance of cluster 1 and cluster 0 and their positional proximity in the UMAP plot led us to speculate whether they share similar transcriptional features. The results of pySCENIC analysis showed that 35 of the top 50 TFs were the same between cluster 0 and cluster 1 including core transcription factors such as Spi1, Cebpb and Irf8 that are associated with macrophage lineage commitment and differentiation [45] (Fig. 5H). We also discovered cluster 0 exhibits some characteristic LAMs (lipid-associated macrophages) transcripts [46,47]. Meanwhile, the expression of cluster 0 was characterized like the previously reported SAMs, which highly express GPNMB, FABP5, SPP1 and CD63 [7]. Thus cluster 0 represents the Gpnmb+Fabp5+ macrophage possessing both LAMs and SAMs related characteristics. Considering the overlap of marker genes between LAMs and SAMs, further investigations are warranted to understand whether LAMs and SAMs represent the same subsets of macrophages in tissues.

Fig. 5

Gpnmb+Fabp5+ macrophage possess characteristics of SAMs. A UMAP plot of macrophages from the integrated mouse scRNA-seq dataset. A total of 8 clusters were identified. B Heatmap showing the expression profiles of top 5 genes ranked by LogFC of each cluster. C Gene expression of Lyve1, Spp1, Cd9 and Trem2. D UMAP plots of macrophages from 0, 1, 3, 7, 14, and 28 days after MI colored by cluster. E Violin plots showing the expression of Trem2, Spp1, Cd9, Gpnmb, Fabp5 and Cd63 in macrophage subsets. F The volcano plots showing the DEGs between cluster 0 vs. cluster 1 (mouse). G The volcano plots showing the DEGs between cluster 0 vs. cluster 3 (mouse). H Heatmap showing the top 50 TFs activity for cluster 0 and cluster 1 and cluster 10 macrophages subsets (mouse). Venn diagram showing TFs common to both subsets. I The volcano plots showing the DEGs between cluster 6 vs. cluster 0 (human). J The volcano plots showing the DEGs between cluster 6 vs. cluster 1 (human). K Heatmap showing the TFs activity for cluster 1 and cluster 6 macrophages subsets (human)

Gpnmb+Fabp5+ macrophage was dramatically reduced in 14 days after MI, making us wonder if it plays a role mainly in the early stages after MI injury. To validate this hypothesis we additionally integrated three mouse scRNA-seq datasets that contained CD45+ cells isolated from the infarcted area of the heart (Figure S5A). We extracted macrophages from the integrated scRNA-seq dataset, sub-clustering of macrophages led to the identification of 10 clusters (Figure S5B). To compare macrophage subpopulations in the newly integrated scRNA-seq dataset, we constructed a gene expression signature containing the top 30 genes of the Gpnmb+Fabp5+ macrophage and then mapped it to the new macrophage dataset. This signature was highly expressed in cluster 2 (Figure S5C). Similar to previous results, cluster 2 appeared mainly in the 3 to 7 days period of MI and decreased significantly in the later stages (14 days post MI) (Figure S5D). At the same time, compared to other macrophage subpopulations cluster 2 highly expressed Spp1 and Gpnmb (Figure S5E). These results suggested that Gpnmb+Fabp5+ macrophage expands rapidly during the acute phase of MI and disappears in the later stages of the disease.

To explore whether Gpnmb+Fabp5+ macrophage transcriptional signature was conserved in the human diseased heart, we extracted macrophages from the above integrated human scRNA-seq dataset and performed unbiased clustering (Figure S6A). Cluster 0 high expression of LYVE1 represented cardiac resident macrophages, and the SPP1 was mainly expressed in cluster 6 (Figure S6B). CD9 and TREM2 were highly expressed in cluster 1 (Figure S6B). We then mapped the Gpnmb+Fabp5+ macrophage signature to the human macrophage dataset and found that these features were mainly expressed in cluster 1 and cluster 0 (Figure S6C). At the same time, we found that CD63 and FABP5 were also predominantly expressed in these two clusters (Figure S6D). Interestingly, we found that SPP1 was predominantly expressed in cluster 6, but not in cluster 1 and cluster 0, which is not consistent with the expression profile in the mouse dataset. We then divided the cells into three different groups based on their origin: control, MI and ICM. Cluster 6 (SPP1hi macrophages) were mainly present in the MI group, cluster 0 (LYVE1hi macrophages) was predominantly found in healthy cardiac tissue, whereas cluster 1 (CD9hiTREMhi macrophages) was predominantly found in ICM patient cardiac tissue. (Figure S6E). GO enrichment analysis of cluster 6 showed that the DEGs were mainly enriched in the regulation of cell migration, extracellular structure organization and external encapsulating structure organization (Figure S6F). The gene expression profiles between cluster 6 and cluster 1, and between cluster 6 and cluster 0 were highly differentiated (Fig. 5I, J). In concordance with this, we found that the TFs were also very different between cluster 6 and cluster 1 (Fig. 5K). Overall, the results of our analyses indicated that, unlike the mouse dataset, macrophages possessing the Gpnmb+Fabp5+ macrophage signature in the human dataset were mainly present in the ICM but not in the MI cardiac tissue. However, in both mouse and human cardiac tissues, macrophages that highly express SPP1 (SPP1hi macrophages) were almost exclusively present in the early stages of MI, not in the later stages of the disease (post-MI > 14 days or ICM group).

Dynamic transition of macrophage subsets in the acute phase of MI

We introduced RNA velocity analysis to profile the dynamics of different subpopulations of mouse macrophages at different stages after MI injury. RNA velocity analysis showed that the majority of macrophage subpopulations in the control group flowed toward cluster 2 and cluster 3 (Fig. 6A). The PAGA plot abstraction obtained similar results, with both cluster 4 and cluster 1 shifting toward cluster 3, and cluster 0 and cluster 1 shifting toward cluster 2 (Fig. 6C). The dynamic flow of macrophages during the subacute phase was similar to the results of the control group (Fig. 6D). The results of the PAGA plot abstraction show that cluster 6 and cluster 0 flow to cluster 3, and cluster 0 and cluster 5 flow to cluster 2 (Fig. 6F). In both the control and subacute groups, cluster 1 showed a higher rate of differentiation (Fig. 6B and 6E). Next, we focused on analyzing the dynamics of macrophages in the acute phase, which possesses a larger number of macrophages. We use the generalized dynamical model to understand the transcriptional dynamics of macrophages. RNA velocity predicted that cluster 0, cluster 5 and cluster 1 may be an intermediate state in macrophage dynamics (Fig. 6G). This is consistent with our above observation that these subpopulations appeared in the acute phase but disappeared in the subacute phase. The results of velocity length and velocity confidence suggested that cluster 1 has a faster speed of differentiation compared to other macrophage subpopulations, which was consistent with the fact that cluster 1 presented in large numbers on the first day after infarction but then decreased dramatically on the fifth day (Fig. 6H). The PAGA velocity graph also showed the transition from cluster 2, cluster 3, cluster 5 and cluster 6 to cluster 0, suggesting that at the acute stage cluster 0 represents a terminally differentiated subpopulation (Fig. 6I). We also used CytoTRACE to assess the differentiation capacity of different macrophage subpopulations, and similar to the results of our RNA velocity analyses, cluster 1 had a high differentiation capacity and represented a rapidly increasing population after MI (Fig. 6J).

Fig. 6

Dynamic transition of macrophage subsets in acute phase of MI. A UMAP plots of RNA velocity results (left) and the recovered latent time (right) in the single cell data sets at control group. B The velocity length (left) and velocity confidence (right) of the macrophages at control group. C UMAP analysis of all macrophages clusters at control group embedded with PAGA connectivities for trajectory inference. D UMAP plots of RNA velocity results (left) and the recovered latent time (right) in the single cell data sets at subacute group. E The velocity length (left) and velocity confidence (right) of the macrophages at subacute group. F UMAP analysis of all macrophages clusters at subacute group embedded with PAGA connectivities for trajectory inference. G UMAP plots of RNA velocity results (left) and the recovered latent time (right) in the single cell data sets at acute group. H The velocity length (left) and velocity confidence (right) of the macrophages at acute group. I UMAP analysis of all macrophages clusters at acute group embedded with PAGA connectivities for trajectory inference. J UMAP plots of differentiation levels and the distribution of macrophages clusters (left). Boxplot showing the comparison of CytoTRACE score between different macrophages clusters (right). K UMAP plots of cardiac macrophages and monocytes in acute stage, macrophages contain 8 clusters. L RNA velocity analysis showing the transition potential among macrophages clusters and monocytes. M UMAP plots of macrophages clusters and monocytes showing time distribution. N Slingshot differentiation trajectory analyses of macrophages-monocytes in the lineage 4. O Slingshot differentiation trajectory analyses of macrophages-monocytes in the lineage 3. P Violin plots showing the expression of Ccr2 and Spp1 in the macrophages-monocytes dataset

To explore whether cluster 0 (Gpnmb+Fabp5+ macrophages) originated from monocytes, we integrated macrophages and monocytes into one dataset for analysis (Fig. 6K). RNA velocity analysis did not give a definitive indication of the direction of cellular differentiation (Fig. 6L), and we subsequently observed cellular composition at different time points and found that Gpnmb+Fabp5+ macrophages appeared on days 5–7 post-MI, whereas most monocytes and cluster 1 macrophages appeared on days 1–3 post-MI (Fig. 6M). We then used Slingshot to infer the differentiation trajectory. Combining the results of the PAGA analysis and the time point at which Gpnmb+Fabp5+ macrophages appeared, we set it as the end point of the trajectory and set the monocyte as the start point. The results of Slingshot analysis showed four different differentiation trajectories. Lineage 3 depicts a trajectory from monocytes to cluster 1 and finally to cluster 4, whereas lineage 4 depicts a trajectory from monocytes to cluster 1 and later to cluster 0 (Fig.