Identification of Independent and Communal Differentially Expressed Genes as Well as Potential Therapeutic Targets in Ischemic Heart Failure and Non-Ischemic Heart Failure

Introduction

Heart failure (HF) is the inability of the heart to deliver sufficient blood to meet the demands of our body under normal filling pressure, which ultimately leads to a complex and severe disease syndrome.^1,2 HF is a rapidly growing public health issue with an estimated prevalence of 37.7 million individuals globally, which confers a substantial burden to the health-care system.^2,3 Although the classification system for HF causes remain debated, it’s undeniable that non-ischemic heart failure (NIHF) and ischemic heart failure (IHF) are the two main etiological categories of HF.⁴ IHF, the most common type, is defined as left ventricular (LV) systolic dysfunction due to the coronary artery disease (CAD).⁵ NIHF is the other heart failure that excludes CAD factors, including all remaining heterogeneous HF etiologies, like valvular diseases, toxic damage, metabolic conditions and genetic cardiomyopathies.^6,7 Recent studies have shown that the expression levels of cytokines, relating to HF, are significantly different between IHF and NIHF.⁷ Meanwhile, large-scale clinical trials have reported that IHF and NIHF respond differently to interventional drug therapies.^4,8 However, the exact biological mechanism of above is still unclear. Interestingly, although the etiology is different between IHF and NIHF, the pathology is very similar, especially in patients with end stage of HF.⁴ So, there could be a common biological mechanism in the development of IHF and NIHF.

Therefore, the purpose of this study is to use Bioinformatic analysis to characterize distinct and independent pathways for IHF and NIHF, and more importantly, to identify the common biological mechanisms of IHF and NIHF, so as to reveal potential biomarkers and therapeutic directions of HF. In this study, two mRNA microarray data sets were downloaded from the Gene Expression Omnibus (GEO) for screening the differentially expressed genes (DEGs) associated with IHF and NIHF. Then, gene ontology and pathway enrichment analysis and protein–protein interaction (PPI) network analysis were performed to help us understand the independent and shared molecular mechanisms of DM and HF. In conclusion, a total of 185 DEGs and 4 hub genes, which might play an important role in the common biological mechanisms between NIHF and IHF.

Materials and Methods

Data Source

GEO (https://www.ncbi.nlm.nih.gov/geo/) is a public gene expression database created by NCBI, which contains high throughput sequencing and microarray data sets.⁹ Among the inclusion criteria were (1) a. NIHF databases: HF patients was diagnosed dilated cardiomyopathy or no obvious evidence of ischemic heart disease; b. IHF database: HF patient was diagnosed with ischemic cardiomyopathy; (2) detection of gene level in left ventricular free wall heart tissue samples; (3) Datasets that included patients and healthy controls. Exclusion criteria included: Patients had participated in a clinical trial for drugs or other treatments. Finally, two datasets were selected from GEO: accession numbers GSE84796 (8 NIHF patients and 10 controls) and GSE26887 (12 IHF patients and 5 controls).

Identification of DEGs

GEO2R (https://www.ncbi.nlm.nih.gov/geo/geo2r/) is used to identify differentially expressed genes (DEGs).¹⁰ GEO2R is based on limma package that comes with the GEO databases. Differentially expressed genes (DEGs) were determined by |logFC| > 0.5 and adj.P-value < 0.05. Next, the online Venn software (http://bioinformatics.psb.ugent.be/webtools/Venn/) was applied to detect the overlap DEGs among three datasets.

Enrichment Analyses of DEGs

DAVID 6.8 (https://david.ncifcrf.gov/) was used for enrichment analyses to elucidate the biological functions of the overlapping DEGs.¹¹ DAVID is a comprehensive bioinformatics analysis tool, providing a set of functional annotation tools for researchers to analyze the biological functions of massive genes. Further evaluation of the pathway enrichment analyses of DEGs was implemented by KOBAS 3.0 (http://kobas.cbi.pku.edu.cn), which annotates the input gene set with putative pathways by mapping to genes with known annotations from 5 pathway databases (KEGG PATHWAY, PID, BioCyc, Reactome and Panther).¹² P-value < 0.05 was considered significant. The genes in modules were also analyzed in the same way.

Protein–Protein Interaction Network Construction and Module Analysis

Protein–protein interaction (PPI) network reveals the specific and unspecific interactions of proteins, and identifies the core protein genes. STRING (http://string-db.org, version 11.0) database was used to predict the PPI network of DEGs and analyze the interactions between proteins.¹³ An interaction with a combined score > 0.4 was recognized as statistical significance. The molecular interaction networks were visualized using the Cytoscape (version 3.7.0).¹⁴ Subsequently, we used the MCODE plugin to identify densely connected modules from the PPI network with the criteria of K-core = 2, degree cutoff = 2, max depth = 100, and node score cutoff = 0.2.¹⁵

Hub Genes Selection and Analyses

The hub genes were selected using the cytoHubba plugin, a Cytoscape plugin, was used to determine the hub proteins or genes in the PPI network.¹⁶ We randomly select 5 of the 12 algorithms in cytoHubba plugin, and take the intersection of the 5 algorithms results to determine hub gene. Subsequently, a network of genes and their co-expression genes was analyzed via GeneMANIA (http://www.genemania.org/), which is a convenient web portal for analyzing gene lists and predicting gene function.¹⁷ Finally, Drug-Gene Interaction database (DGIdb) 3.0 (http://www.dgidb.org/), which helps to predict drug–gene interaction networks, was adopted here to predict drugs based on the module genes.¹⁸ After the prediction of drug-gene pairs associated with the module genes, the network map was then formed by Cytoscape.

Construction of mRNA-miRNA Regulatory Network

We used Mirwalk to predict corresponding miRNA of hub gene. Mirwalk is a publicly available database that mainly focuses on miRNA-target interactions.¹⁹ In order to construct further accurate regulatory network, the screening condition was that the Predicted miRNA could be verified by experiments or other databases. After the prediction of mRNA-miRNA by Mirwalk, the regulatory network were visualized by Cytoscape.

Results

Identification of DEGs

After standardization and identification of the microarray results, DEGs were selected (Figure 1). There were 2780 DEGs in the GSE84796 dataset, including 1597 upregulated genes and 1183 downregulated genes (Figure 2A) (Supplementary Table 1). There were 1331 DEGs in the GSE26887 dataset, including 705 upregulated and 626 downregulated genes (Figure 2B) (Supplementary Table 2). A Venn diagram was generated to show the overlap between GSE26887 and GSE84796 datasets; these include 93 upregulated and 92 downregulated genes (DEGs1) (Figure 2C and D). Furthermore, 1146 DEGs (DEGs2) and 2595 DEGs (DEGs3) were identified independently from the DEGs in GSE26687 and GSE84796, respectively.

Figure 1 Flow diagram of the study design.

Figure 2 Identification of gene expression profiles in the two datasets. (A) Volcano plot of NIHF microarray data. (B) Volcano plot of IHF microarray data. (C) Venn diagram of the 93 communal upregulated DEGs between NIHF and IHF. (D) Venn diagram of the 92 communal downregulated DEGs between NIHF and IHF.

GO Enrichment Analyses of Independent DEGs in IHF and NIHF

To determine the biological functions of DEGs2 and DEGs3, GO enrichment analysis was performed (Table 1). Results were divided into three functional categories, including biological processes (BP), cell component (CC), and molecular function (MF). GO analysis of DEGs2 indicated that changes in BP were significantly enriched in the positive regulation of pri-miRNA transcription from RNA polymerase II promoter (GO:1902895), amino acid transport (GO:0006865), regulation of cell shape (GO:0008360) and leukocyte migration (GO:0050900). As for CC, DEGs2 were particularly enriched in extracellular exosome (GO:0070062), extracellular matrix (GO:0031012), focal adhesion (GO:0005925) and cytosol (GO:0005829). Changes in MF were mostly enriched in low-density lipoprotein particle binding (GO:0030169), amino acid transmembrane transporter activity (GO:0015171), RAGE receptor binding (GO:0050786) and glycoprotein binding (GO:0001948). The GO analysis of DEGs3 returned that the terms response to immune response (GO:0006955), adaptive immune response (GO:0002250), regulation of immune response (GO:0050776) and innate immune response (GO:0045087) under BP were mainly enriched. The terms enriched under CC were external side of plasma membrane (GO:0009897), T cell receptor complex (GO:0042101), integral component of plasma membrane (GO:0005887) and plasma membrane (GO:0005886). Moreover, the terms enriched under MF were receptor activity (GO:0004872), non-membrane spanning protein tyrosine kinase activity (GO:0004715), protein binding (GO:0005515) and receptor binding (GO:0005102).

Table 1 The GO Enrichment Analysis of DEGs2 and DEGs3 (Top 4 Terms According to p.adjust)

Protein–Protein Interaction Network Construction and Module Analysis

The PPI network of DEGs with combined scores greater than 0.4 was generated by Cytoscape, which contained 109 nodes and 119 edges (Figure 3) (Supplementary Table 3). The MCODE plugin identified three densely connected modules in which 12 DEGs were among DEGs1 (Figure 4A). Then, we use online database DAVID 6.8 and KOBAS 3.0 to analyze the enrichment analyses of 12 genes. GO analysis of 12 genes were completed by DAVID 6.8, genes were significantly enriched in ribosomal large subunit biogenesis (GO:0042273) and G-protein coupled purinergic nucleotide receptor signaling pathway (GO:0035589) in BP. As for CC, genes were particularly enriched in nucleolus (GO:0005730) and integral component of plasma membrane (GO:0005887). Changes in MF were mostly enriched in poly(A) RNA binding and G-protein coupled purinergic nucleotide receptor activity (GO:0045028) (Figure 4B and C). Pathway analysis which was completed by KOBAS 3.0 revealed genes to be mainly involved in Class A/1 (Rhodopsin-like receptors) (R-HSA-373076), G alpha (i) signalling events (R-HSA-418594), P2Y receptors (R-HSA-417957) and GPCR ligand binding (R-HSA-500792) (Figure 4D).

Figure 3 Based on database STRING and Cytoscape software, PPI networks of the DEGs were constructed. The pink point represents upregulated genes, and purple point represents downregulated genes.

Figure 4 (A) Top modules from the protein–protein interaction network. (B) The biological process in functional enrichment of the DEGs in Modules was performed using the online biological tool DAVID between HF and T2DM with P-value and (C) gene count. (D) The pathway analysis of the DEGs in Modules by KOBAS 3.0. The abscissa represents the P-value, and the ordinate represents the terms. The size of the circle represents the number of genes involved, and the color represents the frequency of the genes involved in the term total genes.

Hub Gene Selection and Analysis

In the present study, we used cytoHubba to choose hub genes. According to the five classification methods (BottleNeck, Closeness, DMNC, EPC, Stress) in cytoHubba, the top 15 hub genes selected by these ranked methods in cytoHubba are shown in Table 2. Finally, three central genes were identified by overlapping the first 15 genes (Figure 5). CCL5, C5, TLR3 were the three central genes, and three central genes were all upregulated genes. A network of the hub genes and their co-expression genes was analyzed by GeneMANIA online platform (Figure 6). The three genes showed the complex PPI network with the Physical interactions of 67.64%, Co-expression of 13.50%, Co-localization of 6.17%, Predicted of 6.35%, Pathway of 4.35%, Genetic Interactions of 1.40%, and Shared protein domains of 0.59%. Finally, based on the DGIdb predictions of the hub genes, we obtained 14 drug–gene interaction pairs, including two hub genes (C5, TLR3) and 14 drugs (FDA-listed + antitumor drugs) (Figure 7). These results may reveal the therapeutic targets related to HF.

Table 2 The Top 15 Hub Genes Rank in cytoHubba

Figure 5 Three hub genes were identified by overlapping the first 15 genes in the five classification methods of cytoHubba.

Figure 6 Hub genes and their co-expression genes were analyzed using GeneMANIA.

Figure 7 Based on the DGIdb predictions of the module genes, we obtained 14 drug–gene interaction pairs, including two hub genes (C5, TLR3) and 14 drugs. Yellow circle indicates the differentially expressed gene and blank square indicates the drug.

mRNA-miRNA Regulation Network Construction

We used MiRwalk databases to predict the miRNAs of three hub genes and found a total of 45 miRNAs through screening condition that the predicted miRNA could be verified by experiments or other databases. The data of these hub genes and their miRNAs were integrated into a regulatory network, and visualized using Cytoscape software (Figure 8).

Figure 8 mRNA-miRNA regulation network of hub genes was constructed by MiRwalk.

Discussion

In this study, we identified 185 overlapping DEGs (DEGs1) in both NIHF and IHF, of which 93 were upregulated and 92 were downregulated genes. Independent DEGs included 612 upregulated and 534 downregulated genes in NIHF (DEGs2), and 1504 upregulated and 1091 downregulated genes in IHF (DEGs3). GO analysis revealed that DEGs2 were mainly enriched in transcription-related, while DEGs3 were mainly enriched in immune-related. Enrichment analyses of the genes in the key modules of the constructed PPI network revealed that they were mainly enriched in some GPCR pathways and biological processes. Finally, three hub genes, CCL5, C5 and TLR3 were found in the PPI network. Then, we analyzed a network of the hub genes and their co-expression genes, predicted drug–gene interactions and construct mRNA-miRNA Regulation Network.

G protein-coupled receptors (GPCRs), as we know, are a major target of therapeutic intervention in most physiological processes including vision, smell, taste in addition to neurologic, cardiovascular, endocrine and reproductive function. GPCRs play a prominent role in the regulation of cardiovascular function.²⁰ β-adrenergic receptors (Ars), as a member of a family of GPCRs, is widely confirmed closely related to the development of heart failure. Therefore, beta-blockers which can block Ars become one of the classic heart failure drugs.^21,22 At present, most studies focused on the GRK, especially the GRK2 and GRK5.²⁰ In the model of IHF, Phosphorylated GRK2 modified and desensitized AdipoR1 in failing cardiomyocytes, leading to post-MI remodeling and HF progression.²³ Meanwhile, there is also some evidence indicate that GRK2 plays an important role in the development of NIHF. Rockman et al uncovered that upregulating level of GRK2 was closely related to the onset of dilated cardiomyopathy (DCM), and the effect could be reversed by GRK2-inhibitor peptide.^21,24 Similarly, an experiment was also designed to use paroxetine to specifically inhibit GRK2 in mice and found an improvement in heart failure.²⁵ Additional, GRK5 is a regulator of fibroblast activation and cardiac fibrosis.²⁶ In conclusion, the GRCP-related field, especially CRK, remains a potential direction for various heart failure mechanisms and treatments.

CCL5 (C-C motif chemokine ligand 5) is one of several chemokine genes clustered on the q-arm of chromosome 17. Chemokines form a superfamily of secreted proteins involved in immunoregulatory and inflammatory processes.^27,28 Moreover, CCL5 functions as one of the natural ligands for the chemokine receptor 5 (CCR5), which is a chemokine receptor belonging to the GPCR superfamily.^28,29 CCL5 has been suggested to drive cell migration toward the heart tissue of patients with HF, which plays an important role in the occurrence and development of non-ischemic cardiomyopathy.^30,31 The important role of CCL5 in IHF has also been revealed in numerous experiments. Stevenson et al found that CCL5 level increased significantly in the human hearts with ischemic cardiomyopathy compared to control nonfailing hearts.³² What’s more, Montecucco et al reported that the incidence of postinfarct heart failure was significantly reduced by treatment with anti-CCL5 mAb in chronic cardiac ischaemia mouse models.²⁷ These findings show that CCL5 play an important role in the process of IHF and NIHF. So, we speculate that it could be a potential target for the common treatment of all types of HF.

Complement C5 (C5) encodes a component of the complement system, a part of the innate immune system that plays an important role in inflammation, host homeostasis, and host defense against pathogens.³³ Moreover, C5 could be cleaved into the anaphylatoxin C5a and fragment C5b by convertases. During activation and amplification, C5a are constantly released and trigger proinflammatory signaling via their corresponding GPCR.³⁴ Egerstedt et al observed that C5a and other complement levels were significantly elevated in the blood of patients with heart failure.³³ Meanwhile, Lindsey et al reported that C5 gene levels were significantly higher in the extreme ventricular dilatation group than the moderate ventricular dilation group after MI.³⁵ The strong association of C5a with HF and enrichment of complement activation pathways suggests activation of the complement cascade in subjects at risk of HF. It is exciting that recent experiment has pointed to a central role of C5a in myocardial repair and regeneration.³⁶ Therefore, C5 is likely to be a potential marker and provide new therapeutic approaches.

TLR3 is a member of the Toll-like receptor (TLR) family which plays a fundamental role in pathogen recognition and activation of innate immunity. TLR3 primarily protects the heart against viral infection, also mediates inflammatory effects that may exacerbate heart damage.^37–39 TLR3 has been widely confirmed as a potential target for the treatment of HF, Gao et al reported that germline knockout of TLR3 could attenuate HF and improve survival in mouse model of chronic myocardial infarction.⁴⁰ However, the mechanisms and therapeutic direction is still unclear. Loniewski et al reported that the levels of GRK2 protein are markedly increased by TLR2, 3, 4 and 7, which uncover potential cross-talk mechanisms between TLRs and GPCRs.⁴¹ As previous researches mentioned, CRK2-related mechanisms have been confirmed to be closely related to the development of heart failure. So, we believe that potential cross-talk mechanisms between TLRs and GPCRs may be an important HF therapeutic direction.

Finally, based on the DGIdb predictions of the hub genes, we obtained 14 drugs (FDA-listed + antitumor drugs). These results may show potential treatment options to HF, but more experiments should be carried to verify and explore these possibilities. We used MiRwalk databases to predict a total of 45 miRNAs of three hub genes, which could reveal the potential mutual regulation between genes, allowing us to better understand the relationship and the potential regulation between hub genes.

Our study has some limitations.: The identified hub nodes need to be validated in future studies and the sample size in this study was relatively small. Hence, large sample size and further mechanism experiments are still needed to confirm our conclusion.

Conclusion

In conclusion, the independent DEGs provide us new perspectives on different mechanisms between NIHF and IHF, and communal DEGs identified in our study reveal potential common mechanism of NIHF and IHF. Besides, a total of three hub genes (CCL5, C5, TLR3) have been identified, which can be used as biomarkers for HF or as therapy targets.