Role of competitive endogenous RNA networks in the pathogenesis of coronary artery disease.

BACKGROUND: The present study aimed to construct a network of competitive endogenous RNAs (ceRNAs) related to the pathogenesis of coronary artery disease (CAD), to provide a novel rationale for CAD treatment.
METHODS: Bioinformatics methods were applied to screen for differentially expressed long non-coding RNAs (DElncRNAs), microRNAs (DEmiRNAs), and mRNAs (DEmRNAs) from the GSE68506, GSE59421, and GSE20129 datasets of the Gene Expression Omnibus (GEO) database. The miRcode database was used to predict lncRNA-binding miRNAs. The miRTarBase, miRDB, and TargetScan databases were used to predict the target genes of these miRNAs. An mRNA-miRNA-lncRNA ceRNA network of CAD was established.
RESULTS: Between the CAD and normal control groups there were 264 DElncRNAs, 106 DEmiRNAs, and 1,879 DEmRNAs. We screened these differentially expressed gens (DEGs) respectively. There were 21 DElncRNAs, 13 DEmiRNAs, and 143 DEmRNAs in the ceRNA network by using Cytoscape application. The DEmRNAs were involved in the PI3K-Akt signaling pathway and the NF-κB signaling pathway. The key genes in the protein-protein interaction (PPI) network were HSP90AA1, CDKN1A, MCL1, MDM2, MAPK1, ABL1, LYN, CRK, CDK9, and FAS.
CONCLUSIONS: The ceRNA network constructed in this study identified new candidate molecules for the treatment of CAD, providing some more comprehensive and higher-quality choices for the target treatment of CAD. 2021 Annals of Translational Medicine. All rights reserved.

Introduction

Coronary artery disease (CAD) is a commonly occurring type of cardiovascular disease. It can result in angina pectoris, myocardial infarction, heart failure, and arrhythmia. In certain cases, the occurrence of atherosclerosis results in coronary artery stenosis and insufficient blood supply to the myocardium, leading to death. Current treatments for CAD include percutaneous coronary intervention, drug therapy, and coronary artery bypass grafting. The age and sex standardized incidence of CAD was 436 per 100,000 in 2015 (1). Although patient quality of life has improved through advancements in medical treatment and secondary prevention, it still remains that 35% of CAD patients suffer relapse (2).

Recent studies have suggested new molecules involved in the progression of CAD. For example, ANRIL is a long non-coding RNA (lncRNA) expressed at low levels in the serum of CAD patients, and high expression of ANRIL predicts poor prognosis in CAD patients (3). microRNA (miR)-128 negatively regulates the expression of IRS1, which promotes the viability and migration of rat cardiac microvascular endothelial cells and inhibits cell apoptosis (4). The identification of new molecules involved in CAD progression is essential to better understand its pathogenesis and to provide new targets for the treatment of CAD.

A recent study showed that competitive endogenous RNA (ceRNA) regulation networks play an important role in heart diseases. For example, an endogenous competitive relationship between the lncRNA MEG3 and miR-145 was identified. The overexpression of MEG3 decreased the expression of miR-145, which in turn increased the expression of the target gene PDCD4 and promoted cardiomyocyte apoptosis (5). Recently, despite a ceRNA literature report on CAD, it revealed 11 pathways and 15 key genes related to CAD, which provided options for the treatment of CAD (6). However, our study used different datasets from the previous ones. These datasets correspond to lncRNA, miRNA and mRNA chip analysis results respectively, so that we can integrate and analyze data from a wider dimension. From another new perspective, we constructed a ceRNA network to reveal the molecular mechanisms related to CAD, and combined the results reported in previous articles to provide some more comprehensive and higher-quality choices for the target treatment of CAD.

We present the following article in accordance with the MDAR reporting checklist (available at https://dx.doi.org/10.21037/atm-21-2737).

Methods

CAD data

The lncRNA, miRNA, and mRNA expression profiles were downloaded from GEO (https://www.ncbi.nlm.nih.gov/geo/). These expression profiles were from plasma samples of patients with CAD. The lncRNA microarray data were obtained from the GSE68506 (comprising three CAD patients and three normal controls). The miRNA expression data were obtained from the GSE59421 (comprising 33 CAD patients and 63 normal controls), and the mRNA expression data were obtained from the GSE20129 (comprising 48 CAD patients and 71 normal controls). The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

Identification of differentially expressed lncRNAs (DElncRNAs), miRNAs (DEmiRNAs), and mRNAs (DEmRNAs)

The Bioconductor Limma (7) package and Perl were used to identify DElncRNAs, DEmiRNAs, and DEmRNAs in the CAD patients and normal controls. DElncRNAs, DEmiRNAs, and DEmRNAs were screened by thresholds of P<0.05. After the DE analysis (), we visualized the DElncRNAs, DEmiRNAs, and DEmRNAs between CAD patients and normal controls. Clustering heat maps and volcano maps were made using the R package “pheatmap”.

Figure 1

Flowchart of ceRNA network analysis. ceRNA, competitive endogenous RNA; lncRNA, long non-coding RNA; miRNA, microRNA; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; DE, differentially expressed; PPI, protein-protein interaction.

Construction of the ceRNA network

To better comprehend the relationships between the DE mRNAs, miRNAs, and lncRNAs, the lncRNA-mediated ceRNA network of CAD was constructed as follows. First, we used the miRCode database (http://www.mircode.org/) (8) to predict relationships between the lncRNAs and miRNAs. Next, the miRTarBase (http://mirtarbase.mbc.nctu.edu.tw/), miRDB (http://www.mirdb.org/), and TargetScan (http://www.targetscan.org/) (9-11) databases were used to obtain the miRNA-targeted mRNAs. To improve the effectiveness of our results, we showed miRNA-targeted mRNA both in the miRTarBase, miRDB, and TargetScan databases to establish a lncRNA-miRNA-mRNA network. Finally, Cytoscape (http://www.cytoscape.org/) 3.8.1 (12) software was used to visualize the results.

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG)

GO is a popular bioinformatics tool used to analyze the biological functions involved in target genes (13,14). The KEGG is a large-scale molecular dataset, generated using high-throughput experimental methods, that is used to understand the biological signaling pathways involved in genes (15). In the Database for Annotation, Visualization and Integrated Discovery (DAVID; https://david.ncifcrf.gov/) (16), we used GO annotations and KEGG to analyze the biological functions and signaling mechanisms involved in DEmRNAs. P<0.05 was considered statistically significant.

Construction of protein-protein interaction (PPI) network and identification of key genes

The Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) (http://string-db.org/) (17) was used to predict the PPI network based on the gene symbols (18).The PPI network of the DEmRNAs was constructed using the STRING database, using a combination score of >0.4, and the differences were statistically significant. Next, we visualized the molecular interaction network using Cytoscape. The key genes with the highest scores were screened out using the maximum clique centrality (MCC) method in the cytoHubba plug-in of Cytoscape (19).

Statistical analysis

All statistical analyses were performed using R (v.4.0.3) software, Perl (v.5.28.1) and GraphPad Prism 9 software. The P value <0.05 was considered statistically significant.

Results

Identification of DElncRNAs, DEmiRNAs, and DEmRNAs

Based on the screening criterion of P<0.05, a total of 264 DElncRNAs (179 downregulated and 85 upregulated lncRNAs), 106 DEmiRNAs (73 downregulated and 33 upregulated miRNAs), and 1,879 DEmRNAs (1,066 downregulated and 813 upregulated mRNAs) were identified between the CAD and normal control groups. Heatmap clustering indicated that the DElncRNAs, DEmiRNAs, and DEmRNAs had clearly defined differences in expression between the two groups ().

Figure 2

Volcano maps and heat maps of differential expression of lncRNAs (A,B), miRNAs (C,D), and mRNAs (E,F). lncRNA, long non-coding RNA; miRNA, microRNA; CAD, coronary artery disease.

Biological functions and signaling mechanisms related to the DEmRNAs

Through the GO annotations, we found that the DEmRNAs in the GSE20129 were enriched in protein complex assembly, nitric oxide biosynthesis, innate immune response, cytoplasm, plasma membrane, cytosol, etc. (, ). The KEGG pathway analysis showed that the DEmRNAs in the GSE20129 were mainly involved in the PI3K-Akt signaling pathway, tuberculosis, cancer pathways, etc. (, ).

Figure 3

(A) GO terms (BP, CC, MF) and (B) KEGG analysis pathways of DEmRNAs involved in the GSE20129 dataset. GO, Gene Ontology; BP, biological processes; CC, cellular components; MF, molecular functions; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Table 1

Functional enrichment analysis of the DE RNAs from the GSE20129 dataset

GO	ID	Description	P value	Count
BP	GO:0006461	Protein complex assembly	7.20E−04	22
	GO:0045429	Nitric oxide biosynthetic process	7.27E−04	12
	GO:0045087	Innate immune response	8.87E−04	57
	GO:0042493	Response to drug	0.002149315	42
	GO:0050870	Positive regulation of T cell activation	0.002838657	7
	GO:0006919	Activation of cysteine-type endopeptidase activity	0.003999244	16
	GO:0031663	Lipopolysaccharide-mediated signaling pathway	0.004432578	9
	GO:0042110	T cell activation	0.005389073	11
CC	GO:0005737	Cytoplasm	7.06E−06	503
	GO:0005886	Plasma membrane	2.58E−05	403
	GO:0005829	Cytosol	3.84E−05	331
	GO:0005730	Nucleolus	2.91E−04	100
	GO:0005622	Intracellular	3.03E−04	145
	GO:0045121	Membrane raft	3.31E−04	33
	GO:0005759	Mitochondrial matrix	4.08E−04	46
	GO:0030136	Clathrin-coated vesicle	6.83E−04	14
MF	GO:0005515	Protein binding	2.11E−05	819
	GO:0044325	Ion channel binding	4.73E−04	22
	GO:0042803	Protein homodimerization activity	6.77E−04	88
	GO:0005524	ATP binding	9.29E−04	161
	GO:0042802	Identical protein binding	0.002192893	87
	GO:0003824	Catalytic activity	0.004631555	28
	GO:0008266	Poly(U) RNA binding	0.006302894	6
	GO:0019899	Enzyme binding	0.009620338	29

DE, differentially expressed; BP, biological processes; CC, cellular components; MF, molecular functions; GO, Gene Ontology.

Table 2

Pathway enrichment analysis of the DE RNAs from the GSE20129 dataset

ID	Description	P value	Count
hsa05200	Pathways in cancer	4.14E−05	64
hsa04151	PI3K-Akt signaling pathway	0.00157308	52
hsa05152	Tuberculosis	4.19E−04	33
hsa04380	Osteoclast differentiation	1.83E−05	30
hsa04650	Natural killer cell mediated cytotoxicity	0.00147428	24
hsa05215	Prostate cancer	6.94E−04	20
hsa04064	NF-kappa B signaling pathway	0.00159901	19
hsa05140	Leishmaniasis	3.71E−04	18

DE, differentially expressed; NF, nuclear factor.

ceRNA network

Among the DElncRNAs, DEmiRNAs, and DEmRNAs, 21 lncRNAs (18 downregulated and 3 upregulated), 13 miRNAs (13 downregulated), and 143 mRNAs (86 downregulated and 57 upregulated) were involved in the proposed ceRNA network (, ).

Figure 4

ceRNA network of lncRNA-miRNA-mRNA in CAD. Diamonds represent lncRNAs, squares represent miRNAs, and circles represent mRNAs. The red nodes are upregulated RNAs, and the purple nodes are downregulated RNAs. ceRNA, competitive endogenous RNA; lncRNA, long non-coding RNA; miRNA, microRNA; CAD, coronary artery disease.

Table 3

RNAs involved in the lncRNA-miRNA-mRNA ceRNA network

RNA	Upregulated	Downregulated
lncRNA	AC011374.1, NEAT1, KCNQ1OT1	IGF2-AS, H19, HCG27, CLDN10-AS1, SNHG14, HOTAIR, EMX2OS, AC012313.1, SNHG7, AC015987.1, CR381653.1, ZEB1-AS1, ADAMTS9-AS2, DNAJB8-AS1, HOTTIP, LINC00491, NAV2-AS5, DIO3OS
miRNA	–	hsa-mir-17, hsa-mir-93, hsa-mir-141, hsa-mir-152, hsa-mir-18a, hsa-mir-196a, hsa-mir-29a, hsa-mir-23a, hsa-mir-221, hsa-mir-32, hsa-mir-92a, hsa-mir-363, hsa-mir-425
mRNA	FAM57A, ATG16L1, ATL3, BTG3, CAMTA1, CDKN1A, EPHA4, FBXO31, FNBP1L, KATNAL1, KIAA1191, KLHL20, PFKP, PTPDC1, SCAMP5, SESN3, SNTB2, ZBTB4, ZBTB9, GTPBP10, PDE4D, TAOK1, ARID1A, NRP2, PCBD2, ARHGAP28, CALM1, PAFAH1B2, TUB, CCNT1, PABPC1, MCFD2, MYH2, TNPO1, TOP1, MYCN, RPL22, ZBTB5, ABL1, YWHAE, BCAT2, DOCK9, MOAP1, NSF, PAX9, RHPN2, RRN3, TPM3, ACACA, CDK9, CELSR2, HOXC6, SLC7A1, ZNF507, ELAVL1, GRSF1, SLC29A2	ELAVL2, OGT, HSP90AA1, NUP50, ZFYVE9, MXD1, SOS2, TGFA, LPGAT1, NUFIP2, RAB11FIP1, SOD2, ACSL4, CRK, FAM126B, IFNAR1, IRAK4, KLF10, LIMK1, MAPK1, MCL1, MIDN, MTMR3, PDLIM5, PPP1R3B, RORA, RPS6KA5, RRAGD, TMEM127, TRIP10, USP32, WAC, HM13, IFNAR2, NEDD9, CPD, BACH1, GATA6, CAMKK1, CREBZF, MAPK10, PCDHA11, PCDHAC1, CELF1, FAS, FNIP1, IL6R, ZNF257, COL4A2, CPEB3, DNMT3A, MDM2, NASP, RNASEL, SLC16A1, TMOD3, IFRD1, CYB5R4, LYN, ADAM10, BAK1, DUSP10, GRAMD1B, SERTAD3, SNX10, SSFA2, TRIM36, YIPF4, GNG4, MAP3K5, THRB, RCN3, TMEM50A, ZWINT, DNAJB12, ERGIC1, ICAM1, KCNC1, HIPK3, PDIK1L, PLXNA4, CCNG2, IRF1, LYST, TBC1D20, ZNF148

lncRNA, long non-coding RNA; miRNA, microRNA; ceRNA, competitive endogenous RNA.

GO and KEGG analysis of the ceRNA network

DEmRNAs in the ceRNA network were enriched in protein phosphorylation, drug responses, viral processes, cytoplasm, nucleus, cytosol, protein binding, ATP binding, protein serine kinase activity, etc. (, ). The KEGG pathway analysis showed that the DEmRNAs in the ceRNA network were involved in the PI3K-Akt signaling pathway, neurotrophin signaling pathway, cancer pathways, etc. (, ).

Figure 5

GO terms (A) BP, (B) CC, (C) MF and (D) KEGG analysis pathways of DEmRNAs involved in the ceRNA network. GO, Gene Ontology; BP, biological processes; CC, cellular components; MF, molecular functions; KEGG, Kyoto Encyclopedia of Genes and Genomes; DEmRNAs, differentially expressed mRNAs; ceRNA, competitive endogenous RNA.

Table 4

Functional enrichment analysis of the DE RNAs from the ceRNA network

GO	ID	Description	P value	Count
BP	GO:0006468	Protein phosphorylation	6.30E−05	14
	GO:0042493	Response to drug	5.02E−06	13
	GO:0016032	Viral process	0.00252124	9
	GO:0006974	Cellular response to DNA damage stimulus	0.001291559	8
	GO:0009636	Response to toxic substance	5.54E−05	7
	GO:0038096	Fc-gamma receptor signaling pathway	5.00E−04	7
	GO:0008630	Intrinsic apoptotic signaling pathway	5.01E−04	5
	GO:0097192	Extrinsic apoptotic signaling pathway	0.002377636	4
CC	GO:0005737	Cytoplasm	4.66E−06	67
	GO:0005634	Nucleus	2.81E−04	63
	GO:0005829	Cytosol	2.45E−06	50
	GO:0005654	Nucleoplasm	0.002316159	36
	GO:0016020	Membrane	0.001486372	31
	GO:0005925	Focal adhesion	0.01164298	9
	GO:0014069	Postsynaptic density	0.014673212	6
	GO:0033116	Endoplasmic reticulum	0.013650085	4
MF	GO:0005515	Protein binding	1.16E−07	25
	GO:0005524	ATP binding	0.0077614	22
	GO:0004674	Protein serine kinase activity	8.91E−04	11
	GO:0004672	Protein kinase activity	6.25E−04	11
	GO:0044325	Ion channel binding	2.94E−04	7
	GO:0003730	mRNA 3'-UTR binding	0.00704648	4
	GO:0019962	Type I interferon binding	0.01593083	2
	GO:0004905	Type I interferon receptor activity	0.01593083	2

DE, differentially expressed; ceRNA, competitive endogenous RNA; BP, biological processes; CC, cellular components; MF, molecular functions; GO, Gene Ontology.

Table 5

KEGG pathway enrichment analysis of the DE RNAs from the ceRNA network

ID	Description	P value	Count
hsa05200	Pathways in cancer	3.49E−04	13
hsa04151	PI3K-Akt signaling pathway	4.41E−04	12
hsa04722	Neurotrophin signaling pathway	2.02E−06	10
hsa05160	Hepatitis C	2.86E−04	8
hsa04910	Insulin signaling pathway	3.58E−04	8
hsa04012	ErbB signaling pathway	1.85E−04	7
hsa05214	Glioma	3.96E−04	6
hsa05220	Chronic myeloid leukemia	6.37E−04	6

KEGG, Kyoto Encyclopedia of Genes and Genomes; DE, differentially expressed; ceRNA, competitive endogenous RNA.

Key genes in the PPI network

The PPI network was constructed based on STRING in Cytoscape (). The MCC method from the cytoHubba app in Cytoscape was used to screen for genes with higher scores, which were considered key genes. The top 10 key genes were HSP90AA1, CDKN1A, MCL1, MDM2, MAPK1, ABL1, LYN, CRK, CDK9, and FAS (, ).

Figure 6

Identification of hub genes from the PPI network using the MCC method. (A) Eighty-five genes in the PPI network. (B) Top 10 key genes screened by the MCC method; red denotes the highest scores calculated by the MCC method, followed by orange, and lastly yellow. PPI, protein-protein interaction; MCC, maximum clique centrality.

Table 6

Key genes and their scores in the PPI network using the MCC method

Rank	Gene symbol	Description	Score
1	HSP90AA1	Intrinsic ATPase activity	162
2	CDKN1A	Ubiquitin protein ligase binding and cyclin binding	145
3	MCL1	Protein homodimerization activity and BH3 domain binding	128
4	MDM2	Protein binding and ligase activity	114
5	MAPK1	Mediate intracellular signaling	77
6	ABL1	Transferase activity	75
7	LYN	Protein tyrosine kinase activity	75
8	CRK	Protein domain specific binding	32
9	CDK9	Transferase activity, transferring	25
10	FAS	Identical protein binding	15

PPI, protein-protein interaction; MCC, maximum clique centrality.

Discussion

Because of its high risk for emergencies, CAD is the leading disease-related cause of human death (20,21). The World Health Organization estimates that 7.4 million people die of CAD every year (22). Although some progress has been made in the diagnosis and treatment of CAD, its molecular mechanisms are still unclear. Therefore, there is a pressing need for further research to identify potential targets for CAD treatment.

The focus of this study was to screen lncRNA, miRNA and mRNA differential genes related to CAD through GEO database, and then construct lncRNA-miRNA-mRNA network. Finally, 10 key genes and some signaling pathways were identified, which provided a better entry point for the basic research on the pathological mechanism of CAD in the future.

Increasing evidence indicates that lncRNAs can competitively bind miRNAs through sponge adsorption to modulate cell proliferation, metastasis, differentiation, and apoptosis to regulate the initiation and progression of diseases (23). For example, the lncRNA KCNQ1OT1 mediates miR-466i-5p downregulation, inducing high expression of the target gene Tead1 and leading to cardiomyocyte damage (24). In addition, the lncRNAs SNHG14 and SNHG7 competitively sponge miR-322-5p and miR-34-5p, respectively, increasing the expression levels of PCDH17 and ROCK1, leading to cardiomyocyte hypertrophy and fibrosis (25,26).

Downregulation of downstream miR-125a-5p via the lncRNA NEAT1 leads to overexpression of the target gene BCL2L12 and results in cardiomyocyte apoptosis (27). The lncRNA HOTAIR downregulates miR-545 to increase the expression of EGFR and p-ERK, significantly improving cardiomyocyte activity and inhibiting cell apoptosis (28). Thus, these lncRNAs from the ceRNA network play important roles in CAD, aging, and apoptosis. In this study, a lncRNA-miRNA-mRNA ceRNA network was constructed through bioinformatics to identify candidate molecules for the treatment of CAD.

We found the key genes in the constructed PPI network to be HSP90AA1, CDKN1A, MCL1, MDM2, MAPK1, ABL1, LYN, CRK, CDK9, and FAS. HSP90AA1 is the most extensively studied member of the heat shock protein (Hsp) family, whose main role is to maintain protein homeostasis and cell protection. HSP90AA1 overexpression reduced the apoptosis of neonatal rat ventricular cells induced by oxygen glucose deprivation (29). CDKN1A encodes a potent cyclin-dependent kinase inhibitor, which plays a crucial regulative role in cell-cycle progression. Knockdown of CDKN1A can inhibit cardiomyocyte hypertrophy and fibrosis while protecting myocardium in mice (30). MCL1 encodes an anti-apoptotic protein, which is a member of the Bcl-2 family. Knockout of MCL1 gene can cause mitochondrial dysfunction, which impairs the development of autophagy and heart failure (31). In the mouse model of atherosclerosis, the combination of lncRNA-p21 and MDM2 leads to the proliferation of vascular smooth muscle cells (VSMCs), reduces the apoptosis of VSMCs, and participates in the pathogenesis of atherosclerosis (32).

MAPK1 is a protein-coding gene with transferase activity and tyrosine kinase activity that is involved in the transfer of phosphorus-containing groups in signaling pathways. Studies have shown that MAPK1 is upregulated by miR-140-3p and inhibits CAD cell apoptosis (33). Knockout of ABL1 gene inhibits c-Abl activity and significantly reduces apoptosis of VSMCs and synthetic phenotypic transformation induced by Ang II both in vivo and in vitro (34). CRK plays a key role in Rac1-induced membrane ruffling and Rap1-mediated nascent focal complex stabilization, which contributed to ephrin-B1-induced human aortic endothelial cells migration (35).CDK9 has been shown to regulate cardiomyocyte hypertrophy, and recent evidence suggests that it is involved in cardiomyocyte proliferation (36). FAS encoded by this gene is a member of the TNF-receptor superfamily, which contains a death domain. Fas and FasL show interdependence with inflammatory markers in the process of apoptosis in patients with ischemic heart disease (37).

The PI3K-Akt signaling pathway is aberrantly activated during the progression of heart disease. Overexpression of IGF-1 can activate the PI3K-Akt pathway, inducing physiological myocardial hypertrophy and myocardial infarction (38). The NF-κB signaling pathway is also frequently involved in the pathogenesis of heart diseases. For example, miR-21 protects cardiomyocytes from apoptosis that is induced by palmitate through the caspase-3/NF-κB signal pathways (39). Consistent with these results, we found that the PI3K-Akt signaling pathway and the NF-κB signaling pathway were enriched by the DEmRNAs, suggesting that these pathways play an important role in the pathology of CAD. However, further tissue and cell studies still need to be carried out to validate the expression differences of the predicted key genes and determine their roles in the relevant pathways.

Conclusions

The ceRNA network constructed in this study identified new candidate molecules involved in the pathogenesis of CAD and may lead to improved treatment of CAD patients.

The article’s supplementary files as