Differential Expression and Bioinformatic Analysis of the circRNA Expression in Migraine Patients

Background CircRNAs are noncoding RNA molecules that have recently been described and shown to regulate miRNA functionality. While recent studies have suggested such circRNAs to be associated with pain related diseases in humans, no comprehensive migraine-related circRNA profiles have been generated, and there is currently no clear understanding of whether they can serve as regulators of migraine pathology. Methods We initially conducted a circRNA microarray analysis of the plasma of migraine patients and healthy controls. Based upon these data, we then selected 8 differentially expressed circRNAs and confirmed their expression in more migraine patient plasma samples via real-time PCR. We then performed functional and pathway enrichment analyses. Lastly, using a robust rank aggregation approach, we constructed a ceRNA network according to predicted circRNA–miRNA and miRNA–mRNA pairs in these migraine patient samples. Results We were able to detect 2039 circRNAs in our patient samples, with 794 of 1245 these circRNAs being up- and downregulated in migraine patients relative to controls, respectively (fold change ≥ 1.5, p < 0.01). A qRT-PCR analysis confirmed that the expression of hsa_circRNA_100236, hsa_circRNA_102413, and hsa_circRNA_000367 was significantly enhanced in migraine patients, whereas the expression of hsa_circRNA_103809, hsa_circRNA_103670, and hsa_circRNA_101833 was significantly reduced in these individuals relative to healthy controls. We found these differentially regulated circRNAs to be associated with numerous predicted biological processes, with enrichment analyses suggesting that they may modulate the PI3K-Akt signaling so as to promote inflammation to drive migraine development. However, further research will be needed to formally test these mechanistic possibilities and to validate these circRNAs as potential biomarkers of migraine patients. Conclusions Our results offer new potential insights into the mechanistic basis of this condition and suggest that hsa_circRNA_000367 and hsa_circRNA_102413 may offer value as regulators of migraine pathology.


Introduction
Migraines are a very common form of primary headache, typically presenting as unilateral headaches of moderate to severe intensity. Migraine patients frequently suffer from frequent and potentially debilitating attacks, which can also be accompanied by photophobia, nausea, and sensitivity to noise. Migraines are estimated to affect 10-20% of individuals, with incidence peaking in young adults and middleaged individuals, in whom these headaches can markedly reduce productivity and quality of life [1]. The World Health Organization has classified migraine headaches as one of the most serious chronic functional disorders globally, highlighting the clear unmet medical need pertaining to its management and treatment [2]. As migraine pathogenesis is complex and not fully understood, at present, it can be very difficult to accurately diagnose, thus limiting the ability of clinicians to accurately treat patients or to make appropriate prognostic determinations regarding their condition. In an effort to identify more reliable biomarkers of migraine, recent studies have determined that migraine patients exhibit patterns of differential microRNA (miRNA) expression [3][4][5][6], with miRNAs also serving to modulate epigenetics in the context of chronic neuralgia [7,8] and migraine [9,10]. These miRNAs have also been found to provide functional insights into patient drug exposure [11], allowing for the identification of specific miRNA biomarkers of a given patient population or condition [12].
Recent studies have identified circular RNAs (circRNAs) as novel posttranscriptional regulators of the miRNA activity that can control diverse processes including transcription, mRNA splicing, translation, and RNA degradation. These circRNAs are generated through specific and selective splicing events, and they are expressed at high levels within eukaryotic cells. Although first detected within RNA viruses [13], RNAsequencing studies have more recently highlighted the widespread expression of circRNAs associated with a wide range of human genes [14,15]. Importantly, these circRNAs can serve as sponges for specific target miRNAs, binding, and sequestering them and thereby preventing them from executing their normal regulatory functions within cells. In addition, as they are closed loops of RNA, circRNAs are highly resistant to the RNA exonuclease activity. Given their stability, cir-cRNAs are thus potentially ideal as novel diagnostic biomarkers of a range of different conditions. No studies to date, however, have conducted exhaustive circRNA profiling in order to understand whether the circRNA expression and/or activity are associated with migraine pathogenesis.
In the present study, we therefore utilized a highthroughput sequencing approach in an effort to identify differentially expressed circRNAs in migraine patients, and we then performed appropriate bioinformatic analyses on the identified circRNAs. Based on their differential expression profiles, we were able to identify specific migraine-associated cir-cRNAs, with qRT-PCR being used to formally confirm their migraine-specific expression patterns in patient plasma samples. In addition, we conducted GO and KEGG enrichment analyses in order to understand the functionality of these cir-cRNAs, and we constructed a ceRNA network in order to begin to elucidate the mechanisms whereby these circRNAs may influence the migraine onset. Together, the results of this study will serve as a foundation for future research regarding biomarkers and mechanisms of migraine pathology. 2.2. Differential circRNA Differential Expression Analysis. We have completed the Arraystar Human circRNA Array V2 analysis of the 7 samples. Total RNA from each sample was quantified using the NanoDrop ND-1000. The sample preparation and microarray hybridization were performed based on the Arraystar's standard protocols. Briefly, total RNAs were digested with Rnase R (Epicentre, Inc.) to remove linear RNAs and enrich circular RNAs. Then, the enriched circular RNAs were amplified and transcribed into fluorescent cRNA utilizing a random priming method (Arraystar Super RNA Labeling Kit; Arraystar). The labeled cRNAs were hybridized onto the Arraystar Human circRNA Array V2 (8x15K, Arraystar). After having washed the slides, the arrays were scanned by the Agilent Scanner G2505C. Agilent Feature Extraction software (version 11.0.1.1) was used to analyze acquired array images. Quantile normalization and subsequent data processing were performed using the R software limma package. Differentially expressed circRNAs with statistical significance between two groups were identified through Volcano Plot filtering. Differentially expressed circRNAs between two samples were identified through fold change filtering. Hierarchical clustering was performed to show the distinguishable circRNA expression pattern among samples.

qRT-PCR.
After total RNA was isolated from individual samples, SuperScript III Reverse Transcriptase (Thermo Fisher Scientific) was used to prepare cDNA samples via reverse transcription. We then confirmed the differential expression of 7 different circRNAs in these donor samples using FastStart Universal SYBR Green Master (Rox) with the specific primers detailed in Table 1. Thermocycler settings were 95°C for 10 min, 40 cycles of 95°C for 10 s, and 60°C for 1 min. GAPDH served as a reference gene, and all samples were assayed in triplicate. The 2 − ΔΔCt method was used to assess the relative circRNA expression, with the circRNA expression in migraine patients being represented as fold change over control.

Pathway and Functional Enrichment Analyses. The Gene
Ontology project provides a controlled vocabulary to describe gene and gene product attributes in any organism (http://www .geneontology.org). The ontology covers three domains: biological process, cellular component, and molecular function. Fisher's exact test in Bioconductor's top GO is used to find if there is more overlap between the DE list and the GO annotation list than would be expected by chance. The p value produced by top GO denotes the significance of GO terms enrichment in the DE genes. The lower the p value, the more significant the GO Term (p value≤0.05 is recommended). The pathway analysis is a functional analysis mapping gene to KEGG pathways. The p value (EASE score, Fisher p value, or Hypergeometric p value) denotes the significance of the pathway correlated to the conditions. The lower the p value, more significant is the pathway. (the recommend p value cut-off is 0.05).

CircRNA-Related ceRNA Regulatory Network Construction.
CeRNA hypothesis RNA transcripts can crosstalk by competing for common microRNAs, with microRNA response elements (MREs) as the foundation of this interaction [16]. These RNA transcripts have been termed as competing endogenous RNAs-ceRNAs [17]. Any RNA transcripts with MREs might act as ceRNA, and ceRNAs include pseudogene transcripts, lncRNAs, circRNAs, and mRNAs; these transcripts can compete for the same microRNA response elements (MERs) to regulate mutually. To find the potential target of microRNAs, the target/microRNAs is predicted with homemade miRNA target prediction software based on TargetScan Hsa-miR-508-5p Hsa-miR-500a-5p Hsa-miR-7-5p Hsa-miR-29a-5p  [18][19][20].Through merging the common targeted miRNAs, we constructed the ceRNA network. There are three conditions that must exist for the ceRNA network to occur [16]. First, the relative concentration of the ceRNAs and their microRNAs is clearly important; second, the effectiveness of a ceRNA would depend on the number of microRNAs that it can "sponge"; third, not all of the MREs on ceRNAs are equal. So, we only accept these ceRNA-pair relations passing some measures filtering. Besides a measure with the number of common microRNAs, a hypergeometric test is executed for each ceRNA pair separately, which is defined by four parameters: (i) N is the total number of miRNAs used to predict targets; (ii) K is the number of miRNAs that interact with the chosen gene of interest; (iii) n is the number of miRNAs that interact with the candidate ceRNA of the chosen gene; and (iv) i is the common miRNA number between the two genes. The test calculates the p value by using the following formula: 3. Results

Identification of circRNAs Differentially Expressed in
Migraine Patients. We first undertook a high-throughput circRNA microarray approach in order to determine which circRNAs were differentially expressed in migraine patients, analyzing blood samples from 4 migraine patients and 3 healthy donors. This analysis identified 13,209 total cir-cRNAs in these samples distributed across all chromosomes, with the greatest proportion coming from chromosome 1 (Figure 1(d)). From a classification perspective, the majority of these circRNAs were derived from protein-coding exons, with others coming from introns, intergenic regions, and antisense regions (Figure 1(e)). An EdgeR analysis revealed that 2,039 circRNAs were differentially expressed in migraine patients (FC ≥ 1:5, p < 0:05), including 794 and 1245 that were up-and downregulated, respectively. Clear differences in the circRNA expression between control and migraine patients were evident in volcano plots (Figure 1(a)), scatter plots ( Figure 1(b)), and hierarchical clustering analyses (Figure 1(c)), with the top 10 differentially regulated cir-cRNAs shown in Table 1.

GO and KEGG Pathways Analysis of Differentially
Expressed circRNAs. We choose hsa_circRNA_103670, hsa_ circRNA_103809, hsa_circRNA_000367, and hsa_circRNA_ 102413 for analysis. We conducted GO and KEGG analyses for the mRNA molecules corresponding to differentially expressed circRNAs in migraine patients. This analysis revealed these circRNAs to exhibit significant enrichment for 103 GO biological processes, 38 GO cellular components, 29 GO molecular functions, and 3 KEGG pathways. The most significantly enriched biological process-(BP-) related terms included "nerve growth factor signaling pathway", "regulation of nonmotile cilium assembly", and "pro-B cell differentiation" (Figure 2(a)). The most enriched cellular component-(CC-) related terms included "synaptic cleft", "preribosome and small subunit precursor", and "aminoacyl-tRNA synthetase multienzyme complex" (Figure 2(b)). The most enriched molecular function-(MF-) related terms included "oxidoreductase activity, acting on CH or CH2 groups", "interleukin-1 receptor binding", and "proton-exporting ATPase activity and phosphorylative mechanism" (Figure 2(c)). The KEGG pathway analyses revealed the PI3K-Akt signaling to be among those significantly migraine-related ( Figure 3).

Construction of a ceRNA Network.
We choose hsa_cir-cRNA_103670, hsa_circRNA_103809, hsa_circRNA_000367, and hsa_circRNA_102413 for analysis. CircRNAs have been found to be capable of regulating the gene expression at least in part via functioning as sponges that sequester specific target miR-NAs. In order to explore the potential regulatory role of the identified differentially expressed circRNAs in migraine patients, we therefore used the Targetscan and miRanda applications in order to predict the miRNA targets for these circRNAs as well as the mRNA targets downstream of these miRNAs. The final miRNA-gene relationship included 57 miRNAs and 66 genes. We then constructed a ceRNA network based upon this circRNA-miRNA-mRNA coregulatory model ( Figure 4). The results of this network suggested that these circRNAs have many miRNA binding sites and thus have the potential to broadly regulate the gene expression. Additional data pertaining to this network is shown in Supplementary Table 2.

Discussion
There has been extensive research regarding circRNAs in recent years, with regulatory roles for these RNA molecules having been detected in the context of tumor formation [21], Parkinson's disease [22], Alzheimer's disease [23], and ischemic cerebrovascular disease [24]. No studies to date, however, have examined how circRNAs are associated with migraine pathogenesis. To that end, we conducted a microarray analysis to identify those cir-cRNAs which were differentially abundant in circulation in migraine patients. In total, we detected 13,209 circRNAs in our patient samples, including 794 and 1245 that were up-and downregulated in migraine patients, respectively. We then further selected 8 differentially expressed circRNAs in order to validate their differential expression between migraine patients and healthy controls. This analysis confirmed that the expression patterns for hsa_circRNA_103670, hsa_circRNA_103809, hsa_ 4 BioMed Research International    BioMed Research International circRNA_000367, and hsa_circRNA_102413 were consistent with our microarray findings. We additionally conducted GO and KEGG analyses in order to understand the potential functional implications of these differentially expressed circRNAs. We found these cir-cRNAs to be associated with GO terms including "pro-B cell differentiation", "regulation of nonmotile cilium assembly", "nerve growth factor signaling pathway", "positive regulation of inositol phosphate biosynthetic process", "synaptic cleft", "aminoacyl-tRNA synthetase multienzyme complex", "vacuolar proton-transporting V-type ATPase complex", and "interleukin-1 receptor binding". Furthermore, these circRNAs were significantly enriched for the PI3K/Akt signaling pathway.
At present, the mechanistic basis of migraine development is incompletely understood, with trigeminal nerve, vascular reflex, and cortical diffuse inhibition activities all being postulated to be involved in this process [25]. Neurogenic inflammation can occur when meningeal notional receptor peripheral terminals depolarize, leading to the trigeminal nerve and vascular system activation and the release of vasoactive peptides including calcitonin gene-related peptide (CGRP), which can cause dural vasodilation, increased intracranial blood flow, plasma protein extravasation, mast cell degranulation, and inflammation [26]. This, in turn, leads to increased primary neuron excitability in the trigeminal ganglion, as well as secondary and tertiary neuron sensitization, resulting in pain. Both inflammatory factors such as IL-1, IL-6, TNF alpha, and TGF beta, as well as anti-inflammatory factors including IL-4, IL-10, and IL-13 are related to this process [27]. Deng et al. previously found CGRP to induce the upregulation of Ithe L-6 mRNA expression via mmu_circRNA_007893, which served as a sponge for the miRNA mmu-mir-485-5p [28]. The IL-1 cytokine family is composed of 11 proteins, among which IL-1 beta has a strong proinflammatory effect and is an important mediator mediating the immune response [29]. IL-1 beta levels have been shown to be increased in the circulation of migraine patients and in the dural membrane of migraine model rats [30,31]. Han et al. also observed significant increases in circulating IL-6, IL-1 beta, and TNF levels in migraine patients relative to healthy controls [32]. This raises the possibility that circRNAs may regulate IL-6, IL-1 beta, and other inflammatory factors in order to influence migraine development. The PI3K/Akt signaling has also been shown to be active in the brain tissue of migraine model rats. When tensin homolog proteins upstream of the PI3K/Akt signaling pathway, this can lead to enhanced PIP 3, 4, and 5 activation, leading to glycogen synthase kinase-3 beta inhibition [33]. We hypothesize that certain differentially expressed circRNAs  Figure 1: Identification of migraine-associated circRNA expression profiles and qRT-PCR confirmation of differential circRNA expression profiles. Differentially expressed circRNAs are shown in a volcano plot (a) and a scatter plot (b). Those circRNAs that were up-and downregulated by at least twofold in migraine patients relative to controls are shown in red, respectively (p < 0:05). Gray dots correspond to circRNAs for which no differential expression was observed. (c) Hierarchical clustering of these circRNAs. CircRNAs that were differentially expressed in the plasma of migraine patients relative to healthy controls are shown with corresponding chromosomal locations. While circRNAs were identified corresponding to all chromosomes, many were concentrated on chromosomes 1 and 2 (d).
CircRNA types are shown, with may being derived from protein-coding exons, some being derived from introns, and a few from intergenic regions (e).The relative expression of the 8 indicated circRNAs in migraine and control patient samples, as assessed via microarray (f) and qRT-PCR (g). The expression was normalized to control samples. n = 3 samples per group. Data are means with SD. * * p < 0:01, Student's t-test.     11 BioMed Research International identified in our analysis (hsa_circRNA_103670, hsa_cir-cRNA_103809, hsa_circRNA_000367, hsa_circRNA_ 102413) may modulate the PI3K/Akt signaling and consequent production of IL-1 beta, CGRP, or other factors leading to migraine development.

BioMed Research International
Many recent studies have explored the functionality of ceRNAs, revealing them to be capable of functioning by binding to MREs, thereby releasing miRNAs from the RISC complex and preventing the degradation of their target mRNAs, leading to the corresponding gene upregulation [34]. There are multiple classes of ceRNAs, including lncRNAs, pseudogenes, and circRNAs. CircRNAs contain miRNA binding sites, allowing these circRNAs to serve as molecular sponges for specific miRNAs, thereby influencing the expression of their downstream target mRNA molecules, indirectly influencing cellular biology [14,35]. Such circRNA regulatory functions have been shown to influence the incidence of Alzheimer's disease and ischemic cerebrovascular disease [32,36]. Their role in the context of migraine development, however, has not been previously explored. The hsa_circRNA_000367 circRNA has been found to regulate the TGIF1 (inhibiting transformed growth factor-beta signal), IFITM1 (interferon-induced transmembrane protein 1), and PRKRIR (PRKRIR protein kinase, interferon-induced double-stranded RNA-dependent inhibitor) genes via a ceRNA mechanism, with hsa_circRNA_102413 also regulating IFITM1. The specific expression patterns of MCP-1, IL-10, TGF-beta, and other cytokines have also been linked with migraine disease development [37,38]. Furthermore, Saygi et al. suggested that TGF-beta 1 genotypes are associated with pediatric migraine development [39]. This raises the possibility that hsa_circRNA_000367 and hsa_cir-cRNA_102413 may modulate migraine pathogenesis via regulating the expression of these and other inflammatory factors, although further experiments will be needed to confirm this possibility.

Conclusions
In conclusion, our study identified a number of circRNAs that were differentially expressed in migraine patients. We further confirmed the differential expression of specific cir-cRNAs (hsa_circRNA_103670, hsa_circRNA_103809, hsa_ circRNA_ 000367, hsa_circRNA_102413), and further bioinformatics analyses revealed these circRNAs to be associated with migraine development. However, further experiments will be necessary in order to confirm the relationship between these circRNAs and migraine development. Even so, our findings can serve as a foundation for future studies of migraine biomarkers and functional regulators. Shown is a circRNA-miRNA gene network diagram. Yellow, red, and blue nodes correspond to circRNAs, miRNAs, and genes, respectively, with lines indicating interactions between these entities.