3. 1 snRNA-seq data integration and clustering
snRNA-seq data from 17 mice hearts, including sham (n = 3), MI 1 h (n = 2), MI 1 d (n = 6), MI 3 d (n = 3), and MI 7 d (n = 3) heart samples18, were integrated by CCA for a follow-up analysis. After quality control and filtration, 162 047 cells were obtained (Supplementary Fig. S1). After defining principal components (nPCs = 20) and resolution (resolution = 0.2), a total of 11 cell clusters were identified by performing non-linear dimensionality reduction using the uniform manifold approximation and projection (UMAP) method (Fig. 1a,b). Besides, four cell types from 17 hearts were eventually separated based on specific markers (Fig. 1c). CMs exhibited a significant expression of marker genes commonly associated with cardiomyopathy19–21, including Actn2, Rbm20, and Fhl2 markers. Fibroblasts were identified by the presence of markers like Postn, Col1a1, and Col1a2 markers that were closely linked to cardiac fibrosis22. Myeloid cells were characterized by the expression of Mrc1, C1qa, and C1qb markers that were associated with immune response23. Endothelial cells were identified by the presence of Pecam1 and Cdh5 markers that maintain the vascular permeability barrier of the endothelial cell junction24 (Fig. 1d). Concurrently, the proportion of CMs and endothelial cells decreased, while the proportion of fibroblasts and myeloid cells relatively increased as MI progressed, which was consistent with the conventional understanding that cardiac fibrosis gradually intensified with the increased duration of MI. Notably, individual outliers in the data might be attributed to variations between samples (Fig. 1e,f).
3.2 snRNA-seq revealed the heterogeneity and subclusters of MI CMs
To investigate the gene expression profile underlying the heterogeneity of CMs under hypoxic conditions, our study specifically examined CMs derived from the sham and MI heart. Total CMs were clustered into 14 clusters by performing non-linear dimensionality reduction using UMAP analysis (Fig. 2a). Subsequently, 5 unique subclusters (C1-C5) were identified based on gene expression similarity (Fig. 2b). Each subcluster exhibited marker genes that characterized their distinct cell type. For instance, C1 showed high expression levels of PGC1α, Myh6, and Camk1d. C2 demonstrated high expression levels of Dmd, Trdn, Kbtbd12, Chrm2, and Slc25a13. C3 displayed increased expression levels of Pde8b, Lin7a, Pdgfrb, Trpc3, and Myo1b. C4 exhibited elevated expression levels of Kif26b, Itgbl1, Nox4, Pcdh9, and Xkr4. C5 showed elevated expression levels of mitochondrial transcripts, such as mt-Atp6, mt-Cytb, mt-Co3, and mt-Nd2, indicating a notable mitochondrial malfunction that could trigger CM apoptosis during this stage25–28 (Fig. 2c). Therefore, the UMAPs, which included both control and MI CMs, indicated that the highly expressed genes enriched in C5 of MI were potentially specific to MI.
We then conducted GSVA to investigate the possibly biological functions of the 5 CM subclusters. The results demonstrated that different CM subclusters were enriched in distinct biological processes (Fig. 2d). For example, the specific genes enriched in C1 were primarily associated with atrial cardiac muscle tissue morphogenesis, the response to corticotropin releasing hormone, and mitochondrial DNA metabolic process. It is possible that C1 was also involved in the repair of myocardial tissue injury. On the other hand, the specific genes enriched in C2 were predominantly correlated with the cardiac muscle myoblast proliferation, AV node cell to bundle of HIS cell signaling, regulation of voltage gated sodium channel activity, and cardiac muscle cell adhesion, indicating that C2 was possibly involved in cardiac electrical conduction and fibrosis. Furthermore, the specific genes enriched in C3 were mainly associated with the positive regulation of cardiac muscle adaption, suggesting indicating a potential relationship between C3 and post myocardial injury repair. Similarly, the specific genes enriched in C4 were principally correlated with positive regulation of cardiac muscle contraction, immune complex clearance, and regulation of cardiac vascular smooth muscle cell differentiation, indicating that C4 was plausible related to physical function of CMs. Lastly, the specific genes enriched in C5 were predominantly associated with the phosphagen metabolic process, mitochondrial ATP synthesis coupled protein transport, mitochondrial electron transport cytochrome C to oxygen, mitochondrial ADP transmembrane transport, NADH dehydrogenase complex assembly, and oxygen transport, indicating a close correlation between C5 and mitochondrial function, further investigation on this CM subcluster may prove beneficial in understanding the mechanisms of MI and exploring novel therapeutic targets.
3.3 C5 was the main subcluster related to mitochondrial malfunctional and cell apoptosis after MI
To investigate the cellular subgroups associated with functional deterioration and the plausible pathological mechanisms leading to cardiac dysfunction, the Seurat function AddModuleScore using 3 separate mitochondrial malfunction groups and 1 apoptosis group showed that C5 exhibited a significantly higher MD. Score for mitochondrial dysfunction (Fig. 3a,b) and a markedly higher AP. Score for cell apoptosis in comparison to C1-C4 subclusters (Fig. 3c,d), indicating that C5 might be the primary subcluster responsible for ventricular remodeling caused by mitochondrial malfunction and cell apoptosis.
3.4 Trajectory of CM subcluster remodeling after MI
To explore the transition among these CM subclusters, a pseudotime trajectory analysis was conducted on all CM subclusters. The results demonstrated that CMs were separated into distinct lineages by positioning normal CMs (C4) at the initiation point of the trajectory, as the duration of MI escalated (Fig. 4a). Meanwhile, it was observed that the endpoints of pseudotime trajectories, namely C2 and C3, represented distinct end states of CMs. In light of this, two primary transition lineages were identified based on the properties of various subclusters of CMs (Fig. 4b,c).
3.5 GO and KEGG analysis of C1-C5 subclusters
To better understand the transcriptional dynamics that occurred during the trajectories of CM remodeling, we simulated functional changes between C1-C5 subclusters of CMs by GO analysis based on the previously identified two lineages. The C4 subcluster in both lineages expressed specific genes such as Col8a1 and Itgbl1, which were primarily involved in cell-substrate adhesion. Interestingly, a set of genes (C5 cluster) and C5 cell subcluster exhibited significant consistency and specificity, with GO biological processes and functions focusing on mitochondrial transcripts and mitochondrial-related genomic genes including mt-Nd4, mt-Nd5, mt-Nd1, mt-Nd2, Ndufb8, Ndufb9, mt-Nd6, mt-Nd3, mt-Nd4l, mt-Co2, Ndufb11, Atp5b, Atp5a1, Atp5e, mt-Atp6, mt-Atp8, Atp5g1, Atp5g3, Pink1, and mt-Cytb. Meanwhile, C5 cell subcluster demonstrated GO terms related to oxidative phosphorylation, aerobic respiration, cellular respiration, energy derivation by oxidation of organic compounds, and ATP synthesis coupled electron transport. Concomitantly, these genes were predominantly enriched in pathways related to mitochondrial triphosphate electron, purine ATP mitochondrion proton permeability, coupled ribonucleotide reactive species chain, oxygen transport necrotic nucleoside, and permeabilization ribonucleoside synthesis (Fig. 5).
In addition, we simulated functional changes between C1-C5 clusters of CMs by KEGG analysis as well. The C4 subcluster in both lineages includes specific genes such as Col8a1, which were primarily involved in protein absorption. While the C5 subcluster exhibited markable consistency and specificity with C5 gene cluster as well, with KEGG pathological processes focusing on mitochondrial transcripts and mitochondrial-related genomic genes including Cox6a2, Cox8b, Cyc1, Uqcr11, Uqcrfs1, Uqcrq, mt-Co1, mt-Co2, mt-Co3, mt-Cytb, Ndufa4, Ndufb11, Ndufb8, Ndufb9, mt-Nd1, mt-Nd2, mt-Nd3, mt-Nd4, mt-Nd4I, and mt-Nd5. Concomitantly, C5 cell subcluster showed KEGG terms associated with oxidative phosphorylation, diabetic cardiomyopathy, chemical carcinogenesis, and thermogenesis diseases (Fig. 6). These findings were in line with the result that C5 exhibited the most significant MD. Score and AP. Score, aiding in identifying pathological mechanisms associated with the CM mitochondrial malfunction and apoptosis by exploring the gene expression profile of C5.
3.6 The cell-cell communication between C1-C5 subclusters and other cell types
In order to investigate the interactions between C1-C5 and other cell types, such as endothelial cells, smooth muscle cells, and myeloid cells, we conducted cell-cell communication analysis. The results revealed a close communication between C1, C2, C4, and endothelial cells, suggesting that C1, C2, and C4 likely belonged to subclusters that promoted angiogenesis. While C5 exhibited limited communication with other cells, indicating that C5 might lost its ability to interact with other cells due to mitochondrial malfunction and apoptosis (Fig. 7a). We further exampled the receptor-ligand interactions of various pathways. Our observations revealed that C1, C2, and C4 potentially influenced the VEGFR of endothelial cells. While C5 exhibited minimal communication with other cells as well (Fig. 7b), Additionally, our findings indicated that C4 interacted with myeloid cells in the CD86 signaling pathway, which aligned with the immune complex clearance function of C4 (Fig. 7c). C2 and C3 enhanced signal levels in the laminin pathway due to their interaction with fibroblasts, which corresponded to their position at the end of the trajectories and suggested their potential pro-fibrotic effects. Meanwhile, C1, C2, and C4 actively interacted with endothelial cells in the laminin pathway, indicating that C1, C2, and C4 might have a higher potential to promote angiogenesis (Fig. 7d). Besides, C5 displayed little communication with other cells in both the CD86 and CSF signaling pathway (Fig. 7c,e). These findings indicate that C1-C4 have close communication with other cells, while C5 exhibits limited communication with other cells due to its potentially dysfunctional state.
3.7 Slc25a4 was increased in CMs under hypoxic conditions
A set of genes showing increased expression levels in the C5 subcluster were identified after conducting DEA between the C4 subcluster and C5 subcluster, among which Scl25a4 was upregulated in C5 and reported to correlated with mitochondrial malfunction and cell apoptosis29,30 (Fig. 8a,b). The mitochondrial permeability transition pore (MPTP), composed of the Slc25a4 protein, is a large channel located on the inner membrane of mitochondria. When activated, the MPTP leads to a significant increase in osmotic pressure on the inner mitochondrial membrane, causing mitochondrial swelling, membrane rupture, and ultimately resulting in mitochondrial apoptosis during ischemic injury31,32. In our study, immunofluorescence staining demonstrated the co-localization of Slc25a4 with mitochondria in mice cardiomyocyte line (HL-1) (Fig. 8c). To further validate the effect of Scl25a4 on mitochondrial malfunction and apoptosis in hypoxic CMs, 10-week-old C57BL/6 wild type mice underwent sham or MI for 7 days. RT-PCR and western blotting analysis confirmed a high expression level of Scl25a4 in the hearts of MI mice compared to sham mice (Fig. 8d-f). Immunohistochemical and immunofluorescence staining for Scl25a4 in heart slices of sham and MI mice further validated these results (Fig. 8g-i), suggesting that Scl25a4 might play a crucial role in the heart's response to hypoxic stress.
3.8 Slc25a4 promoted mitochondrial malfunction and cell apoptosis in HL-1
To validate the function of Scl25a4 as a regulator of mitochondrial function and apoptosis, we conducted an experiment using siRNA targeting Scl25a4 to reduce its expression in HL-1 cells and using Cocl2 to trigger hypoxia in vitro. Western blotting analysis revealed that the impaired mitochondrial biogenesis and fusion as well as cell apoptosis caused by hypoxia were significantly improved after transfection with Scl25a4 siRNA (Fig. 9a,b). Mito-Tracker staining showed that Scl25a4 knockdown led to an increase in the total number of mitochondria and a reduction in the accumulation of damaged mitochondria in CMs (Fig. 9c). Additionally, the percentage of CM apoptosis under hypoxic conditions was measured using both TUNEL staining and flow cytometry analysis. Our findings indicated that Scl25a4 silence significantly reduced the percentage of CM apoptosis (Fig. 9d-g). Overall, these findings indicate that the knockdown of Scl25a4 greatly improves the mitochondrial function and cell viability of CMs under hypoxic conditions in vitro.
3.9 Slc25a4 was increased in AC16 and promoted mitochondrial malfunction and cell apoptosis
To further validate the role of Slc25a4 in regulating mitochondrial function and apoptosis in human hearts during hypoxia, the bulk RNA-seq dataset (GSE21610) was utilized. The analysis revealed a notably elevated expression of Slc25a4 in failing human hearts compared to non-failing hearts (Fig. 10a). Subsequently, human CM line (AC16) was subjected to Cocl2-induced hypoxic conditions, where various assays including RT-PCR, western blotting, and immunofluorescence demonstrated a significant upregulation of Slc25a4 in hypoxic CMs (Fig. 10b-e). Moreover, when AC16 cells were transfected with Slc25a4 siRNA prior to Cocl2 treatment, Mito-Tracker and TUNEL staining indicated an improvement in mitochondrial function and a reduction in cell apoptosis under hypoxic conditions, respectively (Fig. 10f-h). Overall, these results provide additional confirmation of the involvement of Slc25a4 in the pathogenesis of human failing hearts.