Microarray Analysis of Differentially Expressed LncRNAs with Associated Co-Expression and CeRNA Networks in Coronary Heart Disease

Objectives: The present study aims to explore the expression profiles and biological functions of long-chain noncoding RNA (lncRNA) in coronary heart disease (CHD) and to construct a lncRNA/microRNA (miRNA)/messenger RNA (mRNA) network for mechanism exploration. Methods: Ten human blood samples (five each from the CHD and control groups) were included. The differentially expressed lncRNAs, miRNAs, and mRNAs were evaluated using microarray. Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis for differentially expressed lncRNAs–mRNAs was performed. The miRNAs that interact with differentially expressed lncRNAs were predicted using the miRanda package, whereas the mRNAs with clear biological functions and regulated by miRNAs were predicted using miRWalk2.0. Expression networks, including coding/noncoding gene co-expression, cis-regulation, and lncRNAs-transcription factors (TFs), were constructed using bioinformatics methods and were then combined to form a lncRNA-miRNA-mRNA network. Results: A total of 320 lncRNAs, 25 miRNAs, and 953 mRNAs were differentially expressed in CHD and healthy control cases. GO and KEGG pathway analysis reveal the role of these differentially regulated lncRNAs. Co-expression analysis showed that hundreds of lncRNAs, including lncRNA HIT000066433, ENST00000450016.1, and NR_028272.1, were correlated with CHD. Cis-regulation indicated that uc003qsf.4 had maximum adjacent coding genes called SOD2 and MYADM. The lncRNA-TF diagram indicated that lncRNA HIT000066433 can regulate mRNA expression through Nkx2-5 TF. Finally, the constructed lncRNA-miRNA-mRNA network confirmed that the RNA interactions show a novel perspective for the mechanisms of CHD. Conclusion: LncRNAs harbor miRNA response elements and are involved in the pathogenesis and development of CHD.


Introduction
The incidence of coronary heart disease (CHD) has increased with the improvement in the economic and living standards of people. CHD is a major threat to human health and a complex disease caused by multi-factors. The mortality rate of CHD has increased in 2002-2014. In 2015, the mortality rate of CHD in China was 110.67/100 000 in urban areas and 110.91/100 000 in rural populations, which were both higher than the cases in 2014 [1]. With the changing living standards, dietary habits, modes of work, exercise levels, and epigenetic factors, the onset age for CHD has gradually decreased [2]. Epigenetics are called "bridges" that connect the environment, genes, and diseases. This field connects environmental factors and the genomes affecting gene expression in the absence of genomic variation, thus affecting the sensitivity of the disease. Cardiovascular diseases are complex diseases regulated by epigenomics. The current epigenetic modifications associated with CHD mainly include DNA methylation, microRNA (miRNA), long-chain noncoding RNA (lncRNA), and circular noncoding RNA (circRNA) [3,4]. Approximately 5%-10% of sequences are transcribed in a human genome. Among these transcripts, 10%-20% are protein-coding RNAs, whereas the remaining 80%-90% are non-protein-coding RNAs.
Many studies have revealed the role of noncoding RNAs in cardiovascular diseases [5,6]. LncRNAs are non-coding RNAs that are longer than 200 nt [7] and can be classified as sense lncRNAs, antisense lncRNAs, bidirectional lncRNAs, intronic lncRNAs, and intergenic lncRNAs [8], lncRNAs have the following mechanisms: regulating histone modifications; modulating DNA methylation; regulating chromatin remodeling; interacting with transcription factors (TFs); regulating genome organization by looping of enhancers; and regulating the post-transcriptional level [9] to participate in the regulation of smooth muscle cell and epithelial cell function, macrophage, immunity, and lipid metabolism [10], which are closely related to the occurrence of CHD. Many lncRNAs are molecular markers and can be used as potential treatment targets of CHD. miRNAs form a class of small noncoding RNAs that regulate post-transcriptional gene expression under physiological and pathological conditions.
MiRNAs are highly conserved during evolution and has a significant role in cell proliferation, differentiation, metabolism, apoptosis, development, aging, and in the pathophysiology of many diseases such as cancer, cardiovascular diseases, and neurological disorders [11]. Salmena and colleagues presented a competing endogenous RNA (ceRNA) hypothesis [12] that describes a miRNAmediated complex post-transcriptional regulatory network, which shares miRNA response elements (MREs), competes with proteincoding and noncoding RNAs for binding to miRNAs, and regulate the expression of each other. LncRNAs play an important role in CHD [13]; however, their functions and mechanisms in CHD are poorly understood. The present study investigated the different expression patterns of lncRNAs, miRNAs, and messenger RNAs (mRNAs) in the blood samples of the CHD case group and control group by using microarray.
Several lncRNAs, miRNAs, and mRNAs with significantly different expression levels were identified between the case and control groups (FC > 2, p < 0.05). Gene Ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways were detected, and comprehensive functional landscapes of codingnoncoding gene co-expression (CNC), cis-regulation, lncRNA-TF, and ceRNAs were delineated in the two groups using bioinformatics approaches. The expression profiles of lncRNAs, mRNAs, and miRNAs among the five CHD and five healthy control cases were investigated by using microarray analysis. LncRNA-miRNA-mRNA regulatory networks were constructed, and their functions in CHD were explored to provide reference for further study on the regulating mechanisms of lncRNA in CHD.

Patients and Samples
In this study, 10 individuals were classified into two cohorts (their demographic characteristics are presented in (Table 1). Whole blood samples from five CHD and five healthy control cases in The First Affiliated Hospital of Fujian Medical University and Fujian Medical University Union Hospital were collected from July 2016 to August 2016. Based on previous studies, 10 blood samples were considered to be representative [14]. CHD was diagnosed and defined as coronary stenosis ≥ 50% in at least one of the coronary arteries according to the American College of Cardiology/American Heart Association guidelines [15]. The diagnosis was performed independently by two experienced interventional cardiologists through visual observation. Exclusion criteria were as follows: The control subjects did not show signs of coronary atherosclerosis or Microvascular diseases. This study was approved by the Ethics Committee of the two hospitals. All selected patients provided informed consent, and the study protocol conformed to the ethical guidelines of the Declaration of Helsinki (1975).

RNA Extraction
Total RNA was extracted from the patients by using a fast total RNA extraction kit (Biotech, Beijing, China). This method has been used in previous studies [16]. RNase-free water was used to dissolve the RNA. The OD260/280 reading was measured with a spectrophotometer (NanoDrop ND-1000) to determine the purity and concentration of RNA. RNA integrity was determined by 1% formaldehyde denaturing gel electrophoresis.

Microarray Imaging and Analysis
RNAs of the peripheral blood from the patients were extracted for microarray analysis. The extracted RNAs were digested, dephosphorylated, denatured, amplified, and labeled with Cy3-dCTP according to the manufacturer specifications. Data summarization, normalization, and quality control were performed using Gene Spring software V13.0 (Agilent) to analyze lncRNA + mRNA array data. To select the differentially expressed genes, we used threshold values of > 2 and < −2-fold change and a Benjamini-Hochberg corrected P-value of < 0.05. Adjust Data function of CLUSTER 3.0 software was used to obtain and adjust the mean values of the data for the log2 transform, which were further analyzed with hierarchical clustering with average linkage.

Gene Ontology and KEGG Pathway Analysis
GO analysis (http://www.geneongoloty.org/) was conducted to construct meaningful annotation of genes and gene products in various organisms. The ontology covers areas such as biological processes, cellular components, and molecular functions. The log10 (p-value) denotes the enrichment score indicating the significance of GO term enrichment in the differentially expressed genes. KEGG pathway analysis (http://www.genome.jp/kegg/) was also conducted to capture the clusters of pathways involved in the molecular interaction and reaction networks in differential regulatory gene analysis. Moreover, the -log10 (p-value) denotes the enrichment score showing the significance of the pathway correlations. Both analyses were performed by DAVID.

Construction of the Coding Noncoding Gene Co-Expression Network
The CNC network was constructed based on the correlation analysis between the differentially expressed lncRNAs and mRNAs. For each pair of genes, Pearson correlation was calculated, and the significant correlation pairs were selected to construct the network. LncRNAs and mRNAs with Pearson correlation > 0.99 or < −0.99 were selected to draw the network using the open-source bioinformatics software Cytoscape [17]. In the network analysis, a degree of centralization is defined as the number of nodes that are linked to other nodes. A degree is the simplest and most important measure of the centrality of a gene and can be used to determine the relative importance of a network [18].

Cis-and Trans-Regulation Prediction and Transcription Factor Prediction
According to lncRNA-mRNA co-expression analysis, cis-acting lncRNA was predicted by using their tight correlation (Pearson's correlation coefficient minimum of 0.99) to a group of expressed protein-coding genes. The lncRNAs resided at genomic loci where a protein-coding gene and an lncRNA gene were within 10 kb of each other along the genome [19]. "Cis" therefore refers to the same locus (not necessarily same allele) regulatory mechanisms, which include antisense-mediated regulation by lncRNAs of proteincoding genes that are encoded in the same locus. Trans-prediction was conducted using BLAT tools with default parameter settings to compare the full sequence of the lncRNAs with the 3′ UTR of their co-expression mRNAs. The TFs were predicted using Match 1.0 Public TF prediction tool on the basis of the co-expression results. The initial sites of each lncRNA upstream 2000 bp and downstream 500 bp region were included.

Competing Endogenous RNA Network Analysis
LncRNAs, miRNAs, and mRNAs whose expression levels were meaningfully correlated were analyzed. Potential MREs were searched on the lncRNA and mRNA sequences, and the overlapping of the same miRNA seed sequence binding site on the lncRNA and mRNA was used to predict the lncRNA-miRNA-mRNA interaction. The miRNA-mRNA interactions were predicted by miRWalk2.0, whereas the lncRNA-miRNA interactions were predicted by miRanda. The lncRNA-miRNA-mRNA networks were constructed using Cytoscape.

Baseline Demographic Characteristics of Five CHD and Five Healthy Control Cases
Five CHD and five healthy control cases were enrolled in our study from July 2016 to August 2016 ( Table 1). The risk factors for CHD and control groups were not significantly different (p > 0.05).

Identification of Dysregulated lncRNA, miRNA, and mRNA Expression Profiles in the Two Groups
Microarray is an effective approach for investigating the biological function of RNAs. According to the expression profiles, 320 lncRNAs, 25 miRNAs, and 953 mRNAs were differentially expressed (FC > 2, p < 0.05) between the CHD and healthy control cases. Among these RNAs, 88 lncRNAs, 22 miRNAs, and 70 mRNAs were significantly up-regulated (FC > 2, p < 0.05) in the CHD blood sample, whereas 232 lncRNAs, 3 miRNAs, and 883 mRNAs were down-regulated in the CHD blood sample as compared with those in the healthy control cases (Figure 1). The results showed good clustering between the samples, indicating that the RNAs in the CHD group have abnormal expression.

Classification of Differentially Expressed lncRNAs
The 88 over-expressed lncRNAs and 232 down-regulated lncRNAs were further classified according to their different features, including genome location and context, exerted effect on DNA, and targeting mechanisms ( Table 2) and ( Figure 2). Three up-regulated and three down-regulated sense lncRNAs and 136 (22 up-regulated and 114 down-regulated) intronic antisense lncRNAs were found. For divergent lncRNAs, three were up-regulated and two were down-regulated; for intergenic lncRNAs, 46 were upregulated and 64 were down-regulated; and for intronic lncRNAs, one was up-regulated and 10 were down-regulated.

Screening of Differentially Expressed RNAs
According to the different expression fold changes of lncRNAs and miRNAs, 10 differentially expressed lncRNAs and 10 differentially expressed miRNAs were selected, most of which were significantly up-regulated or down-regulated. In the chart, ENST00000602339.1 (FC: 4.25) was the most down-expressed lncRNA, whereas miR-6130 (FC: 28.43) was the most overexpressed miRNA in CHD cases when compared with their normal counterparts (Tables 3 & 4).

Results of Gene Ontology and KEGG Pathway Analysis
LncRNAs can regulate the expression levels of neighboring and overlapping coding genes. Thus, their functions might be based on the annotations of their co-expressed mRNAs [20]. GO analysis among all differentially expressed mRNAs can reveal the role of the differentially regulated lncRNAs. The analysis showed that the mRNAs associated with biological processes were involved in celltype-specific apoptotic process (GO: 0097285) and endothelium development (GO: 0003158). The mRNAs associated with cellular component were tumor necrosis factor receptor super family complex (GO: 0002947) and membrane-bounded vesicle (GO: 0031988). The mRNAs associated with molecular function were most relevant to BH domain binding (GO: 0051400) and BH4 domain binding (GO: 0051435) ( Figure 3A). The mRNAs and lncRNAs were enriched in GO terms with regard to the regulation of endothelial cell differentiation and endothelium development with the down-regulated expression levels of BTG1, MYADM, CCM2, and STK4. This finding indicated that the differences in endothelial cell differentiation and development may regulate the expression of genes in the blood of patients with CHD. KEGG pathway analysis was further performed to investigate the pathways and molecular interaction related to genes. Most of the significant pathways are viral carcinogenesis (hsa05203), staphylococcus aureus infection (hsa05150), and FoxO signaling (hsa04068) ( Figure 3B). Some identified KEGG pathways include FoxO signaling pathway (hsa04068), GnRH signaling pathway (hsa04912), estrogen signaling pathway (hsa04915), and MAPK signaling pathway (hsa04010). The results suggest that these signaling pathways occupy a significant position in the development of CHD, and lncRNAs may play an important role through these signaling pathways.

MREs of Associated lncRNAs Predicted by miRanda
In order to found what kind of miRNAs is targeted on above lncRNAs, Given that miRNAs interact with lncRNAs through their MREs within the ceRNA network, the potential MREs in lncRNAs were searched. A microarray for miRNAs was performed, and the top 10 key miRNAs were selected for analysis (Table 5). According to microarray analysis, the most significant expression of two miRNAs, namely, miR-6130 and miR-4726-5p, was selected for further study.

mRNAs Targeted by miRNAs
Searching for miRNA-mRNA targets is the next step to establish the lncRNA-miRNA-mRNA networks (ceRNA networks). MiR-6130 and miR-4726-5p interacted with some lncRNAs and were the most significantly expressed according to microarray analysis (Table 6). These two miRNAs were negatively correlated with mRNA.

Construction of lncRNA and mRNA Co-Expression Network and Function Prediction
The co-expression network was constructed based on correlation analysis to investigate the correlation between the differentially expressed lncRNAs and mRNAs. A total of 37548 lncRNA-mRNA pairs were identified with an absolute value of Pearson correlation coefficients. Correlation > 0.90 or < −0.90 was included. The lncRNAs and mRNAs with Pearson correlation coefficients > 0.99 or < −0.99 were selected for constructing the network. One mRNA can be correlated with several lncRNAs, which is in accordance with previous studies [21]. In particular, the most under-expressed lncRNA, HIT000066433, was positively correlated with 133 mRNAs. Among these mRNAs, HIT000066433 was positively correlated with PTEN, SOD2, and QKI and was inversely associated with SNAPC1, TRIOBP, and IFNLR1. Correlation analysis found that 183 mRNAs were correlated with lncRNA ENST00000450016.1, which was significantly down-regulated in myocardium of hypertrophic cardiomyopathy (FC: 7.14) [22]. Among these mRNAs, ENST00000450016.1 was positively correlated with PTEN, lnc-PPIA-1, COTL1, PBXIP1, and BTG1. NR_028272.1, one of the most significantly differentially expressed lncRNA, was correlated with QKI, CALM2, USP3, DUSP11, MXD1, and DUSP1. The top 10 lncRNAs and some of their targeted mRNAs are shown in (Table 7).

Cis-And Trans-Regulating Function Prediction of Lncrnas (Lncrnas Nearby Coding Genes and Lncrna-Tf)
Dysregulated lncRNAs that play a cis-or trans-regulatory role in mRNA genes were further detected according to lncRNA-mRNA coexpression analysis results (correlation > 0.99 or < −0.99, p-value < 0.05). The co-expressed protein-coding genes were defined as cisregulated genes with one differentially expressed lncRNA within 10 kb. Fifty-nine pairs of lncRNAs and their neighboring coding genes were found (Figure 4). For example, NR_037929.1 had one nearby coding gene, namely, ATP5E, whereas uc003qsf.4 had four adjacent coding genes, including SOD2 and MYADM. More valuable clues for these lncRNAs with nearby coding genes could be furnished throughout the networks. The TFs correlated with lncRNAs were searched based on the coexpression results to explore the roles of lncRNAs in CHD. The initial sites of each lncRNA upstream 2000 bp and downstream 500 bp region were included. The trans-regulatory functions of lncRNAs were predicted by the TFs that could regulate their expression. Some lncRNAs may be involved in the specific pathways regulated by TFs. Given that lncRNAs could have trans-regulatory functions, each lncRNA could connect with one to more than a dozen TFs, and each pair of lncRNA-TF was the result of several gene enrichments ( Figure 5). The lncRNAs regulated by TFs were Ap-1, Comp1, Elk-1, Evi-1, FOXD3, HNF-1, HNF-4, Nkx2-5, Oct-1, Pax-6, and Pax-4, suggesting that these TFs can be correlated with the development of CHD. Construction of lnc RNA-MiRNA-mRNA ceRNA Networks Figure 6: ceRNA network in CHD. The ceRNA network has been based on lncRNA/miRNA, miRNA/mRNA, and lncRNA/ mRNA interactions. In this network, blue represents miRNA, yellow represents lncRNA, and green represents mRNA.

Discussion
Obesity, hyperglycemia, hyperlipidemia, hypertension, smoking, and hyperuricemia are closely related to CHD. The role of Epigenetics in CHD development is currently being revealed. The traditional view of gene regulation is focused on proteincoding genes until the discovery of numerous noncoding RNAs, including lncRNAs and circRNAs. In particular, several studies on the deregulated lncRNA expression levels in CHD suggest that the abnormal lncRNA expression might contribute to the occurrence and development of CHD [24,25]. Recently, circRNAs were found to play essential roles in CHD [26]. However, few studies report the coexpressed profiles of lncRNA-miRNA-mRNA in CHD. To explore the functions of lncRNA-miRNA-mRNA regulatory networks in CHD, we enrolled five CHD and five healthy control cases and attempted to construct genome-wide lncRNA-miRNA-mRNA regulatory networks by using microarray.
The expression profiles of lncRNAs, miRNAs, and mRNAs were obtained by microarray. A total of 320 lncRNAs, 25 miRNAs, and 953 mRNAs were differentially expressed (FC > 2, p < 0.05). Among these RNAs, 88 lncRNAs, 22 miRNAs, and 70 mRNAs were significantly up-regulated (FC > 2, p < 0.05), whereas 232 lncRNAs, 3 miRNAs, and 883 mRNAs were down-regulated. The top 10 lncRNAs and miRNAs were identified according to the fold change of differential expression. Considering the large number of differentially expressed lncRNAs, miRNAs, and mRNAs, the occurrence and progression of CHD might be related to the altered expression of lncRNAs, miRNAs, and mRNAs. GO analysis for lncRNAs suggests that some terms in the biological process were related to the occurrence and progression of CHD.
The mRNAs and lncRNAs were enriched in GO terms related to the regulation of endothelial cell differentiation and endothelium development, which was consistent with previous reports [27,28]. KEGG pathway analysis for lncRNAs revealed many pathways, such as FoxO, GnRH, estrogen, and MAPK signaling pathways, which play pivotal roles in CHD. These results provided a new idea for the study of CHD. However, further research is needed to determine these findings. Most lncRNA functions are still poorly understood.
Constructing CNC co-expression network allowed the prediction of lncRNA functions [29]. The results showed that hundreds of lncRNAs were significantly associated with many mRNAs. TLR2 [30] and HSPB1 [31], which are major connecting mRNAs, are associated with CHD.
Therefore, we deduced that these lncRNAs can regulate coexpression genes, which might be related to the development of CHD. LncRNAs can influence the transcription of neighboring or remote genes through cis-and trans-regulatory mechanisms [32]. Several Dysregulated lncRNAs that play cis-or trans-regulatory roles in mRNA genes have been found in this research. NR_037929.1 had one nearby coding gene, namely, ATP5E, whereas uc003qsf.4 had four adjacent coding genes, including SOD2 and MYADM. The four nearby coding genes corresponding to down-regulated uc003qsf.4 were also down-regulated, indicating their positive regulatory role at the transcriptional level. These findings suggest that uc003qsf.4 plays a vital role in CHD. Some lncRNAs were regulated by TFs [33]. The TFs combining cis-elements can also regulate the expression of target genes at promoter locations [34].
The analysis of lncRNA-TF prediction revealed that Ap-1, Elk-1, Comp1, Evi-1, FOXD3, HNF-1, HNF-4, Nkx2-5, Pax-6, and Pax-4 were the most relevant TFs. A previous research has revealed that some of these TFs, such as Nkx2-5, were associated with CHD [35]. LncRNA HIT000066433 can regulate mRNA expression through Nkx2-5 TF. Therefore, we hypothesized that these TFs can play vital roles in the occurrence and development of CHD through regulating the expression of lncRNAs and mRNAs. Further investigation is necessary to understand the relationship of lncRNAs-TFs. LncRNAs were recently proposed to harbor miRNAs and have important effects in regulating the gene expression at the post-transcriptional level. A lncRNA-miRNA-mRNA ceRNA network of CHD was constructed based on microarray data. Several reports revealed that some lncRNAs might be implicated as ceRNAs in CHD. However, the functions of lncRNA-associated ceRNAs in CHD are poorly understood. MiR-6130 and miR-4726-5p were selected for further study because of their highest fold changes in their expression. Some ceRNA networks, including the top 10  built. For instance, lncRNA NR_028272.1 is a ceRNA of miR-4726-5p targeting CECL2, TMBIM6, IDS, FAM65B, and DUSP11; lncRNA ENST00000506718.1 is a ceRNA of miR-6130 targeting MAX, UQCRB, CANX, TPM3, and CELF2; lncRNA TCONS_00029336 is a ceRNA of miR-6130 targeting MAX, CANX, TPM3, CELF2, and QKI; and lncRNA ENST00000429901.1 is a ceRNA of miR-6130 targeting MAX, CANX, TPM3, CELF2, and QKI.
Based on the above results, these RNA interactions would provide a novel perspective for the mechanism of CHD. Given our limited understanding of miRNAs and lncRNAs, more studies should be conducted about their interactions. Among the ceRNAs listed above, QKI, a target gene of miR-6130, was accidentally discovered and has been identified as a novel locus that may serve as a predictor of incident CHD in prospective studies [36]. A previous study proposed that QKI is a central regulator of VSMC phenotypic plasticity, and an intervention in QKI activity can ameliorate pathogenic, fibro proliferative responses to vascular injury [37]. These findings about lncRNAs function in CHD were combined via lncRNA TCONS_00029336-miR-6130-QKI and lncRNA ENST00000429901.1-miR-6130-QKI ceRNA networks. However, the functions of these potential miRNAs, lncRNAs, and mRNAs in CHD are still unclear.

Conclusion
Several differentially expressed lncRNAs were identified by using microarray analysis and might be associated with the occurrence and development of CHD. Our results suggested that lncRNAs harbor MREs and are involved in the pathogenesis and development of CHD. Our findings lay a foundation for the future research of the potential roles, functions, and mechanisms of lncRNAs in CHD.