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

BACKGROUND: Radicular pain, caused by a lesion or autologous nucleus pulposus (NP) implantation, is associated with alteration in gene expression of the pain-signaling pathways. lncRNAs have been shown to play critical roles in neuropathic pain. However, the mechanistic function of lncRNAs in lumbar disc herniation (LDH) remains largely unknown. Identifying different lncRNA expression under sham and NP-implantation conditions in the spinal cord is important for understanding the molecular mechanisms of radicular pain.
METHODS: LDH was induced by implantation of autologous nucleus pulposus (NP), harvested from rat tail, in lumbar 5 and 6 spinal nerve roots. The mRNA and lncRNA microarray analyses demonstrated that the expression profiles of lncRNAs and mRNAs between the LDH and sham groups were markedly altered at 7 days post operation. The expression patterns of several mRNAs and lncRNAs were further proved by qPCR.
RESULTS: LDH produced persistent mechanical and thermal hyperalgesia. A total of 19 lncRNAs was differentially expressed (>1.5-folds), of which 13 was upregulated and 6 was downregulated. In addition, a total of 103 mRNAs was markedly altered (>1.5-folds), of which 40 was upregulated and 63 downregulated. Biological analyses of these mRNAs further demonstrated that the most significantly upregulated genes in LDH included chemotaxis, immune response, and positive regulation of inflammatory responses, which might be important mechanisms underlying radicular neuropathic pain. These 19 differentially expressed lncRNAs have overlapping mRNAs in the genome, which are related to glutamatergic synapse, cytokine-cytokine receptor interaction, and the oxytocin-signalling pathway.
CONCLUSION: Our findings revealed the alteration of expression patterns of mRNAs and lncRNAs in the spinal cord of rats in a radicular pain model of LDH. These mRNAs and lncRNAs might be potential therapeutic targets for the treatment of radicular pain.

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 peripheral1–3 and central nervous systems4 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.7 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 byproducts.8,9 lncRNAs modulate gene transcription by affecting the binding ability of transcription factors and rearranging the chromatin.10 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 micro-arrays. 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.1 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.1 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′-CCTGGTCAACTCC-GTTAATTC-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, >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, T C O N S _ 0 0 0 1 4 2 3 8, T C O N S _ 0 0 0 8 9 3 7 5 , 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

Source lD	P-value	FC (abs)	Regulation
Upregulated lncRNAs
TCONS 00043124	0.043	2.661	Up
TCONS 00036317	0.029	1.768	Up
TCONS 00014238	0.006	1.752	Up
TCONS 00089375	3.31E–04	1.709	Up
TCONS 00031848	0.003	1.645	Up
TCONS 00095748	0.009	1.635	Up
TCONS 00122118	0.002	2.602	Up
TCONS 00050883	8.14E–04	1.544	Up
TCONS 00106609	0.033	1.541	Up
TCONS 00113822	0.031	1.541	Up
TCONS 00121610	0.034	1.536	Up
TCONS 00091440	0.014	1.530	Up
uc.253+	0.009	1.504	Up
Downregulated lncRNAs
TCONS 00040686	0.010	1.777	Down
XR 086254_1	0.032	1.664	Down
TCONS 00122314	0.036	1.655	Down
TCONS 00135380	0.050	1.593	Down
TCONS 00121583	0.005	1.549	Down
TCONS 00036854	0.049	1.548	Down

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,18 Ccr5,19 Cacna1s, P2ry2, and Tnnc2.20 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

Probe name	Gene symbol	FC (abs)	P-value	Regulation	Description
Upregulated genes
A_64_P013437	Cxcl13	11.130	0.012	Up	C-X-C motif chemokine ligand 13
A_64_P050430	LOC679608	5.000	0.010	Up	Hypothetical protein LOC679608
A_64_P052265	Klral	3.548	0.031	Up	Killer cell lectin-like receptor, subfamily A, member 1
A_64_P013753	LOC68 5411	3.228	0.048	Up	Similar to spermatogenesis-associated glutamate (E)-rich protein 4b
A_64_P019790	CYP2d1	2.812	0.036	Up	Cytochrome P450, family 2, subfamily d, polypeptide 1
A_44_P102353S	C3	2.709	0.048	Up	Complement component 3
A_64_P150326	LOC688649	2.364	0.040	Up	Similar to spermatogenesis-associated glutamate (E)-rich protein 4d
A_64_P140936	Zfp39	2.219	0.402	Up	Zinc finger protein 39
A_64_P083999	RGD1564933	2.268	0.026	Up	Similar to olfactory receptor O1r271
A_44_P807861	Six4	1.926	0.012	Up	SIX homeobox 4
A_64_P107991	Gpr31	1182	0.040	Up	G-protein-coupled receptor 31
A_64_P004S54	Ptafr	1.834	0.029	Up	Platelet-activating factor receptor
A_64_P012557	Fam64a	1.808	0.022	Up	Family with sequence similarity 64, member A
A_44_P294914	Sec1	1.214	0.009	Up	Secretory blood group 1
A_44_P34S999	Asb15	1.511	0.045	Up	Ankyrin repeat and SOCS box containing 15
A_64_P135843	Grb7	1.708	0.041	Up	Growth factor receptor bound protein 7
A_PO_47636	S1c17a8	1.619	0.011	Up	Solute carrier family 17 member 8
A_42_P797117	RGD1311429	1.667	0.025	Up	KAT8 regulatory NSL complex subunit 1
A_44_P113974	CCr5	1.666	0.044	Up	Chemokine (C-C motif) receptor 5
A_64_P021440	Slc10a6	1.645	0.041	Up	Solute carrier family 10 member 6
Downregulated genes
A_44_P236197		6gi42	0.011	Down	IL-induced protein with tetratricopeptide repeats 1
A_64_P057196	Cldn14	4.658	0.022	Down	Claudin 14
A_64_P011679	Caznals	2.995	0.035	Down	Calcium voltage-gated channel subunit alpha l S
A_64_P029118	Loxl4	2.368	0.039	Down	Lysyl oxidase-like 4
A_44_P314087	LOC681186	2.346	0.011	Down	Hypothetical protein LOC681186
A_44_P520159	P2ry2	2.312	0.001	Down	Purinergic receptor P2Y2
A_64_P158g22	Ptgdrl	2.309	0.021	Down	Prostaglandin D2 receptor-like
A_64_P115212	I117rb	2.183	0.015	Down	IL 17 receptor B
A_64_P003051	01r5	2.159	0.023	Down	Olfactory receptor 5
A_64_P0930SS	01r1206	2.151	0.011	Down	Olfactory receptor 1206
A_64_P074305	01r440	2.144	0.007	Down	Olfactory receptor 440
A_43_P16403	Clca1	2.032	3.50E–04	Down	Chloride channel accessory 1
A_44_P374618	Cdknlc	1.915	0.031	Down	Cyclin-dependent kinase inhibitor 1C
A_64_P1650S5	Cdknlc	1.955	0.018	Down	Cyclin-dependent kinase inhibitor 1C
A_64_P029486	Tnnc2	1.940	0.045	Down	Troponin C type 2
A_44_P325599	Vwa7	1.132	0.003	Down	von Willebrand factor A domain containing 7
A_64_P096782	Ier5l	1.397	0.030	Down	Immediate early response 5-like
A_44_P231955	Reml	1.271	6.09E–04	Down	RRAD and GEM like GTPase 1
A_43_P11922	Crem	1.746	0.035	Down	cAMP-responsive element modulator
A_44_P401146	Npm2	1.133	0.015	Down	Nucleophosmin/nucleoplasmin 2

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.25 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.26 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.27 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,30 and immune response.31 These findings are consistent with previous studies showing that neuroinflammation32 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.35 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.8 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.