Abstract

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.

1. 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 biomarkers [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 [46]. 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.

2. Materials and Methods

2.1. Microarray Data

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), a public functional genomics data repository. The GSE113079 dataset, deposited by Li L et al., was an expression profile of mRNA-lncRNA based on the GPL20115 Agilent-067406 Human CBC lncRNA+mRNA microarray V4.0. The experiment contained 141 samples of peripheral blood mononuclear cells, including 93 subjects with coronary artery disease (CAD group) and 48 healthy controls (control group). In addition, the GSE31568 dataset, deposited by Keller A et al., was a miRNA expression profile based on the GPL9040 febit Homo Sapiens miRBase 13.0. A total of 40 samples were collected comprising 20 samples of peripheral blood mononuclear cells from patients with coronary artery disease (CAD group) and 20 samples from healthy people (control group).

2.2. 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 and the adjusted 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.

2.3. 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/) [1315]. 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 value < 0.05.

2.4. Construction of the Protein-Protein Interaction (PPI) Network and Hub Gene Identification in CAD

To construct 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.

2.5. Construction of the lncRNA-miRNA-mRNA Competitive Endogenous RNA (ceRNA) Regulation Network in CAD

We constructed a competitive endogenous RNA regulation network among lncRNAs, miRNAs, and mRNAs based on the hypothesis that miRNAs can regulate the activity of target mRNAs, and that lncRNAs can interact with the activity of miRNAs [7]. In this hypothesis, lncRNAs act as miRNA sponges. The miRcode (http://www.mircode.org) was used to predict DElncRNA-DEmiRNA interactions [19]. In addition, we predicted the possible target DEmRNAs of DEmiRNAs by using three data online tools, Targetscan (http://www.targetscan.org) [20], miRTarBase (http://mirtarbase.mbc.nctu.edu.tw) [21], and miRDB (http://www.mirdb.org) [22]. The possible DEmiRNA-DEmRNA interactions were searched in all the three databases. Finally, we established and plotted the ceRNA network from the differentially expressed lncRNAs, miRNAs, and mRNAs by using Cytoscape software (version 3.7.1).

2.6. 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 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.

3. Results

3.1. 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 25, respectively, presents the top 20 upregulated and downregulated mRNAs, lncRNAs, and miRNAs.

3.2. 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.

3.3. 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 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 fatty acid receptor 3 (FFAR3), and opsin 4 (OPN4). The details are shown in 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.

3.5. RT-PCR Verification Results

The expressions of 5 miRNAs in the ceRNA network were verified with RT-PCR. Expressions of miRNA-373-3p, miRNA-146b-5p, and miRNA-132-5p in CAD were 2.70, 1.88, and 2.26 times, respectively. Expressions of miRNA-497-5p and miRNA-124-3p in CAD was 0.57 and 0.76 times, respectively (Figure 5). The trends were consistent with the results of array data.

4. 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 inflammatory response associated with the rejection of human heart transplants [41].

In the present study, the ceRNA network had 5 miRNAs, 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 studies 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 [4851], 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.

5. 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).

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Authors’ Contributions

All authors contributed to the work shown in this paper. Jiabei He designed the calculation process, analyzed the original data, and drafted the manuscript. Yao Zhang and Qingqing Zhang completed the collection of background information and analyzed the original data. Xionglong Li prepared the figures. Lianhong Li was responsible for supervision and full correspondence. The final manuscript was reviewed and approved by all authors.