Integrated Analysis of LncRNA-mRNA Co-Expression Profiles in Patients with Moyamoya Disease

Moyamoya disease (MMD) is an idiopathic disease associated with recurrent stroke. However, the pathogenesis of MMD remains unknown. Therefore, we performed long noncoding RNA (lncRNA) and messenger RNA (mRNA) expression profiles in blood samples from MMD patients (N = 15) and healthy controls (N = 10). A total of 880 differentially expressed lncRNAs (3649 probes) and 2624 differentially expressed mRNAs (2880 probes) were obtained from the microarrays of MMD patients and healthy controls (P < 0.05; Fold Change >2.0). Gene ontology (GO) and pathway analyses showed that upregulated mRNAs were enriched for inflammatory response, Toll-like receptor signaling pathway, chemokine signaling pathway and mitogen-activated protein kinase (MAPK) signaling pathway among others, while the downregulated mRNAs were enriched for neurological system process, digestion, drug metabolism, retinol metabolism and others. Our results showed that the integrated analysis of lncRNA-mRNA co-expression networks were linked to inflammatory response, Toll-like signaling pathway, cytokine-cytokine receptor interaction and MAPK signaling pathway. These findings may elucidate the pathogenesis of MMD, and the differentially expressed genes could provide clues to find key components in the MMD pathway.

Identification of differentially expressed lncRNAs and mRNAs. The list of lncRNAs and their expression profiles were extracted in our previous study 18 . In brief, we identified 880 differentially expressed lncRNAs (3649 probes) and 2624 mRNAs (2880 probes) from the microarrays of MMD patients and healthy controls (P < 0.05; Fold Change > 2.0) (Supplementary Table S1). Of those, 1746 upregulated mRNAs and 878 downregulated mRNAs were identified. A volcano plot was created and scatter analyses were conducted to identify differences among mRNAs (Fig. 1a,b). We further created a heat map of differentially expressed mRNAs (Fig. 2).
We also selected several differentially expressed genes randomly and further performed quantitative real-time polymerase chain reaction (qRT-PCR) to examine their expression levels. The qRT-PCR results suggest that the fold changes observed in the microarray analysis were robust (Fig. 3).
Examination of the function of differentially expressed mRNAs. Gene ontology (GO) and KEGG pathway analyses were conducted to explore the function of the 2624 differentially expressed mRNAs using DAVID (The Database for Annotation, Visualization and Integrated Discovery). The results showed that upregulated genes were enriched for inflammatory response, response to wounding and defense response, etc. (Fig. 4a), while the downregulated genes were enriched for neurological system process, digestion and positive regulation  of transcription from RNA polymerase II promoter, etc. (Fig. 4b). Moreover, the KEGG pathway analysis showed that the upregulated genes were enriched for Toll-like receptor signaling pathway, chemokine signaling pathway and MAPK signaling pathway, etc. (Fig. 4c), while the downregulated genes were enriched for drug metabolism, retinol metabolism and olfactory transduction, etc. (Fig. 4d).
LncRNA-mRNA co-expression networks. We performed lncRNA-mRNA co-expression network analysis including 3649 lncRNAs probes and 2880 mRNAs probes. Our results showed that the co-expression networks were linked to inflammatory response, Toll-like signaling pathway, cytokine-cytokine receptor interaction and MAPK signaling pathway (Fig. 5). Thirty-two lncRNAs interacted with 11 mRNAs in the GO term of inflammatory response (Fig. 5a), 26 lncRNAs interacted with 2 mRNAs in the Toll-like signaling pathway (Fig. 5b), 41 lncRNAs interacted with 6 mRNAs in the cytokine-cytokine receptor interaction (Fig. 5c), and 15 lncRNAs interacted with 6 mRNAs in the MAPK signaling pathway (Fig. 5d).

Discussion
MMD is a chronic occlusive cerebrovascular disease of unknown etiology 21 , and it is usually diagnosed by radiological findings, such as computed tomography (CT) perfusion and magnetic resonance imaging (MRI) 22 . It has been shown that MMD has a high prevalence in Asian countries, such as China, Japan and South Korea [23][24][25][26][27] . However, there are fewer studies focused on lncRNA-mRNA co-expression in MMD, and the molecular mechanisms behind MMDs remain poorly understood. It has been reported that lncRNAs play important roles in a wide range of functional activities 6,7 . The dysregulation of lncRNAs might contribute towards MMD. Therefore, the integrated analysis of the differentially expressed lncRNAs and mRNAs could help to reveal the pathogenesis of MMD. Many studies have associated single nucleotide polymorphisms (SNPs) of genes with MMD [28][29][30][31][32][33][34] . RNF213 and MMP3 were proposed to be susceptibility genes for MMD 28,30 . It was reported that RNF213 is associated with immune response and that it might act cooperatively with other molecules under inflammatory signals based on bioinformatics data. IFNG and TNFA synergistically activated transcription of RNF213 both in vitro and in vivo 35 .  However, IFNG and RNF213 were not co-expressed based on the mRNA microarray. The presence of a heterozygous genotype in TIMP2 promoter could be a genetic factor for familial MMD. Moreover, it was reported that TGFB1 could be involved in vascular growth and transformation processes and may play an important role in the development of MMD 34,36 . SNPs may affect the expression of genes for MMD. Based on the microarray database, RNF213 was not differentially expressed between MMDs and controls. However, TGFB1 and TIMP2 were found to be differentially expressed. These findings need to be validated in further studies including more samples.
The co-expression results were based on the expression of lncRNAs and mRNAs for MMD. It was reported that HIF1A could directly bind to the promoter of HOTAIR, which has been identified in a variety of carcinomas 37 . CXCR2 is involved in migration and activation of leukocytes and plays a key role in several inflammatory diseases 38,39 . MALAT1 could downregulate the expression of CXCR2 via miR-22-3p 40 . These were consistent with the results that HOTAIR and MALAT1 were highly associated with HIF1A and CXCR2, respectively (correlation coefficient − 0.76, 0.87; P < 0.01), based on our microarray. To our knowledge, this is the first study of lncRNA-mRNA co-expression network analysis for patients with MMD. Dai et al. have conducted a serum miRNA signature in MMD, and it identified 94 differential expressed miRNAs 41 . Non-coding RNAs and mRNAs compete for binding to miRNAs by sharing one or more miRNA response elements (MREs) to regulate gene  42 , and the differentially expressed lncRNAs and mRNAs in our study partly correlate with these miR-NAs. Therefore, future studies could focus on finding a competing endogenous RNAs (ceRNA) network with a key role in the pathogenesis of MMDs Next, we investigated whether mRNAs tend to be neighbors with lncRNAs in MMD. We have performed all lncRNA-mRNA (63431 lncRNAs, 39887 mRNAs) co-expression analyses and lncRNA/mRNA reciprocal expression pattern analyses for MMD patients. Surprisingly, the genomic position and orientation of lncRNAs and mRNAs do have relativity with their correlation coefficient.
We further performed GO and KEGG pathway analyses, which could help us understand the pathogenesis of MMD and could provide new potential therapeutic targets. Moreover, the lncRNA-mRNA co-expression analysis showed that 32 lncRNAs interacted with 11 mRNAs in the GO of inflammatory response and 15 lncRNAs interacted with 6 mRNAs in the MAPK signaling pathway (Fig. 5a,d). Association of differentially expressed lncRNAs and mRNAs with pathways relevant to MMD pathogenesis may partly explain the etiology of MMD. Inflammatory response leads to the hyperplasia of intimal vascular smooth muscle cells (VSMCs), which causes lumen stenosis of MMDs 43,44 . The MAPK signaling pathway played important roles in vascular pathological processes [45][46][47] and vascular inflammation 48 . It was reported that many cytokines induce the proliferation of VSMCs via MAPK signaling pathway [49][50][51] . MAPK signaling pathway inhibitors have been successfully used to treat atherosclerosis in vivo 52 . Our results can elucidate key lncRNAs and provide leads to further understand the pathogenesis of MMD.
There are limitations to our study. There were only 15 MMD patients and 10 healthy controls included in our analysis, and we could still identify valuable genes and pathways from our database. Further studies would enlarge the sample size. Moreover, the results were obtained from the bioinformatic analysis and microarray analysis. Therefore, further studies are needed to confirm these differentially expressed genes and pathway mechanisms with experiments. This manuscript is our preliminary work, and further work remains to be done.
In conclusion, we identified 880 differentially expressed lncRNAs (3649 probes) and 2624 differentially expressed mRNAs (2880 probes) from the microarray of MMD patients and healthy controls. The upregulated differentially expressed mRNAs were enriched for inflammatory response, the Toll-like receptor signaling pathway, chemokine signaling pathway and the MAPK signaling pathway. The integrated analysis also indicated interregulation in MMD patients. These findings may reveal the pathogenesis of MMD, and future studies should focus on the inflammatory response and MAPK signaling pathway for MMD patients.

Methods
Patients and Samples. We enrolled 15 MMD patients who had been diagnosed with MMD according to the characteristic angiographic findings with strict inclusion criteria in Beijing Tiantan Hospital and 10 healthy controls in our study. Informed consent was obtained at enrolment and the basic characteristics of all MMD patients and healthy controls are summarized in Table 1. Our study was approved by the Ethics Committee in Beijing Tiantan hospital. The study was carried out in accordance with the Declaration of Helsinki, and all methods were performed in accordance with the relevant guidelines and regulations.
We obtained the peripheral blood anticoagulated with ethylene diamine tetraacetic acid (EDTA) from all patients and healthy controls. Further, we extracted the RNA by using TRIzol reagent (Invitrogen, Grand Island, NY, USA) according to the manufacturer's instructions and quality evaluations were performed using Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA).
LncRNA and mRNA Microarrays. The Agilent Human 4 × 180 K lncRNA and mRNA Microarrays (Agilent, Santa Clara, CA, USA) were performed using a Gene Expression Hybridization Kit (Agilent, Santa Clara, CA, US) according to the manufacturer's instructions. Slides were washed in staining dishes with a Gene Expression Wash Buffer Kit (Agilent, Santa Clara, CA, USA) and scanned by an Agilent Microarray Scanner (Agilent, Santa Clara, CA, USA) with default settings according to the manufacturer's instructions. Raw data were normalized by Quantile algorithm using Gene Spring Software 12.6 (Agilent Technologies). The differentially expressed lncRNAs and mRNAs were identified by using R software (version 3.2.3) with the samr package 53 .