Integrated analysis of lncRNA-miRNA-mRNA ceRNA network in human aortic dissection

Background Many studies on long chain non-coding RNAs (lncRNAs) are published in recent years. But the roles of lncRNAs in aortic dissection (AD) are still unclear and should be further examined. The present work focused on determining the molecular mechanisms underlying lncRNAs regulation in aortic dissection on the basis of the lncRNA-miRNA-mRNA competing endogenous RNA (ceRNA) network. Methods This study collected the lncRNAs (GSE52093), mRNAs (GSE52093) and miRNAs (GSE92427) expression data within human tissue samples with aortic dissection group and normal group based on Gene Expression Omnibus (GEO) database. Results This study identified three differentially expressed lncRNAs (DELs), 19 differentially expressed miRNAs (DEmiRs) and 1046 differentially expressed mRNAs (DEGs) identified regarding aortic dissection. Furthermore, we constructed a lncRNA-miRNA-mRNA network through three lncRNAs (including two with up-regulation and one with down-regulation), five miRNAs (five with up-regulation), as well as 211 mRNAs (including 103 with up-regulation and 108 with down-regulation). Simultaneously, we conducted functional enrichment and pathway analyses on genes within the as-constructed ceRNA network. According to our PPI/ceRNA network and functional enrichment analysis results, four critical genes were found (E2F2, IGF1R, BDNF and PPP2R1B). In addition, E2F2 level was possibly modulated via lncRNA FAM87A-hsa-miR-31-5p/hsa-miR-7-5p or lncRNA C9orf106-hsa-miR-7-5p. The expression of IGF1R may be regulated by lncRNA FAM87A-hsa-miR-16-5p/hsa-miR-7-5p or lncRNA C9orf106-hsa-miR-7-5p. Conclusion In conclusion, the ceRNA interaction axis we identified is a potentially critical target for treating AD. Our results shed more lights on the possible pathogenic mechanism in AD using a lncRNA-associated ceRNA network.


Introduction
Aortic dissection (AD) is a severe aortic disorder involving destruction of aortic wall medial layer, which can induce separation between intima and adventitia, and track blood within the dissection plate in the medial layer, thereby inducing the formation of true and false lumens in the aortic wall. Although the research on aortic dissection has gradually increased and many innovations have been made in research programs, the treatment of thoracic aortic dissection is still extremely challenging [1,2]. Aortic dissections are more common in the non-whites population of elderly men [3], and the incidence rate increases sharply over the age of 50. People between the ages of 50 and 70 are at the highest risk, and increasing age was the important variable associated with AD long-term mortality [4]. In the meanwhile, because it is often misdiagnosed when it appears, it is difficult to assess the exact incidence [5]. Although studies have demonstrated that the diagnostic techniques and treatments for thoracic AD cases are improving, the mortality and morbidity are high, so early diagnosis and treatment are very necessary [6].
The current research always draws their attention to the technical means and pathogenesis of aortic dissection. More and more recent studies have suggested that non-coding RNAs (ncRNAs), including long noncoding RNAs (lncRNAs) and microRNAs (miRNAs), exert decisive roles in the development of cardiovascular diseases (CVDs) [7,8]. lncRNAs are the ncRNAs that are over 200 bp in length, which can be used as an important type of regulatory molecule in the human genome to perform its biological functions in various ways. Many studies have shown that IncRNA can also be used to be the competitive endogenous RNA (ceRNA) as the miRNA sponges and participate in regulating the expression of target genes [9,10]. miRNAs shows negative regulation on protein-coding genes through combining with complementary sequences [11]. As a result, the lncRNA-miRNA-mRNA interaction has a vital effect on CVD development [12,13].
Ren and colleagues revealed that miR-193b-3p and lncRNA H19 had certain effect on vascular smooth muscle cell (VSMC) migration and proliferation, which facilitated to generate new thoughts in AD management [14]. Recently, Zhang et al. carried out a study in order to prove lncRNA XIST's effect on AD pathogenic mechanism and identify its corresponding pathway, with the findings showing that lncRNA XIST regulated smooth muscle cell proliferation and apoptosis through the sponge of miR-17 and the regulation of the subsequent downstream PTEN to affect the development of aortic dissection in mice [15]. Simultaneously, Zhao and colleagues suggested that lncRNA CDKN2B-AS1 regulated STAT3 level through suppressing miR-320d to regulate human VSMC apoptosis and proliferation [16]. To sum up, lncRNA-miRNA-mRNA regulatory network exerts a vital part in AD genesis and progression while the lncRNA targets, roles and underlying mechanisms within different tissues and cells of AD have not been reported. At the same time, there are few reports about the ceRNA regulation mechanism of lncRNA-miRNA related to AD and the interaction between ncRNAs.
In this study, lncRNAs together with the corresponding mechanisms of action within human tissue from AD cases were explored by means of bioinformatic analysis. Firstly, this study applied Gene Expression Omnibus (GEO) database for obtaining the AD-related lncRNA, mRNA and miRNA expression data. Thereafter, we discovered differentially expressed genes (DEGs), differentially expressed miRNAs (DEmiRs), and differentially expressed lncRNAs (DELs) using RStudio. Cytoscape 3.7.2 was utilized to construct a lncRNA-miRNA-mRNA network, followed by the construction of a PPI network.
The present study aims to further screen the key lncRNA-miRNA-mRNA ceRNA axis within AD with microarray data collected from public databases and bioinformatics methods. Besides, our results shed more lights on the AD molecular mechanism, thus providing a novel direction for targeted therapy of AD. To further examine the key gene regulatory axis, we established a lncRNA-miRNA-mRNA regulatory subnetwork.

Materials and methods
Microarray data and data processing Two datasets were obtained by setting the screening criteria for the species type as "Homo sapiens", from GEO database of National Center for Biotechnology Information (NCBI) (https://www.ncbi.nlm.nih.gov/geo/), with such keywords being searched as "Aortic dissection" (AD). Then select the study type as " Expression profiling by array", so that we get the series matrix files and platform file of non-coding RNA. as a result, a total of two datasets were included in this study, namely GSE52093 and GSE92427. As for lncRNA and mRNA expression data (GSE52093 dataset), relative to normal aortic tissue in control group, the experimental group collected the ascending aorta of patients with acute Stanford type A aortic dissection to identify differentially expressed genes. The microRNA microarray is used for expression profiling analysis from subjects of 24 plasma samples (including acute aortic dissection and healthy subjects). Each dataset contains aorta dissected group and aorta normal group. Simultaneously, we applied GEO database for downloading series matrix files and expressive data. The data flow mechanism diagram is displayed in (Figure 1).

DEGs analysis
DEGs were analyzed by R package "Linear Models for Microarray Data (limma)" function for datasets and online analysis software SPSSAU. SPSSAU (https://spssau.com/ index.html) is an online data analysis and visualization software. We inputted the expression matrix data of two datasets into SPSSAU and got the difference between aorta dissected group and aorta normal group of two datasets. In addition, we also obtained the median, extremum and outliers between each sample.
For each dataset, the sva R package was used to remove batch effects and the GeneSoring GX software package was adopted to conduct the quantile normalization for preprocessing, and the annotation information of lncRNA, miRNA and mRNA was downloaded by series matrix files of GSE52093 and GSE92427 datasets. The selection criteria for GSE52093 dataset including log|FoldChange|>1 and adjustp<0.1 were regarded as threshold values. The selection criteria for GSE92427 dataset containing log|FoldChange|>1 and adjustp<0.05 were regarded as threshold values. There existed statistical significance for the selection of this threshold, and those genes that were up-and downregulated can also be selected for performing the subsequent analysis. This indicates statistical significance. R software package was also adopted as a heat map and volcano map drawing of DELs, DEmiRs and DEmRNAs. Subsequently, these DELs, DEmiRs and DEmRNAs between two databases were classified as the up-regulation or down-regulation group. These data would be used in the following ceRNA network construction and protein interaction network construction.
Prediction of miRNA-mRNA and lncRNA-miRNA pairs Interaction relationships between lncRNAs and miRNAs were predicted with R package and miRcode database. By using the R language program to output all the predicted results of DElncRNAs, and intersecting the output results with DEmiRs, the interaction pair of lncRNAs and miRNAs can be obtained. Cytoscape (version 3.7.2) was utilized for visualizing forecast results.
Similarly, we use the R package to predict the miRNAs in the lncRNA-miRNA interaction pair, and intersect the prediction result with the DEmRNAs. As a result, the miRNA-mRNA interaction pair can be obtained. The difference is that we use three databases when predicting mRNAs, and we extract the results with 2 occurrences in the three databases as our prediction result.

Establishment of the ceRNA network of lncRNA-miRNA-mRNA
This study established a ceRNA network with related DELs by using commonly interactive miRNAs with DELs and DEmRNAs. Cytoscape (version 3.7.2) was later utilized to visualize the ceRNA network of lncRNA/ miRNA/mRNA.

PPI network analysis and key gene identification
The plug-in CytoHubba in the Cytoscape software is a visualization software that obtains the dense relationship through the degree, closeness centrality and betweenness centrality algorithms. Those hub genes in ceRNA network were identified by CytoHubba.
The present study adopted the Search Tool for the Retrieval of interacting Genes/Proteins (STRING; version 11.0) for retrieving protein interactions between DEGs identified in GSE52093. After the points without interaction were hidden, the data were imported into Cytoscape (version 3.7.2) for visualizing the protein-protein interaction (PPI) network. CytoHubba was used to identifiy the hub genes of PPI network.
The key genes were determined by taking the intersection of hub genes in the ceRNA network and PPI network. Then, we obtained the key genes for follow-up analysis.
Functional enrichment (DEGs) and differential expression analysis (key genes) Gene ontology (GO) functional annotation and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis were done for analyzing DEGs by using R language software package (Bioconductor and pathview) and DAVID online tool (version 6.8). P<0.05 indicated that the pathways or GO biological process terms were significantly enriched.
The present work adopted GraphPad Prism 8 (version 8.0.1) to visualize key DEGs in different datasets.

Differential expression analysis
In this study, we imported the corrected expression matrix data into the SPSSAU online software, and the resulting box plot was shown in (Figure 2A and B). According to the cutoff criteria of GSE52093 (log|Fold-Change|>1, adjustp<0.1) and the cutoff criteria of GSE92427 (log|FoldChange|>1, adjustp<0.05), volcano map and heat map for differential expressed lncRNAs, miRNAs and mRNAs were obtained ( Figure 2C, D, E and F). The GSE52093 dataset included lncRNA and mRNA data, in which three DELs (two with upregulation and one with down-regulation) and 1046 DEmRNAs (540 with up-regulation and 506 wtih downregulation) were discovered through the comparison of gene expression data between AD and control samples. Altogether 19 DEmiRs of GSE92427 dataset (15 upregulated and four downregulated) were identified from plasma samples including acute aortic dissection and healthy subjects. The most significant differentially expressed genes were shown in ( Table 1).

Forecast of miRNA-mRNA and lncRNA-miRNA pairs
Interactive relationships between DELs and miRNAs were predicted with lncRNA mircode and R language script. We acquired the common miRNAs through the intersection of estimated DELs with DemiRs, while interactions of LncRNAs with miRNAs were predicted by LncBase Predicted v.2. LncRNA FAM87A can bind to hsa-miR-338-3p and hsa-miR-193b-3p. LncRNA UCA1 can bind to hsa-miR-1279 and hsa-miR-455-5p ( Figure 3A). Several DemiRs binding sites that can bind to lncRNA were presented in ( Figure 3B).
Compared with lncRNA prediction, the interaction relationships between DemiRs of binding DELs and DEGs  Table 2). MicroT-CDS was employed for predicting the associations of miRNAs with mRNAs. Several DEGs binding sites that can bind to DemiRs were illuminated in ( Figure 3B).

PPI network analysis and key gene prediction
The top 30 DEGs of ceRNA network acquired based on degree, closeness centrality and betweenness centrality algorithms using cytoHubba were visualized in ( Figure  4A). When uploading the DEGs identified in ceRNA network to STRING website, 395 interaction relationships were found to be present in 211 DEGs incorporated in constructing a PPI network. In the as-constructed network, E2F2(E2F Transcription Factor 2) could interact with CCNE1(Cyclin E1). IGF1R (Insulin Like Growth Factor 1 Receptor) could interact with BNDF (Brain Derived Neurotrophic Factor). BDNF could interact with IGF1R and PPP1R9B (Protein Phosphatase 1 Regulatory Subunit 9B) ( Figure 4B). The top 30 DEGs of PPI network obtained by cytoHubba were visualized in ( Figure 4C). There were four key gens (E2F2, IGF1R, BDNF and PPP1R9B) which were determined by Venn diagram of ceRNA network and PPI network ( Figure 4D).

Functional enrichment and differential expression analysis
In the present study, the DEGs used to analyze GO/ KEGG were imported into DAVID website and R language script. The results of DAVID showed the functions and pathways involved in DEGs: GO (such as plasma membrane and protein binding) and KEGG (like pathways in cancer and PI3K-Akt signaling pathway) ( Table 3). The results of R language software showed the functions and pathways involved in DEGs: GO (such as nuclear receptor activity and steroid binding) and KEGG (such as proteoglycans in cancer and PI3K-Akt signaling pathway) ( Figure 5).
To accurately explore the differential expression of key genes in the three datasets, four key genes were visualized in (Figure 6). In the GSE52093 datasets, E2F2 and PPP1R9B were up-regulated while IGF1R and BDNF were down-regulated, further revealing the underlying mechanisms for key genes.

Discussion
In normal circumstances, the 3 layers of aortic wall (e.g. intima, media, and adventitia) should maintain the complete structures and physiological functions for maintaining the aortic wall stability and managing the great influence resulting from blood flow [17]. The AD pathogenic mechanisms are vascular inflammation, matrix metalloproteinases (MMPs) activity and change in vascular smooth muscle cells (VSMCs) phenotype [18]. Endothelial cell (EC) injury promotes the occurrence of vascular inflammation modulated via immune response. Afterwards, it will activate MMP, the extracellular matrix (ECM) degrading enzyme [19,20]. Meanwhile, AD is also associated with broad alterations of VSMC phenotype and vascular wall at molecular level [21,22].
LncRNAs, with over 200 nucleotides in length, constitute the transcript family without protein coding functions [23]. They are suggested to play roles of ceRNAs for miRNAs; in other words, lncRNAs may play roles of miRNA "decoys" for regulating gene levels [12,24]. As suggetsed by the theory, lncRNAs act as natural sponge for the competition adsorption of certain miRNAs and reduction of miRNA binding to corresponding target genes, thus resulting in alterations of miRNA target gene expression [25]. However, it remains unclear about whether the abnormal lncRNAs play roles of ceRNA for certain miRNAs and have certain influence on arterial wall by the indirect regulation of target mRNA expression in the process of AD occurrence. Recently, lncRNAs have been discovered to be related to CVD genesis and progression. Recent study revealed a pathogenic H19 induce aneurysm by inflammatory pathway related to the formation of AAA, and this has offered a novel treatment for AAA [26]. Meanwhile, Kumar and colleagues investigated the effects of ncRNAs together with the target genes from the perspective of their functions within AAA. They also discussed those animal models adopted to mechanically understand AAA and the possible effects of miRNAs and lncRNAs as diagnostic biomarker and therapeutic targets [27]. Several lncRNAs are suggested to have important functions in aorta-related disease pathogenesis, such as AD [28]. Few lncRNAs show obvious spatiotemporal expression and specificity in the process of tissue growth and differentiation; As a result, they were the favorable diagnostic biomarkers for AD.
E2F belongs to the transcriptional factor family, which functions in controlling G(1)/S transition. In addition, certain E2F components have been recently suggested to regulate functions in addition to cell cycle, like apoptosis induction [29]. Fujita et al. found that E2F5 was induced by fetal bovine serum in VSMC, and E2F5 inducibility has unique function in VSMC, which is related to the feedback modulation for some cell activities in cell proliferation [30]. Nevertheless, Angiotensin II (ANG II) increased E2F3 protein expression, rather than E2F-5 and     FZD1, HMGA1, BRCA1, MRPL12, HIF1A, MLLT11, ZXDA, LBH, TFAP4, MED13, CCNE1, BAMBI,  MYC, CKS2, TASP1,  E2F-1, but it did not increase the mRNA expression. The above alterations are related to the hyperplastic or hypertrophic responses to diverse stimuli or growth factors of VSMCs. Several studies suggested that modulation of E2F had certain effect on treating hypertension, atherosclerosis, restenosis following vascular damage and bypass graft failure [31,32]. In precent study, E2F2 expression was possibly modulated via lncRNA FAM87A-hsa-miR-31-5p/hsa-miR-7-5p in human aortic dissection tissue. Recent study investigated the implication overexpression miR-7-5p, finding that miR-7-5p agomir transfection markedly suppressed mineralization of pulmonary artery smooth muscle cell matrix under hypoxic conditions [33].
It is still unclear about the AD pathogenesis, but more and more studies have suggested the important functions of ncRNAs within AD. In this study, the underlying molecular mechanisms of AD have been investigated based on the lncRNA-miRNA-mRNA ceRNA regulatory network. lncRNAs, mRNAs and miR-NAs expression data are obtained from human tissues of dissection and normal groups. Wang et al. discovered the abnormally expressed miRNAs and lncRNAs within AD samples and identified the role of lncRNA OIP5-  in exacerbating injuries to all the three layer of aortic wall in AD genesis and development by increasing TUB expression through the sponge of miR-143-3p [34]. As discovered by another study, lncRNA PVT1 level increased, yet miR-27b-3p level decreased within AD tissue. They reduced the lncRNA PVT1 level to suppress the migration, phenotype switch and viability of human aortic smooth muscle cells treated with growth factor-BB (PDGF-BB) through targeting miR-27b-3p [35]. For understanding the pathophysiological and biological mechanisms related to AD genesis and development, cir-cRNAs, miRNAs, and lncRNAs with abnormal expression should be discovered and verified in relevant human AD samples or relevant animal models. Extraction of total RNA from critical types of cells, like immune cells, ECs, and SMCs, can shed more lights on alterations specific to cell type in the process of AD progression. When the differential expression is verified, it is necessary to test those mechanisms of ncRNAs in regulating AD progression in vitro and in vivo at molecular level. Certain limitations should be noted in the present work. At first, although it has been shown that E2F2, IGF1R, and BDNF can cause arterial-related diseases, there is no research to prove the relationship between aortic dissection and these target gene. Secondly, the way of action has been predicted based on the measured RNA network, however, which has not been confirmed (dual luciferase reporter gene analysis, gene overexpression or gene knockout). Although several related genes have been screened out in the present study for the first time, further in vitro clinical research and in vivo experiments should be carried out to confirm its expression and functional mechanism in terms of AD.
Currently, little research is conducted to explore the lncRNA mechanism in AD. The present work has its certain strengths. For instance, it is the first study to construct the lncRNA-miRNA-mRNA network based on GEO database. Nonetheless, our findings were just obtained from bioinformatics analysis. Therefore, it is of necessity to conduct a thorough study for verifying the potential effects of those 7 axes within AD. In summary, we confirmed that the ceRNA networks, including the regulated networks, lncRNA FAM87A-hsa-miR-31-5p/ hsa-miR-7-5p-E2F2, lncRNA C9orf106-hsa-miR-7-5p-IGF1R, and lncRNA UCA1-hsa-miR-107-BDNF, might be associated with the pathogenesis of and development of AD. Moreover, our study shed novel lights on the AD pathogenic mechanism. Nevertheless, ceRNA networks and their associations with AD should be validated.

Conclusions
In summary, the ceRNA interaction axis we identified is a potentially critical target for treating AD. According to the difference significance, PPI interaction correlation and regulation relationship, E2F2 axis (LncFAM87Ahsa-miR-7-5p/hsa-miR-31-5p-E2F2) may play a key role in the treatment of AD. Our results shed more lights on the possible pathogenic mechanism in AD using a lncRNA-associated ceRNA network.