This article is mentioned in:

Introduction

Although the medical approaches have developed rapidly, acute myocardial infarction (AMI) remains a principal cause of mortality worldwide (1). AMI induces a complex pathological processes including the death of cardiomyocytes, inflammation, fibrosis and hypertrophic growth of cardiomyocytes, which may progress to heart failure (2). Currently, the most effective therapies for AMI are vascular reflow to the ischemic region by thrombolysis, percutaneous coronary intervention and coronary artery bypass surgery (3).

As the ischemic area remains hypoxic, the heart continues to lose cardiomyocytes, which decreases its contractile function. The condition can be modeled by the permanent ligation of coronary artery in rodent animals. When the myocardium is reperfused for the reoxygenation of ischemic myocardium, the reoxygenation can induce reperfusion injury through cell injury pathways, including free radical production, mitochondrial injury, activation of pro-apoptotic pathway and calcium dysregulation (4). Reoxygenation also induces cardiac remodeling by the transdifferentiation of cardiac fibroblasts into myofibroblasts (4–6). These side effects of reperfusion may increase the risk of mortality (7,8). Therefore, a comprehensive understanding of myocardial biology in ischemia-reperfusion (IR) is important.

Weighted gene co-expression network analysis (WGCNA) has been successfully used in systems biology studies to investigate the intrinsic organization of transcriptomes (9). This approach provides an effective method to examine the gene expression patterns, calculate the adjacency of genes, construct gene networks and assess the importance of genes within the network (10). It has been applied for construction of gene co-expression network modules, the identification of hub genes and identifying the connection of diagnostic genes and diseases (10–12). However, WGCNA has not been applied to identify the network and hub genes associated with IR myocardium from days 2 to 7, to the best of the authors' knowledge.

In the present study, the microarray data on IR remodeling were downloaded from National Center for Biotechnology Information (NCBI) (13) and differentially expressed genes (DEGs) identified by pairwise comparison. Subsequently, WGCNA was performed to investigate the expression modules and key genes for each module. The functional annotation of each module revealed the biological associations with dynamic features of myocardium from days 2 to 7.

Materials and methods

Analysis of differentially expressed genes

The microarray data of rat IR were downloaded from the NCBI Gene Expression Omnibus database (accession number: GSE4105) and normalized by MAS5 (13). All the samples were divided into 4 groups, including day 2 (2d)-sham, 2d-IR, day 7 (7d)-sham and 7d-IR. DEGs were analyzed between two subsequent groups. DEGs between any two consecutive stages were determined by the limma 3.34.9 package in R (version 3.1.1) (14). The threshold of the significance of DEG was set as fold change 31.5 and P<0.05. A total of 3,321 DEGs were obtained with the cutoff set as P<0.05.

Construction of gene co-expression networks

To study the functional organization of the rat IR heart, DEGs between any two consecutive stages were selected and weighted co-expression network analysis performed by using the WGCNA R package (15). First, missing values and outlier samples were checked by an unsupervised hierarchical clustering analysis. Then, a weighted adjacency matrix was created and the soft threshold power (β) set at 10 to analyze scale-free topology. The power β determines the co-expression similarity for the adjacency matrix. The power 10 was chosen, for which the scale-free topology fitting index R2 and mean connectivity may reach a plateau (Fig. 1). It suggested that β=10 was the minimum value for the scale-free topology criterion. Modules were identified by the following parameters: ‘power =10, minModuleSize =30, mergeCutHeight =0.25’. A total of 6 modules containing all the DEGs were identified. Genes of each module are demonstrated in the topological overlap heatmap (Fig. 2).

Figure 1. Different soft-thresholding powers (b) were analyzed for (A) the topology scale-free fitting and (B) mean connectivity as a function of the power evaluated from ischemia reperfusion myocardium data.

Figure 2. Heatmap of the gene coexpression network. The heatmap descriped adjacencies among genes in the analysis. A single gene correspond to each row and column of the heatmap. Low adjacencies are indicated by yellow colour and higher adjacencies are indicated by progressively red colour. The genes in the analysis were clustered in modules labeled by a color code (blue, yellow, brown, green and turquoise), which are shown by squares under gene dendrogram along the diagonal.

Modules were labeled using different colors and shown as a hierarchical clustering dendrogram. The weighted network was visualized using a heatmap plot. Gene expression patterns of DEGs in each module were analyzed by hierarchical clustering analysis.

Gene Ontology (GO) and pathway analysis

To identify the biological process and signaling pathways in each module, GO (16) and Kyoto Encyclopedia of Genes and Genomes (KEGG) (17) pathways were enriched by David online tools (david.ncifcrf.gov) (18,19). The threshold of significant enrichment was set as P<0.05.

Protein-protein interaction (PPI) network analysis

To determine the interaction between DEGs in each module, they were mapped to Search Tool for the Retrieval of Interacting Genes (STRING; version 10.5), an online tool for evaluating protein-protein interaction (20). Interactions with a combined score >0.4 were considered significant. The protein-protein networks were constructed by Cytoscape software (version 3.2.0; http://www.cytoscape.org/) (21). Modules of protein-protein interaction were identified by Molecular Complex Detection (MCODE; version 1.5) (22) with the following cutoff: Nodes of cluster >4, Score of MCODE >3.

Results

Functional annotation of DEGs in each module

Significantly enriched biological process and KEGG pathways in each module are shown in Fig. 3. Different modules represented distinct gene expression patterns following IR. The blue module denoted the DEGs enriched between 7d-IR and other groups (Fig. 3A). ‘Extracellular matrix organization’, ‘wound healing’, ‘angiogenesis’, ‘canonical Wnt signaling’ and ‘PI3K-AKT signaling’ were significantly enriched in this module (Figs. 4A and 5A). In the brown module, DEGs were enriched between IR and sham group (Fig. 3B). ‘Regulation of cell proliferation’, ‘regulation of cell shape’, ‘intergrin mediated signaling pathway’, ‘endocytosis’, ‘Fc gamma R-mediated phagocytosis’, ‘cardiac contraction’, ‘hypertrophic cardiomyopathy’, ‘adrenergic signaling’ and ‘sphingolipid signaling pathway’ were significantly enriched in this module (Figs. 4C and 5C). In the green module, DEGs were enriched between 2d-IR and other groups (Fig. 3C). ‘Response to organic cyclic compound’, ‘epithelial cell proliferation’, ‘response to wounding’, and ‘amino acid metabolism’ were enriched significantly in this module (Figs. 4D and 5D). In the red module, DEGs were enriched between days 2 and 7 (Fig. 3D). ‘Chromatin remodeling’, ‘neural tube closure’, and ‘neural tube formation’ were enriched significantly in this module (Fig. 4F). In the turquoise module, the dynamic pattern of DEGs is shown in Fig. 3E. ‘Fatty acid β-oxidation’, ‘tricarboxylic acid cycle’, ‘ATP metabolic process’ and ‘oxidative phosphorylation’ were enriched significantly in this module (Figs. 4E and 5E). In the yellow module, DEGs were enriched between sham-7 day and other groups (Fig. 3F). ‘Ion transmembrane transport’, ‘regulation of body fluid levels’, ‘adrenergic signaling’ and ‘cGMP-PKG signaling’ pathways were significantly enriched (Figs. 4B and 5B).

Figure 3. Hierarchical clustering analysis of all differentially expressed genes (DEGs) in each module. Blue color represents downregulation of genes. Red indicates upregulation of genes. Different gene expression patterns following IR were shown as different color modules. (A) The blue module enriched the DEGs between 7d-IR and other groups. (B) The brown module enriched DEGs between IR and sham group. (C) The green module enriched DEGs between 2d-IR and other groups. (D) The red module enriched DEGs between days 2 and 7. (E) The turquoise module enriched the dynamic expression of DEGs. (F) The yellow module enriched DEGs between sham-7 day and other groups. IR, ischemia reperfusion.

Figure 4. Gene Ontology enrichment analysis of all differentially expressed genes (DEGs) in each module. DEGs in each module were subjected into biological process enrichment. Top 10 significantly enriched biological process were shown here (P<0.05) for (A) blue, (B) yellow, (C) brown, (D) green, (E) turquoise, and (F) red modules.

Figure 5. Pathway enrichment of all DEGs in each module. DEGs in each module were performed KEGG enrichment. Significantly enriched KEGGs (P<0.05) are shown here for (A) blue, (B) yellow, (C) brown, (D) green and (E) turquoise modules. DEGs, differentially expressed genes; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Identification of the key genes in each module

To elucidate the interaction between the hub genes in each module, PPI networks of hub genes were constructed using Cytoscape software. To explore core modules of protein-protein interactions, PPI modules were identified by MCODE.

The protein interacting with maximum number of DEGs was defined as the key gene in each module. As shown in Fig. 6, the key genes of each module are SRC proto-oncogene, non-receptor tyrosine kinase (Scr) for the blue module, discs large MAGUK scaffold protein 1 (Dlg1) for the yellow module, ATP citrate lyase (Acly) for the brown module, RAN, member RAS oncogene family (Ran) for the green module, tumor protein p53 (TP53) for the turquoise module and polo like kinase 2 (Plk2) for the red module. The most significant cluster of each network was highlighted for yellow color in each module.

Figure 6. Protein-protein interaction network of the differentially expressed genes in each module. The nodes indicate proteins. The edges indicate interactions between proteins. Highlighted yellow nodes represent molecular complexes identified by molecular complex detection algorithm of cytoscape software. Core protein-protein interaction networks were shown here for (A) blue, (B) yellow, (C) brown, (D) green, (E) turquoise, and (F) red modules.

Discussion

The aim of the present study was to identify the genes and networks associated with the heart remodeling process following IR. Distinct modules, biological process, KEGG pathways and co-expression networks during heart remodeling were identified and the results of the present study provided insights into the regulation of heart remodeling following IR.

A total of 6 modules were identified by constructing co-expression networks. Hierarchical clustering analysis revealed that each module represented a different dynamic gene expression pattern and a different biological function. The blue module represented the specific physiological features of the heart following 7 days' IR, associated with ‘canonical Wnt signaling’, ‘TGFβ receptor signaling’ and ‘PI3K-AKT signaling pathway’. It has been reported that canonical Wnt signaling, TGFβ and PI3K-AKT signaling pathways are involved in regulation of cardiac fibroblasts growth and extracellular matrix (ECM) remodeling (23–26).

The brown module represented the common physiological features of the heart following IR, which are associated with phagocytosis. The green module represented the specific physiology of the heart after 2 days' IR, associated with amino acid metabolism. Apoptosis-associated genes were significantly upregulated in the day 2 IR heart, including heme oxygenase 1, FOS like 1, AP-1 transcription factor subunit and eukaryotic translation elongation factor 1 epsilon 1. The red module represented the physiological state of the heart between days 2 and 7, which is associated with neural development. The turquoise module represented the gradual process of heart remodeling which is closely associated with fatty acid metabolism. It has been documented that during heart remodeling and heart failure, the energy metabolism shifts from the fatty acid metabolism to glucose metabolism (27). In this module, fatty acid metabolism genes including acyl-CoA oxidase 3, pristanoyl, carnitine palmitoyltransferase 1B, acyl-CoA synthetase long chain family member 6, hydroxyacyl-CoA dehydrogenase and acyl-CoA dehydrogenase very long chain were gradually and significantly downregulated. The yellow module represented the specific physiological state of the normal heart at day 7, which is closely associated with the cGMP-PKG signaling pathway. The cGMP/PKG pathway is known to inhibit cardiac structural remodeling (28). All these results suggested heart remodeling following IR is a complex process, involving ECM, neural development, apoptosis and energy metabolism shift.

Hub genes were also identified in each module, including Src for the blue module, Dlg1 for the yellow module, Acly for the brown module, Ran for the green module, TP53 for turquoise module and Plk2 for red module. Src, upregulated in the heart 7 days following IR, is a non-receptor tyrosine kinase that serves roles in numerous biological process including cell adhesion, cell cycle and cell migration (29,30). In the heart, Src is activated by AngII and contributes to pathophysiology of cardiac remodeling and hypertrophy (23). It has been suggested that Src may be a candidate target for heart remodeling diseases (31). Dlg1, upregulated in the day 7 normal heart, encodes a multi-domain scaffolding protein, taking part in skeletogenesis, cytoskeleton organization and endothelia proliferation (32,33). Although there have not been any studies on its function in the heart, to the best of the authors' knowledge, it may be a target to improve heart function based on its function in other cell types. Acly, upregulated in the IR heart, is an ATP citrate lyase which is responsible for the synthesis of cytosolic acetyl-CoA (34). Acly-mediated dimethyl α-ketoglutarate prevents autophagy of heart and reduces heart contractile performance (35). Acly may be an adaptive response in the IR heart and may be used for the inhibition of autophagy in the heart. Ran, a member of RAS superfamily, was upregulated in the day 2 IR heart. Although there have not been any studies on its function in the heart, to the best of the authors' knowledge, previous studies have suggested that Ran controls microtubule organization during the apoptotic process (36,37). Ran can perform a similar function in cardiomyocytes (38), which should to be studied in the future. TP53, gradually upregulated in the IR heart, is key for heart disease progression. A number of studies have suggested that TP53 serves an important role in heart remodeling and heart failure, and is responsible for suppression of hypoxia inducible factor 1-induced angiogenesis, apoptosis of cardiomyocytes and fatty acid metabolism (39–41). Plk2, upregulated in the day 7 heart, is a member of the serine/threonine protein kinases. Plk2 mediates antioxidant signaling and is essential for preventing p53-dependent necrotic cell death in cancer cells and so may possess therapeutic implications for heart diseases (42,43).

Taken together, the results of the current study revealed the functional modules for heart remodeling following IR using a WGCNA method and identified key genes in each module by protein-protein interaction network. The results demonstrated that heart remodeling following IR is a complex process. A number of biological processes and signaling pathways are involved. The results of the present study may provide useful targets for heart disease treatment.

Acknowledgements

Not applicable.

Funding

No funding was received.

Availability of data and materials

The analyzed datasets generated during the study are available from the corresponding author on reasonable request.

Authors' contributions

Study design: NG, NZ, LY and ZL; data analysis: JW, FL, YW and XC; manuscript preparation: NG. All authors read and approved the final manuscript.

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Glossary

Abbreviations

Abbreviations:

References