Identification and integrated analysis of neuropathic pain-related circular RNA signature

: Background: More and more evidences show that non-coding RNAs are involved in neuropathic pain ， however, there are few reports on the regulatory mechanism of competitive endogenous RNA (ceRNA) in neuropathic pain. The purpose of this study is to explore the possible molecular mechanisms of neuropathic pain. Methods: We collected neuropathic pain-related microarray datasets providing expression profile of circular RNAs (circRNAs) and mRNAs from the Gene Expression Omnibus (GEO) and then performed bioinformatics analysis on them. Results: The present study has identified that up-regulated circRNAs primarily regulate the activity of focal adhesion-associated biological processes and down-regulated primarily regulate the activity of metabolic-associated biological processes by means of ceRNAs. Conclusions: Our data suggest that circRNAs may be candidates for pathogenesis in neuropathic pain and may be considered as promising therapeutic targets in the future.

Neuropathic pain is the most common symptom of lumbar disc herniation. The main pathogenesis is that spinal nerve roots produce a series of inflammatory reactions based on the compression of the protruding disc [1]. Currently, surgical treatments such as discectomy and lumbar interbody fusion are used, as well as conservative treatments such as drugs and massage acupuncture; however, the overall efficacy is not good, because the pathogenesis of neuropathic pain is still unclear. Therefore, exploring the underlying mechanism of the development of neuropathic pain, and then looking for early intervention methods and effective drug intervention targets, is a major clinical problem that needs to be solved urgently.
Circular RNAs (circRNAs) are a class of noncoding RNAs with closed continuous loop structure that are ubiquitously expressed in mammalian tissues. Due to the lack of terminal 5′ and 3′ ends, circRNAs are more stable than linear RNAs [2]. CircRNAs play pivotal roles in the pathogenesis of neuropathic pain. Cao et al first detected circRNA in the chronic constriction injury (CCI) model of the sciatic nerve. It was found that there were 469 circRNA differential expression between CCI and sham-operated rats [3]. Zhou et al detected the non-coding RNA expression profile involved in neuropathic pain after retention nerve injury by sequencing. It found that the expression of microRNA (miRNA) and 188 circRNA is significantly changed [4]. Wang et al. found that circHIPK3 and miRNA-124 bind to each other and regulate inflammatory factors IL-1β, IL-6, IL-12. And TNF-α expression, while intrathecal injection of circHIPK3 shRNA can be used to treat neuropathic pain in diabetic rats [5]. Although these studies highlighted the important role of circRNAs in neuropathic pain, the expression and function of most circRNAs in neuropathic pain are still largely unknown. In addition, circRNAs are aberrantly expressed in various pathophysiological states. Through the construction of ceRNA networks, we could have a better understanding of the pathogenic mechanism of neuropathic pain.
In this study，we employed a combinative strategy of gene chip and computational biology to investigate novel circRNAs and their potential action mechanisms in neuropathic pain. We collected neuropathic pain-related microarray datasets providing expression profile of circRNAs and mRNAs from the GEO. And then we performed Bioinformatics analysis on these datasets. We aimed to find new molecular targets for therapy of neuropathic pain patients through microarray analysis and provide a solid foundation and effective tool for gene therapy of neuropathic pain.

GEO datasets
Raw data was downloaded from NCBI SRA with accession number of PRJNA558403 and GSE30691。For identification and quantification of circRNAs, reads were mapped to the reference genome Rattus norvegicus(UCSC rn6)by STAR. CIRCexplprer2 were used to quantify and annotate circRNA. For the mRNA part, data were firstly normalized by 'RMA', and low expressed probes were filtered out. Only samples of Spinal Nerve Ligation group were used for the following analysis.

Differential expression analysis
Differential expression analysis of circRNA was performed using NOIseq with parameters, replicates = 'no' and lc=1 [6]. CircRNA with a cutoff of q = 0.8 were considered as differential expressed. Considering circRNA with high expression is meaningful, circRNAs with an absolute ranking value at top10 high and q>=0.6 were also considered differential expressed. Differential expression analysis of mRNA was performed using R 'limma' package.
Probes with a cutoff of p-value<0.05 & absolute fold-changes>=1.5 were considered as differential expressed.

CircRNA-miRNA-mRNA network
The potential binding sites of miRNA to all the differential expression circRNA were predicted by miRanda with the default parameters. The miRNA targets were obtained from targetscan (http://www.targetscan.org/vert_71/) and miRanda. Then only differential expressed targets were retained. We select circRNA acts as ceRNA on the basis of the ceRNA hypothesis. Any circRNA and mRNA pair with the same regulation direction (both up or down regulation between two group of samples) and bound by the same miRNA was considered as ceRNA interaction. From the ceRNA network, with target genes of selected pathway related were retained and played out by Cytoscape 3.6.1 software.

Animals and Spinal nerve ligation (SNL)
Male Sprague-Dawley rats, weighing 200-220g, from Southeast University Laboratory Animal Center. In rats under isoflurane (4%) anesthesia, the left L5 spinal nerve was isolated adjacent to the vertebral column and tightly ligated with 6-0 silk sutures distal to the dorsal root ganglion and proximal to the formation of the sciatic nerve. In sham-operated rats, the L5 spinal nerves were identically exposed without ligation [7]. Animals were maintained under a 12-h light/dark cycle and allowed free activity in their cages. All experiments were reviewed and approved by the Institutional Animal Care and Use Committee in Southeast University School of Medicine, and were conducted according to the committees' guidelines.

Quantitative real-time polymerase chain reaction (qRT-PCR) analysis
Total RNA was extracted from pooled Sham and SNL-treated group samples using Trizol (Thermo Fisher Scientific, Inc.), and 1 μg of total RNA was reverse transcribed into first-strand cDNA using a PrimeScript RT Reagent kit (Takara Bio, Inc., Otsu, Japan), according to the manufacturer's protocols. qPCR was performed with a SYBR-Green realtime PCR kit (Thermo Fisher Scientific, Inc.) using the ABI StepOnePlus Real-Time PCR system (Applied Biosystems; Thermo Fisher Scientific, Inc.). CircRNAs were analyzed with GAPDH as the internal standard. The reactions were prepared as follows: 7.5 μl SYBR Premixm Ex Taq II, 0.25 μl ROX Reference dye II, 0.125 μl forward primer, 0.125 μl reverse primer, 5 μl RNase-free water, and 2 μl cDNA. The thermocycling conditions were: one step at 95 C for 30 sec, followed by 40 cycles of 95 C for 5 sec and 60 C for 30 sec, and a final step of 95 C for 15 sec, 60 C for 15 sec and 95 C for 15 sec. The relative expression of circRNAs was quantified using the 2 -ΔΔcq method [8]. Primer sequences are listed in Table 1.

Statistical analysis
All data are expressed as mean ± standard deviation (SD). Statistical significance was determined using Student's t-test with a p value of <0.05.

CircRNA expression profiling and enrichment analysis of maternal gene function.
We collected neuropathic pain-related microarray datasets from PRJNA558403 of the GEO and conducted data pre-processing ( Fig S1). And then we performed bioinformatics analysis on them. As shown in the heatmaps, each circRNA is presented with the name of its corresponding maternal gene. Two methods were used to screen for significant DECs.
Considering fold-change，the top 20 of the up-regulated and down-regulated sort results are shown separately (Fig 1A, 1C). Considering probability of statistical test, the top 20 of the up-regulated and down-regulated results are shown separately (Fig 1B, 1D). The scatter plot displayed a total of 111 up-regulated and 95 down-regulated DECs (Fig 1E). From the results of KEGG and GO enrichment, these DECs related maternal genes are mainly involved in neural projections, signal transduction, and adhesion (Fig 1F, 1G). This part of the results indirectly suggests that DECs are likely involved in these biological functions.

mRNA expression profiling and enrichment analysis
Considering that circRNAs are likely to be involved in the transcriptional regulation and stability of disease-associated genes, we collected neuropathic pain-related microarray datasets from GSE30691 of the GEO and conducted data pre-processing (Fig S2). And then we performed bioinformatics analysis on them. First, the transcription factors differentially expressed with fold-changes >2.0 over the time course are shown in the heatmap (Fig 2A).
Afterwards, from the results of GO enrichment and KEGG pathway analysis, part of differentially expressed mRNAs involved in axonal, signaling, metabolic and other biological functions (Fig 2B, 2C).

ceRNA analysis of mRNA and circRNA
After the basic analysis of circRNAs in the above section, we then predicted a number of circRNAs that may play a role in disease-related biological functions by means of ceRNA analysis. As shown in the histograms, number of target genes corresponding to up-and down-regulated DECs. And the higher the number is, the more important these DECs are indirectly. (Fig 3A, 3B). We then performed pathway enrichment of target genes corresponding to DECs. The KEGG pathway analysis of up-and down-regulated target genes is shown in histograms (Fig 3C，3D). The results suggest that the two biological functions of adhesion and metabolism are likely to be regulated by DECs via ceRNAs.

Construction of ceRNA network and qRT-PCR validation for the selected circRNAs
We screened TOP1 pathway genes to construct the ceRNA network, with DEGs for Focal adhesion selected in the up-regulation section and DEGs for Metabolic pathway selected in the down-regulation section (Fig S3A, S4A). We construct the ceRNA network with the miRNA part removed (Fig 4A, 4B). GO analysis of selected important up-and downregulated circRNAs is shown in the results (Fig S3B, S3C, S4B, S4C). To validate the bioanalytical results, we collected spinal dorsal horn tissues from rats in the Sham and SNLtreated groups (Day 14), qRT-PCR was used to detect the predicted two up-regulated and two down-regulated circRNAs, which were selected based on the above criteria (Fig 4C,   4D).

Discussion
Numerous studies have been conducted to reveal the pathogenic mechanisms of neuropathic pain. However, the progression of neuropathic pain still remains elusive. The function of circRNAs predicted with KEGG analysis in the mechanisms of neuropathic pain should be studied more in-depth in future work.
Lines of evidence show that focal adhesion signaling modulates axon regeneration of peripheral neuron [10]. Advillin is involved in somatosensory neuron subtype-specific axonal regeneration and neuropathic pain, and has been shown to be highly correlated with neonatal focal adhesion proteins [11]. bioinformatic studies on neuropathic pain [12]. Academic studies have shown that endocannabinoid metabolism and uptake are new targets for neuropathic pain [13]. The effects of endocannabinoids are terminated by the reuptake and metabolism of various enzymes, including fatty acid amide hydrolase (FAAH), monoacylglycerol lipase (MAGL), and cyclooxygenase type 2 (COX2), preventing the metabolism or uptake of endocannabinoids, elevating tissue levels of these lipid compounds, and producing behavioral analgesia in models of acute pain [14][15] [16]. The present ceRNA network analysis demonstrated that a circRNA /miRNA axis may have important roles in metabolic-mediated neuropathic pain. Competing endogenous RNAs (ceRNAs) are transcripts that act as miRNA sponges, modulating each other at post-transcriptional level via competitively binding to shared miRNAs [17]. Related circRNAs may affect the expression of downstream adhesion or metabolism-related genes through competitive combining to miRNAs.
Although only DEGs of the Top1 signaling pathway were selected to construct the ceRNA network, the top2 pathways are likely to be associated with pain. As a traditional Chinese Medicine, paeonia lactiflora pallas has been used to treat pain more than 1000 years in China, through PI3K-Akt signaling pathway [18]. Regulation of PI3K-Akt signaling can be an attractive alternative strategy for neuropathic pain [19] [20]. Previously, epiregulin -mediated activation of EGFR activates dorsal root ganglion (DRG) neurons, producing pain behaviors through a mechanism that involves PI3K/AKT/mTOR signaling pathway [21].
Simultaneously，early HBO therapy could significantly improve symptoms of hyperalgesia of neuropathic pain in rats, possibly via activation of the NO-cGMP-PKG signaling transduction pathway [22]. Previous study has reported that the peripheral activation of A1R plays a role in the regulation of inflammatory hyperalgesia by a mechanism involving the NO/cGMP/PKG/KATP intracellular signaling pathway [23]. Huang et al discovered that cGMP-PKG signaling pathway can be activated by the in vivo chronic compression of DRG (CCD) or in vitro acute dissociation of DRG (ADD) treatment, and that continuing activation of cGMP-PKG pathway is necessary after these two dissimilar forms of injury-related stress [24].
At the end of this study, we selected several important predicted circRNAs and used qRT-PCR to validate the results of ceRNA analysis. Atp2a2, the down-regulated circRNA.
From the results of functional enrichment analysis of the target gene, the effect is mainly on steroid metabolism and also on signaling pathways. Esrrg, the up-regulated circRNA.
From the results of functional enrichment analysis of the corresponding target genes, the main effects are cytoskeleton, ECM, cell adhesion, and so on.

Conclusions
In conclusions, by adopting a comprehensive strategy of the mining of circRNA and mRNA data in the GEO database and computational biology, we constructed a circRNA-miRNA-mRNA network and we found that circRNAs may function as ceRNAs to exert important roles in NP. Up-regulated circRNAs primarily regulate the activity of focal adhesion-associated and down-regulated circRNAs regulate the activity of metabolicassociated biological processes. In addition, qRT-PCR was performed to confirm the ceRNA network results. Our study provides new insights into the pathogenesis and treatment of neuropathic pain from a circRNA-miRNA-mRNA perspective.