Comprehensive Analysis of ceRNA Regulation Network Involved in the Development of Coronary Artery Disease

Background Coronary artery disease (CAD) is one of the most common causes of sudden death with high morbidity in recent years. This paper is aimed at exploring the early peripheral blood biomarkers of sudden death and providing a new perspective for clinical diagnosis and forensic pathology identification by integrated bioinformatics analysis. Methods Two microarray expression profiling datasets (GSE113079 and GSE31568) were downloaded from the Gene Expression Omnibus (GEO) database, and we identified differentially expressed lncRNAs, miRNAs, and mRNAs in CAD. Gene Ontology (GO) and KEGG pathway analyses of DEmRNAs were executed. A protein-protein interaction (PPI) network was constructed, and hub genes were identified. Finally, we constructed a competitive endogenous RNA (ceRNA) regulation network among lncRNAs, miRNAs, and mRNAs. Also, the 5 miRNAs of the ceRNA network were verified by RT-PCR. Results In total, 86 DElncRNAs, 148 DEmiRNAs, and 294 DEmRNAs were dysregulated in CAD. We received 12 GO terms and 5 pathways of DEmRNAs. 31 nodes and 78 edges were revealed in the PPI network. The top 10 genes calculated by degree method were identified as hub genes. Moreover, there were a total of 26 DElncRNAs, 5 DEmiRNAs, and 13 DEmRNAs in the ceRNA regulation network. We validated the 5 miRNAs of the ceRNA network by RT-PCR, which were consistent with the results of the microarray. Conclusions In this paper, a CAD-specific ceRNA network was successfully constructed, contributing to the understanding of the relationship among lncRNAs, miRNAs, and mRNAs. We identified potential peripheral blood biomarkers in CAD and provided novel insights into the development and progress of CAD.


Introduction
Worldwide, coronary artery disease (CAD) is one of the most common causes of sudden death with high morbidity in recent years, severely depriving patients' quality of life and posing a heavy social-economic burden [1]. The earlier the intervention in patients with CAD is, the lesser the occurrence of sudden death. Despite the advances made in the treatment of CAD, its early diagnosis and pathogenesis remain difficult and complex, respectively. The etiology and pathogenesis of CAD are currently not well described. Therefore, the main objective of this study is to explore the mechanisms of CAD progression and seek potential new therapeutic targets for CAD treatment.
High-throughput technology has recently gained traction, providing perfect opportunities for identifying bio-markers [2]. Biomarkers and therapeutic targets of noncoding RNAs play a critical role in CAD [3]. The synergistic action of multiple genes or RNAs may result in complex diseases. The pathogenesis of CAD has received more and more attention, including epigenetic modification and ceRNA [4][5][6]. The mRNAs and proteins they encoded are expressed abnormally in this complex regulatory process, and the involved noncoding RNAs exhibit spatiotemporal specificity. Salmena et al. [7] presented the ceRNA hypothesis as a new pattern of regulatory mechanism for noncoding RNA and mRNA. Besides playing a critical role in the pathological process of various diseases, lncRNAs act as miRNA response elements in ceRNA [8]. There are relatively many studies on the function of miRNA in CAD. Recent research demonstrated that downregulated miRNA-26a-5p induces endothelial cell apoptosis in CAD [9]. In addition, the increased expression levels of miR-24, miR-33a, miR-103a, and miR-122 may be associated with the risk of CAD [10]. Recently, research considers lncRNAs to be noncoding RNAs involved in gene and protein expression regulation at the epigenetic level [11]. lncRNA regulates mRNA by acting as an endogenous sponge. The regulatory role of lncRNA in miRNA has attracted significant attention [12]. Nonetheless, very little information about the ceRNA network in CAD is available. A disorder of lncRNA expression, which breaks the ceRNA interaction of lncRNA/mRNA mediated by miRNA, is important in the ceRNA network. Therefore, a better understanding of the development process of the ceRNA network in CAD may contribute to a better understanding of the mechanism in CAD.
The focus of basic medical experiments is mostly cardiomyocytes or endothelial cells. Organisms respond to damages and change the blood composition during the development and progression of disease due to their interaction and dynamic nature. Peripheral blood may, therefore, be used as an effective sample for CAD analysis and as an alternative to biopsy.
In our study, two peripheral blood microarray datasets of patients with CAD and normal groups were, respectively,  downloaded from GEO, including the information on the expression of lncRNA, miRNA, and mRNA. In addition, we obtained the PPI network and HUB genes. Finally, the ceRNA network of CAD was constructed through comprehensive and integration analysis. This study will help to explore potential molecular targets for new intervention strategies and provide new insights into clinical diagnosis and treatment of CAD.

Data Preprocessing and Differential Expression Analysis.
In addition, we downloaded the GPL20115 and GPL9040 annotation files from GEO. We transformed the probes into the respective gene symbols based on the annotation information. When several samples corresponded to the same gene, we selected the maximum expression. GEO2R, an online analysis tool in GEO, was used to evaluate the differential expression in CAD and control groups from the two profile datasets. We defined the threshold for significant difference as | log 2 fold − change| > 1:2 and the adjusted P value < 0.01. The volcano plot generated in the GraphPad Prism 8.0 software was used to show all the significant differences in mRNAs, lncRNAs, and miRNAs.

Gene Ontology and KEGG Pathway Enrichment Analyses.
The enrichment analyses of differentially expressed mRNAs, including GO and KEGG pathway analyses, were performed using DAVID 6.8 (https://david.ncifcrf.gov/) [13][14][15]. Further structures in GO analysis including biological process (BP), cellular component (CC), and molecular function (MF) were performed to annotate protein functions. Statistical significance was set at P value < 0.05.

Construction of the Protein-Protein Interaction (PPI)
Network and Hub Gene Identification in CAD. To construct 3 BioMed Research International the PPI network, the differentially expressed mRNAs were mapped into the STRING online database (https://string-db .org/) which is used to predict the protein-protein interactions (PPI) and describe the protein functional associations [16]. The interaction pair score of 0.4 was set as the threshold value. The Cytoscape software (version 3.7.1) was used to visualize and construct the PPI network intuitively [17].
In addition, the CytoHubba plugin in Cytoscape used to identify the hubba nodes was selected to measure the core genes in the PPI network [18]. The top 10 nodes were ranked by degree.

Reverse Transcription and Real-Time Quantification of miRNA Expression.
A total of 20 patients were recruited from the Second Affiliated Hospital of Dalian Medical University, China. Ten subjects who had no cardiovascular disease were recruited as the control group. Ten subjects who had been diagnosed with coronary artery disease by DSA were recruited as the CAD group. Peripheral blood was collected by venipuncture in the morning. Peripheral blood mononuclear cells (PBMC) were isolated by the Ficoll reagent (Solarbio, China). Total RNA was extracted from PBMC using the TRIzol Reagent (Thermo Fisher Scientific, USA) according to the procedure supplied by the manufacturer. Reverse transcriptions were performed using EasyScript® One-Step gDNA Removal and cDNA Synthesis SuperMix (TransGen Biotech Co., Ltd., China). Real-time PCR was performed under the iCycler™ Real Time System by mixing it with the SYBR Premix EX Tag Master mixture kit (TransGen Biotech Co., Ltd., China). The reaction conditions of RT-PCR are as follows: 94°C for 30 sec and 45 cycles at 94°C for 5 sec, 60°C for 30 sec, then dissociation stage. The relative expression levels of miRNAs were calculated by the 2 −ΔΔCT method and normalized against U6 small nuclear RNA. The stem-loop RT and RT-PCR primers were designed, and the sequences of the primers can be found in Table 1. This research was approved by the local ethics committee, and each subject gave written informed consent.

Identification of Differentially Expressed mRNAs, lncRNAs, and miRNAs in CAD.
We used the online analysis tool GEO2R to identify the differentially expressed mRNAs and lncRNAs in peripheral blood mononuclear cells of CAD patients and healthy controls in GSE113079. In addition, we identified the differentially expressed miRNAs in GSE31568. In our findings, the GSE113079 dataset was comprised of 294 differentially expressed mRNAs (127 upregulated and 167 downregulated, Figure 1(a)) and 86 differentially expressed lncRNAs (64 upregulated and 22 downregulated, Figure 1(b)), whereas the GSE31568 dataset was comprised of 148 differentially expressed miRNAs (49 upregulated and 99 downregulated, Figure 1(c)). Tables 2-5, respectively, presents the top 20 upregulated and downregulated mRNAs, lncRNAs, and miRNAs.

GO and KEGG Pathway Enrichment Analyses of
Differentially Expressed mRNAs in CAD. The enrichment analysis results revealed 12 GO terms and 5 pathways of differentially expressed mRNAs, including 7 biological processes (BPs), 4 cellular components (CCs), and 1 molecular function (MF). The results indicated that the differentially expressed mRNAs in BPs were mainly enriched in the negative regulation of the apoptotic process, signal transduction, G-protein-coupled receptor signaling pathway, inflammatory response, immune response, positive regulation of transcription from the RNA polymerase II promoter, and apoptotic process. The differentially expressed mRNAs in CCs were enriched primarily in the integral component of the plasma membrane, in the integral component of the membrane, and in the extracellular space. The differentially expressed mRNAs in MF were mainly enriched in RNA polymerase II core promoter proximal region sequence-specific DNA binding. In addition, we obtained the following pathways: Salmonella infection, natural killer cell-mediated cytotoxicity, inflammatory bowel disease (IBD), TNF signaling pathway, and antigen processing and presentation. All the results are shown in the bar graph of Figures 2(a) and 2(b) and summarized in Table 6.

Construction of PPI Network and Hub Gene
Identification in CAD. A PPI network was constructed to elucidate the interactions of the differentially expressed genes. Data analysis was based on the STRING database with the interaction pair score > 0:4 and was visualized in Cytoscape (Figure 3(a)). The PPI network revealed 31 nodes and 78 edges. There were 14 upregulated and 17 downregulated genes, and the identification of hub genes in the PPI network was done using the CytoHubba plugin in Cytoscape. The top 10 genes determined by different-color degree are shown as hub genes (Figure 3(b)). These genes contained phosphoinositide-3-kinase regulatory subunit 1 (PIK3R1), prostaglandin E receptor 1 (PTGER1), C-C motif chemokine ligand 4 (CCL4), sphingosine-1-phosphate receptor 1 (S1PR1), adrenoceptor alpha 1A (ADRA1A), arginine vasopressin receptor 1B (AVPR1B), coagulation factor II thrombin receptor (F2R), bombesin receptor subtype 3 (BRS3), free   Table 7.
3.4. Construction of the lncRNA-miRNA-mRNA ceRNA Regulation Network in CAD. In the ceRNA network, there were a total of 26 lncRNAs, 5 miRNAs, and 13 mRNAs which were visualized using Cytoscape (Figure 4). First, we derived lncRNA-miRNA interactions from the miRcode database. Out of these, 44 of the 86 DElncRNAs and 26 of the 148 DEmiRNAs formed 98 lncRNA-miRNA pairs (Table 8). We only considered the interactions between miRNAs and mRNAs present in all three databases, namely, Targetscan,  miRTarBase, and miRDB databases. As a result, 5 DEmRNAs and 13 DEmiRNAs formed 13 miRNA-mRNA pairs were obtained ( Table 9). The information on the ceRNA network for lncRNA-miRNA-mRNA ceRNA are shown in Table 10.

Discussion
CAD is an international public health problem that has a significant impact on the quality of human life [1]. The most prevalent pathogenesis of sudden death is considered to be CAD, which continues to rise year after year [33]. However, the lack of biomarkers for early diagnosis of CAD poses a serious challenge in making an effective early diagnosis of CAD clinically. High-throughput biological techniques have recently been developed and are widely used in basic researches [2]. From the perspective of genetics, it can demonstrate the difference in thousands of genes in the disease development process, which may provide a biological basis for early accurate diagnosis of CAD [34,35]. We obtained the hub genes of CAD by integrated analysis and constructed a ceRNA network which provides a foundation for further exploration of CAD progression mechanisms.
In this study, we identified a total of 294 DEmRNAs, 86 DElncRNAs, and 148 DEmiRNAs. Based on the results of the enrichment analysis, we received 12 GO terms and 5 pathways of DEmRNAs. The role is primarily linked to the apoptotic process and signal transduction. Moreover, pathways are mainly enriched in inflammation and apoptosis. In this study, we successfully constructed a ceRNA network and established the underlying mechanism between lncRNA and mRNA that can provide a new understanding of CAD therapeutic targets. Among the three types of differentially expressed RNAs in the ceRNA network, there were 26 lncRNAs, 5 miRNAs, and 13 mRNAs.
In the present study, the ceRNA network had 13 mRNAs including KLF3, MYBL1, IRAK1, HNRNPD, TRAF6, MMP16, RARB, NOVA1, SLC10A3, ZNRF3, PHF20L1, RASEF, and AIF1L. Some mRNA genes have been reported to be prominent in the ceRNA network. KLF3 belongs to the KLF family, and it is associated with adipogenesis and inflammation [36,37]. In our present study, the expression of KLF3 was low in CAD. This result could be attributed to the epigenetic modification in CAD. This study also revealed low levels of MYBL1 in CAD. Previous reports implicated MYBL1 in adenoid cystic carcinoma and cutaneous adenocystic carcinoma [38,39]. In addition, a high expression level of RASEF was observed in CAD. RASEF is a member of the Rab family of GTPases which are linked to the regulation of membrane traffic. Overexpression of RASEF increases P38 phosphorylation, leading to apoptosis and proliferation inhibition [40]. The expression of AIF1L, associated with calcium ion binding and actin filament binding, was found to be high in CAD. Previous studies reported low expression of AIF1L in undamaged tissues. Besides, it has been implicated in the  PTGER1 Prostaglandin E receptor 1 9 PTGER1 is one of four receptors identified for PGE2. Hypoxia regulates PGE2 release and promotes PTGER1 expression [24] CCL4 C-C motif chemokine ligand 4 9 The inflammatory microenvironment influences cell recruitment and activation [25] S1PR1 Sphingosine-1-phosphate receptor 1 9 The loss of S1PR1 exacerbates post-MI cardiac remodeling and worsened cardiac dysfunction [ Figure 3: (a) PPI network of the differentially expressed mRNAs are constructed. The interaction between two genes is shown at the edge. Red represents upregulation, and green represents downregulation, respectively. (b) Hub genes in the PPI network were identified with the CytoHubba plugin ranked by degree method. The higher score has a dark red color, and the lower score has a light yellow color. 8 BioMed Research International inflammatory response associated with the rejection of human heart transplants [41].
In the present study, the ceRNA network had 5 miR-NAs, including miRNA-373-3p, miRNA-146b-5p, miRNA-132-5p, miRNA-497-5p, and miRNA-124-3p. Among them, the first three (miRNA-373-3p, miRNA-146b-5p, and miRNA-132-5p) were upregulated in CAD, whereas the last two (miRNA-497-5p, miRNA-124-3p) were downregulated in CAD. Previous studies suggest that miRNA-373-3p has a strong antifibrosis effect which is significant in preventing reconstruction after infarction [42]. In previous studies, miRNA-146b-5p's expression was found to be high in atrial fibrillation in previous stud-ies and associated with fibrosis [43]. According to the report, miRNA-132-5p may cause liver steatosis and hyperlipidemia [44]. We suspect that it might be involved in lipogenesis. Smoking reduces the ability of blood to carry oxygen, affects the role of lipid metabolism, promotes the formation of arterial plaque, and increases the incidence of CAD. The upregulation of miR-497-5p had a protective effect on human bronchial epithelial cells that had been damaged by cigarette smoke [45]. Inflammatory responses may contribute to atherosclerosis [46]. Studies have shown that miRNA-124-3p regulates the inflammatory response induced by ischemia-reperfusion injury of human myocardial cells [47]. In the present study, there were 26 lncRNAs in the ceRNA network, including 16 upregulated and 10 downregulated lncRNAs. Inflammation is not the only etiology of CAD and cannot explain the whole pathogenesis. However, it is an important factor in explaining a possible mechanism of CAD. Among the 26 lncRNAs, GAS5, CRNDE, HCG22, and XIST were associated with inflammation [48][49][50][51], while the rest have not been reported previously. GAS5 and XIST have been extensively researched in previous reports. The role of lncRNA GAS5 in atherogenesis to regulate the apoptosis of macrophages and endothelial cells through exosomes suggests that inhibition of lncRNA GAS5 can be an effective way to treat atherosclerosis [52]. GAS5 could be a successful biomarker for CAD according to the study [53]. The inhibition of XIST could improve myocardial I/R injury by regulating the miRNA-133a/SOCS2 axis and inhibiting autophagy [54]. Bioinformatics-based lncRNA analysis will be helpful for future experimental research.
While our current findings are of crucial clinical significance and may provide a basis for future studies on the mechanism of CAD, we must also pay attention to some limitations. First, based on the analysis of the databases, we should use qPCR for the verification of hub genes in clinical blood samples for further study. Secondly, our sample size of 181 was relatively small; hence, future studies should consider a larger sample size to verify the accuracy of our results. Third, the specific mechanism of the lncRNA-mRNA-miRNA network of CAD needs to be further studied in vivo and in vitro for verification.    -146b-5p  IRAK1, HNRNPD, TRAF6, MMP16,  RARB, NOVA1, SLC10A3,

Conclusion
In summary, we have constructed a CAD-specific ceRNA network to aid in further understanding the relationship among lncRNAs, miRNAs, and mRNAs. Furthermore, we identify potential therapeutic targets through public database integrated analysis and provided novel insights into the development and progress of CAD.

Data Availability
Two microarray expression profiling datasets (GSE113079 and GSE31568) were retrieved and downloaded from the Gene Expression Omnibus (GEO) database of the National Center for Biotechnology Information (NCBI) (http://www .ncbi.nlm.nih.gov/geo).