Coronary artery disease: differential expression of ceRNAs and interaction analyses.

BACKGROUND: Previous studies have demonstrated associations between cardiovascular disease and the expression of various messenger RNAs (mRNAs), microRNAs (miRNAs), and long non-coding RNAs (lncRNAs). This study aimed to investigate the differential expression of mRNAs, lncRNAs, and miRNAs between tissues from patients with coronary artery disease (CAD) and healthy controls, and to determine the interactions between these molecules in CAD.
METHODS: We investigated the differential expression of competitive endogenous RNAs (ceRNAs) between patients with CAD and healthy controls by collecting data from Gene Expression Omnibus (GEO) microarrays. We also investigated the biological function of these differentially expressed ceRNAs by performing Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses. We then created a protein-protein interaction (PPI) network to identify the hub genes. Biosystems and literature searches were also carried out to identify relevant signaling pathways and the potential function of the differentially expressed ceRNAs.
RESULTS: We identified 456 expression profiles for miRNAs, 16,325 mRNA expression profiles, and 2,869 lncRNA expression profiles. With regards to connectivity, GO and KEGG analyses (count ≥9) identified the top 11 PPI network nodes in rank order. We also identified the top 15 significant nodes for the ceRNAs identified according to degree centrality (DC) (P<0.05). Collectively, our analyses confirmed that the differential expression of certain ceRNAs, and their respective signaling pathways were associated with CAD.
CONCLUSIONS: Data arising from 11 GO and KEGG pathways, the top 15 PPI network nodes with the best connectivity rank, and the top 15 ceRNA network nodes, as determined by DC rank in CAD population, indicated that the differential expression of these ceRNAs plays a key role in the CAD. Our findings highlight new molecular mechanisms for CAD and provide new options for the development of therapeutic targets. 2021 Annals of Translational Medicine. All rights reserved.

Introduction

Coronary artery disease (CAD) is a complex phenotype that is driven by a range of genetic and environmental factors. However, current therapies for CAD focus solely on addressing the role of cholesterol and lifestyle. Despite advances in lipid-lowering therapies, clinical trials have shown that the recommended strategies for therapeutic intervention at present can leave patients with a substantial risk of cardiovascular disease (1). It is important to stratify the risk of subsequent coronary events in patients with CAD because this may be able to help guide the application of secondary preventive therapies. In a recent study, eight microRNAs (miRNAs) were identified and used to facilitate the diagnosis of the acute coronary syndrome (2). Another study showed that targeting Angptl3 messenger RNA (mRNA) was able to retard the progression of atherosclerosis and reduce the levels of atherogenic lipoproteins (3). The expressing 9p21.3-associated long non-coding RNA (lncRNA) ANRIL has also been shown to increase the risk of CAD phenotypes in non-risk vascular smooth muscle cells (4). However, the mechanisms underlying these RNA interactions have not yet been elucidated.

Previous research has also reported the existence of different types of RNA molecules that complete with other forms of RNA to bind to miRNA, thus reducing the inhibitory effect of miRNA targeting to mRNA (5). This group of competitive endogenous RNAs (ceRNAs) included various types of RNA transcripts, including lncRNAs, pseudogenes, and protein-encoded mRNAs; these forms of RNA competed for miRNA via mechanisms that are mediated by the miRNA response element (MRE) (6). More recently, researchers used bioinformatic techniques to predict regulatory networks involving ceRNAs. The effects of such ceRNAs on target genes, and the dependence of ceRNAs on miRNAs, can be verified by in vitro experiments involving proteins and RNAs, although functional verification can only be acquired by performing in vivo experiments in cells and animal models.

In the present study, we investigated the differential expression of mRNAs, lncRNAs and miRNAs between tissues from patients with CAD and healthy controls. We also investigated the interactions between these different forms of RNA and constructed a ceRNA regulatory network. Collectively, our findings will provide useful tools to investigate the molecular mechanisms underlying CAD. We present the following article in accordance with the MDAR reporting checklist (available at http://dx.doi.org/10.21037/atm-20-3487).

Methods

Data preprocessing

First, we used public databases to search for mRNA, lncRNA, and miRNA data relating to CAD, however, we were unable to find such data in public databases. Next, we focused on the relationship between lncRNA and mRNA. Finally, we selected a set of miRNA data that we could be used to accurately investigate the regulatory relationships between lncRNA, miRNA, and mRNA. The miRNA expression profile data that we used in this study (dataset serial number: GSE59421) were acquired from the NCBI (Gene Expression Omnibus, GEO, http://www.ncbi.nlm.nih.gov/geo/) Database (7); these data were normalized and then downloaded. A total of 96 blood samples (63 samples from healthy controls and 33 samples from patients with CAD) were acquired and detected using an Agilent 021827 Human miRNA MicroArray (V3) platform (miRBase release 12.0 miRNA ID version).

We also acquired mRNA and lncRNA data from the NCBI GEO database; dataset serial number GSE42148 were firstly normalized and then downloaded for analysis. A total of 24 blood samples (11 controls and 13 patients with CAD) were detected using an Agilent 028004 SurePrintG3 Human GE 8x60K Microarray platform (Feature Number version).

Because the CAD chip data were acquired from a public database, we ensured that we used the same inclusion and exclusion criteria for CAD and controls, as described in these data which were derived from dataset serial numbers GSE59421 and GSE42148 in the NCBI GEO database. The study was approved by the ethics board of Shanghai East hospital (No. 2020161), and granted a waiver of consent for the data analysis. The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

Annotation of mRNAs and lncRNAs

The sequences that matched with probes used on the Agilent-028004 were acquired from the platform annotation file and human reference genome sequences (GRCh38) were downloaded from the GENCODE database (https://www.gencodegenes.org/releases/current.html) (8). Probe sequences were then aligned onto the reference genome using seqmap software (9). We retained the uniquely aligned (uniquely mapped) probes, and also referred their relative positions on the chromosome using positive and negative strand information. The gene that matched with each probe was then identified according to the human gene annotation file (Release 25) provided by GENCODE.

The probe with “protein_coding” annotation was used as the matching probe for mRNA, while probes with “antisense”, “sense_intronic”, “lincRNA”, “sense_overlapping” or “processed_transcript” annotation were used as matching probes for lncRNA.

Finally, the probe numbers and mRNAs/lncRNAs (by gene symbol) were matched one-by-one to remove probes that did not match the gene symbol. When different probes mapped to the same gene, we used an average of the data derived from different probes as the final expression values for individual mRNAs and lncRNAs.

Screening for differentially expressed mRNAs, lncRNAs and miRNAs

We used the limma package in R software with classical Bayesian methodology (version 3.10.3, http://www.bioconductor.org/packages/2.9/bioc/html/limma.html) (10). We then performed differential analysis between CAD and control samples. When analyzing miRNAs, mRNAs, and lncRNAs, we considered both the fold difference of logFC values and the statistical significance (P values). With regards to P values, several papers have reported different thresholds for differential analyses. We compared a variety of thresholds to screen differentially expressed genes; we found that results were optimized when we selected a significance level of P<0.05, as described previously (11). In our analyses, the threshold for differential expression was to a miRNA P value <0.05 and a |logFC| >0.263 (>1.2-fold).

Functional enrichment and pathway analysis of differentially expressed mRNAs

Enrichment analysis was performed using DAVID, a powerful tool for the analysis of enrichment (version 6.8, https://david-d.ncifcrf.gov/); this allowed us to analyze up- and down-regulated genes (12), and map their association with several aspects of Gene Ontology (GO): biological process (BP), cellular component (CC), and molecular function (MF) (13). We also performed Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses (14). Significant enrichment was considered when P<0.05 and at least 2 of enrichment number (count).

Construction of a protein- protein interaction (PPI) network and node connectivity analysis

Next, we used the STRING database (version: 10.0, http://www.string-db.org/) (15) to predict whether there were any specific interactions/relationships between the proteins encoded by the differentially expressed genes. The differentially expressed mRNAs were analyzed as gene sets with the species set to ‘homo’. The PPI score was set to 0.9 (highest confidence) and we looked for protein nodes (interactions) for genes that were up- or down-regulated. Once we identified relationships in the PPI network, we used Cytoscape software (version 3.4.0, http://chianti.ucsd.edu/cytoscape-3.4.0/) (16) to create a network map. Next, node connectivity analysis was performed using non-weighted parameters and the CytoNCA plugin (Version 2.1.6, http://apps.cytoscape.org/apps/cytonca) (17). Data highlighted the important nodes within the PPI network that were involved in key protein interactions; in other words, this allowed us to identify the connectivity of hub proteins by ranking the degree centrality (DC) of individual nodes.

lncRNA and mRNA co-expression analysis

Next, we performed correlation analyses and calculated Pearson’s correlation coefficients for the differentially expressed mRNAs and lncRNAs; we did this by considering matching samples of mRNA and lncRNA. Pair with an r>0.7 (coordinate expression) and a P value <0.05 were subsequently used to screen the ceRNA network. In addition, these mRNAs were considered to be significantly correlated with lncRNAs; mRNAs were considered as potential target genes of the lncRNAs.

Identification of the target genes for miRNA and predicting the lncRNAs upstream

Next, the differentially expressed miRNAs were analyzed by the miRWalk2.0 database (http://zmf.umm.uni-heidelberg.de/apps/zmf/mirwalk2/) database (18), which featured four major components: miRWalk, miRanda, RNA22, and TargetScan. If the predicted target genes were present in all four of these key databases, then we considered the matching mRNA to be regulated by the miRNA. Following the prediction of miRNA-mRNA relationship pairs, the mRNAs were further intersected with the differential mRNAs in order to produce differentially expressed miRNA-mRNA relationship pairs.

Next, we used miRanda software (v3.3a) (19) to predict differentially expressed miRNA-lncRNA relationship pairs using a screen score (sc) ≥140 and an energy (en) ≤−20.

Pathway enrichment analysis for lncRNAs and miRNAs

Based on the identified lncRNA-mRNA co-expression relationship pairs, and the miRNA-mRNA relationship pairs, we then used mRNAs as potential target genes for matching lncRNAs and miRNAs, respectively. KEGG pathway enrichment analysis was then performed using the R package clusterProfiler (version: 3.8.1, http://bioconductor.org/packages/release/bioc/html/clusterProfiler.html) (20); the results of this analysis indirectly predicted the functional roles of the lncRNAs and miRNAs. The threshold for significance was set at P<0.05.

Construction of a ceRNA network

Using the identified mRNA-miRNA and lncRNA-miRNA relationship pairs, we then screened the miRNA-lncRNA-mRNA relationship pairs that were regulated by the same miRNA and then combined the positive co-expression relationship between mRNA and lncRNA (where the correlation coefficient >0.7). Next, we carried out further screening using the miRNA-lncRNA-mRNA relationship pairs in order to construct a ceRNA network. Consequently, the lncRNAs and mRNAs that were regulated by the same miRNAs and exhibited a positive co-expression relationship in the ceRNA network, these lncRNAs and mRNAs were all considered to be ceRNAs.

Finally, we carried out node connectivity (degree) analysis using the Cytoscape plugin CytoNCA using non-weighted parameters. The higher the connectivity, the higher the importance of a specific node within the network.

Results

Basic data relating to differential expression analyses

We identified a total of 456 miRNA expression profiles (Dataset 1), 16,325 mRNA expression profiles (Dataset 2), and 2,869 lncRNA expression profiles (see Dataset 3).

Eighteen differentially expressed miRNAs were identified according to our strict thresholds, including 16 that were down-regulated and 2 that were up-regulated (Dataset 4). A total of 92 differentially expressed lncRNAs were identified, including 46 that were down-regulated and 46 that were up-regulated (Dataset 5). A total of 610 differentially expressed mRNAs were identified, including 244 that were down-regulated and 366 that were up-regulated (Dataset 6).

Based on the identification of these differentially expressed miRNAs, lncRNAs, and mRNAs, we next created a heat map (Figure S1) and a volcano map (Figure S2).

Functional and pathway enrichment analysis for up- and down-regulated mRNAs

Next, we used the differentially expressed mRNAs in GO enrichment analyses with regards to BP, MF, and CC. This was then followed by KEGG pathway enrichment analyses. GO analyses identified enrichment in 36 BPs, 3 CCs, 15 MFs, and 2 KEGG pathways (Dataset 7). Figure S3 shows the top 10 enrichments, along with the corresponding P values. Key results arising from the GO and KEGG pathway enrichment analyses are shown in Table S1.

PPI network construction and module analysis

We identified a total of 388 protein interaction pairs; network construction was performed using Cytoscape software, as shown in . The resultant network included a total of 171 nodes.

Figure 1

Diagram showing a protein-protein interaction (PPI) network. Red indicates up-regulated proteins while green indicates down-regulated proteins. The gray line indicates protein interactions. The size of the nodes indicates the degree of connectivity.

The network was analyzed for node connectivity according to specific parameters; the top 15 nodes, as ranked by DC, are shown in . CXCL8, FPR2, IL6, and PPBP, were all ranked in the top 15 and may represent the hub proteins in our network (Dataset 8). The top 15 nodes arising from connectivity ranking within the PPI network are shown in .

Table 1

PPI network node connectivity rankings (top 15 nodes)

Node	Degree	P value	Effect on regulation	Name
GNG13	22	0.017945	Down	G protein subunit gamma 13
GNG11	22	0.04065	Down	G protein subunit gamma 11
CXCL1	18	0.000525	Up	C-X-C motif chemokine ligand 1
CXCR4	14	0.000134	Up	C-X-C motif chemokine receptor 4
CXCL8	14	0.001358	Up	C-X-C motif chemokine ligand 8
SNAP23	14	0.006229	Up	synaptosome associated protein 23
FPR3	13	0.003603	Down	formyl peptide receptor 3
RAB44	13	0.028049	Up	RAB44, member RAS oncogene family
SKP2	12	0.000943	Up	S-phase kinase associated protein 2
GPR18	12	0.001281	Down	G protein-coupled receptor 18
GALR3	12	0.001499	Up	galanin receptor 3
CXCL6	12	0.002672	Up	C-X-C motif chemokine ligand 6
ITCH	12	0.007857	Up	itchy E3 ubiquitin protein ligase
S1PR3	12	0.010344	Up	sphingosine-1-phosphate receptor 3
GRM2	12	0.010472	Up	glutamate metabotropic receptor 2

Up indicates an upregulation in CAD; down indicates a downregulation in CAD. CAD, coronary artery disease; PPI, protein-protein interaction.

Co-expression analysis of mRNAs and lncRNAs

We performed co-expression analysis of differentially expressed mRNAs and lncRNAs. According to specific thresholds, we screened a total of 1,487 significant relationship pairs that showed the coordinated expression, including 381 mRNAs and 74 lncRNAs (Dataset 9).

miRNA target genes and upstream lncRNA prediction analysis

Based on the differentially expressed miRNAs and lncRNAs, we were able to predict a total of 452 lncRNA-miRNA relationship pairs (Dataset 10), including 18 miRNAs and 72 lncRNAs.

Furthermore, based on our differentially expressed miRNAs, we were able to perform target gene prediction using mirWalk. After considering the intersection formed by differentially expressed mRNAs, we were able to identify 276 miRNA-mRNA relationship pairs, including 17 miRNAs and 170 mRNAs (Dataset 11).

KEGG pathway analysis identified a total of 27 lncRNAs and 12 miRNAs that showed enrichment (Datasets 12 and 13, respectively). These data are presented in .

Figure 2

Results arising from lncRNA and miRNA pathway enrichment analysis. (A) KEGG pathway analysis identified a total of 27 lncRNAs that were enriched. (B) KEGG pathway analysis identified a total of 12 miRNAs that were enriched. Top: lncRNA; bottom: miRNA. A reduction in the significance of the P value is indicated by a color change (from blue to red). The size of the bubbles indicates the proportion of enriched genes (the proportion of the number of involved term genes in the number of input genes). LncRNA, long non-coding RNA; miRNA, micro RNA.

ceRNA network analysis

Next, we screened the identified miRNA-lncRNA-mRNA relationship pairs that were regulated by the same miRNAs, along with the positive co-expression relationships between mRNAs and lncRNAs (with a correlation coefficient >0.7). We also screened the lncRNA-miRNA-mRNA relationship pairs by constructing a ceRNA network (). The network featured a total of 87 lncRNA-miRNA relationship pairs, 88 miRNA-mRNA relationship pairs, and 137 lncRNA-mRNA co-expression relationships (Dataset 14); in total, there were 36 lncRNAs, 64 mRNAs, and 15 miRNAs. Connectivity analysis was performed on each node of the ceRNA network to identify connectivity between mRNAs, miRNAs, and lncRNAs ().

Figure 3

A diagram showing the ceRNA network. The red circles represent up-regulated mRNAs while the green circles represent down-regulated mRNAs. Yellow triangles represent up-regulated miRNAs while gray triangles represent down-regulated miRNAs. Blue diamonds represent down-regulated lncRNAs while pink diamonds represent up-regulated lncRNAs. The blue T-type lines represent miRNA-lncRNA regulatory relationships; the yellow arrows represent miRNA-mRNA regulatory relationships, and the green dotted lines represent co-expression relationships between mRNA and lncRNA. LncRNA, long non-coding RNA; miRNA, micro RNA; mRNA, messenger RNA.

Table 2

Top 15 nodes in a ceRNA network as determined by degree centrality ranking

Node	Degree	Type	logFC	P value
hsa-miR-539-5p	30	mi_down	−0.31167	0.002388
hsa-miR-654-5p	27	mi_down	−0.38709	0.001381
hsa-miR-432-5p	22	mi_down	−0.29313	0.011733
PSMA3-AS1	20	lnc_up	0.629594	0.032669
hsa-miR-136-5p	19	mi_down	−0.34012	0.010986
MIR29A	18	lnc_up	0.849701	0.001608
RP11-420L9.4	17	lnc_down	−0.69127	0.00077
LINC01619	14	lnc_up	0.610105	0.005996
NACC2	14	m_up	1.011445	7.23E-05
hsa-miR-329-3p	13	mi_down	−0.28995	0.006895
hsa-miR-1246	12	mi_up	0.489247	0.000567
hsa-miR-133b	11	mi_down	−0.40401	0.045067
hsa-miR-431-5p	11	mi_down	−0.37006	0.001259
ST20-AS1	10	lnc_up	1.167706	0.000122
STK35	10	m_up	1.147951	0.000764

Up indicates an upregulation in CAD; down indicates a downregulation in CAD. mi, microRNA; lnc, lncRNA; m, mRNA. CAD, coronary artery disease; ceRNA, competitive endogenous RNA; lncRNA, long non-coding RNA; mRNA, messenger RNA.

Discussion

The study aimed to identify statistically significant and differentially expressed genes between patients with CAD and healthy controls. With regards to connectivity, GO and KEGG analyses (count ≥9) identified the top 11 PPI network nodes in rank order. We also identified the top 15 significant nodes for ceRNAs that were identified according to degree centrality (P<0.05). Collectively, our analyses confirmed that the differential expression of certain ceRNAs, and their respective signaling pathways, were all associated with CAD.

GO and KEGG pathway analyses identified enrichment in 36 BPs, 3 CCs, 15 MFs, and 2 KEGG pathways (Table S1). These enrichments were all statistically significant and were probably associated with the development and progression of CAD. These enrichments should be investigated further in future research as they could form targets for diagnostic or therapeutic intervention.

By ranking the nodes in our PPI network by connectivity, we were able to identify the top 15 nodes and identified statistically significant protein interactions in CAD () that were also associated with specific signaling pathways (). Some of these nodes and pathways were identified to play a known role in cardiovascular diseases. For example, GNG11 is a member of the gamma subunit family of heteromeric G-proteins. The over-expression of GNG11 has been shown to activate ERK1/2 in the MAP kinase family, but not Ras (21). These previous findings provide clinically relevant biological insight into the heritable variation of vagal heart rhythm regulation and demonstrate that genetic variants (GNG11, RGS6) play a key role by influencing the action of G-protein heterotrimers during GIRK channel- induced pacemaker membrane hyperpolarization (see ) (22).

Table 3

Differentially expressed genes in a PPI network as determined by bioinformatics analysis

Differentially expressed genes	Possible signaling pathways	BSID
GNG13	Calcium signaling; endothelin pathways; G protein signaling pathways; myometrial relaxation and contraction pathways	198906; 672453; 198849; 198759
GNG11	Calcium regulation; G protein signaling pathways; myometrial relaxation and contraction pathways	198906; 198849; 198759
CXCL1	Cytokines and inflammatory response; IL23-mediated signaling events; senescence and autophagy in cancer; spinal cord injury	198794; 138000; 198780; 739007
CXCR4	Amplification and expansion of oncogenic pathways as metastatic traits; cardiac progenitor differentiation; CXCR4-mediated signaling events; ephrin B reverse signaling; GPCRs, Class A Rhodopsin-like; hematopoietic stem cell differentiation; HIF-1-alpha transcription factor network; peptide GPCRs S1P3 pathway; syndecan-4-mediated signaling events	1458200; 712094; 137910; 138052; 198886; 1458271; 138045; 198897; 138063; 138064
CXCL8	Allograft rejection; ATF-2 transcription factor network; bladder cancer; calcineurin-regulated NFAT-dependent transcription in lymphocytes; corticotropin-releasing hormone signaling pathway; EBV LMP1 signaling; glucocorticoid receptor regulatory network; hepatitis C and hepatocellular carcinoma; IL-3 signaling pathway; LPA receptor mediated events; lung fibrosis; overview of nanoparticle effects; photodynamic therapy-induced NF-kB survival signaling; regulation of toll-like receptor signaling pathway; regulation of Wnt-mediated beta catenin signaling and target gene transcription; RIG-I-like receptor signaling; senescence and autophagy in cancer; spinal cord injury; SIDS susceptibility pathways; syndecan-2-mediated signaling events; syndecan-3-mediated signaling events; TSLP signaling pathway; TLR4 signaling and tolerance; toll-like receptor signaling pathway; validated transcriptional targets of AP1 family members Fra1 and Fra2	920963; 138006; 1458294; 137993; 920957; 198790; 138014; 1458278; 198881; 137928; 1458233;1458277; 1458260; 920955; 169352;1458257; 198780; 739007; 198901; 138031; 137952; 672451; 1458263; 198855; 169349
SNAP23	Insulin signaling	198845
FPR3	GPCRs, class A rhodopsin-like; peptide GPCRs	198886; 198897
RAB44	NA	NA
SKP2	Cell cycle; C-MYC pathway; FOXM1 transcription factor network; FoxO family signaling; notch signaling pathway	198811; 169344; 137935; 138036; 169345
GPR18	GPCRs, class A rhodopsin-like; GPCRs, other	198886; 198765
GALR3	GPCRs, class A rhodopsin-like; peptide GPCRs	198886; 198897
CXCL6	NA	NA
ITCH	Calcineurin-regulated NFAT-dependent transcription in lymphocytes; CXCR4-mediated signaling events; EGF/EGFR signaling pathway; ErbB4 signaling events; notch signaling	137993; 137910; 198782; 137991; 198879
S1PR3	S1P3 pathway; signal transduction of S1P receptor; small ligand GPCRs; sphingosine 1-phosphate (S1P) pathway	138063; 198819; 198892; 137931
GRM2	GPCRs, class C metabotropic glutamate, pheromone	198839

BSID indicates the ID number of the identified signaling pathway, as determined by bioinformatics analysis while NA indicates that there were no data available. AP1, APETALA1; ATF-2, activating transcription factor-2; CXCR4, C-X-C chemokine receptor type 4; EBV LMP1, EBV virus latent membrane protein 1; EGF/EGFR, Epidermal growth factor/Epidermal growth factor receptor; GPCRs, G-protein-coupled receptors; NFAT, nuclear factor of activated T cells; NF-kB, nuclear factor-kappa B; PPI, protein-protein interaction.

CXCL1 is known to be produced mainly by TNF-stimulated endothelial cells (ECs) and pericytes; these support luminal and sub-EC neutrophil crawling. CXCL1 and CXCL2 act in a sequential manner to guide neutrophils through venular walls as governed by their distinct cellular sources (23). The angiotensin II-induced infiltration of monocytes in the heart is largely mediated by CXCL1-CXCR2 signaling which initiates and aggravates cardiac remodeling. Consequently, the inhibition of CXCL1 and/or CXCR2 may represent new therapeutic targets for treating hypertensive heart diseases (see ) (24).

Wnt-Cxcr4 (C-X-C motif chemokine receptor 4) signaling is known to regulate the interactions between oligodendrocyte precursor cells (OPCs) and endothelial cells and coordinates OPC migration and differentiation (25). It has been shown that a number of neutrophils re-enter the vasculature and follow a pre-programmed journey that entails transportation to the lungs to up-regulate CXCR4 before entering the bone marrow, where they subsequently undergo apoptosis (26). Vascular CXCR4 limits atherosclerosis by maintaining arterial integrity, and by preserving endothelial barrier function and a normal contractile smooth muscle cell phenotype. Enhancing these beneficial functions of arterial CXCR4 by selective modulators might provide new therapeutic options for atherosclerosis () (27).

Oleic acid treatment has been shown to reduce the insulin sensitivity of heart muscle cells and that such sensitivity can be completely restored by transfection with SNAP23. Thus, SNAP23 might represent a link between insulin sensitivity and the flow of fatty acids into a cell () (28).

Co-activator-associated arginine methyltransferase 1 (CARM1) is a crucial component of autophagy in mammals. CARM1-dependent histone arginine methylation is a crucial nuclear event in autophagy and takes part in the AMPK-SKP2-CARM1 signaling axis in the regulation of autophagy induction following nutrient starvation () (29).

GPR18 is a cannabinoid-activated orphan G protein-coupled receptor (GPCR) that is selectively expressed on immune cells (30). GPR18 is known to play a role in the cardiovascular system, at least in part, by mediating the elevation of adiponectin levels () (31).

By ranking the top 15 nodes in our ceRNA network by degree centrality, we identified a range of significant ceRNA interactions in CAD (see ); these findings concurred with previous literature (see ) (32-47). It was evident from our results that several ceRNAs were related to cardiovascular diseases. For example, by regulating CDKN2A and inhibiting G1- to S-phase transition, STK35L1 may act as a central kinase that links the cell cycle and migration of endothelial cells. Specifically, the interaction between STK35L1 and nuclear actin may be critical in the regulation of these fundamental endothelial functions (45). Serine/threonine kinase 35 (STK35) is a recently identified human kinase that plays a key role in autophosphorylation, which is functionally linked with actin stress fibers, cell cycle progression, and survival (see ) (46). Microarray analysis previously showed that nuclear-retained importin α2 binds with DNase I-sensitive nuclear component(s) and controls the selective upregulation of mRNA encoding STK35. Chromatin immunoprecipitation and promoter analysis further demonstrated that importin α2 can access the promoter region of STK35 to accelerate transcription in response to hydrogen peroxide exposure. Furthermore, the constitutive overexpression of STK35 proteins enhances caspase-independent cell death under oxidative stress conditions (see ) (47).

Table 4

Differentially expressed genes in the ceRNA network that were derived from literature searches

Differentially expressed genes	Description	Reference
hsa-miR-539-5p	LOC100129148 enhanced the KLF12 expression through functioning as a competitive ‘sponge’ for miR-539-5p	(31)
hsa-miR-654-5p	hsa-miR-654-5p can directly repress the mRNA and protein expressions of BMP2 by binding to a specific target site	(32)
hsa-miR-432-5p	hsa_circ_0008039/miR-432-5p/E2F3 axis might be a potential therapeutic target for breast cancer	(33)
PSMA3-AS1	a unique role of exosomal PSMA3 and PSMA3-AS1 in transmitting proteasome inhibitors’ (PI) through a novel exosomal PSMA3-AS1/PSMA3 signaling pathway	(34)
hsa-miR-136-5p	hsa-miR-136-5p, an ERα36 3’UTR-targeting microRNA, the expression of which inversely correlated to the ERα36 one in breast cancer cells	(35)
MIR29A	the mechanism of microRNA29a (miR29a)-mediated modulation of early collagen 3 expression that facilitate tissue remodeling in the tendon after injury	(36)
RP11-420L9.4	GENE NAME: ENSG00000238158	NA
LINC01619	LINC01619 functions as a ceRNA and regulates miR-27a/FOXO1-mediated endoplasmic reticulum stress and podocyte injury in diabetic nephropathy	(37)
NACC2	NACC2 was found to associated with survival	(38)
hsa-miR-329-3p	hsa-miR-329-3p be antiviral therapeutics against MERS-CoV infection	(39)
hsa-miR-1246	miRNA-1246 are useful as biomarkers for tracking lung cancer progression and as therapeutic targets	(40)
hsa-miR-133b	hsa-miR-133b showed highly significant differential expression in Parkinson’s disease	(41)
hsa-miR-431-5p	hsa-miR-431-5p may play a role in diffuse large B-cell lymphoma by regulating tyrosine protein kinase Fyn	(42)
ST20-AS1	The immune-related lncRNA ST20-AS1 signature was identified in patients with anaplastic gliomas	(43)
STK35	STK35 may act as a central kinase linking the actin stress fibers, cell cycle progression, migration of endothelial cells, and survival. Its protein constitutive overexpression enhances caspase-independent cell death under oxidative stress conditions	(44-46)

NA indicates no data available. BMP2, bone morphogenetic protein 2; ceRNA, competitive endogenous RNA; KLF12, Krüppel like factor 12; NACC2, nucleus accumbens associated 2; PSMA3, proteasome subunit alpha type-3; STK35, serine/threonine kinase 35.

Our hypothesis was that the analysis of ceRNAs, which compete with miRNA for the MRE, might help us to distinguish between patients with CAD and controls. Our findings identified a significant number of differentially expressed genes and identified enrichment of these genes in CAD. We also identified key GO terms, pathways, and several key interactions by degree centrality ranking.

There are some limitations associated with our study that should be taken into consideration. For example, our bioinformatics analysis was based on GEO and could have been affected by variability between different subjects. Furthermore, the lncRNA molecules identified from specific tissues would acquire the reliable conclusion, thus the magnitude of connectivity and the matching lncRNA-miRNA-mRNA relationships should be combined in future experimental validations.

Conclusions

Our research identified 11 GO terms and KEGG pathways; based on these, we identified the top 15 nodes in a PPI network by connectivity ranking. We also identified the top 15 nodes in a ceRNA network by degree centrality ranking in a population of CAD subjects. Our findings are likely to facilitate the further investigation of molecular mechanisms underlying CAD.

The article’s supplementary files as