Bioinformatics Analysis and Identification of Genes and Pathways in Ischemic Cardiomyopathy

Introduction

Heart failure (HF), a leading cause of high morbidity and mortality of cardiovascular disease and a global health care issue, has an estimated global prevalence of approximately 38 million.^1,2 Ischemic cardiomyopathy (ICM) accounts for more than 40–60% of systolic HF, which is characterized by prior myocardial revascularization or myocardial infarction (MI) and severe coronary heart disease.^3,4 The severity and duration of ischemia that contributes to myocardium injury and irreversible myocardium necrosis is common in the earlier phase of ICM pathogenesis, which results in cardiac fibrosis and remodeling and eventual HF.^5,6 To date, the detailed pathogenic mechanism in the progression of ICM has remained poorly understood, and efforts in drug development for ICM are currently focused on the symptomatic treatment and prognosis, such as angiotensin-converting enzyme inhibitor (ACEI), β receptor blocker, and so on.^7,8

Dyregulated hub genes, proteins and key pathways represent the genomic mechanism of the process of diseases.^9,10 Therefore, the mechanism of genes and proteins correlation and key pathways in the pathogenesis of ICM are necessary to improve the prevention, diagnosis, and treatment. Recently, with the continuous development of genomic technologies, microarray technology and RNA sequence have increased the ability of experts to investigate the regulation of genes of cardiovascular disease.¹¹ Several studies have reported that unique genes contribute to the pathogenesis of heart failure. An analysis of myocardial tissues in eight ICM patients and eight non-failing humans revealed 160 different expressed microRNA and 679 long non-coding RNA, which suggest an essential role for the expression signature of miRNA and lncRNA in the pathogenesis of ICM.¹² A previous study of 196 failing and 16 non-failing human hearts demonstrated that MEF2, NKX, GATA, and NF-AT transcription factors play a crucial role in the progression of HF.¹³ Although these databases have revealed differentially expressed genes (DEGs) in ICM, co-expressed DEGs by the integrated analyses of multiple datasets have not been thoroughly investigated.

In the present study, original data from microarray analyses and RNA sequences conducted on ICM patients’ heart samples were downloaded from the Gene Expression Omnibus (GEO) database, and a comprehensive analysis was implemented.

Then, the gene ontology (GO) and Kyoto Encyclopedia of Genes/Genomes (KEGG) pathway enrichment analyses and protein–protein interactions (PPI) were performed to explore key genes and the molecular mechanism involved in the pathogenesis of clinical ICM, which may provide therapeutic targets for treatment and prevention.

Materials and Methods

Acquisition of Data of Gene Expression Profiles

The datasets of gene expression profiles with the sequence numbers of GSE116250, GSE46224 and GSE5406 were downloaded from GEO, which performed microarray or RNA sequence and matched normal tissues. Totally 14 non-failing donors and 14 ICM patients’ left ventricular tissues consisted of GSE116250; in GSE46224, there were eight non-failing control and eight ICM patients' heart tissues; in GSE5406, the number of cases from ICM patients and normal humans was 108 and 16, respectively.

Screening of Differential Expressed Genes

R language was applied to analyze the initial data to identify DEGs. We utilized the Affy package to perform the data normalization, including conversion to raw data and background correction. The data were then subjected to DEGs analysis using the Limma package. P<0.05 and |log₂(FC)|>1 were set as the threshold, and the genes that met the criteria were considered DEGs. Finally, the “heatmap” package of “R” software was applied to construct heat maps, and the regions in which the differential genes were mainly concentrated were highlighted.

Gene Ontology (GO) Enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) Pathway Analysis

Annotations of cellular components, biological processes and molecular functions of DEGs were determined by using Gene Ontology (GO) enrichment analysis. The route of the gene cluster and related functions was determined using Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis. R cluster profile package was used to explore GO enrichment and KEGG pathway analysis with a cut-off criterion of adjusted P<0.05.

Protein–Protein Interaction (PPI) Network Analysis

In our study, to determine functional interactions between DEGs, protein–protein interaction (PPI) networks were created using the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) online. PPI networks were visualized using Cytoscape software. In a co-expressed network, Cytohubba, a plugin in the Cytoscape software, was adopted to calculate the degree of each protein node. In our study, the top ten genes were selected as hub genes.

RNA Extraction and cDNA Synthesis

Total RNA was extracted from heart tissues of ICM patients and the control using RNA isolation plus reagent (Takara, Osaka, Japan). Then, the concentration and quality of total RNA were examined by Thermo Fisher NanoDrop, the PrimerScript^TM RT Master Mix (Takara, Osaka, Japan) was used to reversely transcript the RNA into cDNA.

Real-Time Quantitative PCR

All the protocols and the use of human heart tissues were approved by the Xiangya Hospital of Central South University Institutional Review Board. Four patients diagnosed with ICM and two controls were enrolled in this study. The failing hearts were harvested from patients who were diagnosed with ICM and accepted heart transplantation. Normal hearts were obtained from the gift of hope and were used as controls. The detailed characteristics of the patients were listed in Supplementary Table 5. cDNA was generated with oligo (dT) primers and the Superscript first-strand synthesis kit (Invitrogen) from 1ug of total RNA as the manufacturer’s instructions. Then, cDNA was diluted and used for RT-PCR using SYBR Green following the application instructions. The primer sequences are summarized as in Supplementary Table 6.

Statistical Analysis

The results are expressed as Mean±SD. T-tests or nonparametric tests were used to analyze the scores of two different groups. A two-tailed value of P < 0.05 was considered statistically significant. Statistics were done using GraphPad Prism 8.0 and SPSS 20.0.

Results

Identification of DEGs in ICM Samples Compared with Not-Failing Samples

The basic information of the datasets related to ICM was shown in Supplementary Table 1. A total of 844 differential genes were screened out between the ICM and not-failing groups from GSE116250, of which 447 were up-regulated genes and 397 were down-regulated genes, respectively (P<0.05, fold change >2). The volcano plots and heat map displayed the DEGs from GSE116250 are shown in Figure 1A and D. A total of 99 DEGs were screened out between the ICM and non-failing groups from GSE46224, of which 58 were up-regulated genes and 41 were down-regulated genes, respectively (P<0.05, fold change >2). As shown in Figure 1B and E, the volcano plots and heat map displayed the DEGs from GSE46224. A total of 30 differential genes were screened out between the ICM and control group from GSE5406, of which 10 were up-regulated genes and 20 were down-regulated genes, respectively (P<0.05, fold change >2). The volcano plots and heat map displayed the DEGs from GSE5406 are shown in Figure 1C and F. The top 10 DEGs of the three datasets, respectively, are shown in Table 1. Subsequently, we mapped the Venn diagram of the differentially expressed genes in the three datasets, and identified eight co-expressed differentially expressed genes (co-DEGs) (Figure 2A). A total of eight co-DEGs are shown in Table 2.

Table 1 The Top 10 Differential Expressed Genes (DEGs) in the Left Ventricular Myocardial Tissue of Patients with Ischemic Cardiomyopathy in GSE116250, GSE46224, GSE5406

Table 2 Co-Expressed Differential Expressed Genes (co-DEGs)

Figure 1 Volcano plot and heat map of all the differentially expressed genes (DEGs) between Non-failure patients and heart failure patients in cardiac tissues. (A), volcano plot of DEGs in GSE116250 dataset; (B), volcano plot of DEGs in GSE46224 dataset; (C), volcano plot of DEGs in GSE5406 dataset; Green dot represent DEGs with fold change <-1, red dot represents DEGs with fold change >1, adjust P <0.05. (D), heat map of DEGs in GSE116250 dataset; (E), heat map of DEGs in GSE46224 dataset; (F), heat map of DEGs in GSE5406 dataset; The gradient colour change from blue to red represents the changing process from downregulation to upregulation.

Figure 2 Functional enrichment analysis of common differentially expressed genes (co-DEGs). (A), Venn diagram of co-DEGs; (B), Gene Ontology (GO) analysis of co-DEGs regard to biological process (BP), cellular component (CC), and molecular function (MF).

Gene Ontology Function Enrichment Analysis

GO analyses were performed on the detected eight co-DEGs to examine their biological functions in detail. The enrichment analysis related to biological processes were primarily enriched in collagen fibril organization (P<0.0001), response to transforming growth factor beta (P<0.0001), and platelet activation (P<0.0001), whereas those relating to cellular components, were mainly enriched in collagen-containing extracellular matrix (P<0.0001), extracellular matrix (P<0.0001), and the fibrillar collagen trimer (P<0.0001). In addition, molecular function analysis showed that co-DEGs were involved in extracellular matrix structural constituent (P<0.0001), platelet-derived growth factor binding (P<0.0001), and extracellular matrix structural constituent conferring compression resistance (P<0.0001) (Figure 2B). In the next step, we conducted separate GO analyses on these three databases. The DEGs of GSE116250 datasets related to biological processes were primarily enriched in extracellular structure organization (P<0.0001), extracellular matrix organization (P<0.0001), and collagen fibril organization (P<0.0001). Homoplastically, the enrichment analysis related to cellular components mainly were enriched in collagen-containing extracellular matrix (P<0.0001), extracellular matrix (P<0.0001), and basement membrane (P<0.0001), whereas those relating to molecular function analysis were mostly enriched in extracellular matrix structural constituent (P<0.0001), glycosaminoglycan binding (P<0.0001) and heparin binding (P<0.0001) (Figure 3A). The DEGs of GSE46224 datasets related to biological processes were primarily enriched in extracellular structure organization (P<0.0001), extracellular matrix organization (P<0.0001), and collagen fibril organization (P<0.0001). The enrichment analysis related to cellular components showed that DEGs were involved in collagen-containing extracellular matrix (P<0.0001), extracellular matrix (P<0.0001), and fibrillar collagen trimer (P<0.0001). In addition, the molecular function analysis was mainly enriched in extracellular matrix structural constituent (P<0.0001), heparin binding (P<0.0001) and glycosaminoglycan binding (P<0.0001) (Figure 3B). The DEGs of GSE46224 datasets with regard to biological processes were primarily enriched in extracellular structure organization (P<0.0001), collagen fibril organization (P<0.0001) and muscle system process (P<0.0001), and related to cellular components mainly were enriched in extracellular matrix (P<0.0001), collagen-containing extracellular matrix (P<0.0001), and fibrillar collagen trimer (P<0.0001). The molecular function analysis of DEGs was mainly focused on extracellular matrix structural constituent (P<0.0001), extracellular matrix structural constituent conferring compression resistance (P<0.0001), and platelet-derived growth factor binding (P<0.0001) (Figure 3C).

Figure 3 Gene Ontology (GO) analysis of all DEGs. (A) the top 10 enriched GO terms in biological process (BP), cellular component (CC), and molecular function (MF) of GSE116250 dataset; (B), the top 10 enriched GO terms in biological process (BP), cellular component (CC), and molecular function (MF) of GSE462224 dataset; (C), the top 10 enriched GO terms in biological process (BP), cellular component (CC), and molecular function (MF) of GSE5406 dataset.

Kyoto Encyclopedia of Genes and Genomes Pathway Analysis

To get in more information about the noteworthy pathways of DEGs, a Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis was exerted, and the top 12 significant pathways are shown in Figure 4.

Figure 4 Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways analysis of all DEGs. (A), the enriched KEGG pathways in co-DEGs; (B), the enriched KEGG pathways of DEGs in GSE116250 dataset; (C), the enriched KEGG pathways of DEG in GSE46224 dataset; (D), the enriched KEGG pathways in GSE5406 dataset.

A total of ten key pathways were found through the KEGG pathway analysis of co-DEGs, which were primarily enriched in protein digestion and absorption (P=0.0008), AGE-RAGE signaling pathway in diabetic complications (P=0.0009), and amoebiasis (P=0.0009) (Figure 4A). Respectively, the top 12 significant pathways of GSE116250 were mainly enriched in Rheumatoid arthritis (P<0.0001), Graft-versus-host disease (P<0.0001), and Human infection (P<0.0001) (Figure 4B). And the top 12 significant pathways of GSE46224 were mainly enriched in Malaria (P<0.0001), ECM–receptor interaction (P<0.0001), and Phagosome (P=0.0001). Only one pathway was enhanced to the GSE5406, which was related to the AGE-RAGE signaling pathway in diabetic complications (P=0.0008) (Figure 4D).

PPI Network Analysis

A PPI network analysis of DEGs was conducted with the strings database. As shown in Figure 5A, protein Collagen alpha-1(III) chain (COL3A1), Collagen alpha-2(I) chain (COL1A2), Lumican (LUM), Asporin (ASPN) and Matrix-remodeling-associated protein 5 (MXRA5) are interconnected among the eight co-DEGs, while the remaining three proteins had no remarkable interactions. 449 nodes, 69 nodes and 16 nodes were found to be interconnected in GSE11250, GSE46224 and GSE5406 severally. Collagen alpha-1(I) chain (COL1A1), COL3A1, COL1A2, matrix metalloproteinase-2 (MMP2), biglycan (BGN), collagen alpha-1(V) chain (COL5A1), connective tissue growth factor (CTGF), LUM, fibromodulin (FMOD) and proenkephalin A (PENK) were the top 10 hub genes in the GSE116250 dataset (Figure 5B, Supplementary Table 2). And COL1A1, COL3A1, COL1A2, BGN, LUM, periostin (POSTN), thrombospondin-2 (THBS2), FMOD, CTGF and Collagen alpha-1(XVI) chain (COL16A1) were the top 10 hub genes in the GSE46224 dataset in which there were seven hub genes the same as in GSE116250 (Figure 5C, Supplementary Table 3). In the GSE5406 dataset, COL3A1, COL1A2, osteoglycin (OGN), LUM, ASPN, Filamin-C (FLNC), myotilin (MYOT), myosin-6 (MYH6), plasminogen activator inhibitor 1 (SERPINE1) and pentraxin 3 (PTX3) were the top 10 hub genes in the GSE5406 dataset, which have three hub genes the same as both GSE116250 and GSE46224 datasets (Figure 5D, Supplementary Table 4).

Figure 5 Protein–protein interaction (PPI) network analysis. (A), PPI network of co-DEGs identified by Cytoscape; (B and C), PPI network and hub gene of DEGs in GSE116250; (D and E), PPI network and hub gene of DEGs in GSE46224; (G and H), PPI network and hub gene of DEGs in GSE5406.

Validation of the Eight co-DEGs in Myocardial Tissue of ICM Patients

Eight co-DEGs validated in GSE116250, GSE46224 and GSE5406 were chosen for qRT-PCR analysis (Table 2). As shown in Figure 6, 5 co-DEGs were down-regulated: SERPINA3 (NC 0.82±0.25 vs ICM 0.05±0.06, P<0.01), FCN3 (NC 1.35±0.49 vs ICM 0.25±0.16, P<0.05), LUM (NC 0.86±0.20 vs ICM 0.31±0.18, P<0.05), COL1A2 (NC 0.74±0.37 vs ICM 0.30±0.26, P=0.158) and ASPN (NC 0.63±0.52 vs ICM 0.47±0.32, P=0.533). In addition, gene HBB (NC 0.55±0.64 vs ICM 9.83±12.27, P=0.133) and COL3A1 (NC 0.81±0.27 vs ICM 0.97±0.61, P=0.752) were up-regulated, while MXRA5 (NC 0.99±0.21 vs ICM 1.01±0.58, P=0.800) was no significant difference between the two groups.

Figure 6 Relative gene expression verification of eight co-DEGs. The expression of SERPINA3 (A), FCN3 (B), LUM (C), COL1A2 (D), ASPN (E), HBB (F), COL3A1 (G), MXRA5 (H) were measured by real time-qPCR. *P < 0.05; **P < 0.01.

Discussion

Ischemic cardiomyopathy (ICM) is the main cause of HF with reduced ejection fraction (EF), accounting for approximately two-thirds of all cases, which is correlated with a worse prognosis.¹⁴ Although some new therapies, such as stem cell transplantation and new drugs, have emerged, ICM remains high in incidence and mortality.¹⁵ Elucidating the underlying mechanism of ICM may contribute to the early diagnosis and effective treatment of HF patients. Bioinformatics analysis provided information for the identification of different genes and proteins for the progress of ICM. In the present study, eight co-DEGs were identified in three datasets (GSE116250, GSE46224 and GSE5406). Analysis of GO enrichment and KEGG pathway were performed, and a PPI network was constructed to analyze further molecular processes underlying ICM progression. We selected SERPINA3, FCN3, LUM, COL1A2, COL3A1, ASPN, HBB, MXRA5 to verify their expression in ICM. Expression of five DEGs (SERPINA3, FCN3, COL3A1, HBB, MXRA5) in heart tissues using qRT-PCR results was consistent with our GEO analysis, while expression of three DEGs (ASPN, LUM, COL1A2) was opposite with GEO analysis. These results may contribute to a better understanding of the mechanism of ICM development.

Collagen types I and III (Col I and Col III) are the major fibrillar collagens produced fibroblasts, including cardiac fibroblasts of the heart.¹⁶ The triple-helical collagen molecules are assembled by independently coded α1(I), α2(I), and α1(III) chains.¹⁷ Previous studies demonstrated that Col I and Col III are involved in extracellular matrix protein (ECM) during the progression of cardiac remodeling in the failing heart.^18,19 Several lines of evidence indicate that expression of COL3A1 and COL1A2 coding for the two collagens mainly increased in cardiovascular diseases, such as ICM and dilated cardiomyopathy.^20–22 In our study, our results showed that COL3A1 was up-regulated in heart tissues of ICM compared with normal tissues, which showed a correlation with the cardiac remodeling with ICM.

Lumican (LUM) and Asporin (ASPN) belong to the small leucine-rich proteoglycan (SLRP) family consisting of a core protein with leucine-rich repeats and one or more attached glycosaminoglycan chains, which are considered to control the assembly of collagen fibers in the ECM.^23–25 LUM is an ECM localized proteoglycan associated with inflammatory conditions and known to bind collagens.^26,27 Previous study has indicated that LUM protein and mRNA levels are increased in cardiac tissues of patients and mice with HF compared with control hearts.²⁸ Moreover, a moderate lack of LUM attenuated cardiac fibrosis and improved diastolic function following pressure overload in mice.²⁹ ASPN, firstly identified as an extracellular secreted protein in 2001, is distinct from other class 1 SLRP family members because of its unique aspartate residues named the D-repeat.³⁰ It has been demonstrated that ASPN in the myocardium of the db/db mice was significantly higher than that in the control mice.³¹ Thus, ASPN could play an essential role in pathological cardiac remodeling. The findings suggest that LUM and ASPN may contribute to the fibrinogenesis of cardiac remodeling of ICM.

MXRA5, a member of the MXRA protein family, plays a vital role in cell adhesion and ECM remodeling.³² The protein is found in humans to be encoded by the MXRA5 gene, which contains seven leucine-rich repeats and 13 immunoglobulin-like C2-type domains related to perlecan with a molecular weight of 312KD and is reported to express in primates but not in rats or mice.^33,34 MXRA5 mRNA is increased in human chronic ischemic myocardium and primary fibroblast isolated from aging patients. Our results further confirmed that MXRA5 mRNA levels are enhanced in heart tissues of ICM patients, although there is no statistical significance compared with normal hearts. Recently, MXRA5 is found to be closely associated with the mitogen-activated protein kinases (MAPKs) pathways, which play a key role in the regulation of cell proliferation, survival, differentiation and apoptosis.³⁵ However, the role of MXRA5 in the progression of ICM remains unclear.

As with most bioinformatics network analysis studies of human diseases, this study had several limitations. The decreased expression of COL1A2, LUM and ASPN and no statistical difference of up-regulated expression of COL3A1 and MXRA5 were not consistent with the GEO analysis, which may be due to the limited number of heart tissues (n=4 ICM, n=2 normal hearts) for qRT-PCR confirmation and demands further investigation to increase the number of each group to confirm these results. In addition, the sample of GSE46224 was relatively small, and the sampling method did not eliminate the effect of gender and medications that may alter the gene expression. Although we detected the hub genes expression in cardiac tissues, the related pathways and mechanisms were not confirmed in in vitro studies or other functional studies, which need further investigation.

Conclusion

The present study was intended to identify DEGs and pathways with comprehensive bioinformatic analysis to find the potential biomarkers and predict the progression of ICM. We further confirmed COL3A1, COL1A2, LUM, ASPN and MXRA5 might be exploited as diagnostic and prognostic indicators for ICM. Finally, AGE-RAGE signaling pathway, ECM–receptor interaction and protein digestion and absorption may be essential signaling pathways in the development of ICM.

Therefore, our findings provided a potential basis to understand the cause and mechanism of ICM and therapeutic targets for clinic treatment.