Differential ceRNA Expression and Interaction Analysis in Coronary Artery Disease

Previous studies had shown that mRNA, miRNA and lncRNA wer e associated with cardiovascular diseases. The study was aimed to explore the differential expressions of mRNA, lncRNA and miRNA between coronary artery disease (CAD) and healthy control, and their interaction in CAD. We investigated the differential expression of ceRNA between CAD and healthy control through data collected from Gene Expression Omnibus (GEO) microarrays. Furthermore, we investigated the biological function of these differential expressions of ceRNAs by Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses. Protein-protein interaction (PPI) network was created to identify the hub genes. Biosystems and literature search were performed for signaling pathways and their function of the included differential expression ceRNAs. A total of 456 miRNA expression profiles, 16,325 mRNA expression profiles, and 2,869 lncRNA expression profiles were obtained . Eleven Go and KEGG pathways (count ≥9), top 15 of PPI network node connectivity rank, and top 15 of ceRNA network node degree centrality rank were achieved at the statistical significance level (P<0.05). We further identified that several differential expressions of ceRNAs and their signaling pathways were associated with CAD through biosystems and literature search. Based on eleven Go and KEGG pathways, top 15 of PPI network node connectivity rank, and top 15 of ceRNA network node degree centrality rank in CAD population, our findings would contribute to further exploration for the molecular mechanism of CAD.


Introduction
Coronary artery disease (CAD) is a complex phenotype driven by genetic and environmental factors. However current therapies focus on addressing the role of cholesterol and lifestyle in CAD. Despite advances in the development of lipid-lowering therapies, clinical trials have shown that a substantial risk of cardiovascular disease persists af ter currently recommended medical therapy. 1 Stratification for subsequent coronary events among patients with CAD is of considerable interest because of the potential to guide secondary preventive therapies. Recently, eight microRNAs (miRNAs) were identifi ed to facilitate acute coronary syndrome diagnosis. 2 Targeting Angptl3 messenger RNA (mRNA) retarded the progression of atherosclerosis and reduced levels of atherogenic lipoproteins. 3 The expressing 9p21.3-associated long non-coding RNA ANRIL induces risk CAD phenotypes in non-risk vascular smooth muscle cells. 4 So far, it is not clear the mechanism of these RNAs in CAD and their interaction.
Noticeably, the different types of RNA molecule compe ted to bind to miRNA, which reduced the inhibitory effect o f miRNA targeting on its mRNA. 5 These competitive endogenous RNA (ceRNA) included various types of RNA transcripts, such as circular RNA (circRNA), long -chain non-coding RNA (lncRNA), pseudogenes a nd protein-encoded mRNA, which competed for miRNA through the "language" mediated by the miRNA response element (MRE) .6 After that, researchers used bioinformatics methods to predict ceRNA regulatory networks. The effect of ceRNA on the target gene and the dependence of ceRNA on miRNA would be verified at the experiments of proteins and RNAs, but the functional verification would be pe rformed at the experiments of cells and animal models.
Thus, the study was aimed to explore the d ifferential expressions of mRNA, lncRNA and miRNA between CAD and healthy control, and the interaction of them, including constructure of ceRNA regulatory netw orks, which would contribute to the molecular mechanism of CAD.

Basic Information Statistics of Differential Expression Analysis
As described in the Methods, a total of 45 6 miRNA expression profiles (supplementary Table1), 16,325 mRNA expression profiles (supplementary Table 2), and 2,869 lncRNA expression profiles were obtained (supplementary Table 3).
According to the set threshold, 18 differentially expres sed miRNAs were finally obtained, including 16 down-regulated and 2 up-regulated (supplementary Table 4). a total of 92 differential lncRNAs were obtained, including 46 down-regulated and 46 up-regulated (supplementary Table 5).
A total of 610 differential mRNAs were obtained, including 244 down-regulated and 366 up-regulated (supplementary Table 6).
Based on the obtained differential miRNA, lncRNA and mRNA, the heat map was shown in Figure 1 Table 1. Go and KEGG pathway enrichment analysis of differential genes

. Protein Interaction Network Construction (PPI) and Module Analysis
As described in the methods, we achieved a total of 388 protein interaction relationship pairs, and the network construction was performed using Cytoscape software as shown in Figure 4. A total of 171 nodes were included in the network. The network was analyzed for node connectivity according to the parameters set by the Method, the top15 of Degree Centrality (DC) of each node was ranked in Table 2. Notably， CXCL8, FPR2, IL6, and PPBP were ranked in Top15, which might be hub proteins in the network (supplementary Table 8 Table 2. PPI network node connectivity rank ( TOP15)

. MRNA and lncRNA Co-expression Analysis
We performed co-expression analysis of differentially expressed mRNAs and lncRNAs. According to the threshold set by the Methods, we screened a total of 1487 significantly coordinately expressed relationship pairs, including 381 mRNAs and 74 lncRNAs ( supplementary Table 9)

. MiRNA Target Genes and Upstream lncRNA Prediction Analysis
Based on the differentially expressed miRNAs and differential lncRNAs, a total of 452 lncRNA-miRNA relationship pairs were predicted as described in the Methods (supplementary  (supplementary Table 13), here we showed a part of the results in Figure 5.

. CeRNA Network Analysis
As   Table 14). There were a total of 36 lncRNAs, 64 mRNAs, and 15 miRNAs. Connectivity analysis was performed on each node of the ceRNA network to obtain mRNA, miRNA, and lncRNA connectivity as detailed in Table 3.  Table 3. ceRNA network node degree centrality rank (TOP15)

Discussion
The study analyzed the differential genes at the statistical significance level between CAD and healthy control, and foun d eleven Go and KEGG pathways (count ≥9), top 15 of PPI netw ork node connectivity rank, and top 15 of ceRNA network node degree centrality rank, which would contribute to further exploration for the molecular mechanism of CAD.
Firstly, the Go and KEGG pathways (count ≥9) showed their function in Table 1, and they contained a large number of differential genes at the statistical significance level in CAD, thus these Go and KEGG pathways might play a role of molecular level of CAD. Several further explored experiments related to these pathways would achieve the interesting and important findings in CAD field.
Secondly, in top 15 of PPI network node connectivity rank, we found that the extensive protein interaction relationship pairs at the statistical significance level in CAD (Table 2), which matched with their signaling pathways in the biosystems (supplementary Table 15). Some of them were identified to play a role in cardiovascular diseases. For example, GNG11 was a member of the gamma subunit family of heteromeric G-protein. Overexpression of GNG11 activated ERK1/2 of the MAP kinase family, but did not Ras. 20 These findings provide clinically relevant biological insight into heritable variation in vagal heart rhythm regulation, with a key role for genetic variants (GNG11, RGS6) that influence G-protein heterotrimer action in GIRK-channel induced pacemaker membrane hyperpolarization (supplementary Table 15). 21

CXCL1 was produced mainly by TNF-stimulated endothelial cells (ECs) and
pericytes and supported 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. 22 Angiotensin II-induced infiltration of monocytes in the heart is largely mediated by CXCL1-CXCR2 signalling which initiates and aggravates cardiac remodelling. Inhibition of CXCL1 and/or CXCR2 may represent new therapeutic targets for treating hypertensive heart diseases (supplementary Table 15). 23 Wnt-Cxcr4 (C-X-C motif chemokine receptor 4) signaling in regulation of oligodendrocyte precursor cells (OPCs)-endothelial interactions coordinates OPC migration with differentiation. 24 Many of the neutrophils reenter the vasculature and have a preprogrammed journey that entails a sojourn in the lungs to up-regulate CXCR4 before entering the bone marrow, where they undergo apoptosis. 25 Table 15). 26 Oleic acid treatment decreases the insulin sensitivity of heart muscle cells, and this sensitivity is completely restored by transfection with SNAP23. Thus, SNAP23 might be a link between insulin sensitivity and the inflow of fatty acids to the cell (supplementary Table 15). 27 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 identify a new signalling axis of AMPK-SKP2-CARM1 in the regulation of autophagy induction after nutrient starvation (supplementary Table 15). 28 GPR18 is a cannabinoid-activated orphan G protein-coupled receptor (GPCR) that is selectively expressed on immune cells. 29 A salutary cardiovascular role for GPR18, mediated, at least partly, via elevation in the levels of adiponectin (supplementary Table 15). 30 Thirdly, in top 15 of ceRNA network node degree centrality rank, the extensive ceRNA interaction relationship pairs at the statistical significance level in CAD (Table 3) , which matched with their description in the literature search (supplementary Table 16). Certainly, we also found that several ceRNA was  Table 16). 45 Table 16). 46

. MRNA and lncRNA Annotation
The sequences matched the probes of Agilent-028004 were obtained from the platform annotation file, and the human reference genome (GRCh38) sequences was downloaded from the GENCODE database 8 (https://www.gencodegenes.org/releas es/current.html), and the probe sequences were aligned onto the reference genome using the seqmap software. 9 Firstly, we retained the uniquely aligned (unique map) probes, and secondly, we referred their position to the chromosome with positive and negative strand information, the gene matched each probe was obtained according to the human gene annotation file ( Release 25) provided by GENCODE.
We kept the probe with the annotating information "protein_coding" as the matching probe for mRNA, the probes with the annotating information with "antisense", "sense_intronic", "lincRNA", "sense_overlapping" or "processed_transcript" were considered to the matching lncRNA probe.
Finally, the probe numbers and mRNA/lncRNA (Gene symbol) were matched one by one to remove probes that did not match to Gene symbol.
For different probes mapping to the same gene, we used the average of different probes as the final expression value of the mRNA/lncRNA.

Differential MRNA, lncRNA and MiRNA Screening
We took the R software limma package with the classical Bayesian method 10 (version 3.10.3, http://www.bioconductor.org/packages/2.9/bioc/html/limma.html) . The differential analysis was performed between CAD and CTRL. Importantly, the miRNAs, mRNAs, and lncRNAs were analyzed to obtain their p values and logFC values, which were evaluated at the levels of both fold difference and statistical significance. The threshold of differential expression was set as miRNA: p value < 0.05 and |logFC| > 0.263 (>1.2 times).

Functional Enrichment and Pathway Analysis of Differentially
Expressed MRNA Enrichment analysis was performed with the common enrichment analysis tool DAVID 11 (version 6.8, https://david-d.ncifcrf.gov/) to analyze the up-and down-regulated genes, which were involved in the pathways of Gene Ontology BP (biological process), 12 CC (cellular component), MF (molecular function) and KEGG. 13 Significant enrichment results were considered to a significance threshold p value < 0.05 and an at least 2 of enrichment number (count).

Protein Interaction Network (PPI) Construction and Node Connectivity Analysis
In combination with STRING (version: 10.0, http://www.string-db.org/), 14 the database predicted whether there was an interaction relationship between the proteins encoded by the analyzed genes. The differential mRNAs were inputted into gene sets, and homo was inputted into species. The parameter of PPI score was set to 0.9 (highest confidence), the interactional protein nodes were required in the up-and down-regulated genes. After the PPI relationship pairs were obtained, the data were analyzed using Cytoscape software (version 3.4.0, http://chianti.ucsd.edu/cytoscape-3.4.0/), 15 for which a network map was constructed. The node connectivity analysis was performed with parameters of no weigh by using the CytoNCA plugin (Version 2.1.6, http://apps.cytoscape.org/ apps/cytonca). 16 The results obtained the important nodes in the PPI network that were involved in protein interaction relationships, i.e., hub proteins, through the connectivity Degree Centrality (DC) rank of individual nodes.

lncRNA and MRNA Co-expression Analysis
The correlation test was performed and their pearson correlation coefficients of the differential mRNA and lncRNA were respectively calculated by using the matched sample of mRNA and lncRNA data. The relationship pairs with r > 0.7 (coordinate expression) and p value < 0.05 were focused on screening for the subsequent ceRNA network co nstruction, and these mRNAs were considered to be s ignificantly correlated with lncRNAs, while mRNAs were considered as potential target genes of lncRNAs.

Target Genes MiRNA and The Prediction for Their Upstream lncRNA
Based on the differential miRNAs obt ained from the differential analysis, miRWalk2.0 (http://zmf.umm.uni-heidelberg.de/apps/zmf/mirwalk2/) database 17 was used, which integrated the four typical databases including miRWalk, miRanda, RNA22, and TargetScan. If the predicted target genes were presented in each of four databases, the marching mRNA was considered to be regulated by the miRNA. After the predicted miRNA-mRNA relationship pairs were obtained, the mRNAs were f urther intersected with the differential mRNAs to obtain the differential miRNA-differential mRNA relationship pairs.

Pathway Enrichment Analysis of lncRNAs and MiRNAs
Based on the obtained lncRNA-mRNA co-expression relationship pairs and miRNA-mRNA relationship pairs, mRNAs were used as potential target genes of matching lncRNAs and miRNAs, respectively. KEGG Pathway enrichment analysis was performed by using the R package clusterProfiler (version: 3.8.1, http://bioconductor.org/packages/release/bioc/html/clusterProfiler.html), 19 and its results indirectly predicted the functions of lncRNAs and miRNAs.
The threshold was set at p value < 0.05.

CeRNA Network Construction
Based on the obtained mRNA-miRNA and lncRNA-miRNA relationship pairs, we firstly screened the miRNA-lncRNA-mRNA relationship pairs regulated by the same miRNA, then combined the positive co-expression relationship between mRNA and lncRNA (correlation coefficient > 0.7), and further screened the miRNA-lncRNA-mRNA relationship pairs for network construction, i.e., the ceRNA network . Thus, the lncRNAs and mRNAs with positive co-expression relationship regulated by the same miRNA in the ceRNA network were each other ceRNAs.
Finally, the node connectivity (degree) analysis was also pe rformed using the Cytoscape plugin CytoNCA with the parameter set to no weight.
The higher the connectivity, the higher the import ance of this node in the network.

Conclusions
Based on eleven Go and KEGG pathways, top 15 of PPI network node connectivity rank, and top 15 of ceRNA network node degree centrality rank in CAD population, our findings would contribute to further exploration for the molecular mechanism of CAD.