Identification of the Spinal Expression Profile of Non-coding RNAs Involved in Neuropathic Pain Following Spared Nerve Injury by Sequence Analysis

Neuropathic pain (NP) is caused by damage to the nervous system, resulting in aberrant pain, which is associated with gene expression changes in the sensory pathway. However, the molecular mechanisms are not fully understood. A non-coding Ribose Nucleic Acid (ncRNA) is an RNA molecule that is not translated into a protein. NcRNAs are involved in many cellular processes, and mutations or imbalances of the repertoire within the body can cause a variety of diseases. Although ncRNAs have recently been shown to play a role in NP pathogenesis, the specific effects of ncRNAs in NP remain largely unknown. In this study, sequencing analysis was performed to investigated the expression patterns of ncRNAs in the spinal cord following spared nerve injury-induced NP. A total of 134 long non-coding RNAs (lncRNAs), 12 microRNAs (miRNAs), 188 circular RNAs (circRNAs) and 1066 mRNAs were significantly regulated at 14 days after spared nerve injury (SNI) surgery. Next, quantitative real-time polymerase chain reaction (PCR) was performed to validate the expression of selected lncRNAs, miRNAs, circRNAs, and mRNAs. Bioinformatics tools and databases were employed to explore the potential ncRNA functions and relationships. Our data showed that the most significantly involved pathways in SNI pathogenesis were ribosome, PI3K-Akt signaling pathway, focal adhesion, ECM-receptor interaction, amoebiasis and protein digestion and absorption. In addition, the lncRNA-miRNA-mRNA and circRNA-miRNA-mRNA network of NP was constructed. This is the first study to comprehensively identify regulated ncRNAs of the spinal cord and to demonstrate the involvement of different ncRNA expression patterns in the spinal cord of NP pathogenesis by sequence analysis. This information will enable further research on the pathogenesis of NP and facilitate the development of novel NP therapeutics targeting ncRNAs.

Introduction

Neuropathic pain (NP) is one type of chronic pain and is caused by primary damage and dysfunction of the nervous system, which is characterized by dysesthesia, hyperalgesia, and allodynia (Gaskin and Richard, 2012; Finnerup et al., 2016). The number of patients with this type of pain has increased rapidly in recent years (Langley et al., 2010, 2013). Indeed, the mechanisms of NP remain poorly understood, but it is known to involve nerve injury, inflammation ectopic discharge, anatomical remodeling, N Methyl D Aspartate (NMDA) and P2X receptors, etc. (Chen et al., 2012, 2013; Detloff et al., 2014). Thus, an improved knowledge of the pathogenesis of NP is crucial for the development of an effective therapeutic strategy to prevent NP and improve the curative effect. A better understanding of the genetic and various neurobiological bases could provide the clinician with important diagnosis and treatment tools.

Non-coding Ribose Nucleic Acids (ncRNAs) comprise a class of RNA molecules that typically do not encode proteins but functionally regulate protein expression (Mattick and Makunin, 2006). Based on their size, ncRNAs can be subdivided into small ncRNAs (<200 nucleotides long), including microRNAs (miRNAs) and long ncRNAs (lncRNAs), with a length >200 bp, as well as circular RNAs (circRNA) consisting of a closed continuous loop (Thum, 2014; Thum and Condorelli, 2015). It has been speculated that these ncRNAs are emerging key regulators of gene expression under physiological and pathological conditions (Wang and Chang, 2011; Yang et al., 2014; Li et al., 2016). Moreover, emerging data have shown that ncRNAs are of crucial importance in various types of pain, particularly NP (Chen et al., 2014; Shao et al., 2015; Zhang and Banerjee, 2015; Wang et al., 2016). However, the regulatory functions of ncRNAs in NP and their underlying functional mechanisms have not been systematically described. Thus, a comprehensive forecasting and analysis of the ncRNAs underlying the progression of NP are essential for the development of effective strategies to treat this troublesome disorder and to prevent its progression.

The spared nerve injury (SNI) model induces symptoms of NP, such as mechanical allodynia due to tactile stimuli that do not normally provoke a painful response (Costigan et al., 2009). The advantages of the SNI model are the robustness of the response and the lack of a requirement for expert microsurgical skills (Decosterd and Woolf, 2008). In this study, we predicted and analyzed ncRNAs of the spinal cord in SNI-induced NP using RNA sequencing techniques.

Results

Model Identification of Neuropathic Pain

To determine the effect of NP on ncRNAs in the spinal cord, the SNI model was generated using SD rats. In the SNI model in rats, the common peroneal and tibial nerves are injured, producing consistent and reproducible pain hypersensitivity in the territory of the spared sural nerve (Bourquin et al., 2006). Twelve rats were randomly divided into the CON group (n = 6) and SNI group (n = 6). The rats were evaluated by mechanical allodynia and thermal sensitivity using von Frey and Hargreaves plantar tests at 3 h before and 1, 3, 7, and 14 days after surgery (T₀, T₁, T₂, T₃, and T₄). Three rats were randomly sacrificed and samples of L₄₋₅ spinal cord were collected after the detection of pain threshold at T₄ in each group. All SNI model rats developed mechanical allodynia in the ipsilateral side at 1, 3, 7, and 14 days after SNI surgery compared to the CON group (Figure 1A). The SNI model did not modify the thermal sensitivity of the rats at any time point (Figure 1B).

FIGURE 1

Figure 1. Nociceptive behavior developed in SNI model rats. MWT (mechanical withdrawal threshold) in each time points (A), TWL (thermal withdrawal latency) in each time points (B). n = 6, ***p < 0.001 compared to the CON group.

Differentially Expressed (DE) ncRNAs and mRNAs

To determine if ncRNAs are involved in the pathological process of NP, the L₄₋₅ spinal cord of rats were analyzed using a sequencing technique at T₄. We analyzed DE ncRNAs and mRNAs using significance analysis of the microarrays method with Cuffdiff software, following the criteria q < 0.05. DE ncRNAs and mRNAs in the samples between T₀ and T₄ were shown using a Volcano plot, Venn diagram and clustering map. Information of the top 20 up-regulated and 20 down-regulated lncRNAs, circRNAs, miRNAs, and mRNAs in the SNI group compared with the CON group at 14 days after SNI are listed in Tables 1–4, respectively. All DE miRNAs are listed in Table 3 because only 12 DE miRNAs were detected. Figures 2A–C indicate the Volcano plot, Venn diagram and clustering map of DE lncRNAs, respectively. Figures 3A–C indicate the Volcano plot, Venn diagram and clustering map of DE circRNAs, respectively. Figures 4A–C indicate the Volcano plot, Venn diagram and clustering map of DE miRNAs, respectively. Figures 5A–C indicate the Volcano plot, Venn diagram and clustering map of DE mRNAs, respectively.

TABLE 1

Table 1. The detail information of the top 12 up-regulated and 20 down-regulated lncRNAs.

TABLE 2

Table 2. The detail information of the top 20 up-regulated and 20 down-regulated circRNAs.

TABLE 3

Table 3. The detail information of the up-regulated and down-regulated miRNAs.

TABLE 4

Table 4. The detail information of the top 20 up-regulated and 20 down-regulated mRNAs.

FIGURE 2

Figure 2. The expression profiling changes of lncRNAs in spinal cord of SNI rats Vocalno Plot indicate up and down regulated lncRNAs of rats in group SNI compared with group CON (A); Venn diagram showing the number of overlap lncRNAs in group SNI compared with group CON (B); Heat map of lncRNAs showing hierarchical clustering of changed lncRNAs of rats in group SNI compared with group CON. In clustering analysis, up- and down-regulated genes are colored in red and blue, respectively (C).

FIGURE 3

Figure 3. The expression profiling changes of circRNAs in spinal cord of SNI rats Vocalno Plot indicate up and down regulated circRNAs of rats in group SNI compared with group CON (A); Venn diagram showing the number of overlap circRNAs in group SNI compared with group CON (B); Heat map of circRNAs showing hierarchical clustering of changed circRNAs of rats in group SNI compared with group CON. In clustering analysis, up- and down-regulated genes are colored in red and blue, respectively (C).

FIGURE 4

Figure 4. The expression profiling changes of miRNAs in spinal cord of SNI rats Vocalno Plot indicate up and down regulated miRNAs of rats in group SNI compared with group CON (A); Venn diagram showing the number of overlap miRNAs in group SNI compared with group CON (B); Heat map of miRNAs showing hierarchical clustering of changed miRNAs of rats in group SNI compared with group CON. In clustering analysis, up- and down-regulated genes are colored in red and blue, respectively (C).

FIGURE 5

Figure 5. The expression profiling changes of mRNAs in spinal cord of SNI rats Vocalno Plot indicate up and down regulated mRNAs of rats in group SNI compared with group CON (A); Venn diagram showing the number of overlap mRNAs in group SNI compared with group CON (B); Heat map of mRNAs showing hierarchical clustering of changed mRNAs of rats in group SNI compared with group CON. In clustering analysis, up- and down-regulated genes are colored in red and blue, respectively (C).

The results of the DE ncRNAs are shown as follows: there were 144 DE lncRNAs (15 up-regulation and 129 down-regulation), 188 DE circRNAs (68 up-regulation and 120 down-regulation) and 12 DE miRNAs (6 up-regulation and 6 down-regulation) in the rat spinal cord at 14 days after SNI in the SNI group compared to the CON group, respectively. The results of the DE mRNAs are shown as follows: There were 531 up-regulated mRNAs and 535 down-regulated mRNAs at 14 days after SNI in the SNI group compared to the CON group, respectively.

Figure 6A showed the summary histogram of DE lncRNA, DE circRNA, DE miRNA, and DE mRNA.LncRNAs regulated the expression of target gene (mRNA) expression by co-localization or co-regulation. When target genes of lncRNAs are the same as the DE mRNAs, the DE mRNAs could be more directly or indirectly controlled by lncRNAs. Intersectional analysis between target mRNAs of the co-localization or co-expression with lncRNAs and DE mRNAs are shown in the Venn diagram (Figure 6B). Take the intersection of DE mRNAs with DE circRNAs (Figure 6C) and DE miRNAs (Figure 6D) to obtain the information that whether DE circRNAs and DE miRNAs could reflect the change of corresponding DE mRNAs.

FIGURE 6

Figure 6. Count of relatived ncRNAs and mRNAs in spinal cord of SNI rats. Histogram showing the number of up and down regulated ncRNAs and mRNAs in spinal cord of SNI rats (A); Venn diagram showing the overlap number of targeted mRNA of up-regulated lncRNAs, targeted mRNA of down-regulated lncRNAs, up-regulated mRNAs and down-regulated mRNAs (B); Venn diagram showing the overlap number of targeted mRNA of DE circRNAs and DE mRNAs (C); Venn diagram showing the overlap number of targeted mRNA of DE miRNAs and DE mRNAs (D).

Real-Time Quantitative Polymerase Chain Reaction (qPCR) Validation of ncRNAs and mRNAs Expression

To validate the reliability of the sequencing results and to provide the research basis for further study, the changes of some ncRNAs and mRNAs expression at 14 days after SNI were analyzed. The lncRNA (XLOC_021333 and Rn50_8_0646.1), the circRNAs (rno circ 0004058 and rno circ 0005854), the miRNAs (rno-miR-344b-1-3p and rno-miR-490-3p) and the mRNAs (Slc4a1 and Thrsp) were analyzed by qPCR (Figure 7A). The expression levels of corresponding ncRNAs and mRNAs by sequencing method were showed in Figure 7B. All of the validated ncRNAs and mRNAs were consistent with the results obtained from second generation sequencing.

FIGURE 7

Figure 7. QPCR validations of eight regulated ncRNAs in the spinal cord from SNI rats. The expressions of lncRNAs, circRNAs, miRNAs and mRNAs (A) were significantly regulated at 14 days after SNI. One-way ANOVA followed by Tukey's multiple comparison test. ***P < 0.001. The sequencing results of corresponding lncRNAs, circRNAs, miRNAs and mRNAs with qPCR verification were showed in (B).

Functional Prediction of DE ncRNAs in SNI

To examine the gene function of ncRNAs in NP, we further selected genes in which the absolute value of correlation was >0.95 to predict the function of ncRNAs using GO and KEGG pathway analysis. LncRNA might be play a regulative role with the nearby protein-coding genes. We set the threshold of co-localization as 100 kb with upstream and downstream of lncRNAs, and then major function of DE lncRNAs were predicted with enrichment analysis.

Gene Ontology (GO, http://www.geneontology.org/) is the international standard classification system of gene function (Young et al., 2010). Directed acyclic graph (DAG) is a graphical display of DE genes GO enrichment analysis results. The branch represents the relationship of inclusion, which defines the scope from increasingly small from top to bottom. The top 10 results of GO enrichment analysis are selected as the master node of DAG, and are shown together related by GO term by including the relationship, and systematically GO term shown together, where the color depth represents the enrichment degree. DAGs were drawn in the biological process (BP), cellular component (CC) and molecular function (MF), respectively. The DAG of BP, CC and MF, which GO enrichment analysis by co-localization and co-expression of genes of DE lncRNAs are shown in Figures 8A–C. The DAG of BP, CC and MF, which GO enrichment analysis by target genes of DE miRNAs are shown in Figures 9A–C. The DAG of BP, CC and MF, in which GO enrichment analysis by intersection of co-localization and co-expression of genes of DE lncRNAs, target genes of DE miRNAs and predicted mRNAs are shown in Figures 10A–C. The DAG of BP, CC and MF, in which GO enrichment analysis by intersection of targeted genes of DE circRNAs, target genes of DE miRNAs and predicted mRNAs are shown in Figures 11A–C.

FIGURE 8

Figure 8. GO analysis of lncRNAs in spinal cord of SNI rats. The significant molecular function, biological process and cellular component changed mRNAs targeted lncRNAs in spinal cord of SNI rats. Directed Acyclic Graph (DAG) is the graphical display of GO enrichment results with candidate targeted genes (A–C). The number of genes in GO term were showed in histograph (D).

FIGURE 9

Figure 9. GO analysis of miRNAs in spinal cord of SNI rats. The significant molecular function, biological process and cellular component changed mRNAs targeted miRNAs in spinal cord of SNI rats. Directed Acyclic Graph (DAG) is the graphical display of GO enrichment results with candidate targeted genes (A–C). The number of genes in GO term were showed in histograph (D).

FIGURE 10

Figure 10. GO analysis of lncRNA-micRNA-mRNA in spinal cord of SNI rats. The significant molecular function, biological process and cellular component with changed lncRNA-micRNA-mRNAs in spinal cord of SNI rats. Directed Acyclic Graph (DAG) is the graphical display of GO enrichment results with candidate targeted genes (A–C). The number of genes in GO term were showed in histograph (D).

FIGURE 11

Figure 11. GO analysis of circRNA-micRNA-mRNA in spinal cord of SNI rats. The significant molecular function, biological process and cellular component with changed circRNA-micRNA-mRNAs in spinal cord of SNI rats. Directed Acyclic Graph (DAG) is the graphical display of GO enrichment results with candidate targeted genes (A–C). The number of genes in GO term were showed in histograph (D).

According to the distribution of predicted target genes in the Gene Ontology, the gene function of DE ncRNAs were clarified. The number of genes was statistically analyzed with significant enrichment of each GO term, and displayed in the form of a histogram. On the basis of the GO analysis of the co-localization and co-expression genes of DE lncRNAs, the most significantly enriched BP were developmental process, localization and multicellular organismal development, and the most noteworthy enriched CC were membrane-bounded vesicle, extracellular region and extracellular region part. The most significantly enriched MF were binding, protein binding and receptor binding (Figure 8D). Based on the GO analysis of targeted genes of DE miRNAs, the most significantly enriched BP were localization, transport and establishment localization, and the most noteworthy enriched CC were intracellular membrane-bounded organelle, membrane-bounded organelle and nucleus. The most significantly enriched MF were binding, catalytic activity and protein binding (Figure 9D). Based on the GO analysis of intersection of co-localization and co-expression of genes of DE lncRNAs, target genes of DE miRNAs and predicted mRNAs, the most significantly enriched BP were developmental process, single-organism developmental process and localization, and the most noteworthy enriched CC were membrane, extracellular region and extracellular region part. The most significantly enriched MF were protein binding, ion binding and transporter activity (Figure 10D). Based on the GO analysis of intersection of targeted genes of DE circRNAs, target genes of DE miRNAs and predicted mRNAs, the most significantly enriched were the same as the analysis results obtained from co-localization and co-expression of genes of DE lncRNAs, target genes of DE miRNAs and predicted mRNAs (Figure 11D). The most striking category of gene function will be the focus in a future study.

Different genes coordinate with each other to exercise their biological functions in the organism. The main biochemical pathways and signal transduction pathways involved with the candidate target genes can be determined through the pathway of significant enrichment. KEGG (Kyoko Encyclopedia of Genes and Genomes) is a major public database on pathways, which can adopt KEGG Pathway as a unit, apply the hypergeometric test, and determine the significant enrichment pathway in the candidate target genes compared with the entire genome background (Mao et al., 2005; Kanehisa et al., 2008). An enriched scatter diagram of the candidate target genes is the graphic display of KEGG enrichment analysis. In this graphic, the degree of KEGG enrichment is assessed by the Rich factor, Q-value and number of genes. When the rich factor is greater, the Q-value is closer to zero, and the number of genes is greater, then the enrichment is more significant. The top 20 pathways were displayed in the figure. When the data were KEGG-analyzed with intersection of co-localization and co-expression of genes of DE lncRNAs, target genes of DE miRNAs and predicted mRNAs, the most significantly involved pathways in SNI pathogenesis were ribosome, Phosphatidyl Inositol 3-kinase (PI3K)-Akt signaling pathway, focal adhesion, extracellular matrixc (ECM)-receptor interaction, amoebiasis and protein digestion and absorption (Figure 12A). These results were similar with enriched pathways, which were analyzed with predicted genes in SNI pathogenesis (Figure 12B). The main biochemical pathways and signal transduction pathways determined by KEGG analysis will provide further insight toward future research directions of ncRNAs.

FIGURE 12

Figure 12. ncRNAs and mRNAs enriched KEGG pathway scatterplot in spinal cord of SNI rats. A,B,C,D co-location; E,F,G,H co-expression. LncRNA-micRNA-mRNAs enriched KEGG pathway scatterplot showing statistics of pathway enrichment in spinal cord of SNI rats (A), mRNAs enriched KEGG pathway scatterplot showing statistics of pathway enrichment in spinal cord of SNI rats (B).

Regulatory Network of ncRNA and mRNA

To examine the molecular mechanism in ncRNAs involved in NP pathogenesis, we performed regulatory network analysis of ncRNA and mRNA in SNI pathogenesis. LncRNA has extensive regulatory function, and it can not only directly regulate the structure of DNA and transcription and translation of RNA, but it can also inhibit target gene regulation of miRNA to indirectly regulate gene expression, thereby acting as a miRNA sponge to bind miRNA competitively with its binding sites. Based on the theory of competing endogenous RNA (ceRNA), lncRNA-gene pairs with same binding sites of miRNA were identified. By constructing lncRNA–miRNA-gene pairs with lncRNA as a decoy, miRNA as the center, mRNA as the target, the ceRNA regulation network were built (Figure 13). Because there are binding sites between the circRNA and miRNA, the regulatory role for the miRNA target gene was inhibited and the gene expression were indirectly regulated by circRNA, which competitively bound the miRNA as a miRNA sponge. Based on the theory of ceRNA, circRNA-gene pairs with same binding sites of miRNA were identified. By constructing circRNA–miRNA-gene pairs with circRNA as a decoy, miRNA as center, mRNA as target, the ceRNA regulation network were generated (Figure 14). The mechanism of gene expression regulated by ncRNA was revealed at the full transcriptome level via the regulatory network of the ceRNA. These results illustrate the regulatory relationship between ncRNA and mRNA in the mechanism of NP. Indeed, as previously described above, the regulatory role of ncRNAs in NP pathogenesis was very complicated such that an in-depth study should be implement in the future.

FIGURE 13

Figure 13. lncRNA-micRNA-mRNAs regulatory network analysis of ncRNAs in spinal cord of SNI rats. Figrue 13 is interaction network of lncRNA-micRNA-mRNAs in spinal cord of SNI rats.

FIGURE 14

Figure 14. cirRNA-micRNA-mRNAs regulatory network analysis of ncRNAs in spinal cord of SNI rats. Figrue 14 is interaction network of circRNA-micRNA-mRNAs in spinal cord of SNI rats.

The most common mode of action with ncRNA pairs is sponging effect as ceRNA. Considering the apoptosis of neurons in spinal cord play an important role in NP, and miR-184 had been confirmed to be involved in cell apoptosis (Bi et al., 2016). Therefore, we selected two pairs (miR-184 and LNC_001457, miR-184 and rno_circ_0006928) to verify in neural stem cells of rat with dual-luciferase reporter system. Luciferase assay revealed that miR-184 displayed a sponging effect for LNC_001457 (Figure 15A) and crno_circ_0006928 (Figure 15B) and decreased luciferase activity. Acturally, miR-184 was downregulated in SNI in Table 3 and LNC_001457 and circ_0006928 were upregulated in Figures 13, 14, respectively. These results were consistent with our speculation. On the one hand, the results verified the accuracy of the network interaction of ncRNAs. On the other hand, although the sequencing experiment is not a single cell line, the data combined with references of cell function could effectively compensated for the insufficient in the follow-up study.

FIGURE 15

Figure 15. Confirmation of the pairs relationship. Luciferase assays using reporter constructs lncRNA (LNC_001457) (A) or cirRNA (rno_circ_0006928) (B) were performed in NSCs transfected with miRNA-184 (or a control). n = 5, *p < 0.05 compared to the blank group.

Discussion

In the present study, we determined that some ncRNAs and mRNAs were induced significant changes in the spinal cord in response to SNI-induced NP. We also predicted potential functions of DE ncRNAs using GO and KEGG pathway analysis in the SNI model. In addition, the regulatory networks of ncRNA and mRNA were constructed. These findings prompted the proposal that ncRNAs played a significant role in NP processing, and sequencing analysis revealed a potential therapeutic target of NP.

Neuropathic pain (NP) is a debilitating disease that affects the central and peripheral nervous systems, (Zhuo et al., 2011). In addition to spontaneous pain and hyperalgesia, a debilitating symptom of NP is pain hypersensitivity to normally innocuous stimuli, known as tactile allodynia (Baron, 2006). The development of NP is a complex mechanism and is still poorly understood. Future studies of the NP mechanism may have important clinical implications and combined treatments may provide insight toward new therapeutic targets. The spinal dorsal horn receives sensory information from the peripheral nervous system following nociceptive stimuli (Prescott et al., 2014). The pain pathway of the spinal dorsal horn critically contributes to pathologically enhance NP processing (Molander and Grant, 1986; Beggs and Salter, 2006). However, the mechanisms of the spinal dorsal horn are complex and have implications. Non-coding RNA (ncRNA) is a class of genetic, epigenetic and translational regulators, which consists of two major classes according to their size: short ncRNA (<200 nt) and long ncRNA (>200 nt). The transcribed genomic output, which consists of diverse classes of ncRNAs, plays important roles and/or are affected by many biochemical cellular processes, and dysregulation of ncRNA is associated with a variety of diseases, including NP (Guttman and Rinn, 2012; Li et al., 2012; Yu et al., 2013). However, how ncRNAs play such an important and complicated role in NP pathogenesis remains unknown. Thus, we investigated DE ncRNAs and predicted the function and regulatory interaction among ncRNAs and mRNAs in SNI-induced NP, which is essential to thoroughly and comprehensively determine how ncRNAs are involved in NP in further studies.

The SNI model was chosen in the present study to examine ncRNAs in the NP process. The SNI model presents consistent and reproducible pain hypersensitivity in the territory of the spared sural nerve by ligating the common peroneal and tibial nerves, and appear to be advantageous in some aspects, such as less tissue damage, consisting of simple and intuitive experiment steps, good behavior reproducibility and more similar clinical manifestations compared to other NP models, such as chronic constrictive injury (CCI), partial sciatic nerve ligation (PSL), and spinal nerve ligation (SNL) (Bourquin et al., 2006). In our study, all SNI model rats developed mechanical allodynia on the ipsilateral side at T₁₋₄, and did not show modification of the thermal sensitivity at any time point. We observed a change in hyperalgesia according to the classic behavior characteristics of the SNI model, which was a reliable precondition of the sequencing results. The spinal cord is a crucial site for nociceptive processing, which materializes to receive, process and relay sensory perception and participates in the coordination of sensory-motor output (Guo and Hu, 2014). Numerous cellular and molecular mechanisms, which involve neurons or glia and astrocytes in the spinal cord were proven to play a critical role in the processing of NP (Alfonso Romero-Sandoval and Sweitzer, 2015). Thus, we examined the DE ncRNAs in the spinal cord of rats in the present study, and our results demonstrated that numerous ncRNAs and mRNAs were significantly up- or down-regulated in the spinal cord of rats in response to SNI. These data indicated that the spinal cord, as an important pain signal transmission center, played a critical role in the cellular and molecular mechanisms of ncRNAs in NP.

Findings obtained from a large body of studies have revealed causal roles of microRNAs in chronic pain (Sakai and Suzuki, 2015). Indeed, miRNAs had been frequently reported as a therapeutic potential target of NP in the past few years (Andersen et al., 2014; Heyn et al., 2016). Recent studies have also shown that peripheral noxious stimuli induced changes in the expression of lncRNAs, which are associated with pain hypersensitivity underlying NP (Zhao et al., 2013). Previous research studies performed gene chip studies to screen and predict DE lncRNAs in the spinal cord SNL model mice (Jiang et al., 2015). Although much evidence has been found, no comprehensive analysis on ncRNAs involved in NP has been performed. Thus, we performed second generation sequencing to analyze DE ncRNAs in the spinal cord of rats with SNI induced-NP. Our results showed that a total of 134 lncRNAs, 12 miRNAs, 188 circRNAs, and 1066 mRNAs were significantly dysregulated at 14 days after SNI surgery. These data indicated that many lncRNAs were involved in NP, and these lncRNAs differed from those observed in Bao-Chun Jiang's research study. One reason for this discrepancy is different species used in the study. Another reason is the analytical results from second generation sequencing were more full-scale than the gene chip studies. Because more novel lncRNAs were detected in the present study, this only highlights the lack of previous knowledge regarding the role of lncRNA in NP. Importantly, circRNA are highly stable, circularized long non-coding RNAs. CircRNAs are conserved across species and are specifically enriched in the nervous system (Van Rossum et al., 2016). Numerous circRNAs have been recently identified; however, the biological function of most circRNAs remains unclear. In addition, there are few research studies on circRNAs in NP. Our research study screened 188 DE circRNAs in NP pathogenesis, which provides further insight on the underlying mechanisms of circRNAs in NP.

Although we screened DE ncRNAs in the spinal cord of rats with SNI induced-NP and a few ncRNAs were confirmed in the present study, these underlying mechanisms of ncNRAs in NP are poorly understood. The Gene Ontology (GO) is a major bioinformatics tool that unifies the representation of genes and gene product attributes across all species (Altshuler et al., 2008). GO terms and GO annotations have proven to be good predictors of gene function and trend (Li et al., 2009; Du et al., 2016). KEGG pathway databases store the higher order functional information for systematic analysis of gene functions and are more widely used in current enrichment analysis platforms (Camon et al., 2004). Thus, we analyzed ncRNAs-related gene functions and the corresponding pathways in the spinal cords of SNI induced-NP rats with GO and KEGG term enrichment analyses in the present study. Our data showed that the most significantly involved pathways in SNI pathogenesis were ribosome, PI3K-Akt signaling pathway, focal adhesion, ECM-receptor interaction, amoebiasis and protein digestion and absorption. The function of ncRNAs predicted with GO and KEGG analysis in the mechanisms of NP should be studied more in-depth in future work.

RNA transcripts with microRNA response elements (MREs) might act as ceRNA, which include pseudogene transcripts, lncRNAs, circRNAs and mRNAs, which act as natural miRNA sponges to suppress miRNA function using shared MREs for mutual regulation (Ashwal-Fluss et al., 2014). Although theoretical and experimental studies have partially revealed that lncRNAs and/or miRNAs can regulate gene expression involved in mechanisms of NP (Bali and Kuner, 2014; Xiong et al., 2015), circRNAs appears to be more correlated with pain because they were more enriched in neurons and have been proven to regulate synaptic function (You et al., 2015). However, the molecular mechanisms underlying the interaction of ncRNAs and mRNAs in NP remain largely unclear. Thus, for the first time, the lncRNA-miRNA-mRNA and circRNA-miRNA-mRNA network of NP was constructed based on our second sequencing data. These pioneering discoveries might enrich our understanding on the mechanisms underlying the role of non-coding RNAs in NP pathogenesis.

It is worth mentioning that the spinal cord contains a variety of cell populations, including many kinds of neuron subtypes, astrocytes and microglia. The GO, KEGG and ncRNA regulative network analyses we proposed in this paper just made the spinal cord as a whole. Therefore, lots of experiment in vitro and in vivo need to further verified and analyzed. The single-celled sequencing analysis with neurons of spinal cord are implementing with our team. Related results will be combined with the current data, can be more accurately described ncRNAs function of neurons or gliocytes in spinal cord in mechanism of NP.

Because lncRNA, miRNA, and circRNAs and mRNA are studied together for their predicted role and regulated relationship, our study presented an innovative data integration analysis of ncRNA and mRNA in NP pathogenesis. We also provided a catalog of predicted ncRNAs in NP and identified and confirmed several genes using sequencing, which may be of relevance for other research groups. These abundant data prompted the proposal that ncRNAs may interact and regulate their related protein-gene expression and play a key role in the pathogenesis of NP. The next step is to further study these predicted ncRNAs with regard to their complete proteomic and relevant signaling pathway, which may ultimately enable the full disclosure of the mechanisms underlying NP.

Methods

Animals

Adult male Sprague Dawley rats (body weight 250–280 g, the Laboratory Animal Center of The First People's Hospital of Foshan, Foshan, China) were randomly assigned to standard cages with 4 to 5 animals per cage. All animal procedures were in accordance with national and international animal care and ethical guidelines and have been approved by the institutional animal welfare committee. The environment was maintained at a constant temperature (22 ± 0.581°C) and relative humidity (60–70%) with a 12-h light/dark cycle (lights on at 7 a.m.). All animals were fed with a standard laboratory diet and tap water was provided ad libitum. The animals were placed in the experimental room 24 h before behavioral testing for acclimatization.

SNI Model

The rats were fixed to an operating table after anesthesia, and the femur at the right mid-thigh level was used as the incision landmark. The 3 peripheral branches (the sural, common peroneal, and tibial nerves) of the sciatic nerve were exposed at mid-thigh level distal to the trifurcation and exposed without stretching of the nerve structures and freed of connective tissue under sodium pentobarbital anesthesia (45 mg/kg, intraperitoneally). The tibial and common peroneal nerves were tightly ligated by two knots, 4 mm apart, using 6.0 silk (Ethicon, Johnson & Johnson Inc., Brussels, Belgium) and were completely severed between the knots. The sural nerve was left intact. The sural nerve was carefully preserved by avoiding any nerve stretch or contact with surgical tools. The muscle and skin were closed in 2 distinct layers using silk 5.0 suture (Decosterd and Woolf, 2008).

Nociceptive Behavior

All behavioral tests were executed in a blind manner. Nociceptive responses to mechanical and thermal stimuli were measured 3 h before and 1, 3, 7, and 14 days after surgery. Mechanical withdrawal threshold (MWT) was measured to an applied von Frey hair (Stoelting, Wood Dale, Illinois, USA) to the right hind paw until a positive sign of pain behavior was elicited (Chaplan et al., 1994). The test area was the mid-plantar claw in the distribution area of the sciatic nerve. Von Frey filaments with logarithmically incremental stiffness (0.4–15.1 g) were serially applied to the paw using the up-down method. The hairs were presented in ascending order of strength, perpendicular to the plantar surface, with sufficient force to cause slight bending against the paw and were held for 6–8 s. A positive response was noted if the paw was sharply withdrawn. Flinching immediately upon removal of the hair was also considered a positive response. The 15.1-g hair was selected as the upper cutoff for testing. Animals that exhibited no response to the 15.1-g filament were assigned this cutoff value. The bending force that evoked 50% of paw withdrawal occurrence was set as the MWT. The Hargreaves test was used to measure paw thermal withdrawal latency (TWL) to heat stimuli and to determine the presence of thermal hyperalgesia (Hargreaves et al., 1988). The rat was placed on the surface of a 2-mm-thick glass plate covered with a plastic chamber (20 × 20 × 25 cm). TWLs were measured using a radiant thermal stimulator (BME410A, Institute of Biological Medicine, Academy of Medical Science, Tianjin, and China). Heat was concentrated on the hind paw, which was flush against the glass, and radiant heat stimulation was delivered to the site. The stimulus terminated with hind paw movement, and a 25-s cutoff was imposed on the stimulus duration to prevent tissue damage. The intensity of thermal stimulation remained constant throughout the experiment. Five stimuli were applied to the same site, and the average TWL from three thermal stimulations was obtained. This mean TWL was used as the steady state of the TWL value.

Tissue Collection

We prepared 12 rats for SNI and sham operation. At 14 days after operation, three rats were randomly selected in each group and deeply anesthetized with isoflurane after the behavioral test. The L₄₋₅ spinal cord segments were collected after perfusion. Total RNA was extracted from the spinal cord dorsal horn tissue using Trizol reagent (Invitrogen, Carlsbad) according to the manufacturer's instructions.

RNA Isolation and RNA Quantification and Qualification

RNA degradation and contamination was monitored on 1% agarose gels. RNA purity was confirmed using the NanoPhotometer^® spectrophotometer (IMPLEN, CA, USA). The RNA concentration was measured using the Qubit^® RNA Assay Kit in Qubit^® 2.0 Fluorometer (Life Technologies, CA, USA). RNA integrity was assessed using the RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, CA, USA). The quantitative method of lncRNAs and miRNAs were same as conventional procession of mRNAs. Quantification of circRNAs should be performed with exonuclease to exclude non-circRNAs. Take two copies of the same amount of RNA extraction, one sample was digested the linear RNA with RNase R (Cat. No. RNR07250, Epicentre Company, USA) and leaved only circRNAs, another sample did not treated with RNase R. The two RNA were processed to reverse transcriptase meanwhile. The samples with RNase treatment were used to detect circRNAs, the samples untreated with RNase were used to detect β-actin. Library preparation for lncRNA sequencing.

A total of 3 μg RNA per sample was used as input material for the RNA sample preparations. First, ribosomal RNA was removed using the Epicenter Ribo-zero™ rRNA Removal Kit (Epicenter, USA), and rRNA-free residue was cleaned up by ethanol precipitation. Subsequently, sequencing libraries were generated using the rRNA-depleted RNA by NEBNext^® Ultra™ Directional RNA Library Prep Kit for Illumina^® (NEB, USA) following manufacturer's recommendations. Briefly, fragmentation was performed using divalent cations under elevated temperature in NEBNext, First Strand Synthesis Reaction Buffer (5X). First strand cDNA was synthesized using random hexamer primers and M-MuLV Reverse Transcriptase (RNaseH-). Second strand cDNA synthesis was subsequently performed using DNA Polymerase I and RNase H. In the reaction buffer, dNTPs with dTTP were replaced by dUTP. Remaining overhangs were converted into blunt ends via exonuclease/polymerase activities. After adenylation of 3′ ends of DNA fragments, NEBNext, Adaptor with hairpin loop structure were ligated to prepare for hybridization. To select cDNA fragments of preferentially 150–200 bp in length, the library fragments were purified with AMPure XP system (Beckman Coulter, Beverly, USA). Then, 3 μl USER Enzyme (NEB, USA) was used with size-selected, adaptor-ligated cDNA at 37°C for 15 min followed by 5 min at 95°C before PCR. Then, PCR was performed with Phusion High-Fidelity DNA polymerase, Universal PCR primers and Index (X) Primer. Finally, products were purified (AMPure XP system) and library quality was assessed on the Agilent Bioanalyzer 2100 system.

Library Preparation for Small RNA Sequencing

A total of 3 μg total RNA per sample was used as input material for the small RNA library. Sequencing libraries were generated using NEBNext^® Multiplex Small RNA Library Prep Set for Illumina^® (NEB, USA.) following the manufacturer's recommendations and index codes were added to attribute sequences to each sample. Briefly, the NEB 3′ SR Adaptor was directly and specifically ligated to the 3′ end of miRNA. After the 3′ ligation reaction, the SR RT Primer hybridized to the excess of 3′ SR Adaptor (that remained free after the 3′ ligation reaction) and transformed the single-stranded DNA adaptor into a double-stranded DNA molecule. This step is important to prevent adaptor-dimer formation. In addition, dsDNAs are not substrates for ligation mediated by T4 RNA Ligase 1 and thus do not ligate to the 5' SR Adaptor in the subsequent ligation step. Next, the 5'end adapter was ligated to 5'ends of miRNAs. First strand cDNA was synthesized using M-MuLV Reverse Transcriptase (RNase H–). PCR amplification was performed using LongAmp Taq 2X Master Mix, SR Primer for Illumina and index (X) primer. PCR products were purified on an 8% polyacrylamide gel (100 V, 80 min). DNA fragments corresponding to 140–160 bp (the length of small noncoding RNA plus the 3′ and 5′ adaptors) were recovered and dissolved in 8-μL elution buffer. At last, library quality was assessed on the Agilent Bioanalyzer 2100 system using DNA High Sensitivity Chips.

Clustering and Sequencing of ncRNA

The clustering of the index-coded samples was performed on a cBot Cluster Generation System using TruSeq PE Cluster Kit v3-cBot-HS (Illumina) according to the manufacturer's instructions. After cluster generation, the libraries were sequenced on an Illumina HiSeq 2500 platform, 125 bp paired-end and 50-bp single-end reads were generated. The reads number of lncRNAs, circRNAs and miRNAs for each samples were, respectively 94210880-113312518, 94210880-113312518 and 14658590-18771808. The clean bases of lncRNAs, circRNAs and miRNAs for each samples were, respectively 13.74G-16.43G, 13.74G-16.43G, and 0.733G-0.939G. The error rate of lncRNAs, circRNAs and miRNAs for each samples were respectively <0.02, <0.02, and <0.01%. The transcription with splicing of each sample were combined and screened as lncRNAs with Cuffmerge Software, and the conditions were as follows: the number of exon≥2, length > 200 bp, FPKM ≥0.5 (Cuffquant) and to eliminate overlapping and coding potential transcription with annotation of database at exon region (Cuffcompare Software). CircRNAs were identified base on the data of lncRNAs with find_circ (Memczak et al., 2013). Clean reads were screened the lengh of 21–22 nt as miRNA, and located to reference sequence with bowtie. Combined with miREvo Software (Wen et al., 2012) and mirdeep2 Software (Friedländer et al., 2012) to analysis the funtions of new miRNAs. Adopt DESeq2 with negative binomial distribution (Love et al., 2014) to analyse differentially expression of ncRNAs All sequencing program were performed by Novogene Company (China, Beijing).

Quantitative Real-Time RT-PCR

The remaining three rats in each group were validated using real-time PCR. Spinal cord tissues were separated under anesthetic after the behavioral test. The first-strand cDNA was synthesized using reverse transcription kit (Vazyme) according to the manufacturer's instruction and amplified in triplicate using IQ SYBR green SuperMix reagent (Bio-Rad, Hercules, CA) with an Opticon real-time PCR machine (MJ Research, Waltham, MA), according to the manufacturer's instructions. The specificity of real-time PCR was confirmed via routine agarose gel electrophoresis and melting-curve analysis. The sequences of all primers are shown in Table 5. The comparative Ct method (ΔΔCt) was used to quantify gene expression, and the relative quantification was calculated as 2−ΔΔCt. Calculation of lncRNAs and circRNAs expression were consistent with mRNA, the internal control were β-actin. The expression levels of miRNAs were normalized to U6 level in each sample.

TABLE 5

Table 5. Primers designed for qRT-PCR validation of candidate ncRNAs and mRNAs.

Go Annotations and KEGG Pathways Analysis

Gene Ontology (GO) annotations and KEGG pathway analysis were applied to investigate the roles of all DE ncRNAs. Briefly, GO analysis was applied to elucidate genetic regulatory networks of interest by forming hierarchical categories according to the BP, CC, and MF aspects of the differentially expressed genes (http://www.geneontology.org). Pathway analysis was performed to explore the significant pathways of the differentially expressed genes, according to KEGG (http://www.genome.jp/kegg/).

Analysis of ncRNAs Regulatory Network

To reveal the role and interactions among ncRNAs and mRNAs in NP pathogenesis, we constructed a ncRNAs regulatory network. The interaction network was built and visually displayed using Cytoscape software based on the screening of lncRNA-miRNA-gene pairs and circRNA-miRNA-gene pairs with miRanda software and psRobot software. Different shapes represent different RNA types, and different colors represent the regulated relationship. The size of the node was in direct proportion to the degree. Thus, these significant nodes are in a core position in the regulated network, which are more associated with NP.

Luciferase Assay

A dual-luciferase reporter system E1960 (Promega, Madison, WI, USA) was used to perform luciferase activity assay. In brief, embryonic neural stem cells (NSCs) of rats were cultured on 12-well tissue culture plates at a density of 2 × 10⁵ cells per well. Cells were co-transfected with the luciferase reporter constructs contain lncRNA (LNC_001457) or cirRNA (rno_circ_0006928), miRNA(miR-184) mimics and Renilla luciferase construct for 5 h (Lipofectamine^® MessengerMAX™ Transfection Reagent, Thermo Fisher Scientific). After 3d culture at 37°C, the transfected cells were lysed by 150 μl of passive lysis buffer. In total, 30 μl of lysates were mixed with 50 μl of LAR II, and then firefly luciferase activity was measured by a luminometer. For the internal control, 50 μl of Stop and Glo reagent was added to the sample.

Statistical Analysis

The data are presented as the means ± SEM. The results from the behavioral study were statistically analyzed using one-way or two-way analysis of variance (ANOVA). The RT-PCR results were analyzed using one-way analysis of variance followed by Tukey's multiple comparison test. Significance was established at p < 0.05.

Ethics Statement

The treatment of animals used in the study was approved by the First People's Hospital of Foshan Ethics Committee.

Author Contribution

JZ is responsible for the construction of animal model and collection of results and analysis of ncRNAs regulatory network; QX is responsible for the RNA isolation and RNA quantification and qualification; HC is responsible for the library preparation for lncRNA and Small RNA sequencing and Quantitative Real-Time RT-PCR; CY is responsible for tissues collection and the satistical analysis; YF is responsible for the preparation for lncRNA and Small RNA sequencing and analysis of ncRNAs regulatory network.

Funding

National Natural Science Foundation of China (81300974), the Natural Science Foundation of Guangdong Province (2015A030313899) and the Medical Scientific Research Foundation of Guangdong Province (A2015013).

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

References

Alfonso Romero-Sandoval, E., and Sweitzer, S. (2015). Nonneuronal central mechanisms of pain: glia and immune response. Prog. Mol. Biol. Transl. Sci. 131, 325–358. doi: 10.1016/bs.pmbts.2014.11.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Altshuler, D., Daly, M. J., and Lander, E. S. (2008). Genetic mapping in human disease. Science 322, 881–888. doi: 10.1126/science.1156409

PubMed Abstract | CrossRef Full Text | Google Scholar

Andersen, H. H., Duroux, M., and Gazerani, P. (2014). MicroRNAs as modulators and biomarkers of inflammatory and neuropathic pain conditions. Neurobiol. Dis. 71, 159–168. doi: 10.1016/j.nbd.2014.08.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Ashwal-Fluss, R., Meyer, M., Pamudurti, N. R., Ivanov, A., Bartok, O., Hanan, M., et al. (2014). CircRNA biogenesis competes with pre-mRNA splicing. Mol. Cell. 56, 55–66. doi: 10.1016/j.molcel.2014.08.019

PubMed Abstract | CrossRef Full Text | Google Scholar

Bali, K. K., and Kuner, R. (2014). Non-coding RNAs: key molecules in understanding and treating pain. Trends Mol. Med. 20, 437–448. doi: 10.1016/j.molmed.2014.05.006

CrossRef Full Text | Google Scholar

Baron, R. (2006). Mechanisms of disease: neuropathic pain–a clinical perspective. Nat. Clin. Pract. Neurol. 2, 95–106. doi: 10.1038/ncpneuro0113

PubMed Abstract | CrossRef Full Text | Google Scholar

Beggs, S., and Salter, M. W. (2006). Neuropathic pain: symptoms, models, and mechanisms. Drug Dev. Res. 67, 289–301. doi: 10.1002/ddr.20094

CrossRef Full Text | Google Scholar

Bi, X., Cao, Y., Chen, R., Liu, C., Chen, J., and Min, D. (2016). MicroRNA-184 promotes proliferation and inhibits apoptosis in HaCaT cells: an in vitro study. Med. Sci. Monit. 22, 3056–3061. doi: 10.12659/MSM.897250

PubMed Abstract | CrossRef Full Text | Google Scholar

Bourquin, A. F., Süveges, M., Pertin, M., Gilliard, N., Sardy, S., Davison, A. C., et al. (2006). Assessment and analysis of mechanical allodynia-like behavior induced by spared nerve injury (SNI) in the mouse. Pain 122, 14.e1–14. doi: 10.1016/j.pain.2005.10.036

PubMed Abstract | CrossRef Full Text | Google Scholar

Camon, E., Barrell, D., Lee, V., Dimmer, E., and Apweiler, R. (2004). The Gene Ontology Annotation (GOA) database-an integrated resource of GO annotations to the UniProt knowledgebase. In Silico Biol. 4, 5–6.

PubMed Abstract | Google Scholar

Chaplan, S. R., Bach, F. W., Pogrel, J. W., Chung, J. M., and Yaksh, T. L. (1994). Quantitative assessment of tactile allodynia in the rat paw. J. Neurosci. Methods 53, 55–63. doi: 10.1016/0165-0270(94)90144-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, H. P., Zhou, W., Kang, L. M., Yan, H., Zhang, L., Xu, B. H., et al. (2014). Intrathecal miR-96 inhibits Nav1.3 expression and alleviates neuropathic pain in rat following chronic construction injury. Neurochem. Res. 39, 76–83. doi: 10.1007/s11064-013-1192-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, Y. W., Lin, M. F., Chen, Y. C., Hung, C. H., Tzeng, J. I., and Wang, J. J. (2013). Exercise training attenuates postoperative pain and expression of cytokines and N-methyl-D-aspartate receptor subunit 1 in rats. Reg. Anesth. Pain Med. 38, 282–288. doi: 10.1097/AAP.0b013e31828df3f9

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, Y. W., Li, Y. T., Chen, Y. C., Li, Z. Y., and Hung, C. H. (2012). Exercise training attenuates neuropathic pain and cytokine expression after chronic constriction injury of rat sciatic nerve. Anesth. Analg. 114, 1330–1337. doi: 10.1213/ANE.0b013e31824c4ed4

PubMed Abstract | CrossRef Full Text | Google Scholar

Costigan, M., Scholz, J., and Woolf, C. J. (2009). Neuropathic pain: a maladaptive response of the nervous system to damage. Annu. Rev. Neurosci. 32, 21–32. doi: 10.1146/annurev.neuro.051508.135531

PubMed Abstract | CrossRef Full Text | Google Scholar

Decosterd, I., and Woolf, C. J. (2008). Spared nerve injury: an animal model of persistent peripheral neuropathic pain. Pain 87, 149–158. doi: 10.1016/S0304-3959(00)00276-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Detloff, M. R., Smith, E. J., Quiros Molina, D., Ganzer, P. D., and Houlé, J. D. (2014). Acute exercise prevents the development of neuropathic pain and the sprouting of non-peptidergic (GDNF- and artemin-responsive) c-fibers after spinal cord injury. Exp. Neurol. 255, 38–48. doi: 10.1016/j.expneurol.2014.02.013

PubMed Abstract | CrossRef Full Text | Google Scholar

Du, J., Li, M., Yuan, Z., Guo, M., Song, J., Xie, X., et al. (2016). A decision analysis model for KEGG pathway analysis. BMC Bioinformatics 17, 407. doi: 10.1186/s12859-016-1285-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Finnerup, N. B., Haroutounian, S., Kamerman, P., Baron, R., Bennett, D. L., Bouhassira, D., et al. (2016). Neuropathic pain: an updated grading system for research and clinical practice. Pain 157, 1599–1606. doi: 10.1097/j.pain.0000000000000492

PubMed Abstract | CrossRef Full Text | Google Scholar

Friedländer, M. R., Mackowiak, S. D., Li, N., Chen, W., and Rajewsky, N. (2012). miRDeep2 accurately identifies known and hundreds of novel microRNA genes in seven animal clades. Nucleic Acids Res. 40, 37–52. doi: 10.1093/nar/gkr688

PubMed Abstract | CrossRef Full Text | Google Scholar

Gaskin, D. J., and Richard, P. (2012). The economic costs of pain in the United States. J. Pain 13, 715–724. doi: 10.1016/j.jpain.2012.03.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Guo, D., and Hu, J. (2014). Spinal presynaptic inhibition in pain control. Neuroscience 283, 95–106. doi: 10.1016/j.neuroscience.2014.09.032

PubMed Abstract | CrossRef Full Text | Google Scholar

Guttman, M., and Rinn, J. L. (2012). Modular regulatory principles of large non-coding RNAs. Nature 339, 482. doi: 10.1038/nature10887

CrossRef Full Text | Google Scholar

Hargreaves, K., Dubner, R., Brown, F., Flores, C., and Joris, J. (1988). A new and sensitive method for measuring thermal nociception in cutaneous hyperalgesia. Pain 32, 77–88. doi: 10.1016/0304-3959(88)90026-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Heyn, J., Luchting, B., Hinske, L. C., Hübner, M., Azad, S. C., and Kreth, S. (2016). MiR-124a and miR-155 enhance differentiation of regulatory T cells in patients with neuropathic pain. J Neuroinflammation 13, 248. doi: 10.1186/s12974-016-0712-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Jiang, B. C., Sun, W. X., He, L. N., Cao, D. L., Zhang, Z. J., and Gao, Y. J. (2015). Identification of lncRNA expression profile in the spinal cord of mice following spinal nerve ligation-induced neuropathic pain. Mol Pain 11, 43. doi: 10.1186/s12990-015-0047-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Kanehisa, M., Araki, M., Goto, S., Hattori, M., Hirakawa, M., Itoh, M., et al. (2008). KEGG for linking genomes to life and the environment. Nucleic Acids Res. 36, D480–D484. doi: 10.1093/nar/gkm882

PubMed Abstract | CrossRef Full Text | Google Scholar

Langley, P. C., Van Litsenburg, C., Cappelleri, J. C., and Carroll, D. (2013). The burden associated with neuropathic pain in Western Europe. J. Med. Econ. 16, 85–95. doi: 10.3111/13696998.2012.729548

PubMed Abstract | CrossRef Full Text | Google Scholar

Langley, P., Mller-Schwerfe, G., Nicolaou, A., Liedgens, H., Pergolizzi, J., and Varrassi, G. (2010). The societal impact of pain in the European Union: health-related quality of life and healthcare resource utilization. J. Med. Econ. 13, 571–581. doi: 10.3111/13696998.2010.516709

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, H., He, Z., Gu, Y., Fang, L., and Lv, X. (2016). Prioritization of non-coding disease-causing variants and long non-coding RNAs in liver cancer. Oncol Lett. 12, 3987–3994. doi: 10.3892/ol.2016.5135

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, L., Zhang, K., Lee, J., Cordes, S., Davis, D. P., and Tang, Z. (2009). Discovering cancer genes by integrating network and functional properties. BMC Med. Genomics 2:61. doi: 10.1186/1755-8794-2-61

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, S., Yu, B., Wang, S., Gu, Y., Yao, D., Wang, Y., et al. (2012). Identification and functional analysis of novel micro-RNAs in rat dorsal root ganglia after sciatic nerve resection. J. Neurosci. Res. 90, 791–801. doi: 10.1002/jnr.22814

PubMed Abstract | CrossRef Full Text | Google Scholar

Love, M. I., Huber, W., and Anders, S. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 15, 550. doi: 10.1186/s13059-014-0550-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Mao, X., Cai, T., Olyarchuk, J. G., and Wei, L. (2005). Automated genome annotation and pathway identification using the KEGG orthology (KO) as a controlled vocabulary. Bioinformatics 21, 3787–3793. doi: 10.1093/bioinformatics/bti430

PubMed Abstract | CrossRef Full Text | Google Scholar

Mattick, J. S., and Makunin, I. V. (2006). Non-coding RNA. Hum. Mol. Genet. 15, R17–R29. doi: 10.1093/hmg/ddl046

PubMed Abstract | CrossRef Full Text | Google Scholar

Memczak, S., Marvin, J., and Antigoni, E. (2013). Circular RNAs are a large class of animal RNAs with regulatory potency. Nature 495, 333–338. doi: 10.1038/nature11928

PubMed Abstract | CrossRef Full Text | Google Scholar

Molander, C., and Grant, G. (1986). Laminar distribution and somatotopic organization of primary afferent fibers from hindlimb nerves in the dorsal horn. A study by transganglionic transport of horseradish peroxidase in the rat. Neuroscience 19, 297–312. doi: 10.1016/0306-4522(86)90023-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Prescott, S. A., Ma, Q., and De Koninck, Y. (2014). Normal and abnormal coding of somatosensory stimuli causing pain. Nat. Neurosci. 17, 183–191. doi: 10.1038/nn.3629

PubMed Abstract | CrossRef Full Text | Google Scholar

Sakai, A., and Suzuki, H. (2015). microRNA and Pain. Adv. Exp. Med. Biol. 888, 17–39. doi: 10.1007/978-3-319-22671-2_3

PubMed Abstract | CrossRef Full Text | Google Scholar

Shao, J., Cao, J., Wang, J., Ren, X., Su, S., Li, M., et al. (2015). MicroRNA-30b regulates expression of the sodium channel Nav1.7 in nerve injury-induced neuropathic pain in the rat. Mol. Pain 12, 187–193. doi: 10.1177/1744806916671523

PubMed Abstract | CrossRef Full Text | Google Scholar

Thum, T. (2014). Non-coding RNAs and myocardial fibrosis. Nat. Rev. Cardiol. 11, 655–663. doi: 10.1038/nrcardio.2014.125

CrossRef Full Text | Google Scholar

Thum, T., and Condorelli, G. (2015). Long non-coding RNAs and microRNAs in cardiovascular pathophysiology. Circ. Res. 116, 751–762. doi: 10.1161/CIRCRESAHA.116.303549

CrossRef Full Text | Google Scholar

Van Rossum, D., Verheijen, B. M., and Pasterkamp, R. J. (2016). Circular RNAs: novel regulators of neuronal development. Front. Mol. Neurosci. 26:74. doi: 10.3389/fnmol.2016.00074

CrossRef Full Text | Google Scholar

Wang, K. C., and Chang, H. Y. (2011). Molecular mechanisms of long non-coding RNAs. Mol Cell. 43, 904–914. doi: 10.1016/j.molcel.2011.08.018

CrossRef Full Text | Google Scholar

Wang, S., Xu, H., Zou, L., Xie, J., Wu, H., Wu, B., et al. (2016). LncRNA uc.48+ is involved in diabetic neuropathic pain mediated by the P2X₃ receptor in the dorsal root ganglia. Purinergic Signal. 212, 139–148. doi: 10.1007/s11302-015-9488-x

CrossRef Full Text | Google Scholar

Wen, M., Shen, Y., Shi, S., and Tang, T. (2012). miREvo: an integrative microRNA evolutionary analysis platform for next-generation sequencing experiments. BMC Bioinformatics 13:140. doi: 10.1186/1471-2105-13-140

PubMed Abstract | CrossRef Full Text | Google Scholar

Xiong, W., Jiang, Y. X., Ai, Y. Q., Liu, S., Wu, X. R., Cui, J. G., et al. (2015). Microarray analysis of long non-coding RNA expression profile associated with 5-fluorouracil-based chemoradiation resistance in colorectal cancer cells. Asian Pac. J. Cancer Prev. 16, 3395–3402. doi: 10.7314/APJCP.2015.16.8.3395

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, Q., Yang, K., and Li, A. (2014). MicroRNA-21 protects against ischemia-reperfusion and hypoxia-reperfusion-induced cardiocyte apoptosis via the phosphatase and tensin homolog/Akt-dependent mechanism. Mol. Med. Rep. 9, 2213–2220. doi: 10.3892/mmr.2014.2068

PubMed Abstract | CrossRef Full Text | Google Scholar

You, X., Vlatkovic, I., Babic, A., Will, T., Epstein, I., Tushev, G., et al. (2015). Neural circular RNAs are derived from synaptic genes and regulated by development and plasticity. Nat. Neurosci. 18, 603–610. doi: 10.1038/nn.3975

PubMed Abstract | CrossRef Full Text | Google Scholar

Young, M. D., Wakefield, M. J., Smyth, G. K., and Oshlack, A. (2010). Gene ontology analysis for RNA-seq: accounting for selection bias. Genome Biol. 11:R14. doi: 10.1186/gb-2010-11-2-r14

PubMed Abstract | CrossRef Full Text | Google Scholar

Yu, B., Zhou, S., Hu, W., Qian, T., Gao, R., Ding, G., et al. (2013). Altered long non-coding RNA expressions in dorsal root ganglion after rat sciatic nerve injury. Neurosci. Lett. 534, 117–122. doi: 10.1016/j.neulet.2012.12.014

CrossRef Full Text | Google Scholar

Zhang, J., and Banerjee, B. (2015). Role of MicroRNA in visceral pain. J. Neurogastroenterol. Motil. 21, 159–171. doi: 10.5056/jnm15027

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, X., Tang, Z., and Zhang, H. (2013). A long non-coding RNA contributes to neuropathic pain by silencing Kcna2 in primary afferent neurons. Nat. Neurosci. 16, 1024–1031. doi: 10.1038/nn.3438

CrossRef Full Text | Google Scholar

Zhuo, M., Wu, G., and Wu, L. J. (2011). Neuronal and microglial mechanisms of neuropathic pain. Mol. Brain 4, 31–36. doi: 10.1186/1756-6606-4-31

PubMed Abstract | CrossRef Full Text | Google Scholar