Construction of competitive endogenous RNA network reveals regulatory role of long non-coding RNAs in intracranial aneurysm

Background Rupture of intracranial aneurysm (IA) is the main cause of devastating subarachnoid hemorrhage, which urges our understanding of the pathogenesis and regulatory mechanisms of IA. However, the regulatory roles of long non-coding RNAs (lncRNAs) in IA is less known. Results We processed the raw SRR files of 12 superficial temporal artery (STA) samples and 6 IA samples to count files. Then the differentially expressed (DE) mRNAs, miRNAs, and lncRNAs between STAs and IAs were identified. The enrichment analyses were performed using DEmRNAs. Next, a lncRNA-miRNA-mRNA regulatory network was constructed using integrated bioinformatics analysis. In summary, 341 DElncRNAs, 234 DEmiRNAs, and 2914 DEmRNAs between the STA and IA. The lncRNA-miRNA-mRNA regulatory network of IA contains 91 nodes and 146 edges. The subnetwork of hub lncRNA PVT1 was extracted. The expression level of PVT1 was positively correlated with a majority of the mRNAs in its subnetwork. Moreover, we found that several mRNAs (CCND1, HIF1A, E2F1, CDKN1A, VEGFA, COL1A1 and COL5A2) in the PVT1 subnetwork served as essential components in the PI3K-Akt signaling pathway, and that some of the non-coding RNAs (ncRNAs) (PVT1, HOTAIR, hsa-miR-17, hsa-miR-142, hsa-miR-383 and hsa-miR-193b) interacted with these mRNAs. Conclusion Our annotations noting ncRNA’s role in the pathway may uncover novel regulatory mechanisms of ncRNAs and mRNAs in IA. These findings provide significant insights into the lncRNA regulatory network in IA. Supplementary Information The online version contains supplementary material available at 10.1186/s12868-021-00622-7.


Background
Rupture of intracranial aneurysm (IA) is the main cause of subarachnoid hemorrhage (SAH), leading to exceedingly high mortality and morbidity [1,2]. The prevalence of IA in the general population is reported to be 3.2% [1]. Approximately 85% of non-traumatic SAH is caused by IA rupture [3]. Moreover, although unruptured IAs (UIs) are typically asymptomatic, UIs are being detected more and more, and have become an important health issue [1]. Since surgical clipping and intravascular coiling are invasive procedures with potentially serious complications, the effective management of UIs remains a challenge [4]. A deeper understanding of the pathogenesis of IA helps to find more effective treatment for IA.
Recently, an increasing number of studies have focused on the competitive endogenous RNA (ceRNA) network's impact in cardiac hypertrophy, type 2 diabetes mellitus, and various types of cancer [5][6][7][8]. The hypothesis of ceRNA described a molecular regulatory mechanism for post-transcriptional regulation. Several studies revealed that lncRNAs can function as microRNA (miRNA) sponges in a ceRNA regulatory network, and further regulate mRNA expression levels [9]. LncRNAs are non-coding RNA of greater than 200 nucleotides in length, which were reported to play a crucial role in the biological processes [10]. However, the potential role of the ceRNA regulatory network in IA formation remains unclear.
In this study, we identified differentially expressed lncRNAs, miRNAs, and mRNAs between STAs and IAs. Then, we developed a lncRNA-miRNA-mRNA regulatory network. Previous study found that knockdown of lncRNA PVT1 could reverse the murine angiotensin II-induced abdominal aortic aneurysm (AAA) associated alterations, including inhibition of vascular smooth muscle cell (VSMC) apoptosis and ECM disruption [11]. However, the regulatory mechanism of lncRNA PVT1 in development of IA was not clear. Thus, lncRNA PVT1 was selected for further analysis. The correlations between PVT1 and mRNAs in the PVT1 subnetwork were analyzed. After a functional enrichment analysis, we found that mRNAs in the PVT1 subnetwork were annotated into the PI3K-Akt signaling pathway. In addition, we found that several ncRNAs were connected with these mRNAs, and we annotated these ncRNAs to the PI3K-Akt signaling pathway, which may uncover novel regulatory mechanisms of ncRNAs and mRNAs in IA.

Patients and RNA-seq data processing
The SRR files (SRR1819905-SRR1819922), including 12 superficial temporal artery (STA) samples and 6 intracranial aneurysm (IA) samples, were retrieved from GEO NCBI [12] for single ended data, which were converted to 'fastq' format data using sratoolkit (version 2.8.2). The reference genome of human (GRCh38.p13) was downloaded from Ensembl (https ://asia.ensem bl.org/). STAR (version 2.4.0i) [13], a highly efficient RNAseq alignment tool. This was used to align the 'fastq' format data with the human reference genome, using the custom default parameters. The SAM files from STAR alignment were converted to the BAM format by using SAMtools (version 1.3.1) [14]. The BAM files were converted to counts files by using HTSeq [15].We generated ".csv" files, each consisting of all gene counts for that particular sample. Furthermore, we combined all of the ".csv" files and obtained a single file with sample names depicted as columns, and gene names depicted as rows. Then, the counts files were normalized using the edgER package [16]. LncRNAs and mRNAs were annotated using the Ensembl database [16]. Regarding the microRNA data in GSE66240, the raw probe-level data was preprocessed through the robust multi-array average (RMA) algorithm in the Affy package [17]. For genes corresponding to multiple probes, we used the average probe value as the expression level [18]. The missing data in these gene expression matrices were imputed with the k-Nearest Neighbor (KNN) method (k = 10) [19]. As all the data were retrieved from the GEO database, the approval from the local Ethics Committee was not needed.

Functional enrichment analysis
Both the Gene ontology (GO) enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were performed using the Database for Annotation, Visualization, and Integrated Discovery platform (DAVID 6.8, https ://david .ncifc rf.gov/) [20]. DEmRNAs were utilized for GO and KEGG pathway enrichment analyses. The pathways and GO terms, with corrected P-values < 0.05 using the Benjamini method, were considered significant categories. Figure 1 depicts the flow chart of the ceRNA network construction. Firstly, the lncRNA-miRNA potential interaction file (Highly conserved microRNA families) was downloaded from miRcode [21]. MiRcode provides a whole transcriptome human microRNA target prediction based on comprehensive GENCODE gene annotation, including > 10,000 long non-coding RNA genes. Secondly, DElncRNAs were put into the miRcode database and lncRNA-miRNA pairs were identified. Then, we eliminated miRNAs that expressed no difference between the STA and IA tissues. Thirdly, miRNA-targeted mRNAs were identified using three miRNA reference databases: miRDB, miRTarBase, and TargetScan, which contain predictive or experimental validated miRNAtargeted mRNAs [22][23][24]. Only mRNAs appearing in the three databases were defined as miRNA-targeted mRNAs for increasing the prediction reliability. Furthermore, targeted mRNAs that expressed no difference between STA and IA tissues were filtered out. The ceRNA network was constructed and viewed using Cytoscape (http://www. cytos cape.org/). The relationship between PVT1 and mRNA expression level was analyzed with Pearson correlation. The heatmaps based on DEmRNAs, DEmiRNAs and DElncRNAs included in the ceRNA regulatory network were generated by "pheatmap" package.

Statistical analysis
GraphPad Prism (version 6.0, GraphPad Software, San Diego, CA, USA) and R language (3.4.0) were used for statistical analysis. A P-value < 0.05 was considered statistically significant.

Functional enrichment analysis of DEmRNAs
GO enrichment and KEGG pathway analyses were performed to explore the potential functions of the 2.914 differentially expressed mRNAs. In the GO enrichment analysis, a total of 378 enriched GO terms have been identified in the Biological Process (BP). And the top 10 significantly enriched terms are shown in Fig. 2a. Moreover, the DEmRNAs were primarily enriched in collagen and extracellular matrix-related BPs, such as "extracellular matrix organization", "collagen catabolic process", "collagen fibril organization", and "cell adhesion". Moreover, we also found that these DEmRNAs were also enriched in "angiogenesis", "skeletal system development", "muscle contraction", and "inflammatory response", which indicates that these BPs may play roles in IA formation. In addition, a total 53 enriched pathways were identified after the KEGG pathway analysis. The top 10 enriched pathways are shown in Fig. 2b. And "vascular smooth muscle contraction" is closely associated with IA formation in these pathways. Furthermore, the enriched "focal adhesion", "ECM-receptor interaction", and "PI3K-Akt signaling pathway" indicated that the PI3K-Akt signaling pathway may play an important role in IA, which was further analyzed in this study.

Construction of a ceRNA regulatory network in IA
A lncRNA-miRNA-mRNA regulatory network of IA was constructed (Fig. 3) for exploring the regulatory mechanism of IA. The mRNA-miRNA and lncRNA-miRNA relationship pairs (Additional file 1: Table S1, Additional file 2: Table S2) were combined into the ceRNA network, following a negative regulation pattern. The positively co-expressed mRNA-miRNA and lncRNA-miRNA pairs were excluded. Finally, the network was constructed with 91 nodes (60 mRNAs, 9 miRNAs and 22 lncRNAs) and 146 edges. The degree of a node is defined as the total number of edges connecting it in the network, with an increased degree indicating that the node is highly connected, and likely to be a hub gene. The heatmaps based on DEmRNAs, DEmiRNAs and DElncRNAs included in the ceRNA regulatory network were shown in Fig. 4. Among the DEmiRNAs in the network, hsa-miR-17 has the highest degree (degree = 52), indicating that 52 nodes directly connect to it (Fig. 5). Moreover, other DEmiR-NAs that connect lncRNAs and mRNAs in the network are also shown (Fig. 5). The highest degree DElncRNAs were PVT1, NEAT1 and KCNQ1OT1 (degree = 8).

The components of PVT1 subnetwork
PVT1 was the core lncRNA, and had the highest degree in the network. A previous study found that knockdown of lncRNA PVT1 reversed the murine angiotensin IIinduced AAA associated alterations, such as attenuation of aortic diameter dilation and marked adventitial thickening [11]. To explore the regulatory mechanism of PVT1 in IA, we focused on the PVT1 subnetwork. The components of the PVT1 subnetwork are shown in Fig. 6a. The predicted or validated binding sites of hsa-miR-17, hsa-miR-143, hsa-miR-152, hsa-miR-761 and hsa-miR-383 with PVT1 are shown (Fig. 6b-e).

Interaction between lncRNA and mRNA in the ceRNA network
LncRNA can interact with mRNAs indirectly following post-transcriptional regulatory mechanisms. We analyzed the correlation between PVT1 and mRNA's expression levels in its subnetwork, which showed that there were strong positive correlations. There were 52 mRNAs interacting indirectly with PVT1. The correlation coefficients between 38.5% (20/52) mRNAs and PVT1 were more than 0.6 (Additional file 3: Table S3), and the correlation coefficients between 90.4% (47/52) mRNAs and PVT1 were more than 0.3 (Additional file 3: Table S3). Scatter plots of the top 10 mRNAs with the highest correlation coefficients are shown in Fig. 7.

Discussion
LncRNA are a type of non-coding RNA with more than 200 bp transcript length. Various studies have found that lncRNAs play important roles in different types of diseases. Large amounts of lncRNAs have been annotated,  mellitus, and various types of cancer [5][6][7][8]. However, the potential role of the ceRNA network in IA formation remains unclear.
In this study, all of the DEmiRNAs in the ceRNA network were downregulated in IAs relative to STAs. According to ceRNA theory, the down-regulated lncRNAs and mRNAs were excluded. Previous studies found that downregulation of miR-143/145 cluster could upregulate their targeted gene KLF5, and further promote proliferation and migration of VSMCs [29]. In our study, we found that miR-143 interacted with lncRNA PVT1 and HOTAIR, and that miR-143 could target COL1A1 and COL5A2 and further affect the PI3K-Akt signaling pathway. Thus, we may provide a potential axis PVT1/HOTAIR-miR-143-COL1A1/ COL5A2 as an upstream regulator of the PI3K-Akt signaling pathway.
The PI3K-Akt signaling pathway is activated by various types of cellular stimuli or toxic insults, and regulates fundamental cellular functions such as transcription, translation, proliferation, growth, and survival. Dysregulation in the PI3K-Akt signaling pathway is implicated in various diseases such as cancer, type 2 diabetes mellitus, and IA [7,30,31]. Previous research reported that activation of PI3K-Akt signaling pathway promoted proliferation of vascular smooth muscle cells which were found to be involved in aortic aneurysms [32]. Inhibition of VEGF/PI3K/Akt signaling pathway mediated by miR-195 could suppress formation of abdominal aortic aneurysm [33]. In addition, the PI3K-Akt signaling pathway was less studied in IA formation. Sun et al. found that bone marrow mesenchymal stem cells-derived exosomes supressed the activation of the PI3k/Akt/NF-κB signaling pathway and maintained Th17/Treg balance, which could inhibite development of IA. Recently, various studies revealed that various lncRNAs affect the PI3K-Akt signaling pathway by targeting miRNAs, and further regulating mRNAs in the pathway [34,35].
In our study, we found that PVT1 and HOTAIR were involved in the regulation of PI3K-Akt signaling pathway by targeting different miRNAs and regulating upstream/downstream effectors (COL1A1, COL5A2, VEGFA, CCND1, HIF1A, E2F1 and CDKN1A). Several members of collagen family, such as COL1A1, COL5A2, were overexpressed in IAs, which were consistent with previous studies [36,37]. And these two genes were involved in collagen formation that could result in the ECM remodeling of intracranial blood vessels, which promotes formation of IA [36]. VEGFA, related to aSAH, was reported to be associated with thrombospondin and vascular endothelial growth factor [38]. A previous study revealed that HIF1A may be associated with rupture of human saccular IA wall [39]. Another research found that HIF1A was pivotal for the development of AAA [40]. HIF1A could upregulate expression levels of MMP-2 and MMP-9, and further promote aneurysmal progression [40]. Downregulation of PVT1 could reduce proinflammatory cytokines, MMP-2 and MMP-9, increase TIMP-1, and further reverse angiotensin II-induced AAA-associated alterations in mice [11]. In this study, we annotated these ncRNAs to PI3K-Akt signaling pathway. These findings establish novel connections among lncRNAs, miRNAs, and mRNAs in PI3K-Akt signaling pathway, which play an important role in IA.
In this study, we also found that there were strong correlations between PVT1 expression and several genes' expression levels ( Fig. 7 and Additional file 3: Table S3), including DNMT1, E2F1 et al. Jin et al. reported that PVT1 could recruit DNMT1 via EZH2 to miR-18b-5p DNA promoter and inhibited the miR-18b-5p transcription via DNA methylation [41]. Previous study found that PVT1 could regulate E2F1 expression levels in tumor development [42]. These genes might also play roles in development of IA, which needs more verification in future research. The potential axis PVT1/HOTAIR-miR-143-COL1A1/ COL5A2 needs further verification in future research, and the function of the axis in formation of IA needs to be explored.

Conclusion
First, we identified differentially expressed lncRNAs, miRNAs and mRNAs between STAs and IAs. Functional enrichment analysis of differentially expressed mRNAs was performed. We further developed a lncRNA-miRNA-mRNA regulatory network. Then, the subnetworks of each core miRNAs were extracted and analyzed. Moreover, PVT1, the hub lncRNA in the ceRNA network, was reported to be related to aneurysm. The subnetwork of PVT1 was extracted and analyzed. Furthermore, the correlations of expression levels between PVT1 and mRNAs in its subnetwork were analyzed. In addition, we found that several mRNAs in the PVT1 subnetwork, including CCND1, HIF1A, E2F1, CDKN1A, VEGFA, COL1A1, and COL5A2, were an important part of the PI3K-Akt signaling pathway. We found that several ncRNAs were connected with these mRNAs. Then, we annotated these ncRNAs to the PI3K-Akt signaling pathway. This study may provide a comprehensive view of the underlying mechanisms of gene regulation and interaction in IA. Additionally, this study also revealed that ncR-NAs involved in the PI3K-Akt signaling pathway could play an important role in IA formation. Fig. 9 The ncRNAs involved the PI3K-Akt signaling pathway. The nodes labeled in red triangles were mRNAs in ceRNA network. The lncRNAs were labeled in green, and miRNA were labeled in red. The ncRNAs were annotated to PI3K-Akt signaling pathway according to previous studies and the predicted interactions in ceRNA network. The red and green arrows represented interactions between these RNAs