Construction and Comprehensive Analysis of Dysregulated Long Noncoding RNA-Associated Competing Endogenous RNA Network in Moyamoya Disease

Background Moyamoya disease (MMD) is a rare cerebrovascular disease characterized by chronic progressive stenosis or occlusion of the bilateral internal carotid artery (ICA), the anterior cerebral artery (ACA), and the middle cerebral artery (MCA). MMD is secondary to the formation of an abnormal vascular network at the base of the skull. However, the etiology and pathogenesis of MMD remain poorly understood. Methods A competing endogenous RNA (ceRNA) network was constructed by analyzing sample-matched messenger RNA (mRNA), long non-coding RNA (lncRNA), and microRNA (miRNA) expression profiles from MMD patients and control samples. Then, a protein-protein interaction (PPI) network was constructed to identify crucial genes associated with MMD. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes pathway (KEGG) enrichment analyses were employed with the DAVID database to investigate the underlying functions of differentially expressed mRNAs (DEmRNAs) involved in the ceRNA network. CMap was used to identify potential small drug molecules. Results A total of 94 miRNAs, 3649 lncRNAs, and 2294 mRNAs were differentially expressed between MMD patients and control samples. A synergistic ceRNA lncRNA-miRNA-mRNA regulatory network was constructed. Core regulatory miRNAs (miR-107 and miR-423-5p) and key mRNAs (STAT5B, FOSL2, CEBPB, and CXCL16) involved in the ceRNA network were identified. GO and KEGG analyses indicated that the DEmRNAs were involved in the regulation of the immune system and inflammation in MMD. Finally, two potential small molecule drugs, CAY-10415 and indirubin, were identified by CMap as candidate drugs for treating MMD. Conclusions The present study used bioinformatics analysis of candidate RNAs to identify a series of clearly altered miRNAs, lncRNAs, and mRNAs involved in MMD. Furthermore, a ceRNA lncRNA-miRNA-mRNA regulatory network was constructed, which provides insights into the novel molecular pathogenesis of MMD, thus giving promising clues for clinical therapy.


Introduction
Moyamoya disease (MMD) is a rare cerebrovascular disease characterized by chronic progressive occlusion or stenosis of the bilateral internal carotid artery (ICA), the anterior cerebral artery (ACA), and the middle cerebral artery (MCA) [1,2]. MMD is secondary to the formation of an abnormal vascular network at the base of the skull. Because the abnormal vascular network of the skull base looks like "smoke" on cerebral angiography images, it is called "moyamoya disease" [3]. The MMD incidence rate in Eastern Asian countries is higher [4], and it mainly occurs in children and young adults, peaking at the ages of 5 to 9 and 35 to 45 years [5]. MMD can seriously affect the mental and physical health of patients. However, the etiology and pathogenesis of MMD remain poorly understood; it may be related to genetics, inflammation, immune response, and environmental factors [6][7][8][9][10][11].
Moreover, with the development of microarray and sequencing technology, investigators have begun to explore factors other than direct disease-causing genes, including noncoding RNAs (ncRNAs). Gao et al. revealed the expression profile of lncRNAs and mRNAs in MMD patients in 2016 [9], and Dai et al. analyzed miRNAs in the serum of MMD patients and healthy controls in 2012 [32]. miRNAs can posttranscriptionally regulate gene expression by binding to MREs (miRNA-response elements) of their target transcript. mRNAs, lncRNAs, and other RNA transcripts could act as endogenous miRNA sponges to inhibit miRNA function. These interactions illustrate the famous ceRNA hypothesis presented by Salmena in 2011 [33], which gave us a new "language" in different types of RNA transcripts. After that, the ceRNA hypothesis was applied to many fields [34]. The Linc2GO database was constructed by Liu et al. in 2013 [35]. StarBase v2.0 was published by Li et al. to predict miRNA-ceRNA interactions [36]. Moreover, continued analysis of ceRNA networks would deepen our knowledge about how different subtypes of noncoding RNAs work with each other.
In this study, a comprehensive analysis of the miRNA, mRNA, and lncRNA expression profiles in MMD was done, and then, we constructed MMD-specific ceRNA networks using a large cohort from an online database. As far as we know, this is the first study to establish a ceRNA lncRNA-miRNA-mRNA network in MMD, which provides novel insight into the molecular pathogenesis of MMD, thus giving promising clues for clinical therapy. In addition, core regulatory miRNAs (miR-107 and miR-423-5p) and key mRNAs (STAT5B, FOSL2, CEBPB, and CXCL16) were enriched in immune system/inflammation biological processes, indicating their potential role in MMD.

Materials and Methods
2.1. Data Collection. miRNA microarray data were downloaded from Gene Expression Omnibus (GEO, https://www .ncbi.nlm.nih.gov/geo/) in NCBI (The National Center for Biotechnology Information). GEO is an unrestricted open access repository that provides high-throughput microarray and next-generation sequence datasets that have been submitted by researchers around the world. GSE45737 is a miRNA expression profile of the serum from 10 MMD patients and 10 normal healthy controls [32]. The lncRNA and mRNA expression profiles in blood samples from 15 MMD patients and 10 healthy controls were kindly provided by a collaborating academician, Zhao [9].

Identification of Differentially Expressed RNAs in MMD
Patients Compared to Healthy Controls. R software with packages ggplot2, edgeR, and pheatmap (http://bioconductor.org/ bioclite. R) was adopted to identify differentially expressed RNAs (DERs). In brief, datasets were standardized after conversion of formats, variance normalization, and the addition of missing values as well as statistical testing of differentially expressed probes. The expression levels of all targets, including mRNA, miRNA, and lncRNA, within the datasets were subjected to analysis with R. The threshold was set as a P value < 0.05 and |log 2 FC | >1. According to these criteria, DERs were identified for further analysis.

Gene Ontology and Pathway Enrichment
Analyses. The Database of Annotation, Visualization and Integrated Discovery (DAVID, http://david.ncifcrf.gov) is a public database with comprehensive online tools for functional annotation. Kyoto Encyclopedia of Genes and Genomes (KEGG) is a collection of databases that contain information about genomes, biological pathways, diseases, and chemical substances [37]. Gene Ontology (GO) is an international standardized gene functional classification system that offers a dynamically updated controlled vocabulary and a strictly defined concept to comprehensively describe properties of genes and their products in any organism. GO has three ontologies: molecular function, cellular component, and biological process [38].
In the present study, GO and KEGG pathway enrichment analyses were performed using DAVID. P < 0:05 was considered statistically significant.  [36]. The unique algorithm of ENCORI enables all obtained interactions to be confirmed by at least one other major RNA-RNA prediction website, such as miRanda, Pic-Tar, or TargetScan. In addition to sequence matching, the prediction was approved by multidimensional sequencing data. All these features make ENCORI a reliable source for predicting RNA-RNA interaction, especially the miRNA-mRNA interaction. Two other databases, miRcode (http://www.mircode.org) and DIANA (http://carolina.imis .athena-innovation.gr), were applied in the study for predicting miRNA-lncRNA interactions. Afterwards, all interactions were input into Cytoscape (version 3.7.2, http:// cytoscape.org) to visualize ceRNA regulatory networks. The flow chart can be seen in Figure 1. 2.5. Protein-Protein Interaction (PPI) Network. All DEGs were imported into STRING 10.5, which is a search tool used to identify gene interactions (https://string-db.org/). The PPIs were used to construct a network, which was visualized by using Cytoscape software 3.6 (http://www.cytoscape.org). The color of edges in the network indicate protein-protein associations: light blue and purple indicate known interactions from curated database and experimentally determined, respectively; dark green/red/dark blue indicate predicted interactions by gene neighborhood/gene fusions/gene cooccurrence, respectively; and light green/black/blue indicate text mining/coexpression/protein homology.

Gene Expression Signature Analysis with a Connectivity
Map. The DEGs were used to perform gene expression signature analysis with connectivity maps (CMap, clue.io). The upregulated and downregulated genes were used as tags, changed into probe IDs referred to Affymetrix U133 Gene-Chip and uploaded into the CMap database to calculate their values from other drug-target datasets. According to the similarity of gene expression profiles, pairs of gene expression signatures and targeted drugs were used to obtain a value. If the value was a positive number, the target drug would have an effect that was similar to that of the MMD-induced gene expression signature. If the value was a negative number, the targeted drug would have an effect that was opposite that of the MMD-induced gene expression signature; namely, the targeted drug might have an effect that could be useful in treatment.

Statistical Analysis.
We used SPSS 11.0 (SPSS, Chicago, IL) to analyze the dataset from the microarray experiments. All data are represented as the mean ± SD. Statistical significance was determined at P < 0:05.

Result and Discussion
3.1. Differentially Expressed mRNAs, miRNAs, and lncRNAs between MMD Patients and Healthy Controls. After differential expression analysis, a total of 2294 DEmRNAs were screened between MMD patients and healthy controls, 865 of which were downregulated and 1429 of which were upregulated in MMD patients. (Table S1, Figure 2(a)). Several genes reported in previous studies in MMD, such as HIF1α (log 2 FC = 1:214), SORT1 (log 2 FC = 1:628), and MMP9 (log 2 FC = 2:40), are marked in Figure 2. HIF1α was found to be overexpressed in the intima of the MCA of MMD patients. HIF1α is a master transcriptional regulator of the adaptive response to hypoxia. Under hypoxic conditions, HIF1α translocates to the nucleus, where the HIF1 complex (HIFα/HIFβ) binds to the hypoxia-response element and activates the expression of many genes that can increase oxygen delivery and respond to oxygen deprivation in MMD [7]. MMP9 belongs to a family of zinc-binding proteolytic enzymes that are capable of degrading all the components of the extracellular matrix in a variety of physiologic and pathophysiological conditions. Fujimura et al. inferred that the higher expression of MMP9 in MMD patients may play an integrated role in physiologic and pathologic angiogenesis and to the instability of the cerebral vascular structure [39]. SORT1 is another gene reported to be associated with MMD. Increased expression of SORT1 inhibited endothelial cell tube formation and regulated major angiogenic factors and MMP9 expression, implying that SORT1 participated in the pathogenesis of MMD [28].
In addition, 94 DEmiRNAs and 3649 DElncRNAs from GEO datasets were identified. Representative DERs are shown in Figure S1 (a-d).

Construction of a Competing
Endogenous RNA Regulatory Network. The ENCORI database was employed to screen potential interactions between DERs. A synergistic, competitive module of the ceRNA network was constructed separately according to upregulated or downregulated DEmRNAs, which contained 84 nodes in the upregulated group and 66 nodes in the downregulated group. In addition, there were 68 mRNA-miRNA interactions and 16 lncRNA-miRNA interactions in the upregulated group (Figure 3(a)). In the downregulated group, there were 61 interactions between mRNAs and miRNAs and 35 interactions between lncRNAs and miRNAs ( Figure 3(b)). The ceRNA network was generated using Cytoscape, as previously discussed.
Based on the network organization, we found that miR-107 competed with 16 mRNAs and 4 lncRNAs (LINC02434, AL589642.1, AC003092.1, and AL035425.3) in the module (Figure 3(b)). A previous study showed that miR-107 is upregulated in response to low-oxygen conditions [40]. Subsequently, miR-107 was found to be abnormally expressed in several cancers, such as PDAC. When miR-107 expression was downregulated in PDAC, cell migration and invasion were inhibited, implying the important role of miR-107 in tumor cell activity [41]. Furthermore, they found that the expression of caveolin-1 was upregulated by a miR-107 inhibitor. Caveolin-1 was reported to be associated with negative remodeling in MMD through the inhibition of angiogenesis in endothelial cells and the induction of apoptosis in VSMCs [30,31]. Another study by Meng et al. found that miR-107 can inhibit endothelial progenitor cell (EPC) differentiation via HIF1β [42]. HIF1β is another subunit of HIF1 that generally heterodimerizes with HIF1α. Together, they play key roles during hypoxic conditions, which are similar to the conditions in MMD: low oxygen because of vascular occlusion. EPCs can differentiate into mature endothelial cells and play important roles in the recovery of endothelial function and tissue repair. The role of EPCs reflects the mixed state of vascular obstruction and abnormal angiogenesis in the pathogenesis of MMD [43]. The ceRNA network near miR-107 reveals that FoxC1 is one of the potential 3 Computational and Mathematical Methods in Medicine downstream target genes, and it is necessary in the process of vascular development, involving arterial specification and lymphatic sprouting. Abnormal expression of FoxC1 leads to unusual angiogenesis in many tissues [44,45].
In addition, miR-423-5p competed with 36 mRNAs (CXCL16, FOSL2, etc.) and 4 lncRNAs (NEAT1, HCG18, AL137145.2, and LINC00963) in the module (Figure 3(a)). miR-423-5p was reported to play important roles in the inhibition of the cell proliferation and invasion of cancer cells such as colon cancer and ovarian carcinoma [46,47]. Therefore, the downregulation of miR-423-5p in MMD patients may increase the proliferation of vascular smooth muscle cells, which is one likely reason for vessel occlusion. In addition, numerous studies focusing on NEAT1's role in cancer biology found that this lncRNA plays a crucial role in carcinogenesis [48]. NEAT1 mainly works as a ceRNA by sponging antitumor miRNAs [49]. NEAT1 is also involved in immune system responses, viral diseases, and     Computational and Mathematical Methods in Medicine neurodegeneration disorders [50]. To study FOSL2, also named Fra 2, Maurer et al. created Fra 2 knockout mice and found that the mice developed pulmonary arterial occlusion due to vascular SMC proliferation and inflammation and pulmonary fibrosis [51,52]. All of the above results imply that the ceRNA lncRNA-miRNA-mRNA regulatory network we constructed provides many new clues regarding MMD pathogenesis.

Functional
Annotation of the mRNAs Involved in the ceRNA Network. After the ceRNA network was established with the help of the DAVID database, functional annotation and pathway analysis of this small group of DEmRNAs were performed to identify potential candidate pathways or biological processes related to MMD. As shown in Figure 4, some of pathways require our attention, and processes related to the immune response and inflammatory reaction, including immune system process, T cell aggregation, T cell activation, lymphocyte aggregation, and lymphocyte activation, were significantly enriched. Additionally, another enrichment also occurred in biological processes associated with cell development and differentiation, including paraxial mesoderm development and mesenchymal cell differentiation; these results suggest important roles for these biological activities in MMD. The Kyoto Encyclopedia of Genes and Genomes showed that DEmRNAs were enriched in chemokine signaling, ErbB signaling, axon guidance, and vascular smooth muscle contraction (Table 1).
Recently, many studies have shown that immunological/inflammatory factors are involved in the occurrence and development of MMD. According to IHC staining, there were T cells and macrophages infiltrating in the stenosed and thickened vascular intima of MMD patients [53]. The abnormal deposition of IgG in the elastic layer of the ICA and MCA suggests that the infiltration of immune cells and the damage to the immune functions are related to MMD [54]. Moreover, the overexpression of inflammatory factors in MMD patients, such as MCP-1, IL-1β, and SDF-1α, suggests that inflammation may also affect the progression of MMD [55]. Consistently, in this study, several mRNAs that encode critical inflammatory molecules, such as chemokines and cytokines, were dysregulated and were determined to be DEmRNAs in MMD patients. Nevertheless, although varied mRNAs were clearly enriched in terms of GO analysis, there were few found in the ceRNA network. However, several important genes involved in the regulation of inflammation in MMD were modulated by ceRNAs. CXCL16 is considered to be an important pathogenic mediator of atherosclerosis (clinical severity is graded according to the severity of carotid stenosis) [56]. CXCL16 is a vascular-derived factor that induces angiogenesis [57]. CXCL16 also exists in a soluble form and interacts with its specific chemokine receptor, CXCR6, to recruit the migration of activated T cells into the inflammatory tissue [58]. As shown in Figure 3(b), four potential lncRNAs, including LINC00963, NEAT1, HCG18, and AL137145.2, could act as ceRNAs to regulate CXCL16 through miR-107. The work on this interesting ceRNA network remains to be done in the future. Figure 5, a PPI network for DEmRNA-involved ceRNA  Figure 5: The protein-protein interaction network. The red and green nodes represent the regulated and downregulated genes, respectively. (For interpretation of the references to color in this figure legend, the reader is referred to the online version of this article.)    Figure 7: The potential mechanism of a DElncRNA sponging miR-423-5p. 9 Computational and Mathematical Methods in Medicine analyzed by CMap to identify small molecule drugs. Strong negative correlations were found between MMD and enzastaurin, cyproheptadine, flupirtine, indirubin, and mitoglitazone (CAY-10415); strong positive correlations were found between MMD and flavokavain-b, CGS-20625, vinburnine, apicidin, and cytochalasin-d (Table S5). The drugs that had a strong negative correlation with the pathogenesis of MMD might have therapeutic effects on MMD (Table 2). CAY-10415 and indirubin gained our attention. The structures of the two potential molecular drugs were investigated using the PubChem database ( Figure 6). CAY-10415 is a member of a new class of compounds that modulate mitochondrial pyruvate carrier (MPC), a key controller of cellular metabolism that influences mTOR activation [59]. It is commonly known that CAY-10415 can be used as an insulin sensitizer, and it can play this role without activating PPARᵧ. Therefore, CAY-10415 can avoid negative side effects observed in currently used insulin sensitizers, such as pioglitazone and rosiglitazone. CAY-10415 has been used in Alzheimer's disease patients [60]. It is generally accepted that insulin sensitizers can not only improve diabetes but also improve blood lipid disorders, reduce the level of free fatty acids in plasma, reduce the effect of fat toxicity, and indirectly protect the function of β cells [61]. By inhibiting the proliferation and migration of vascular smooth muscle cells and reducing the intima-media thickness of arteries, it can play a protective role in the intima. Likewise, indirubin, a red isomer of indigo, is the active ingredient of the traditional Chinese drug Danggui Longhui Wan, which was used for the treatment of chronic myelocytic leukemia (CML) [62]. Enzyme-based in vitro studies have observed that indirubin and its derivatives, such as indirubin-3′-monoxime, indirubin-5-sulfonate, and indirubin-3′-monoxime-5-sulphonic acid, are potential inhibitors of CDKs [63]. Furthermore, different indirubin derivatives showed antiangiogenesis activity by blocking VSMC proliferation and endothelial cell function through the inhibition of the STAT signaling pathway and reduction of neointima formation in vivo [64]. All of the above findings suggest that CAY-10415 and indirubin may be used in MMD patients to avoid vascular aberration and occlusion.

Conclusions
In summary, using bioinformatics analysis of candidate RNAs, the present study identified a series of clearly altered lncRNAs, miRNAs, and mRNAs involved in MMD. Furthermore, a ceRNA lncRNA-miRNA-mRNA regulatory network was constructed, which provides a novel insight into the molecular pathogenesis of MMD, thus giving promising clues for clinical therapy. In addition, core regulatory miR-NAs (miR-107 and miR-423-5p) and key mRNAs (STAT5B, FOSL2, CEBPB, and CXCL16) were enriched in immune system/inflammation biological processes, indicating their potential role in MMD (Figure 7). In the future, more attention should be paid to the validation of competing endogenous RNA interactions with experimental techniques.
Finally, two potential small molecule drugs, CAY-10415 and indirubin, were identified by CMap to be candidate drugs for treating MMD.

Data Availability
(1) The miRNA microarray data used to support the findings of this study have been deposited in the Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/) in NCBI (GSE45737). (2) The lncRNA and mRNA microarray data included in this study are available upon request by contact with the corresponding author. The data were kindly provided by a collaborating Prof. Zhao.

Conflicts of Interest
The authors declare no conflicts of interest.