Characterization of novel lnc RNAs in the spinal cord of rats with lumbar disc herniation

Introduction

Radicular pain induced by a lumbar herniated disc is an unbearable symptom for patients and has become a public health problem. More and more research have shown that the inflammatory and immune responses occurring in the peripheral^1–3 and central nervous systems⁴ play a critical role in the progression of sciatica induced by lumbar disc herniation (LDH). Central sensitization in the spinal cord dorsal horn is closely involved in pain process in lumbar radiculopathy.^5,6 Recent researches have been concentrated in the molecular targets in pain-signaling systems, such as G-protein-coupled receptors, ion channels in the glial, and neurons localized within the nociceptive pathways.⁷ However, the molecular mechanism underlying radicular pain is largely unknown.

Recent studies found that lncRNAs play roles in regulation of gene expression via the control of cell cycle distribution, cell differentiation, and transcriptional and posttranscriptional regulation, and they are not simply some transcriptional by-products.^8,9 lncRNAs modulate gene transcription by affecting the binding ability of transcription factors and rearranging the chromatin.¹⁰ lncRNAs also affect miRNA functions by controlling pre-mRNA splicing or acting as miRNA sponges.^11,12 Evidence suggests that aberrant expression of lncRNAs occurs in various types of pain, including diabetic neuropathic pain, peripheral nerve injury pain caused by spinal nerve ligation and chronic constriction injury, cancer-associated pain, and osteoarthritis inflammation-related pain.^13–15 Although the expression signatures of many lncRNAs have been defined in the context of pain, expression levels to association and function based on expression levels of lncRNAs in LDH and pain-signaling pathways remain largely unknown

The spinal cord is a crucial area in processing of pain-related signals.^16,17 Here, we analyzed the expression patterns of mRNAs and lncRNAs in the spinal dorsal horn between a group of rats treated with nucleus pulposus (NP) implantation and a sham group. We identified 103 mRNAs and more than 19 unique lncRNAs that were dramatically different in expression in rats with LDH by using microarrays. Some differentially expressed (DE) mRNAs and lncRNAs were further verified by quantitative PCR (qPCR) in six additional experiments. Then we used gene ontology (GO) and pathway analyses to analyze the biologic roles of these mRNAs and lncRNAs. Our study might provide a novel therapeutic target for the treatment of radicular pain in patients with LDH.

Materials and methods

Animals

Experiments were performed on adult male Sprague-Dawley rats (body weight 220±20 g), which were approved by Soochow University. The total rats used in the present study were 68. These rats were housed under standard conditions (07:00–19:00 lighting, 24°C±2°C) with a unified laboratory diet and water. Experimental procedure of rats was performed according to the guidelines of the International Association for the Study of Pain. Because the paw withdraw threshold and WBD were at the maximum different point of the time–response curve (Figure 1A, C) and in order to minimize the suffering from pain of rats, almost all of the experiments were performed within 7 days after NP treatment.

Figure 1 LDH induces persistent neuropathic pain of rat hindpaw. Notes: NP treatment produced persistent mechanical allodynia (A) and heat hyperalgesia (B). In rats treated with NP, the body weight-bearing difference was dramatically increased between hindlimbs at 3, 7, 14, 21, and 28 days after surgery compared to the sham group (C). NP-treated rats did not produce any effect on the time to fall for LDH rats in the rotarod test (D). Data are expressed as mean ± SEM (n=6 for each group). *P<0.05, two-way repeated-measures ANOVA. Abbreviations: LDH, lumbar disc herniation; NP, nucleus pulposus; PWT, paw withdrawal threshold; WB, weight borne.

Disc herniation model

Surgeries for the LDH model were performed as previously described.¹ Briefly speaking, all experimental rats were deeply anesthetized with sodium pentobarbital (50 mg/kg body weight, intraperitoneally). The hair of the rat’s lower back was shaved, and the skin was regularly disinfected by iodophor. NP was harvested from the coccygeal intervertebral disc level of each tail. Harvested NP (~5 mg) was implanted next to the left L5 and L6 nerve roots just proximal to the corresponding DRG. Nonoperation was did on the right side of the dorsal root in all rats. The surgical procedure for the sham group was similar to the NP-treated group, but NP implantation was omitted.

Pain behavior test

Rats (n=8 for each group) were measured for mechanical pain sensitivity in plantar surface of the hindpaws at 3 days prior to surgery and 3, 7, 14, 21, 28, and 35 days after surgery in a double-blinded manner. Mechanical withdrawal thresholds (PWT) and thermal threshold were tested as previously described.^1,3 In brief, mechanical withdrawal thresholds (PWT) were tested utilizing von Frey filaments (VFF). A series of calibrated VFF were applied to the plantar surface of the rat hind paw with sufficient force to bend the filaments for 10 seconds or until the rat withdrew. A filament of the next lower force was applied when a response occurred. The cutoff strength of VFF was set at 20.30 g to avoid the potential injury of the tissue. A 50% likelihood of stimulus production was determined by the “up-down” calculating assay. The average value of VFF force was determined as a withdrawal response. Each test was duplicated 3 times at a 2-minute interval.

Body weight-bearing test

WBD test was performed by an incapacitance tester (PH-200, Taimeng Chengdu, China) as described previously.^1,3 Results were determined as averages of three trials, with each trial measuring the hindlimb weight over 5 seconds. WBD is presented as weight on the contralateral limb minus weight on the ipsilateral limb. When one hindpaw has pain, the rat will not put more body weight to the leg. Therefore, the other leg will bear more body weight. The greater the body WBD is, the more pain the rat has. Thus, the weight-bearing asymmetry is indicative of hyperalgesia. All rats were tested 3 days before the operation and 3–35 after implantation of NP.

Rotarod test

The rotarod system (ZH-300, Zhenghua, Anhui Province, China) for locomotor assessment was used to measure the time period for an animal to maintain its balance on a moving cylinder as previously described.¹ Thirty minutes after the conditioning, each animal was placed on the rotarod and its latency until falling was determined and expressed in seconds.

Tissue extraction and RNA isolation

Three rats from LDH treated and three from sham group were anesthetized with isoflurane and perfused through the ascending aorta with saline at 7 days after operation. The L5–L6 spinal cord segments were then extracted. Total RNAs were collected from the dorsal horn of spinal cord using Trizol reagent (Invitrogen, Carlsbad, CA, USA) based on the manufacturer’s instructions. The RNA concentration and purity were examined using the NanoDrop 1000 Spectrophotometer (Thermo). Equal mRNA from one rat with the same treatment was mixed as one sample after these testing. Total of six samples (three from LDH and three from sham group) were selected for microarray analysis.

lncRNA and mRNA microarray expression profiling

The microarray profiling was tested in the research laboratory of the OE Biotechnology Company in Shanghai (China). The sample labeling, microarray hybridization, and washing were handled according to the manufacturer’s standard protocols. In brief, mRNA was purified from total RNA by using an mRNA-only Eukaryotic mRNA Isolation Kit (Epicentre Biotechnologies, Madison, WI, USA). The RNA samples were hybridized to mouse lncRNA microarray. After washing, the arrays were scanned with the Agilent Scanner G2505C (Agilent Technologies). Array images and raw data were analyzed by using Extraction software (version 10.7.1.1, Agilent Technologies). Data were then normalized with the quantile algorithm. The probes with at least one condition out of two conditions flagged as “P” were chosen for further data analysis. Then, DE mRNAs and lncRNAs were identified with P-values and FC as calculated with t-test. The FC ≥1.5 and P-value <0.05 were set for the up- and downregulated gene threshold. Hierarchical clustering was performed to determine the different mRNA and lncRNA expression patterns among these samples.

Bioinformatics analysis

Volcano plot filtering identified DE lncRNAs and mRNAs with statistical significance. The threshold used to screen RNAs with an FC >1.5 was P<0.05. Hierarchical clustering was carried out by Cluster 3.0, and the heat maps were generated in Java Treeview. The DE mRNAs that were adjacent to or overlapping with the DE lncRNAs were recognized as DE lncRNAs associated with mRNAs using the UCSC Genome Browser.

We first calculated coexpressed mRNAs for each DE lncRNA, and then we conducted a functional enrichment analysis of this set of coexpressed mRNAs. The enriched functional terms were used as the predicted functional terms of given lncRNA. The coexpressed mRNAs of lncRNAs were identified by calculating Pearson correlation with a correlation P-value <0.05. We then used hypergeometric cumulative distribution function to calculate the enrichment of functional terms in annotation of coexpressed mRNAs. GO analysis and KEGG analysis were applied to determine the biologic roles of these DE mRNAs, based on the latest KEGG (database (http://www.genome.jp/kegg/). The P-value (hypergeometric P-value) denotes the significance of the pathway correlated to the conditions. The recommended P-value cutoff is 0.05.

RT-PCR

The microarray results were further confirmed by RT-PCR. The yield of RNA was determined using a NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA), and the integrity was evaluated using agarose gel electrophoresis stained with ethidium bromide. Quantification was performed with a two-step reaction process as described previously. The primer sequences were designed in the laboratory and synthesized by Generay Biotech (Generay, Shanghai, China) based on the mRNA sequences obtained from the NCBI database as follows: BMP4 (forward 5′-CCTGGTCAACTCCGTTAATTC-3′; reverse 5′-TTCAACACCACCTTGTCG-3′), SLC17a8 (forward 5′-GCTGTGTCATGTGTGTGAG-3′/reverse 5′-AGAGGTTGTGGCTAGACGA-3′), Cacnals (forward 5′-CACCAATAAGATCCGCGTC-3′/reverse 5′-CTGATTCCTCATGGAGTCG-3′).

Statistical analysis

Statistical analyses were carried out by using SPSS 17.0 software (SPSS Inc., Chicago, IL, USA). All data are expressed as the mean ± SEM or proportions where needed. One-way ANOVA followed by Tukey’s multiple comparison tests, paired t-test, or unpaired t-test was used for comparisons. GraphPad Prism 5.0 (GraphPad Software, San Diego, CA, USA) was used to manage all graphs. P-values <0.05 (two-tailed) were considered statistically significant.

Results

Pain hypersensitivity induced by transplantation of autologous NP

Autologous NP was harvested from the tails of rats and implanted in the lumbar 5 and lumbar 6 spinal nerve roots, which caused pain hypersensitivity in the rats when compared with the sham group. In the LDH group, the mechanical paw withdrawal threshold (PWT) markedly reduced at 3 days after surgery, which indicated mechanical pain allodynia. The change returned to the normal level on day 35 after implantation (*P<0.05, Figure 1A). The NP-treated rats exhibited a remarkable decrease in the thermal paw withdrawal latency (PWL) compared to the sham group. The PWL was decreased at 7 days after NP treatment and returned to the baseline level at day 14 (Figure 1B, n=6 for each group, *P<0.05). Compared to the sham group, the LDH group also demonstrated a dramatic increase in body weight-bearing difference (WBD). The WBD was increased on day 3 after the NP implantation and returned to the baseline level on day 35 (Figure 1C, *P<0.05). We performed rotarod tests to determine whether the motor performance was affected by the surgery. There was no significant difference in the time to fall on the rotarod at a speed of 20 rpm (Figure 1D, P>0.05).

mRNAs and lncRNAs expression profiles in LDH

We next examined the expression profiles of the mRNAs and lncRNAs in the L5–L6 spinal cord at 7 days after LDH or sham by microarray. From the hierarchical clustering and scatter plot analysis results, we obtained an overview of the expression signatures of the mRNAs and lncRNAs. The scatter plots demonstrated that large numbers of the lncRNAs and mRNAs were DE between the LDH and the sham group (Figure 2A, B). Hierarchical cluster analysis of the lncRNAs and mRNAs demonstrated that the three LDH or three sham samples were clustered together, and the signal intensity was consistent in the LDH and sham groups (Figure 2C, D). The heatmap of the DE lncRNAs or mRNAs that were upregulated or downregulated by 1.5-fold was magnified (Figure 2E, F), which indicated concordance in the LDH and sham groups. These differential alterations in the mRNAs and lncRNAs in the spinal cord are correlated with radicular pain.

Figure 2 LDH results in the expression profiling changes of mRNA and lncRNA in spinal cord dorsal horn. Notes: Scatter plot comparing global mRNA (A) or lncRNA (B) gene expression profiles in the spinal cord between the LDH and sham groups. Green lines indicate 1.5-fold differences in either direction in mRNA and lncRNA expression. Heat map showing hierarchical clustering of overall lncRNAs (C) or mRNA (D) expression pattern of reliably measured probe sets. Heat map showing hierarchical clustering of lncRNAs (E) or mRNA (F), whose expression changes were more than 1.5-fold. In clustering analysis, upregulated and downregulated genes were colored in red and green, respectively. Abbreviation: LDH, lumbar disc herniation.

DE lncRNAs and mRNAs

To further analyze the DE lncRNAs, we used the significance analysis of microarrays method, following the criteria P-value <0.05 and fold change (FC) >1.5. The results revealed that 19 lncRNAs, which comprised 13 upregulated and 6 downregulated lncRNAs, were greatly altered in the LDH group compared to the sham group. The most enhanced lncRNAs were TCONS_00043124, TCONS_00036317, TCONS_00014238, TCONS_00089375, TCONS_00031848, and TCONS_00095748. The most reduced lncRNAs were TCONS_00040686, XR_086254.1, TCONS_00122314, TCONS_00135380, TCONS_00121583, and TCONS_00036854. Of the significantly altered lncRNAs, TCONS_00043124 was the most upregulated, with an FC of 2.661, whereas TCONS_00040686 was the most downregulated, with an FC of 1.777. The detailed information of the DE lncRNAs is shown in Table 1.

Table 1 The details of the differentially expressed lncRNAs Abbreviation: FC, fold change.

In the DE mRNAs, there were 103 genes that exhibited greater than a 1.5-FC. The number of downregulated mRNAs (63) was greater than the number of upregulated mRNAs (40) in LDH. These DE mRNAs corresponded to many known genes involved in pain processing, including Cxcl13, Slc17a8,¹⁸ Ccr5,¹⁹ Cacna1s, P2ry2, and Tnnc2.²⁰ Detailed information regarding the top 20 upregulated and 20 downregulated mRNAs is shown in Table 2.

Table 2 The details of the top 20 upregulated and 20 downregulated mRNAs

Real-time qPCR confirmation of mRNA and lncRNA expression

To confirm the reliability of the microarray data and analyze the temporal alterations in mRNA and lncRNA expression after LDH, the lncRNAs TCONS_00031848 and TCONS_00122118 and the mRNAs BMP4, SLC17a8, and Cacnals were randomly selected and analyzed by qPCR. The spinal cord tissues were collected from the sham and LDH rats before operation and 3, 7, 14, 21, 28, and 35 days post operation. BMP4 (bone morphogenetic protein 4), a protein-coding gene associated with the TGF-beta-signaling pathway, was significantly decreased at 7 days and persisted until 14 days post operation (Figure 3A). SLC17a8 (solute carrier family 17 member 8, also known as Vglut3) exhibits glutamate transporter activity and may play a role in glutamate transport in the synaptic vesicles of cholinergic and serotoninergic neurons. SLC17A8 was significantly increased at 7 days and until 14 days post operation (Figure 3B). The lncRNAs TCONS_00031848 and TCONS_00122118 were both significantly increased at 7 and 14 days post operation (Figure 3C, D). In addition, the FCs of these lncRNAs and mRNAs detected by qPCR in the LDH tissue at 7 days post operation were consistent with the results from the microarray (Figure 3E), which further supported the reliability of the array data.

Figure 3 Real-time quantitative PCR validations of two mRNAs and two lncRNAs in the spinal cord from NP-treated rats. Notes: The expressions of mRNA BMP4 were significantly downregulated at 7 and 14 days after LDH (A), mRNA SLC17a8 (B), lncRNA TCONS_00031848 (C), and lncRNA TCONS_00122118 (D) were significantly upregulated at 7 and 14 days after LDH. The fold changes of these lncRNAs and mRNAs at 7 days after NP treatment were consistent with the results from the microarray (E). One-way ANOVA followed by Tukey’s multiple-comparison test. *P<0.05. Abbreviations: BMP, bone morphogenetic protein; LDH, lumbar disc herniation; NP, nucleus pulposus.

Functional prediction of DE mRNAs and lncRNA in LDH

To investigate the molecular mechanisms in radicular pain, we conducted GO and pathway analysis of the DE genes in LDH vs sham. The GO results demonstrated that the most significantly enriched biologic processes of the upregulated genes in LDH were chemotaxis, immune response, and inflammatory responses (Figure 4A). The most significantly enriched cellular components of the upregulated genes in LDH were extracellular space, apical dendrite, and extrinsic component of plasma membrane (Figure 4B). The most significantly enriched molecular functions of the upregulated genes in LDH were chemokine activity and CXCR chemokine receptor activity (Figure 4C). The most significantly enriched biologic processes of the downregulated genes in LDH were negative regulation of phosphorylation, membrane depolarization during action potential, positive regulation of the BMP-signaling pathway, and negative regulation of mitogen-activated protein kinase activity (Figure 4D). The most significantly enriched cellular components of the downregulated genes in LDH were the extracellular region, the anchored component of the membrane, and the nucleolar ribonuclease P complex (Figure 4E). The most significantly enriched molecular functions of the downregulated genes in LDH were co-receptor binding, IL-17 receptor activity, glycerophosphodiester phosphodiesterase activity, and calcium-dependent phospholipase A2 activity (Figure 4F).

Figure 4 Biologic functions of differentially expressed mRNAs. Notes: The significant biologic process, cellular component, and molecular function of upregulated mRNAs (A–C). The significant biologic process, cellular component, and molecular function of downregulated mRNAs (D–F). Abbreviations: BMP, bone morphogenetic protein; MAP, mitogen-activated protein.

Similarly, we analyzed the DE genes in Kyoto Encyclopedia of Genes and Genomes (KEGG). The results demonstrated that the upregulated genes in LDH were involved in the chemokine-signaling pathway, Staphylococcus aureus infection, cytokine–cytokine receptor interaction, glycosphingolipid biosynthesis, the toll-like receptor-signaling pathway, and the NF-kappa B-signaling pathway (Figure 5A). The downregulated genes in LDH were involved in ether lipid metabolism, basal cell carcinoma, Hippo-signaling pathway, neuroactive ligand–receptor interaction, and calcium-signaling pathway (Figure 5B).

Figure 5 Pathway analysis for differentially expressed mRNAs with fold changes >1.5. Notes: The significant pathways for upregulated genes in NP-treated group (A). The significant pathways for downregulated genes in NP-treated group (B). Abbreviations: KEGG, Kyoto Encyclopedia of Genes and Genomes; NP, nucleus pulposus.

Function prediction of mRNAs and lncRNAs

As some lncRNAs have been reported to play critical roles in regulation of the expression of their neighboring genes, we conducted additional screening for DE lncRNAs associated with mRNAs. Briefly, we identified and performed a functional enrichment analysis of the coexpressed mRNAs for each different lncRNA. The coexpressed mRNAs of the lncRNAs were identified by calculating the Pearson correlation with a correlation P-value <0.05. We performed the GO and pathway analysis with the top 200 DE mRNAs and lncRNAs to investigate the potential biologic relationships. GO and pathway analyses suggested that some functional pathways were identified. Among these pathways, voltage-gated potassium channel activity, BMP receptor binding, sensitive calcium-release channel activity, and RNA polymerase II regulatory region sequence-specific DNA binding were the most closely associated with LDH (Figure 6A). Additionally, a KEGG pathway analysis demonstrated associations with some pathways, including sphingolipid metabolism, cell cycle, oxytocin-signaling pathway, glutamatergic synapse, and Hippo-signaling pathway (Figure 6B).

Figure 6 Function prediction of lncRNAs and mRNAs. Notes: Molecular function enrichment analysis of differentially expressed lncRNA-related mRNAs. Gene ontology analysis indicated that several functional pathways were enriched (A). KEGG pathway analysis showed associations with some pathways (B). Abbreviations: BMP, bone morphogenetic protein; KEGG, Kyoto Encyclopedia of Genes and Genomes; PPAR, peroxisome proliferator-activated receptor.

Discussion

LDH caused radicular pain that affects the nervous system.^21,22 Evidence has indicated that gene expression plays important roles in the initiation and development of neuropathic pain.^23,24 In recent years, the underlying molecular mechanisms of radicular pain have been widely studied; however, the epigenetic process of pain is still unclear.²⁵ lncRNAs have complicated molecular functions, including regulating protein activities, serving structural or organizational roles, serving as precursors to small RNAs, modulating transcriptional patterns, and altering RNA-processing events.²⁶ The function of lncRNAs in regulating the gene expression of neuropathic pain is unclear. In this study, we identified the DE lncRNAs in the spinal cord under a radicular pain condition and analyzed their characteristics and possible relationships with mRNAs.²⁷ A total of 9,984 lncRNAs were found in the spinal cord of the LDH rats. Among these lncRNAs, 13 lncRNAs were upregulated, and 6 lncRNAs were downregulated at 7 days after the NP treatment. Additionally, qRT-PCR was used to confirm the microarray analysis results, and findings were consistent with the microarrays data of spinal cords from the LDH and sham groups, which suggested that lncRNAs might be involved in radicular pain processing. These findings might provide novel insights into the molecular mechanism of neuropathic pain.

Among the DE mRNAs, the number of upregulated mRNAs did not exceed the number of downregulated mRNAs in the LDH spinal cords, which reveals the new pathways and biologic processes in neuropathic pain. Many pain-related genes, including Cacna1s, Cxcl13, Cxcl12, Ccr5, Mki67, Asb15, and Tnnc2, were significantly enhanced after LDH. In addition, many other mRNAs, such as BMP4, Nkd1, Fzd6, Slc17a8, and Gpr35, with unclear functions in the spinal cord (especially in neuropathic pain), were also screened out. As the DE genes may be related to nerve damage and glial cell activation to immune and inflammatory responses, further investigations are needed to confirm whether they are involved in pain hypersensitivity.

According to the GO and KEGG enrichment analyses of the DE mRNAs, we showed that the significantly enriched biologic processes and molecular functions of the upregulated genes in LDH vs sham were mainly involved in the chemokine signaling pathway,^28,29 inflammation,³⁰ and immune response.³¹ These findings are consistent with previous studies showing that neuroinflammation³² and the inflammatory mediators produced in the peripheral and central nervous system (CNS) play a very important role in hyperalgesia.^33,34 From the significant pathway analyses of the DE genes, we showed some significantly enriched pathways of the upregulated genes in LDH vs sham, such as the cytokine–cytokine receptor interaction, the toll-like receptor-signaling pathway, and the NF-kappa B-signaling pathway. Indeed, Tlr4 and NF-kappa B have been indicated as potential mechanisms in neuropathic pain and other pain models.³⁵ Collectively, the results suggest that anti-inflammation may be an effective therapeutic target for radicular pain.

Although lncRNAs play key roles in the regulation of gene expression, the major functional principles and regulatory mechanisms of lncRNAs are complex and obscure. Unlike miRNAs, there are no common approaches that can be used to predict the target genes and functions of lncRNAs by their sequence information or secondary structure. A large number of studies suggest that a number of lncRNAs can activate or repress the expression of their neighboring or overlapping genes.⁸ Our further analysis showed that an upregulated lncRNA, TCONS_00089375, regulates the expression of a neighboring gene, Cstf3, through a cis transcriptional mechanism, and the downregulated lncRNA, TCONS_00036854, regulates a neighboring gene, Ppp1r15b, through a cis transcriptional mechanism in the LDH model. In the spinal cord after LDH, TCONS_00031848 was found to be associated with the Hippo-signaling pathway and influenced BMP4, Bmp2, Fzd3, Ndk1, Wnt16, and Wnt8a signaling. Of note is that our present study showed that autologous NP application to the spinal nerve produced persistent mechanical pain hypersensitivity. However, the thermal hypersensitivity is not remarkable. It is transient and not as prominent as the mechanical data. The mechanisms remain unknown. Future investigation of roles of DE genes in mechanical and thermal hypersensitivity is definitely helpful.

In summary, our results demonstrated that lncRNA transcripts were DE in the spinal cord after LDH. A total of 19 DE lncRNAs were shown to have neighboring or overlapping mRNAs in the genome. These lncRNAs may regulate the expression of their associated proteins and genes and play important roles in the pathogenesis of radicular pain. Future study of the lncRNAs found in our investigations will likely focus on their functions and their relationship with LDH. Nevertheless, our study might provide important information for exploring potential targets for the treatment of neuropathic pain.

Acknowledgments

This study was supported by grants from the Science and Technology Project of Suzhou (SYS201544), Livelihood Technology of Suzhou Science and Technology Bureau (Basic research on medical and health, SYS201731), and Natural Science Research of Jiangsu Higher Education Institutions (17KJB320012).

Author contributions

QW performed experiments, analyzed data, prepared figures, and drafted the manuscript. HA performed experiments, analyzed data, prepared figures, and drafted the manuscript. JL performed experiments and analyzed data. MX performed experiments and analyzed data. ZZ performed experiments. CQ analyzed data and drafted the manuscript. YX analyzed data and drafted the manuscript. JY designed and supervised the experiments and edited the manuscript. All authors read and approved the final manuscript. All authors contributed to data analysis, drafting or revising the article, gave final approval of the version to be published, and agree to be accountable for all aspects of the work.

Disclosure

The authors report no conflicts of interest in this work.

References

1.		Wang Q, Zhu H, Zou K, et al. Sensitization of P2Χ3 receptors by cystathionine β-synthetase mediates persistent pain hypersensitivity in a rat model of lumbar disc herniation. Mol Pain. 2015;11:15.
2.		Yan J, Hu S, Zou K, et al. Inhibition of cystathionine β-synthetase suppresses sodium channel activities of dorsal root ganglion neurons of rats with lumbar disc herniation. Sci Rep. 2016;6(1):38188.
3.		Yan J, Zou K, Liu X, et al. Hyperexcitability and sensitization of sodium channels of dorsal root ganglion neurons in a rat model of lumber disc herniation. Eur Spine J. 2016;25(1):177–185.
4.		Zhu X, Cao S, Zhu MD, Liu JQ, Chen JJ, Gao YJ. Contribution of chemokine CCL2/CCR2 signaling in the dorsal root ganglion and spinal cord to the maintenance of neuropathic pain in a rat model of lumbar disc herniation. J Pain. 2014;15(5):516–526.
5.		Stemkowski PL, Smith PA. Sensory neurons, ion channels, inflammation and the onset of neuropathic pain. Can J Neurol Sci. 2012;39(4):416–435.
6.		Deliu E, Brailoiu GC, Arterburn JB, et al. Mechanisms of G protein-coupled estrogen receptor-mediated spinal nociception. J Pain. 2012;13(8):742–754.
7.		Kong X, Wei J, Wang D, et al. Upregulation of spinal voltage-dependent anion channel 1 contributes to bone cancer pain hypersensitivity in rats. Neurosci Bull. 2017;33(6):711–721.
8.		Jiang BC, Sun WX, He LN, Cao DL, Zhang ZJ, Gao YJ. Identification of lncRNA expression profile in the spinal cord of mice following spinal nerve ligation-induced neuropathic pain. Mol Pain. 2015;11:43.
9.		Koch L. Functional genomics: screening for lncRNA function. Nat Rev Genet. 2017;18(2):70.
10.		Zhang HY, Yang W, Zheng FS, Wang YB, Lu JB. Long non-coding RNA SNHG1 regulates zinc finger E-box binding homeobox 1 expression by interacting with TAp63 and promotes cell metastasis and invasion in lung squamous cell carcinoma. Biomed Pharmacother. 2017;90:650–658.
11.		Wang GQ, Wang Y, Xiong Y, et al. Sirt1 AS lncRNA interacts with its mRNA to inhibit muscle formation by attenuating function of miR-34a. Sci Rep. 2016;6(1):21865.
12.		Abidin SZ, Leong JW, Mahmoudi M, et al. In silico prediction and validation of Gfap as an miR-3099 target in mouse brain. Neurosci Bull. 2017;33(4):373–382.
13.		Zhang M, Gu H, Xu W, Zhou X. Down-regulation of lncRNA MALAT1 reduces cardiomyocyte apoptosis and improves left ventricular function in diabetic rats. Int J Cardiol. 2016;203:214–216.
14.		Zhao X, Tang Z, Zhang H, et al. A long noncoding RNA contributes to neuropathic pain by silencing Kcna2 in primary afferent neurons. Nat Neurosci. 2013;16(8):1024–1031.
15.		Liu S, Zou L, Xie J, et al. lncRNA NONRATT021972 siRNA regulates neuropathic pain behaviors in type 2 diabetic rats through the P2X7 receptor in dorsal root ganglia. Mol Brain. 2016;9(1):44.
16.		Peng H, Zou L, Xie J, et al. lncRNA NONRATT021972 siRNA decreases diabetic neuropathic pain mediated by the P2X3 receptor in dorsal root ganglia. Mol Neurobiol. 2017;54(1):511–523.
17.		Wang S, Xu H, Zou L, et al. lncRNA uc.48+ is involved in diabetic neuropathic pain mediated by the P2X3 receptor in the dorsal root ganglia. Purinergic Signal. 2016;12(1):139–148.
18.		Ruel J, Emery S, Nouvian R, et al. Impairment of SLC17A8 encoding vesicular glutamate transporter-3, VGLUT3, underlies nonsyndromic deafness DFNA25 and inner hair cell dysfunction in null mice. Am J Hum Genet. 2008;83(2):278–292.
19.		Sun S, Chen D, Lin F, et al. Role of interleukin-4, the chemokine CCL3 and its receptor CCR5 in neuropathic pain. Mol Immunol. 2016;77:184–192.
20.		Akgün E, Javed MI, Lunzer MM, et al. Inhibition of inflammatory and neuropathic pain by targeting a mu opioid receptor/chemokine receptor5 heteromer (MOR-CCR5). J Med Chem. 2015;58(21):8647–8657.
21.		Suri P, Rainville J, Hunter DJ, Li L, Katz JN. Recurrence of radicular pain or back pain after nonsurgical treatment of symptomatic lumbar disk herniation. Arch Phys Med Rehabil. 2012;93(4):690–695.
22.		Imai S, Ikegami D, Yamashita A, et al. Epigenetic transcriptional activation of monocyte chemotactic protein 3 contributes to long-lasting neuropathic pain. Brain. 2013;136(Pt 3):828–843.
23.		Lee SE, Kim JH. Involvement of substance P and calcitonin gene-related peptide in development and maintenance of neuropathic pain from spinal nerve injury model of rat. Neurosci Res. 2007;58(3):245–249.
24.		Liu X, Liu H, Xu S, et al. Spinal translocator protein alleviates chronic neuropathic pain behavior and modulates spinal astrocyte-neuronal function in rats with L5 spinal nerve ligation model. Pain. 2016;157(1):103–116.
25.		Gao YJ, Zhang L, Samad OA, et al. JNK-induced MCP-1 production in spinal cord astrocytes contributes to central sensitization and neuropathic pain. J Neurosci. 2009;29(13):4096–4108.
26.		Tian P, Li ZJ, Fu X, Ma XL. Role of interleukin-17 in chondrocytes of herniated intervertebral lumbar discs. Exp Ther Med. 2015;10(1):81–87.
27.		Jurga AM, Rojewska E, Piotrowska A, et al. Blockade of toll-like receptors (TLR2, TLR4) attenuates pain and potentiates buprenorphine analgesia in a rat neuropathic pain model. Neural Plast. 2016;2016(3):1–12.
28.		Chen S, Xiong J, Zhan Y, Liu W, Wang X. Wogonin inhibits LPS-induced inflammatory responses in rat dorsal root ganglion neurons via inhibiting TLR4-MyD88-TAK1-mediated NF-κB and MAPK signaling pathway. Cell Mol Neurobiol. 2015;35(4):523–531.
29.		Yuan B, Liu D, Liu X. Spinal cord stimulation exerts analgesia effects in chronic constriction injury rats via suppression of the TLR4/NF-κB pathway. Neurosci Lett. 2014;581:63–68.
30.		Fritah S, Niclou SP, Azuaje F. Databases for lncRNAs: a comparative evaluation of emerging tools. RNA. 2014;20(11):1655–1665.
31.		Lai F, Orom UA, Cesaroni M, et al. Activating RNAs associate with mediator to enhance chromatin architecture and transcription. Nature. 2013;494(7438):497–501.
32.		Paraskevopoulou MD, Hatzigeorgiou AG. Analyzing miRNA-lncRNA interactions. Methods Mol Biol. 2016;1402:271–286.
33.		Jalali S, Bhartiya D, Lalwani MK, Sivasubbu S, Scaria V. Systematic transcriptome wide analysis of lncRNA-miRNA interactions. PLoS One. 2013;8(2):e53823.
34.		Kuang L, Lei M, Li C, et al. Identification of long non-coding rnas related to skeletal muscle development in two rabbit breeds with different growth rate. Int J Mol Sci. 2018;19(7):2046.
35.		Li L, Guo J, Chen Y, Chang C, Xu C. Comprehensive circRNA expression profile and selection of key circRNAs during priming phase of rat liver regeneration. BMC Genomics. 2017;18(1):80.