Identifying Pivotal microRNAs and Target Genes Associated With the Pathogenesis of Atrial Fibrillation

Background: Atrial brillation (AF) is the most common arrhythmia. However, specic molecular mechanism of AF remains unclear. Our study aimed to identify pivotal target genes and miRNAs in the process of AF, which help provide the basis for clinical diagnosis and the methods for early intervention. Methods: Three gene expression array datasets (GSE31821, GSE41177 and GSE79768) and a miRNA expression array dataset of AF dataset (GSE68475) were downloaded. Differential expressed genes (DEGs) were identied using the LIMMA package and differential expressed miRNAs (DEMs) were screened from GSE68475. Target genes of DEMs were predicted using the miRTarbase database, the number of the intersection between DEGs and these target genes was 26, named CDEGs. The common DEGs (CDEGs) was subject to following analysis. Results: A total of 264 DEGs and 40 DEMs were identied between the AF and control groups. Functional and pathway enrichment analyses of up-regulated DEGs and down-regulated DEGs were performed. The CDEGs were mainly enriched in PI3K-Akt signaling pathway, negative regulation of cell division and response to hypoxia. Subsequently, the protein-protein interaction (PPI) network, the microRNA ‐ transcription factor ‐ target regulatory network and drug ‐ gene network were also constructed by Cytoscape software. Conclusion: The present study revealed several novel genes and miRNAs involved in AF. We speculated that PI3K-Akt signaling pathway might participate in the pathogenesis of AF with the interaction of MYC proto-oncogene (MYC), heat shock protein 90 kDA alpha, class B, member 1 (HSP90AB1) and DNA damage-inducible transcript 4 (DDIT4), (superoxide 2) could target miR-671-5p,


Background
Atrial brillation (AF) is one of the most prevalent sustained arrhythmias and estimated to affect 34 million people worldwide and is increasing as the population ages [1,2]. However, the pathophysiological mechanism of many cases of AF remains unclear, leading to a lack of effective treatment [3]. Only a small number of patients with AF can have their heart rhythm restored by catheter ablation or heart surgery [4,5]. The higher prevalence and limited treatments of AF lead to substantial public health and economic burdens [6]. Therefore, it is very important to illustrate the molecular mechanism of AF.
In recent years, as a new interdisciplinary subject combining molecular biology and information technology, bioinformatics has been developing rapidly [10]. With the method, Wang, T. has found that the screened miRNAs and target genes, including ZBTB20, YY1, FOXO3, miR-221, miR-101 and so on, may be a target molecule for the development of AF, which may be helpful for the early diagnosis and future treatment of AF [11]. In this article, three gene expression pro le datasets and a miRNA expression array related to AF were collected from Gene Expression Omnibus (GEO) database and identi ed signi cant genes, miRNAs as well as pathways provide a new idea for the study of the occurrence, development and precision treatment of AF at the gene level.

Study design and differentially expressed genes screening
The research was performed according to the experimental work ow in Figure 1. Using robust multiarray average algorithm in R package software version 3.6.2; (http://www.R-project.org/), three datasets (GSE31821, GSE41177 and GSE79768) were analyzed with the Affymetrix platform to convert raw array data into expression value, followed by background correction, quintile normalization and probe summarization [13]. The GSE31821, GSE41177 and GSE79768 datasets were then merged into the integrated dataset via the ComBat algorithm of the Bioconductor sva package [14]. The LIMMA package was then applied to screen differential expressed genes (DEGs) [15]. Differential expressed miRNAs (DEMs) were screened from GSE68475. p-value < 0.05, and |log fold change (FC)| > 0.5 was set as a cutoff point for selecting DEGs. DEMs were selected by using cut-off values of p-value <0.05 and |log Fc|> 0.

Functional enrichment analysis
The online tool Database for Annotation, Visualization and Integrated Discovery (DAVID) (version 6.8; david.abcc.ncifcrf.gov) [16] was used to annotate the Gene Ontology (GO) enrichment analyses [17] of the DEGs. KEGG Orthology Based Annotation System (KOBAS) is a web server for annotating and identifying the KEGG enriched pathways. [18] The signi cant enrichment for GO and KEGG analyses threshold was pvalue <0.05 and count ≥2.

Protein-protein interaction (PPI) network construction
MiRTarbase database has stored more than three hundred and sixty thousand miRNA-target interactions (MTIs), which are validated experimentally. Target genes of 40 miRNAs from GSE68475 were predicted by the miRTarbase database [19]. The Venn Diagram was used to present the intersection of DEGs and target genes of miRNAs, named common differential expressed genes (CDEGs). Search Tool for the Retrieval of Interacting Genes (STRING, https://string-db.org) is a biological resource that provides a critical assessment and integration of protein-protein interactions [20]. In this study, the list of CDEGs was submitted to STRING database to detect signi cant protein interactions with con dence (combined score) >0.4. Then, PPI network was constructed and visualized by Cytoscape 3.7.2 software [21].

MiRNA-transcription factor (TF)-target regulatory network
The transcription factors (TFs) targeted the CDEGs were predicted using the Enrichr database (http://amp.pharm.mssm.edu/Enrichr/) [22]. The results with the p-value < 0.05 were screened out. After miRNA-TFs-target regulatory relationships were obtained, the miRNA-TF-target regulatory network was constructed using Cytoscape software.

Drug-gene network analysis
Drug-Gene Interaction database [23] is developed for consolidating different data sources that involve gene druggability and drug-gene interactions. Using the DGIdb database (http://www.dgidb.org/), druggene pairs were predicted and the drug-gene network was built by Cytoscape software.

Statistical analysis
A difference of p-value <0.05 was considered signi cant.

Functional enrichment analysis
Functional enrichment analysis indicated that the up-regulated DEGs were mainly involved in biological process (BP) terms such as signal transduction, immune response. In the cell component (CC) ontology, the up-regulated DEGs were signi cantly enriched in extracellular exosome, extracellular region. Molecular function (MF) analysis also showed that the up-regulated DEGs were mainly enriched in protein binding, calcium ion binding (Figure 3(a)) (Additional le 3). Additionally, the KEGG pathway of up-regulated DEGs was found to be enriched in metabolic pathways, cytokine-cytokine receptor interaction (Figure 4(a)) (Additional le 4). For down-regulation, DEGs were mainly enriched in 10 GO terms, including 4 BP terms (negative regulation of cell proliferation, etc), 5 CC terms (extracellular region, etc), and 1 MF terms (clathrin binding) (Figure 3(b)) (Additional le 5). In addition, down-regulated DEGs were signi cantly enriched in 5 KEGG pathways such as cytokine-cytokine receptor interaction, Metabolic pathways ( Figure   4(b)) (Additional le 6).
Target gene prediction A total of 2383 target genes of 40 miRNAs were predicted by the miRTarbase database. The Venn Diagram was performed to show the intersection of DEGs and target genes of miRNAs ( Figure 5(a)).

Protein-protein interaction network
To screen out the most important genes, we constructed protein-protein interaction (PPI) network by Cytoscape software. The PPI network including 10 nodes and 9 edges, in which MYC, SOD2 and TXNIP has the highest number of nodes ( Figure 5(b)). The CDEGs were mainly enriched in negative regulation of cell division, response to hypoxia and PI3K-Akt signaling pathway ( Figure 5(c)).

Drug-gene network analysis
For the CDEGs, 98 drug-gene pairs were acquired. In the drug-gene network, there were 93 drugs, 4 upregulated CDEGs (including MMP9, GLUL, EIF2S3 and HSP90AB1) and 1 downregulated CDEG (SOD2) (Figure 7). Our results revealed DIPYRIDAMOLE could interact with HSP90AB1 to be the target of treatment of AF, but the speci c mechanism remains unclear.

Discussion
As a disease with high morbidity, AF is a major public cardiac healthy concern, which causes a high incidence of thrombosis [24]. Although scientists have made prodigious progress in the treatments of AF including the reduction of symptoms, the control of rate, the prevention of thromboembolism and cardiomyopathy [25], the awareness of accurate molecular mechanisms of AF remains suboptimal.
We integrated three publicly available gene-related AF datasets with bioinformatics analysis. We identi ed 264 DEGs, of which expression of 179 was upregulated and 65 was downregulated. Besides, 40 DEGs were uncovered in miRNA-related AF microarray data and 2383 DEGs were sorted using the miRTarbase database. Furthermore, the intersection of DEGs was constructed and 26 CDEGs were uncovered.
Studies have shown that activation of the PI3K-Akt signaling pathway promotes the growth and proliferation of cells, inhibits apoptosis [26], reduces blood glucose levels [27], enhances the in ammatory response, and aggravates the vulnerability of unstable atherosclerotic plaques [28].
Jalife and colleagues discovered that the interaction between AF and atrial remodeling leads to arrhythmia exacerbation [29]. McMullen and collaborators found that inhibition of the PI3K-Akt signaling pathway could increase AF incidence [30]. Xue and coworkers showed that exogenous hydrogen sul de might be helpful for reduction of the atrial remodeling and AF caused by diabetes mellitus through activation of the PI3K/Akt/endothelial nitric oxide pathway [31]. Zhao and collaborators postulated that aliskiren might upregulate expression of the PI3K/Akt pathway to exert cardioprotective effects against rapid atrial pacing [32]. Taken together, these results suggest that regulation of the PI3K-Akt signaling pathway might participate in AF progression. Cardiac brosis occupies an important position in cardiac remodeling, which is consistent with AF [33].
Zhang and coworkers demonstrated that c-MYC expression was upregulated by lncRNA ROR and facilitated the proliferation and differentiation of cardiac broblasts [34]. Moreover, MYC could be an important molecule downstream of the PI3K/Akt signaling pathway in various tumors [35]. DDIT4 has been shown to be activated under stress situations [36]. Li and colleagues demonstrated that DDIT4 mediates methamphetamine-induced autophagy and apoptosis through the mammalian target of rapamycin (mTOR) signaling pathway in cardiomyocytes [37]. HSP90AB1 supports the appropriate folding and maintenance of stability of proteins [38]. García and colleagues veri ed that HSP90AB1-TGFβ receptor I complex is an active participant in collagen production in TGF-β-activated broblasts [39]. Based on our results, we speculate that MYC, DDIT4, and HSP90AB1 might function in AF progression through the PI3K-Akt signaling pathway.
We constructed an miRNA-TFs-target regulatory network. miRNAs and TFs function as regulators of expression of target genes [40,41]. In most cases, miRNA inhibits expression of the target gene, so we were more focused on "reverse-regulated" miRNA-mRNA pairs. With the combination of "hub" genes, four pairs were visualized. Of these, superoxide dismutase 2 (SOD2) dominated the pairs. SOD2 is a mitochondrial antioxidant enzyme. Xu and colleagues found the protein expression of SOD2 was upregulated in AF rats after treatment with the proliferator-activated receptor-γ activator pioglitazone: those data are consistent with our ndings [42]. miRNAs related to SOD2 included miRNA-671-5p, miRNA-4306, miRNA-3125 and miRNA-4298.
Numerous studies have shown that the function of miRNA-671-5p varies in different types of cancer, including inhibition of cell proliferation, transformation and promotion of apoptosis in gastric cancer, breast cancer, osteosarcoma, and esophageal squamous cell carcinoma [43][44][45][46], but the reverse functions in colon cancer [47]. miR-4306 has been shown to be related to triple-negative breast cancer [48]. Wang and collaborators found that miRNA-4298 expression in glioblastoma patients differed from that in healthy cohorts. Little research has been undertaken on miRNA-3125 until now [49].
With regard to TFs, with the combination of genes enriched in the PI3K-Akt signaling pathway, hypoxiainducible factor 1 (HIF1)a and signal transducer and activator of transcription 1 (STAT1) were found to be associated with HSP90AB1 and MYC. HIF1a expression has been shown to increase in the right atrial appendages of AF patients [50]. Tsai and coworkers demonstrated STAT1 to be activated in pigs with AF [51].
Our study had three main limitations. First, we concentrated only on reverse regulation of miRNA-mRNA pairs and ignored the more complicated mechanisms of miRNA-mRNA pairs, Second, the miRNA and mRNA we obtained were not from identical samples. Finally, we concentrated on public databases: in vitro and in vivo studies are required to validate our ndings.

Conclusion
The present study speculated several pivotal genes and miRNAs involved in AF. We conjectured that MYC, HSP90AB1 and DDIT4 might function in the pathogenesis of AF through PI3K-Akt signaling pathway.

Competing interests
None declared.

Funding
None.

Authors' contributions
Conception and design of the research: PDF, XSJ, and ZYF; acquisition of data: PDF and XSJ; analysis and interpretation of data: XSJ and ZYF; drafting the manuscript: PDF; revision of manuscript for important intellectual content: LQZ, ZTT, ZH and ZYF. All authors read and approved the nal manuscript.