Identifying Immune Cell Infiltration and Hub Genes During the Myocardial Remodeling Process After Myocardial Infarction

Introduction

Acute myocardial infarction (MI) is the leading cause of global morbidity and mortality. About 30% of MI patients undergo myocardial remodeling after MI, and these patients have a higher risk of mortality.¹ The myocardial remodeling after MI is a complex reparative process with 3 overlapping phases: the inflammatory phase, proliferative phase, and healing phase.² In the early phase of infarction, the inflammatory response is activated, and immune cells then accumulate in the ischemic infarct area to participate in the inflammatory cascade.³ Meanwhile, resident cardiac fibroblasts are activated by early tissue damage.⁴ During the proliferative and healing phases, these cells proliferate, migrate into wounded tissue, and only differentiate into myofibroblasts, which initiate and promote scar formation and fibrosis.⁴ Moderate inflammation response helps cardiac repair, whereas excessive inflammation worsens cardiac injury and post-MI myocardial remodeling.⁵ Current treatments for MI consist solely of cause-dependent intervention, which is incapable of minimizing myocardial necrosis and maximizing post-MI heart healing. In addition to revascularization, new treatment strategies should consider the molecular mechanisms that regulate post-MI myocardial remodeling.⁶ The discovery of new collagen genes, essential plasma proteins, and a novel macrophage phenotype has enhanced our understanding of the molecular mechanisms underlying post-MI myocardial remodeling in humans.^7–9 Nevertheless, the underlying mechanisms of post-MI myocardial remodeling remain unclear.

Over the years, bioinformatic analysis based on RNA-sequencing (RNA-seq) data has been increasingly applied to the identification of the key genes, yielding crucial insights into underlying mechanisms of disease progression.^10–12 Moreover, CIBERSORTx, a computational algorithm for characterizing relative fractions of diverse cell subsets of gene expression profiles,¹³ is increasingly being used to estimate the immune cell composition of complex tissues from their RNA-seq data.^14–16 However, RNA profiling is limited to ensemble-based approaches that cannot provide biological insights into specific cell types.¹⁷ With the emergence and development of single-cell techniques, researchers can now explore the function and role of distinct cell types in disease progression.^17,18 Nevertheless, the number of detected genes based on single-cell techniques is significantly lower than that of bulk RNA-seq.¹⁸ Thus, the integration of RNA-seq and single-cell techniques is essential for elucidating the biomarkers and molecular mechanisms underlying the development of disease.¹⁷

To investigate hub genes and immune cell landscape during myocardial remodeling after MI, we performed a comprehensive analysis of RNA-seq data of human post-MI myocardial remodeling, which identified differential expression genes (DEGs), explored the underlying biological functions and pathways, and further screened hub genes by constructing protein-protein interaction (PPI) network. The expression of hub genes was then validated in the RNA-seq data of post-MI mice. Next, we used CIBERSORTx to determine the percentage of immune infiltration of human post-MI myocardial remodeling and performed immune infiltration correlation analyses of hub genes. Furthermore, our candidate gene P3H3 was further validated using immunohistochemistry (ICH) and Quantitative Real-Time Polymerase Chain Reaction (qRT-PCR). Finally, single nuclei RNA sequencing (snRNA-seq) data of post-MI mice was used to determine the expression and distribution of hub genes. The workflow of this study is shown in Figure 1.

Figure 1 Workflow chart of the study.

Materials and Methods

Data Acquisition

Two RNA-seq datasets (GSE132143 and GSE151834) and an snRNA-seq dataset (GSE176092) were acquired from the GEO database (https://www.ncbi.nlm.nih.gov/geo/). The GSE132143 dataset contained 20 post-MI myocardial remodeling samples and 12 normal samples from the human ventricle. The GSE151834 dataset, containing 4 injured tissues and 4 uninjured remote tissues from 4 post-MI (2 weeks) mice ventricles, was used as an independent validation set. The snRNA-seq dataset GSE176092, containing filtered cardiomyocytes of 2 injured tissues and 2 uninjured remote tissues from 2 post-MI (2 weeks) mice ventricles, was used to clarify the expression and distribution of hub genes. The uninjured remote tissues from post-MI mice can be considered as the control samples of infarcted tissues.¹⁹

Processing of the Datasets and DEGs Screening

The dataset GSE132143 was translated into the gene symbol by the R package “org.Hs.eg.db”. Then, it was log2 transformed and quantile normalized by the “normalizeBetweenArrays” function of the “limma” R package.^20,21 Principal component analysis (PCA) was calculated using “prcomp” function to verify the repeatability of samples. Next, DEGs in the post-MI myocardial remodeling were screened by using the “limma” R package²¹ with cut-off criteria of the false-discovery rate (FDR, Benjamini-Hochberg correction method) < 0.05 and |log fold change (FC)| > 1.5.

Functional Enrichment Analysis

Gene ontology and signaling pathways associated with DEGs have a significant impact on the development of post-MI myocardial remodeling. GO analysis and KEGG pathway analysis were performed by the R package “clusterProfiler”.²¹ FDR < 0.05 was considered statistically significant.

Identification of Protein-Protein Interaction (PPI) Networks and Hub Genes

The construction and visualization of the PPI network based on DEGs were carried out with the online STRING database (http://string-db.org/).²² The combined score > 0.4 was set as the threshold value. Next, the result of the PPI network was further imported into Cytoscape software (version 3.10.0) to identify critical subnetworks and hub genes through the plug-ins in Cytoscape.²³ Specifically, Molecular Complexity Detection (MCODE),²⁴ a plugin in Cytoscape, was utilized to identify the most important modules based on the PPI network. Degree Cutoff = 2, Node Score Cutoff = 0.2, K-core = 2, Max. depth = 100 and MCODE scores ≥10 were set as filtering parameters. Furthermore, CytoHubba,²⁵ a plugin in Cytoscape, was utilized to calculate the top 10 core genes of the PPI network based on Maximal Clique Centrality (MCC) and Density of Maximum Neighborhood Component (DMNC) algorithms. The hub genes were determined by taking the overlapping genes of MCODE (significant module 1), MCC, and DMNC algorithms.

The Accuracy of the Hub Genes

The ROC curves were generated by using the “pROC” package,²⁶ while the area under the curve (AUC) was calculated to evaluate the predictive accuracy of the hub genes for myocardial remodeling after MI. An AUC > 0.9 indicated that the hub genes had a good fitting effect.

Immune Cell Landscapes and Immune Correlation Analysis

The CIBERSORTx¹³ (https://cibersortx.stanford.edu/) was used to analyze the fractions of infiltrating immune cells in the GSE132143 samples with gene signature file LM22 set at 100 permutations. Spearman correlation analysis performed by “corrplot” package was used to evaluate the association of immune cells with each other. The associations of the hub genes with infiltrating immune cells were evaluated using Pairwise Spearman correlation analysis with family-wise error rate (FWER, Holm correction method). The resulting associations were visualized using the “ggplot2” package. FWER < 0.05 was considered statistically significant.

Animal Models

The infarct tissue from 2 weeks post-MI mice ventricular was defined as being in a late stage of myocardial remodeling.²⁷ Eight-week-old male C57BL/6 mice used in this study were purchased from SiPeiFu (Beijing, China) and randomly assigned to the MI (n=6) or sham group (n=4). As previously described, the MI model was constructed by permanent ligation of the left anterior descending (LAD) artery.²⁸ After chest hair removal, rodents were placed in an anesthesia box for rapid isoflurane (2%) anesthesia and then placed supine on the experimental plate. Following a left-sided thoracotomy, the heart was swiftly removed. The LAD artery was ligated 2 mm from its origin using a silk suture of size 6–0. After returning to the heart and expelling the gas from the chest cavity, the chest was closed with a 5–0 silk suture. The sham group received the same surgical procedure, but coronary artery ligation was not performed. All mice were housed in an environment with a light/dark cycle of 12/12 h and had free access to food and water. All animal experiments were approved by the Institutional Animal Care and Use Committee of Shanghai General Hospital (NO.2023AW004) and followed the standards of the National Institutes of Health Guide for the Care and Use of Laboratory Animals.

Echocardiography

See Supplemental Material.

Immunohistochemical Analysis

See Supplemental Material.

qRT-PCR Assay of the P3h3

See Supplemental Material.

SnRNA-Seq Data Analysis

See Supplemental Material.

Statistical Analysis

The data are expressed as the mean ± SD. Unpaired Student’s T-test and Mann–Whitney U-test were applied for comparisons between the two groups. All statistical analyses were performed via R software 4.2.0 and GraphPad Prism 9. P value < 0.05 was considered statistically significant.

Results

Identification of DEGs in Myocardial Remodeling After MI

To determine any significant differences in the gene expression profile among samples in GSE132143, we analyzed data distribution before and after quantile normalization through box plots. After normalization, the differences among samples were eliminated (Supplementary Figure 1A and B). We further validated the interpretability and scalability of the data based on PCA dimension reduction. PCA demonstrated that post-MI myocardial remodeling samples clustered together (Figure 2A). Principal component-1 (PC1) separated the transcriptional fingerprints of the two groups, suggesting that there were differences in gene expression between the two groups (Figure 2A). A total of 566 up-regulated DEGs and 409 down-regulated DEGs were identified (Figure 2B). The heatmap was generated to visualize the DEGs in the meta-data cohort (Figure 2C).

Figure 2 Identification of DEGs. (A) PCA analysis was performed on all the myocardial samples. (B) Volcano plot for DEGs. The number of significantly down-regulated (blue) and up-regulated (red) DEGs was shown on the top. (C) Heatmap for the top 50 upregulated and 50 downregulated DEGs.

Functional Enrichment Analysis

To investigate the potential molecular mechanisms and pathways of post-MI cardiac remodeling, we performed GO and KEGG analyses on DEGs. GO enrichment analysis revealed that DEGs in the biological process (BP) were mainly enriched in the cellular divalent inorganic cation homeostasis, modulation of chemical synaptic transmission, regulation of trans-synaptic signaling, extracellular matrix organization, extracellular structure organization, external encapsulating structure organization (Figure 3A). The results of cellular component (CC) were significantly enriched in the collagen-containing extracellular matrix, presynapse, postsynaptic specialization, and neuron-to-neuron synapse (Figure 3A). As for molecular function (MF), the DEGs were mainly focused on extracellular matrix structural constituents, glycosaminoglycan binding, and immune receptor activity (Figure 3A). KEGG analysis revealed that the DEGs were mainly enriched in pathways of neutrophil extracellular trap (NET) formation, protein digestion and absorption, phagosome, and extracellular matrix (ECM) receptor interaction (Figure 3B).

Figure 3 Functional enrichment analysis of DEGs. (A) GO analysis on DEGs. (B) KEGG analysis on DEGs.

Identification and Validation of Hub Gene

A total of 975 DEGs were used to construct a PPI network using the STRING online platform. There were 680 nodes and 4740 edges (Figure 4A). To screen the hub genes, the critical subnetworks were extracted with an MCODE plug-in based on the PPI network. Three significant modules (subnetworks) closely related to the PPI network were found and separated (Supplementary Table 1). The most significant module 1 composed of 14 nodes (Figure 4B). Then, the Cytoscape software was also used to analyze the PPI network. The critical subnetworks and hub genes of MCC and DMNC topological algorithms are shown in Figures 4C, D and Supplementary Table 2. Among the results of MCODE (significant module 1), MCC, and DMNC algorithms, the overlapping genes were as follows: P3H3, COL15A1, COL16A1, and COL27A1 (Figure 4E). Furthermore, to demonstrate the accuracy of these hub genes, the dataset GSE151834 was used for external verification. The results revealed that the expression levels of P3h3, Col15a1, Col16a1, and Col27a1 were significantly higher in the injured group than in the uninjured remote group (Figure 5). These 4 hub genes may serve as significant biomarkers in the myocardial remodeling process.

Figure 4 Identification of hub genes. (A) PPI network of DEGs. The no connected dots were hidden. (B) The 14 protein nodes of significant module 1 were generated by the MCODE algorithm. (C) Top 10 protein nodes generated by the MCC algorithm. (D) Top 10 protein nodes generated by the DMNC algorithm. (E) Hub genes identified by Venn diagram. (B-D) The deeper the color of the dot indicates the more advanced ranking of the hub genes.

Figure 5 Validation of the expression of the four hub genes in the GSE151834 dataset. P-values were calculated using an unpaired Student’s T-test. ***P < 0.001; ****P < 0.0001. Abbreviation: RPKM, Reads Per Kilobase of exon model per Million mapped reads.

The Accuracy of Hub Genes in Myocardial Remodeling After MI

The diagnostic ability of the four hub genes to distinguish post-MI myocardial remodeling samples from the normal samples exhibited a favorable fitting effect, with AUCs of 1.000, 0.946, 0.996, and 1.000 for P3H3, COL16A1, COL15A1, and COL27A1 in the meta-data cohort, respectively (Figure 6A). Moreover, to confirm the discrimination ability of those hub genes, the AUCs of all four hub genes, P3h3, Col16a1, Col15a1, and Col27a1, were 1.000 in the GSE151834 dataset (Figure 6B).

Figure 6 The ROC curve of the diagnostic effectiveness of the four hub genes. (A) ROC curve of P3H3, COL15A1, COL16A1 and COL27A1 in the metadata cohort. (B) ROC curve of P3h3, Col15a1, Col16a1 and Col27a1 in the GSE151834 dataset.

Immune Cell Landscapes and Immune Correlation Analysis

The immune microenvironment plays an important role in myocardial remodeling after MI.³ To further explore the relationship between hub genes and immune cells, we studied the subpopulations of infiltrated immune cells in post-MI myocardial remodeling and normal samples. The abundance distribution of 22 infiltrated immune cell types in each sample was presented in a histogram (Supplementary Figure 2A) and heatmap (Supplementary Figure 2B). Post-MI myocardial remodeling samples had significantly higher proportions of CD4+ naive T cells (p < 0.05), regulatory T cells (p < 0.05), and monocytes (p < 0.01) than normal samples. The proportions of M2 macrophages (p < 0.0001) and neutrophils (p < 0.01) were significantly lower in post-MI myocardial remodeling samples than in normal samples (Figure 7A). Interestingly, both post-MI myocardial remodeling and normal samples lacked activated mast cells, CD4+ memory-activated T cells, and gamma delta T cells. Furthermore, both P3H3 and COL16A1 were positively correlated with Regulatory T cells (FWER < 0.05) and monocytes (FWER < 0.05), whereas they were negatively correlated with M2 macrophages (both FWER < 0.01) (Figure 7B). The M2 macrophages showed a significantly negative correlation with monocytes (r = −0.6, p < 0.001) (Figure 7C).

Figure 7 Landscape of immune infiltration between two groups and co-expression patterns among immune cell components. (A) Comparison of immune infiltration between the two groups. The non-significant p-values were hidden. (B) Correlation between the hub genes and 19 infiltrating immune cells. (C) Correlation matrix of 19 immune cells. Red: negative correlation; blue: positive correlation. P-values were calculated using a Mann–Whitney U-test. *P < 0.05; **P < 0.01; ***P < 0.001; ****P < 0.0001.

Overexpression of P3h3 in Post-MI Myocardial Remodeling Samples

After reviewing the literature, three hub genes including Col15a1, Col16a1, and Col27a1 have been reported to increase in infarcted tissue from 14-day post-MI mice via qRT-PCR.¹⁹ Therefore, we choose to validate the expression of P3h3 by IHC and qRT-PCR analyses. Figure 8A depicts echocardiographic images representative of the MI and sham samples. The cardiac function confirmed by echocardiography showed a significant reduction in contractility, as reflected by left ventricle EF (29.7 ± 3.9% versus 55.0 ± 4.6%, p < 0.001) and FS (14.0 ± 2.0% versus 28.3 ± 3.0%, p < 0.001) in the MI group compared with that in the sham group (Figure 8B and C). As depicted in Figure 8D, the cardiac tissue in the MI group showed significant myocardial infarction, but not in the sham group. The positive area of P3h3 in the MI group was significantly higher (p < 0.0001) than that in the sham group (Figure 8E and F). Moreover, the P3h3 mRNA expression level in the MI samples was significantly higher (p < 0.001) than that in the sham samples (Figure 8G).

Figure 8 Cardiac function and validation of P3h3 gene. (A) Representative echocardiographic images. (B) Left ventricle ejection fraction (n = 4). (C) Left ventricle fractional shortening (n = 4). (D) Representative images of hematoxylin and eosin staining. Magnification 20×, scale bar = 100μm (n = 4). (E and F) Representative images and quantification of P3h3 immunohistochemical staining. Magnification 20×, scale bar = 100μm (n = 4). (G) The expression level of P3h3 by qRT-PCR (n=4). P-values were calculated using an unpaired Student’s T-test. ***P < 0.001; ****P < 0.0001.

Expression Prediction of Hub Genes in snRNA-Seq Profiling

To elucidate the expression and distribution of hub genes in post-MI myocardial remodeling, we analyzed the data of filtered cardiomyocytes from a public snRNA database. According to the UMAP algorithm, cells were successfully classified into 6 separate clusters (Figure 9A), which include cardiomyocytes (Ryr2), endothelial cells (Pecam1), fibroblasts (Dcn), pericardial cells (Upk3b), myeloid cells (C1qa), and pericytes (Rgs5) (Supplementary Figure 3). The marker genes in 6 cell clusters are listed in Supplementary Figure 3. We found that the P3h3, Col16a1, Col15a1, and Col27a1 were highly expressed in fibroblasts of post-MI myocardial remodeling (Figure 9B). Moreover, P3h3 and Col16a1 were highly expressed in pericytes and pericardial cells of post-MI myocardial remodeling, respectively (Figure 9B).

Figure 9 The analysis of expression level and distribution of hub genes based on snRNA-seq data. (A) The UMAP algorithm was applied and 6 cell clusters were successfully classified. (B) The expression and distribution of hub genes in 6 cell clusters of two groups.

Discussion

Although significant progress has been made in the treatment of MI, strategies to restore damaged tissue and reconstruct organ function are still required. Therefore, it is necessary to explore the specific mechanisms of post-MI myocardial remodeling to motivate the development of new therapeutic strategies and thereby improve the prognosis of patients with MI. In this study, we identified hub genes and the landscape of immune infiltration in human post-MI myocardial remodeling.

We first identified 975 DEGs of post-MI myocardial remodeling by bioinformatics analysis. The potential biological functions of these DEGs were investigated by GO and KEGG enrichment analyses, which mainly indicated the key roles of inflammatory response and extracellular matrix. The dynamic alterations in ECM composition may contribute to cardiac repair following MI.²⁹ Some approaches that strictly regulate ECM deposition and crosslinking may prevent adverse remodeling and heart failure after MI.³⁰

Subsequently, among the significant DEGs, we identified four hub genes for post-MI cardiac remodeling, including P3H3, COL15A1, COL16A1, and COL27A1, and validated them using the GSE151834 dataset. COL15A1 is essential for the organization of the extracellular matrix and structural and functional microvessels in the heart that can promote myocardial proliferation.^31,32 It can modulate the elasticity of the myocardium, and its deficiency leads to left ventricular dilatation and left ventricular wall thinning.³¹ COL16A1 is a linker protein that contributes to the organization of large fiber networks, thereby regulating the integrity and stability of the ECM. It is positively associated with left ventricular dysfunction in patients with ischemic cardiomyopathy.⁷ COL27A1 is known to be regulated by SOX9.³³ The role of COL27A1 in post-MI myocardial remodeling requires additional study. P3H3, a member of the prolyl 3-hydroxylase family, is involved in collagen post-translational modifications. Database analyses of vertebrate genomes showed that the prolyl 3-hydroxylase family had three isoenzymes named P3H1, P3H2, and P3H3.³⁴ Although the prolyl 3-hydroxylase family was discovered in collagens more than half a century ago, its exact functions remain unclear.^35,36 P3h3 has been reported to form a complex with Plod1 and P3h4 that is necessary for the normal assembly and cross-linking of collagen fibrils.³⁷ P3h3^−/− mice displayed severe under-hydroxylation of lysines, especially at cross-linking sites in type I collagen from the bone, cornea, and aorta.³⁸

In addition, we observed significant differences in the abundance of infiltrating immune cell subpopulations between the two groups, emphasizing the role of immune cells in the development of post-MI myocardial remodeling. In our study, the proportion of CD4+ naive T cells,³⁹ regulatory T cells, and monocytes⁴⁰ increased, whereas the proportion of M2 macrophages and neutrophils decreased. It has been demonstrated that regulatory T cells in the post-MI-damaged myocardium enhance cardiac function, promote scar repair, and inhibit ventricular remodeling.^41,42 CD4+ T cells become activated after MI, thus promoting wound healing and scar formation.⁴³ Nonetheless, persistently increased CD4+ T cells regulated by TNFR1/ TNFR2 promote ventricular remodeling and progressive cardiac dysfunction during ischemic heart failure.^39,44 Moreover, neutrophils can promote M2 macrophage polarization.⁴⁵ M2 macrophages can produce multiple anti-inflammatories, proangiogenic, and pro-reparative factors, as well as engulf and clear dead cells, thereby enhancing infarct healing and reducing adverse cardiac remodeling following acute MI.⁴⁶ Our bioinformatics analysis highlights the importance of these cells in the pathogenesis of myocardial remodeling after MI in humans.

Given that the expression of P3H3 has not been described in other studies, we further verified the elevated expression of P3h3 using immunohistochemical and qRT-PCR analyses. Additionally, snRNA-seq data analysis revealed that P3h3, Col15a1, Col16a1, and Col27a1 were highly expressed in fibroblasts in post-MI myocardial remodeling, among which COL15A1 is expressed in fibroblasts of the human heart.^32,47 A recent study showed that P3H3 may be associated with cardiac fibroblast activation mediated by transforming growth factor beta.⁴⁸ We hypothesize that P3H3 plays a crucial role in the assembly mechanism and molecular packing of collagen fibrils following cardiac injury. However, this may require further in vitro experiment validation.

This research has several limitations. The precise mechanisms of P3H3, COL15A1, COL16A1, and COL27A1 on the development of myocardial remodeling after MI remain unclear and require further investigation. In addition, the expression and distribution of these genes in subpopulations of immune cells have yet to be determined.

In this study, we assess the immune cell landscape in post-MI myocardial remodeling by using RNAseq data from human heart tissue, and identify P3H3 expression after MI. These findings increase our understanding of the role of immuno-microenvironment in cardiac injury repair and emphasize the significance of collagen post-translational modifications in myocardial remodeling. However, the impact of the prolyl 3-hydroxylase family on post-MI myocardial remodeling is unknown and should be the subject of future research.

Conclusion

Our study identified four hub genes, including P3H3, COL15A1, COL16A1, and COL27A1, as well as the immune cell landscapes in post-MI myocardial remodeling. These findings may contribute to a greater understanding of the mechanism underlying the development of myocardial remodeling after MI and provide novel ideas for the drug treatment of post-MI myocardial remodeling.