Construction of miRNA-mRNA network reveals crucial miRNAs and genes in acute myocardial infarction

Kai Wang, Zhongming Li, Wenjie Ma, Yan Sun, Xianling Liu, Lijun Qian, Jian Hong,Dasheng Lu, Jing Zhang, Di Xu,✉

1Department of Cardiology, 2Department of Geriatrics, the First Affiliated Hospital of Nanjing Medical University,Nanjing, Jiangsu 210029, China;

3Department of Cardiology, the Second Affiliated Hospital of Wannan Medical College, Wuhu, Anhui 241000, China.

Abstract Acute myocardial infarction (AMI) is a severe cardiovascular disease. This study aimed to identify crucial microRNAs (miRNAs) and mRNAs in AMI by establishing a miRNA-mRNA network. The microarray datasets GSE31568, GSE148153, and GSE66360 were downloaded from the Gene Expression Omnibus (GEO) database.We identified differentially expressed miRNAs (DE-miRNAs) and mRNAs (DE-mRNAs) in AMI samples compared with normal control samples. The consistently changing miRNAs in both GSE31568 and GSE148153 datasets were selected as candidate DE-miRNAs. The interactions between the candidate DE-miRNAs and DEmRNAs were analyzed, and a miRNA-mRNA network and a protein-protein interaction network were constructed, along with functional enrichment and pathway analyses. A total of 209 DE-miRNAs in the GSE31568 dataset, 857 DE-miRNAs in the GSE148153 dataset, and 351 DE-mRNAs in the GSE66360 dataset were identified. Eighteen candidate DE-miRNAs were selected from both the GSE31568 and GSE148153 datasets. Furthermore, miR-646, miR-127-5p, miR-509-5p, miR-509-3-5p, and miR-767-5p were shown to have a higher degree in the miRNA-mRNA network. THBS-1 as well as FOS was a hub gene in the miRNA-mRNA network and the protein-protein interaction (PPI) network, respectively. CDKN1A was important in both miRNAmRNA network and PPI network. We established a miRNA-mRNA network in AMI and identified five miRNAs and three genes, which might be used as biomarkers and potential therapeutic targets for patients with AMI.

Keywords: acute myocardial infarction, miRNAs, mRNAs, miRNA-mRNA regulatory network

Introduction

Acute myocardial infarction (AMI) is a severe cardiovascular disease[1]. It is one of the most prevalent causes of morbidity and mortality worldwide[2]. Early diagnosis and treatment can significantly reduce mortality and improve the prognosis of patients with AMI[3]. Currently, cardiac troponin I and cardiac troponin T represent the 'goldstandard' for diagnosing AMI. However, false positive results often occur in patients with heart failure,chronic kidney diseases, and sepsis[4]. Therefore, it is critical to identify other novel biomarkers with high sensitivity and specificity for early diagnosis of AMI and thereby to improve clinical outcomes.

MicroRNAs (miRNAs) are reported to have pivotal roles in the regulation of AMI[5]. By using neonatal mice myocytes transfected with agomir or antagomir,Tonyet alrevealed that miR-208a alters apoptosis gene expression and promotes apoptosis in ischemic cardiomyocytes[6]. Yuet alfound that miR-133 can protect cardiomyocytes against myocardial infarction[7]. Qiaoet alshowed that miR-125b could improve AMI in rats by regulating p38/SIRT1/p53 signaling pathway[8]. Therefore, miRNAs could be one of the potential diagnostic markers and therapeutic targets in AMI. There are many studies on the expression and function of miRNAs in AMI.However, studies on miRNA-mRNA regulatory network in AMI are scarce.

In this study, the results of both miRNA and mRNA expression profiles were integrated by reanalyzing the public datasets (GSE31568, GSE148153, and GSE66360) obtained from the Gene Expression Omnibus (GEO) (https://www.ncbi.nlm.nih.gov/geoprofiles). We identified the differentially expressed miRNAs (DE-miRNAs) and mRNAs (DE-mRNAs)between the AMI and normal control groups, and constructed the miRNA-mRNA regulatory network.Based on comprehensive bioinformatics analysis, we expected to find potential diagnostic markers and therapeutic targets for AMI.

Materials and methods

Data collection

In this study, the series of matrix files of the GSE31568, GSE148153, and GSE66360 datasets were downloaded from the GEO. The miRNA dataset GSE31568 based on Platform GPL9040 included 70 normal control samples and 20 AMI samples, whereas the miRNA dataset GSE148153 performed on GPL20712 included 7 normal control samples and 5 AMI samples. Additionally, the mRNA expression profile GSE66360 based on GPL570 consists of 50 normal control samples and 49 AMI samples. Both miRNAs and mRNAs were examined from blood samples. The flow diagram of the study design is shown inFig. 1.

Fig. 1 Flow diagram of constructing the miRNA-mRNA regulatory network associated with AMI. AMI: acute myocardial infarction;GO: Gene Ontology; KEGG: Kyoto Encyclopedia of Genes and Genomes; PPI: protein-protein interaction.

Differential expression analysis

The limma package of R-studio was used to obtain the DE-miRNAs and DE-mRNAs between AMI and control samples. To remove batch effect, Bayesian methods were used. Thet-test was applied to calculate theP-values of significant gene expression differences. The DE-miRNAs and DE-mRNAs were screened by theP-values <0.05 and |log2(fold change)| >1.Meanwhile, DE-miRNAs and DE-mRNAs identified by limma were visualized using heatmaps and volcano plots. The differentially upregulated and downregulated miRNAs identified from GSE31568 and GSE148153 dataset were separately intersected by Venn diagram to obtain the candidate DE-miRNAs that were closely related to MI.

Analysis of transcription factors and DE-miRNADE-mRNA pairs

The candidate DE-miRNAs were uploaded to FunRich (Version 3.1.3), a functional enrichment and interaction network analysis tool, to perform transcription factor enrichment analysis[9]. The target mRNAs of candidate DE-miRNAs were predicted by employing three miRNA-target tools: miRDB,miRTarBase, and TargetScan. Then, the obtained miRNA-target relationship was matched with DEmRNAs obtained by microarray analysis to obtain the interaction between DE-miRNAs and DE-mRNAs(DE-miRNAs-DE-mRNAs relationship pairs).

Construction of the miRNA-mRNA regulatory network

Based on the interaction information of the DEmiRNAs-DE-mRNAs, the miRNA-mRNA network was constructed. Then, the miRNA-mRNA network was visualized by using the Cytoscape 3.7.2 software[10]. The node degrees, betweenness centrality,and closeness centrality of the regulatory network were calculated using the cytoHubba plugin of Cytoscape 3.7.2[11]. The significant miRNAs and genes were selected according to the node degree ranking, betweenness centrality ranking and closeness centrality ranking.

Functional enrichment analysis of genes in the miRNA-mRNA network

In order to investigate the functions of the mRNAs in the network, the Gene Ontology (GO) and the Kyoto Encyclopedia of Genes and Genomes (KEGG)pathway analyses were performed in R-studio. The GO covers three domains: biological process (BP),cellular component (CC), and molecular function(MF). AP-value <0.05 was considered statistically significant.

Construction of protein-protein interaction network

The protein-protein interaction (PPI) network of genes in the miRNA-mRNA network was established by using the STRING database. An interaction score≥0.4 of PPI pairs was considered statistically significant. Subsequently, Cytoscape software was used to visualize the PPI network. After that, hub genes, considered as the key genes, were obtained by using cytoHubba from Cytoscape.

Statistical analysis

All statistical analyses were performed using R software[12]. Differences between groups were compared using the Student'st-test. TheP-values were corrected using the Benjamini & Hochberg method. A value ofP<0.05 or an adjustedP-value<0.05 was considered statistically significant.

Results

Identification of DE-miRNAs and DE-mRNAs

After the data preprocessing, 209 DE-miRNAs (94 upregulated and 115 downregulated) in the GSE31568 dataset (Fig. 2A) and 857 DE-miRNAs (191 upregulated and 666 downregulated) in the GSE148153 dataset (Fig. 2B), as well as 351 DEmRNAs (289 upregulated and 62 downregulated) in the GSE66360 dataset (Fig. 2CandD) were identified, consequently. By using the Venn-diagram,18 miRNAs were obtained in both GSE31568 dataset and GSE148153 dataset (Fig. 2EandF). All intersecting DE-miRNAs are shown inTable 1.

Transcription factor enrichment analysis

Transcription factor enrichment analysis for the candidate DE-miRNAs was performed by FunRich software, which includes three different background databases namely FunRich, UniProt, and Custom[13].TheP-value <0.05 was considered statistically significant. According to the result of FunRich software, there were a total of 15 253 genes in background, while 592 genes were mapped into 196 transcription factors. The top 10 transcription factors which had strong closeness to the candidate DEmiRNAs includedEGR1,SP4,SP1,POU2F1,SOX1,NKX6-1,IRF1,FOXK1,STAT1, andSOAT1(Fig. 3).

The miRNA-mRNA network in AMI

Based on the three miRNA-target tools, we obtained 8468 miRNA-target gene pairs of candidate 18 miRNAs. Then, the target genes were further filtered by matching the DE-mRNAs. After matching with DE-mRNAs, we obtained 95 pairs of DEmiRNAs-DE-mRNAs, which included 1 upregulated and 12 downregulated miRNAs, as well as 65 upregulated and 3 downregulated mRNAs. For the sake of investigating the role of miRNAs in AMI, a miRNA-mRNA regulatory network, consisting of 81 nodes and 95 interactions, was constructed, which is shown inFig. 4. We used the plugin cytoHubba in Cytoscape to find hub miRNAs and genes according to parameters such as degree, closeness centrality, and between closeness,etc.The top five miRNAs in the network weremiR-646,miR-127-5p,miR-509-5p,miR-509-3-5p, andmiR-767-5p, and the top 5 genes wereCDKN1A,CREB5,KCNJ2,EPAS1, andTHBS1(Table 2).

Table 1 The candidate differentially expressed miRNAs in both GSE31568 and GSE148153 datasets

Fig. 2 Identification of DE-miRNAs and DE-mRNAs between the AMI and normal control samples. A–B: Volcano plots of DEmiRNAs in GSE31568 (A) and GSE148153 (B). Red dot indicates upregulated miRNAs and green dot indicates downregulated miRNAs. C:Volcano plot of DE-mRNAs in GSE66360. Red dot represents upregulated mRNAs and green dot represents downregulated mRNAs. D:Heatmap of DE-mRNAs in GSE66360. The horizontal axis represents sample names and the vertical axis represents gene symbol names.The left 50 samples were from the control group, and the right 49 samples were from the AMI group. Red indicates upregulated mRNAs.Green indicates downregulated mRNAs. P-values were determined using the Student's t-test. Adjusted P-value <0.05 and |log2(fold change)| >1 were considered statistically significant. E and F: Venn-diagram of upregulated (E) and downregulated (F) DE-miRNAs in AMI between GSE31568 and GSE148153. AMI: acute myocardial infarction; DE-miRNAs: differentially expressed miRNAs; DE-mRNAs: differentially expressed mRNAs.

Functional enrichment analysis

In this study, the GO analysis of the genes in the regulatory network showed that a total of 806 BP terms, 35 MF terms, and 30 CC terms were significantly enriched. The BP terms were mainly enriched in the positive regulation of leukocyte activation, negative regulation of transferase activity,regulation of smooth muscle cell proliferation, smooth muscle cell proliferation, and muscle cell proliferation(Fig. 5AandB). The significantly enriched CC terms included secretory granule membrane, tertiary granule, ficolin-1-rich granule, ficolin-1-rich granule membrane and tertiary granule membrane (Fig. 5C).The significantly enriched MF terms included 14-3-3 protein binding, protein kinase inhibitor activity,kinase inhibitor activity, SMAD binding, and cAMP response element binding (Fig. 5D). The five most significantly enriched pathways in KEGG included the TNF signaling pathway, FOXO signaling pathway,colorectal cancer, Kaposi sarcoma-associated herpesvirus infection, and p53 signaling pathway (Fig. 6).

Fig. 3 The transcription factor enrichment for candidate differentially expressed miRNAs predicted by Funrich. The vertical axis represents percentage of genes and the horizontal axis represents the top 10 transcription factors. Percentage of genes: the number of genes associated with transcription factor/the number of all genes in the dataset. P-value <0.05 was considered as statistically significant.

PPI network construction

Genes in the miRNA-mRNA network were used to construct a PPI network, which included 50 edges and 39 nodes. The top 10 hub genes ranked by degree method areFOS,LYN,CDKN1A,FPR1,TNFAIP3,CLEC4D,GJA1,KLF4,PTAFR, andSLC11A1(Fig. 7).The degree score is shown inTable 3.

Table 3 Top 10 genes in the PPI network ranked by degree method

Fig. 5 The GO analyses of genes in the miRNA-mRNA network. The GO analysis of genes involved in the miRNA-mRNA network was carried out using R clusterprofiler package. A–D: The enriched GO terms in biological process (A), the relationship between genes and enriched pathways in biological process terms (B), cellular component (C) and molecular function (D). Gene ratio represents the ratio of the number of genes enriched in one GO term to the total number of differentially expressed genes. P-value <0.05 was considered as statistically significant. GO: Gene Ontology.

Fig. 6 The KEGG pathway analyses for genes in the miRNA-mRNA network. The KEGG pathway analysis of genes involved in the miRNA-mRNA network was carried out using R clusterprofiler package. The horizontal axis represents the number of genes enriched in one KEGG pathway. P-value <0.05 was considered statistically significant. KEGG: Kyoto Encyclopedia of Genes and Genomes.

Fig. 7 PPI network of genes in the miRNA-mRNA network by STRING database. A: The PPI network. B: The top 10 hub genes in the PPI network. The higher the degree score, the higher the core level of the gene in the network. PPI: protein-protein interaction.

Table 2 Top five miRNAs and top five genes in the miRNAmRNA network

Fig. 4 The miRNA-mRNA regulatory network. Based on the interaction information of the DE-miRNAs-DE-miRNAs, the construction of miRNA-mRNA regulatory network was performed using the Cytoscape software. Red triangle represents upregulated miRNA, green arrows represent downregulated miRNAs, red circles upregulated genes, and green circles downregulated genes.

Discussion

AMI is an acute and severe cardiovascular disease with increasingly high incidence and prevalence[14].Given the high mortality rate of MI, it is urgently needed to identify effective interventions to prevent and treat MI. Multiple studies proved that miRNAs are involved in the occurrence and development of AMI[15–16]. However, the molecular mechanisms of MI remain unknown and there are still many unidentified miRNAs associated with MI. Therefore, the constructed miRNA-mRNA network, which may unveil the mechanisms of MI, could help identify potential novel biomarkers to improve the accuracy of MI diagnosis.

In this study, we carried out an integrated analysis to identify changes of miRNA and mRNA expression in AMI. A total of 209 DE-miRNAs (94 upregulated and 115 downregulated) in the GSE31568 dataset and 857 DE-miRNAs (191 upregulated and 666 downregulated) in the GSE148153 dataset, as well as 351 DEMs (289 upregulated and 62 downregulated) in the GSE66360 dataset were identified. There were 18 candidate DE-miRNAs in both GSE31568 and GSE148153 datasets. In the miRNA-mRNA network,miR-646,miR-127-5p,miR-509-5p,miR-509-3-5p,andmiR-767-5pwere shown to have a higher degree.Besides, genes in the miRNA-mRNA network were mainly involved in the regulation of TNF signaling pathway and FOXO signaling pathway. In addition,THBS-1served as an important regulator in the miRNA-mRNA network andFOSwas a significant regulator in the PPI network. Also,CDKN1Awas important in both miRNA-mRNA network and PPI network.

Close functional relationships were found between transcription factors and miRNA by some previous studies. They can form regulatory networks and play an important role in the occurrence and development of diseases[17–18]. Thus, transcription factors for candidate DE-miRNAs were predicted. A study showed that EGR-1 could function as a master regulator of remote preconditioning inducing a protective effect against myocardial ischemia reperfusion injury[19]. Additionally, Wanget alrevealed that EGR-1 could be involved in coronary microembolization-induced myocardial injuryviaBim/Beclin-1 pathway-mediated autophagy inhibition and apoptosis activation. Specificity protein 1 (SP1), a ubiquitous transcription factor, plays an important role in apoptosis[20], fibrosis[21], and inflammation[22]. In 2016, Liet alrevealed that miR-7a/b could improve post-myocardial infarction remodeling by inhibiting SP1 and PARP-1 expression[23]. Similarly, previous studies showed that pathological stiffness can upregulate the transcription factor POU2F1 to promote cardiac fibroblast differentiation in mouse myocardial infarction tissue[24]. Therefore, the role of transcription factors is in accord with the importance of the candidate DE-miRNAs in the pathogenesis of AMI.

In addition, we found several miRNAs with a higher degree in the miRNA-mRNA regulatory network, includingmiR-646,miR-127-5p,miR-509-5p,miR-509-3-5p, andmiR-767-5p. Zhaoet alfound that miR-646 serves an important role in AMI[25]. In their study, they constructed circRNA-miRNA interaction networks and identified crucial cirRNAs in AMI patients. Four upregulated circRNAs were predicted to four miRNAs includingmiR-491-3p,miR-646,miR-603, andmiR-922at the same time.Similarly,miR-646was considered as a crucial miRNA in our study. Also, a study showed that the overexpression of miR-646 could downregulate EGFR/Akt pathway and inhibit lung cancer cell proliferation and metastasis[26]. Moreover, EGFR pathway is closely involved in promoting cardiac fibrosis after MI[27]. ThemiR-646network might play its role in myocardial infarctionviaregulating endothelial function. Previous studies showed thatmiR-127-5pwas considered to be potential genetic diagnosis biomarkers for MI[28]. Moreover, TGFβ1/Smad3 signaling pathway is related to left ventricular remodeling after myocardial ischemiareperfusion[29]. Xuet alproved thatmiR-127could aggravate myocardial failure by promoting the TGFβ1/Smad3 signaling pathway[30]. ThemiR-127might be significant in myocardial infarctionviathe TGFβ1/Smad3 signaling pathway. Karolinaet alfound thatmiR-509-5pcould be a potential contributor of dyslipidemia in metabolic syndrome, which is one of the most dangerous initiators responsible for MI[31–32].The miR-509-5p may serve as a key function in the progression of MI. Previous studies showed that miR-509-3-5p could cause aberrant mitosis and antiproliferative effect by suppressing PLK1 in human lung cancer A549 cells. Although MI can result in massive cardiomyocyte loss, several lines of evidence suggest that the heart could be capable of efficient regeneration[33]. PLK1 is an essential component of cardiomyocyte proliferation during heart regeneration[33]. Thus, miR-509-3-5p might be a key regulator for cardiomyocyte proliferation after AMI by targeting PLK1, and miRNA-767-5p is a novel oncogenesis-associated miRNA[34]. However, there is no evidence regarding the role of miR-767-5p in AMI.Nevertheless,miR-646,miR-127-5p,miR-509-5p,miR-509-3-5p, andmiR-767-5pmay be effective biomarkers and potential therapeutic targets for AMI.

CDKN1A, a cell cycle regulator, was a hub in the miRNA-mRNA network and PPI network.CDKN1Ahas been linked to human fibroblast proliferation[35].MI may lead to structural remodeling with fibroblast activation and differentiation[36]. Therefore,CDKN1Amay be related to the progression of MI. However,experiment verification should be performed in future.Thrombospondin-1 (THBS-1) andFOSwere key genes in miRNA-mRNA network and PPI network,respectively.THBS-1, a potent inhibitor of angiogenesis and activator of transforming growth factor,could prevent the expansion of healing myocardial infarcts.FOS, a proto-oncogene, is related to tissue development and cellular stress, as well as malignant transformation and progression in numerous tumors[37]. Isoyamaet alrevealed thatc-FOSexpression was diminished in aged rat hearts compared with young adult hearts against ischemia injury[38],implying its potential role in AMI. Thus,CDKN1A,THBS-1, andFOSmay all be necessary for AMI.

KEGG pathway analysis showed that TNF signaling pathway and FOXO signaling pathway were the two most significant pathways. It is well documented that the activation of TNF signaling pathway is involved in AMI. A study performed by Suganoet alshowed that the inhibition of TNF-α reduced the infarct size in AMI following ischemia and reperfusion[39]. In contrast, administration of lowdose TNF prior to ischemia reperfusion could reduce the infarct size[40]. Likewise, FOXOs have special functions in ischemia-associated myocardial injury.Hsuet alshowed that FOXO1 has a protective effect against cardiac ischemia through its anti-oxidative property[41].FOXO3knockout mice exhibit enlarged myocardial infarct sizes upon ischemia/reperfusion[42].Therefore, the intervention of TNF signaling pathway or FOXO signaling pathway may represent novel therapeutic strategies for AMI.

In conclusion, we constructed the miRNA-mRNA network in AMI and identified several miRNAs (miR-646,miR-127-5p,miR-509-5p,miR-509-3-5p, andmiR-767-5p), genes (CDKN1A,THBS1, andFOS),and pathways (TNF signaling pathway and FOXO signaling pathway), which may play a significant role in the occurrence and development of AMI.Therefore, these miRNAs and genes could be used as candidate biomarkers and potential therapeutic targets for patients with AMI.

Limitations

However, there were still some limitations. These findings were identified with relatively small sample size, which may have restricted the generalizability of the results. Further validation is needed with a large sample study when possible. Besides, experimental validation of these findings was not performed due to resource limitation. Molecular biology experiments are expected to validate the findings of this study in the future.

CLINICAL TRIAL REGISTRATION

TheJournalrequires investigators to register their clinical trials in a public trials registry for publication of reports of clinical trials in theJournal. Information on requirements and acceptable registries is available at https: //clinicaltrials.gov/.

Acknowledgments

This study was supported by the funds from the National Natural Science Foundation of China (Grant No. 81871359 and No. 81800445). The authors acknowledge the Gene Expression Omnibus (GEO)database for providing data of AMI available.