Long Non-Coding RNA and mRNA Profiles in the Spinal Cord of Rats with Resiniferatoxin-Induced Neuropathic Pain

Introduction

Neuropathic pain caused by nerve injury or disorder negatively impacts life quality, and the treatment is still challenging.^1,2 Investigating the potential mechanisms of this disorder may contribute to the development of novel therapeutics for neuropathic pain. In our previous research, we depleted the capsaicin-sensitive afferents in adult rats via systemic approach by using RTX, an ultrapotent transient receptor potential vanilloid 1 (TRPV1) agonist, to generate long-lasting tactile allodynia and impaired thermal sensitivity,^3,4 which is a typical feature of postherpetic neuralgia (PHN). Moreover, tactile allodynia induced by resiniferatoxin (RTX) is probably resulted from aberrant formation of myelinated afferent nerve fibers (MANF) in lamina II of spinal cord.^5–7 Therefore, the spinal dorsal horn (SDH) is a critical site for initiating allodynia in PHN. Although the RTX model has been widely used to explore the mechanisms of neuropathic pain and develop novel analgesic agents,^8,9 little is known about the differential gene expression profile in the SDH induced by RTX treatment.

RNA isoforms and long noncoding RNA (lncRNAs) play an important role in regulating cellular activities and lead to pathological changes in many diseases.^10–12 LncRNAs, containing 200 nucleotides, exist in the cytoplasm and nucleus.^11,13 LncRNAs, as decoys, scaffolds, and guides, may serve as a regulator in the expression and localization of genes.¹² Recently, numerous lncRNAs were found, and altering the expression pattern of lncRNAs might be related to neuropathic pain.^14,15 However, there was one report on the regulation of neuropathic pain by lncRNA. The Kcna2 antisense RNA expression level is elevated after peripheral nerve is injured, which may increase the excitability of neurons, and blocking KCNA 2 antisense RNA attenuates neuropathic pain.¹⁵ It is of necessity to develop in-depth new insight into the molecular events in SDH by analyzing the lncRNA and mRNA expression changes and their relationship in neuropathic pain.

Microarray techniques can be used to analyze expression changes of numerous mRNAs and lncRNAs.^16–18 Previous studies have used mRNA expression profiling to identify the gene spectrum in SDH of animals with neuropathic pain.^19,20 However, it is unclear whether differential expressions of mRNAs play a part in lncRNAs changes in the spinal cord during the development of PHN. Therefore, in this study, we used microarray to determine changes in mRNAs and lncRNAs expression in the spinal dorsal horn of rats with RTX-stimulated PHN. Gene ontology (GO) enrichment analysis was used to investigate the exact biological processes associated with the differential gene expression.

Materials and Methods

In-vivo Models

Sprague–Dawley rats (male, 250–280 g, adult rats) were provided by the Experimental Animal Center of Tongji Medical College of Huazhong University of Science and Technology. Our experiments obtained approval from the Animal Care Committee at Huazhong University of Science and Technology, and abided by the ethical rules of the International Association for the Study of Pain.²¹ The animals were fed in the cage, and allowed free access to water and food. The environment followed a 12-hour light/darkness cycle. In the RTX group, the rats were intraperitoneally injected with RTX (250μg/kg, LC Laboratories, Woburn, MA) after anesthetized by halothane (2% in O₂). RTX solution was prepared by dissolving the agent in normal saline containing ethanol (10%) and Tween 80 (10%), with a final concentration of 0.2μg/μL.²² Normal saline containing 10% ethanol and 10% Tween 80 was given to the rats in the vehicle control group. Before treatment, the baseline data of sensitivity to thermal and mechanical induction were recorded.

Nociceptive Behavioral Tests

Before RTX treatment, the behavioral tests were conducted 3 times. 1) Thermal sensitivity evaluation. The animals stayed in the test environment for 30 min for habituation. As described in the previous research (UgoBasile, Italy), hind paws of the rats were exposed to radiant heat beam via transparent glass to evaluate thermal sensitivity by using a plantar analgesia meter.²³ The average value of withdrawal latency from 2 to 3 trials was calculated. Possible tissue injury was avoided by using a cutoff value of 30s.²² 2) Mechanical allodynia evaluation. The rats were placed on the elevated mesh floor to evaluate mechanical allodynia. The “up-down” method was used to determine the tactile threshold (Chaplan et al, 1994; Chen and Pan, 2002). After 30 min, calibrated von Frey filaments (provided by Stoelting, Wood Dale, IL) were applied vertically to the plantar plane of the two hind paws using appropriate force to make the filament bend for 6s. The withdrawal of the brisk or flinching of the paw was regarded as a positive response. The average value was calculated after repeating the test for 2 to 3 times.

Tissue Preparation and RNA Extraction

Six weeks after injection, the dorsal spinal cord tissues of the rats in RTX or vehicle groups were excised immediately after the animals were sacrificed. Trizol reagent (provided by Invitrogen, Carlsbad, CA, USA) was utilized for extraction of total RNA as per the manual. Six RNA samples for microarrays were obtained (n=3 in vehicle group, and n=3 in RTX group). The purity and quantity of RNA were evaluated by employing a K5500 micro-spectrophotometer (BoKai company, Beijing). RNA purity was regarded as acceptable when A260/A230≥0.5 and A260/A280≥1.5. RNA integrity was regarded as acceptable if RIN value≥7. Gel electrophoresis was used to assess genomic DNA contamination.

Expression Microarray

RiboArray microArray (Ribobio Co., Guangzhou, China) was applied in this study. The RiboArray microArray covering 29,659 unique rat mRNAs and 6713 unique rat lncRNAs (refseq60) was provided by Ribobio (Ribobio Co. Guangzhou, China) via the Combimatrix platform (Washington-Seattle, USA). This experiment was performed as per the protocol. In brief, antisense RNA (aRNA) was acquired via reverse transcription of 1μg of the extracted RNA. Then, Amino Allyl MessageAmpTM IIaRNA Kit (Life Technologies, USA) was used to label 4 μg of aRNA. The microarrays were hybridized with the labeled samples in a 40°C environment overnight and then washed by buffer as per the manual. The Genepix 4000B laser scanner (Molecular Device, USA) was used to scan the slides and Genepix Pro 7.0 software (Molecular Device, USA) was applied to analyze the images. The intensity of the signal weaker than 1.5 times of the background was removed. LncRNA and mRNA microarray, including labeling, hybridization, scanning, normalization, and data analysis, was performed by Ribobio. We selected the lncRNAs and mRNAs with expression levels of at least 1.5-fold difference between RTX-treated samples and control samples.

Microarray Analysis

All files related to genes were input into RiboArray Expression Console (Ribobio Co, Guangzhou, China; version 1.2.1). The quantile method was used for normalization.^24,25 The normalized intensity was adjusted by using Combat Software to eliminate batch effects. The differentially expressed genes (DEGs) between the two groups of rats were discriminated by using the random variance model (RVM) ANOVA which could effectively increase the degree of freedom when the sample size was small. DEGs were selected if the p-value was smaller than 0.05, and the altered expression was more than 1.5-fold. Unsupervised hierarchical clustering analysis was conducted, and RiboBio Co. Ltd performed microarray analysis.

Functional Group Analysis

The functions of the identified DEGs were investigated by using GO analysis. According to molecular functions and biological processes, the genes were grouped into hierarchical categories and their regulatory networks were revealed in GO analysis (http://www.geneontology.org/).^26,27 Fisher’s exact test (two-sided) was applied to the classification of GO category. P-value was adjusted by FDR.²⁸ P-values were computed for GOs enriched in DEGs (the recommended cutoff p-value was 0.05). The pathway analysis of DEGs was conducted by using Reactome (http://www.genome. jp/kegg/), Biocarta, and Kyoto Encyclopedia of Genes Genomes (KEGG). The pathways were recognized by using Fisher’s exact tests, and p-value was used to define the significance threshold. The enrichment was calculated similarly as the GO analysis.^29,30

qPCR Analysis

We extracted the total RNA from SDH specimens by employing Trizol reagent (provided by Invitrogen, Carlsbad, USA). cDNA was harvested after reverse transcription by ReverTra Ace-a-TM (Toyobo, Osaka, Japan). The total volume of PCR reaction system was 10 µL, containing SYBR qPCR mix (5 µL) (Toyobo, Osaka, Japan), diluted cDNA (1 µL), corresponding primer (1 µL), and free water (3 µL). CFX96 system (Bio-Rad) was utilized to perform RT-qPCR. Bio-Rad CFX Manager software was applied to quantify the relative expression. The internal control was β-actin. The primers used in this experiment are shown in Table 1. The threshold cycle (CT) method was used to determine the expression level of the genes. CT value was obtained for ENSRNOG00000022535, ENSRNOG00000042027, NR_027478, NR_030675, Apobec3b, and β-actin. After subtraction of CT value of β-actin by that of ENSRNOG00000022535, ENSRNOG00000042027, NR_027478, NR_030675, and Apobec3b, the delta CT (ΔCt) value was calculated. The 2^−ΔΔCt method was used in the calculation. The relative mRNA levels of 8 independent experiments relative to the vehicle control group were expressed as percentage.

Table 1 Primers Used for Real-Time PCR Analysis of lncRNA and mRNA Levels

Statistical Analysis

Data were expressed as mean ± SD. Unpaired two-tailed Student’s t-test was used to analyze the differences between 2 groups, and P < 0.05 was regarded statistically significant. Differentially expressed lncRNA and mRNA were determined if the fold change ≥2.0 and P value < 0.05. GraphPad Prism software was used for statistical analysis.

Results

Nociceptive Behavioral Data

The nociceptive behavioral tests were conducted prior to RTX injection.⁴ The test results showed the withdrawal thresholds were similar in the rat models at baseline before RTX exposure. 4 days after RTX exposure, the thermal sensitivity was significantly decreased, while the tactile sensitivity was significantly increased.

Microarray Hybridization Data

The baseline withdrawal thresholds of rats treated with RTX and vehicle were quantified by von Frey filaments, which is similar to what we reported previously.⁴ The RNA quantity and purity were validated (Figure 1A and B). The microarray data were stored in the Gene Expression Omnibus (GEO) database (GEO accession GSE46756). Finally, 6713 noncoding RNAs and 29,659 coding genes were screened from RTX and vehicle control groups and compared. The heat map was plotted using the hierarchical clustering results, and distinguishable expression profiling of lncRNAs and mRNAs between the two groups was demonstrated (Figure 1C and D).

Figure 1 The expression profile analysis of the lncRNAs/mRNAs and overview of microarray analysis. The purity and quantity of RNA were evaluated by a K5500 micro-spectrophotometer. RNA purity was acceptable if A260/A230 ≥ 0.5 and A260/A280 ≥ 1.5. RNA integrity was acceptable when RIN value ≥ 7 detected by Agilent 2200 RNA assay. Gel electrophoresis (A) was conducted to assess the genomic DNA contamination. Expression signal distribution after every mRNA array dataset was normalized (B). The microarray heat map of spinal dorsal horn tissues in the PHN model and analysis of differentially expressed lncRNAs (Vehicle group: A1, A2 and A3; PHN model group: B1, B2 and B3) (C and D). Red color represented that the relative expression was high and green color represented that the relative expression was low. For advanced data analysis, all biological replicates were pooled and calculated to identify differentially expressed mRNAs based on the threshold of |Fold change| ≥ 1.5 and P ≤ 0.05. Unsupervised hierarchical clustering analysis was conducted to demonstrate the correlation of expression profiles between treatment and biological replicates.

Analysis of Differentially Expressed lncRNAs

Microarray technology was used to determine lncRNA expressions in spinal dorsal horn samples of adult rats 6 weeks after treatment with RTX or vehicle. RVM ANOVA revealed that differentially expressed lncRNAs included 139 upregulated and 140 down-regulated lncRNAs. The top 25 upregulated and downregulated lncRNAs are listed in the Supplementary Materials. The fold change of 404 lncRNAs was greater compared to vehicle group (p<0.05, fold change ≥1.5), of which the expression levels were increased in the RTX group. As shown in Figures 1C and S1, there were significant differences in these lncRNAs between the two groups. According to our microarray datasets, the expression levels of these lncRNAs were increased in the RTX group. Particularly, ENSRNOG00000032914 (upregulated) and ENSRNOG00000011429 (downregulated) were the most significantly altered lncRNAs. The maximal and minimal fold change was 57.03159 and 1.506366, respectively.

Analysis of Differentially Expressed mRNAs

According to the mRNA expression profiles, we discriminated the differentially expressed mRNAs between RTX and vehicle control groups. Analysis of the microarray data showed that there were 29,659 differentially expressed mRNAs (Figure S2), in which levels of 745 mRNAs were increased and levels of 590 mRNAs were decreased. The top 25 upregulated and downregulated mRNAs were listed in the Supplementary Materials. The maximal and minimal fold change was 57.03159 and 1.506366, respectively. Among these mRNAs, 46 showed greater fold change (p<0.05, ratio >5). The most significantly changed mRNAs were zinc finger protein 90 (zfp90) (upregulated) and LOC100912010 (downregulated).

Figure 2 LncRNAs–coexpressed mRNAs were subjected to gene ontology enrichment and pathway analysis. Three sections were covered in the Gene ontology: cellular component, molecular function, and biological process.

GO Analysis and KEGG Pathway Analysis of Differentially Expressed mRNAs

Figure 2 presented the GO ontology analysis of the differentially expressed mRNAs classified according to biological processes, cellular components, and molecular function. The KEGG pathway analysis indicated that the differentially expressed mRNAs participated in processing of the environment information, diseases in human, metabolism activities, pathways in organismal systems (Figure 3), mRNA KEGG up-thumbnail (Figure 4), and mRNA KEGG down-thumbnail (Figure 5). The network heat map for functional enrichment analysis is detailed in the Supplementary Materials (Figures S1–S4). Among these pathways, 52 genes were related to neuropathic pain, including VEGF, calpain, netrin-1, DCC, TRPV1, NGF, sema3a, slit3, zinc finger protein, caspase-1, Nlrp3, sema3b, BDNF, etc. And there were more than 2 differentially expressed mRNAs acting as regulators (Figures S3 and S4).

Figure 3 KEGG classification map of differential genes: Environmental Information Processing; Human Diseases Metabolism; Organismal Systems.

Figure 4 mRNA KEGG up-thumbnail. Upregulated mRNAs are listed.

Figure 5 mRNA KEGG down-thumbnail. Downregulated mRNAs are listed.

LncRNA-mRNA Interaction Network

All lncRNA target genes were obtained using lncRNA target gene prediction. The interaction between lncRNA and target genes was obtained. LncRNA-mRNA interaction network diagram is shown in Figure 6A. Due to a large amount of data in this project, lncRNAs with more than 20 target genes were screened. The intersection of mRNAs predicted by lncRNAs and differential mRNA results was revealed. The interaction network diagram between the intersection genes and the corresponding lncRNAs is shown in Figure 6B. An overview of the abnormal lncRNAs was obtained after analyzing microarray data (Figure 6C). To validate microarray data, RT-qPCR analysis was conducted by using four differentially expressed lncRNAs (ENSRNOG00000022535, NR_027478, ENSRNOG00000042027, and NR_030675) and the Apobec3b mRNA. ENSRNOG00000022535 was downregulated, and NR_027478, ENSRNOG00000042027, NR_030675, and Apobec3b mRNA levels were upregulated (Figure 7A–E). The expression of these four lncRNAs and one mRNA showed the change was similar to the microarray data. According to the correlation analysis, RT-qPCR results were closely correlated with the microarray data (Figure 7).

Figure 6 LncRNA-mRNA interaction network diagram. The red box represents lncRNAs, the blue circle represents predicted mRNAs, and the line represents the relationship between lncRNAs and mRNAs. All LncRNA target genes were obtained by LncRNA target gene prediction. The interaction between lncRNAs and target genes was obtained. LncRNA-mRNA interaction network diagram was drawn. Due to the large number of results of this project, lncRNAs with more than 20 target genes were screened (A). A network diagram of the relationships between overlapping genes. The intersection of mRNAs predicted by lncRNAs and differential mRNA results was revealed, and the interaction network diagram between the intersection genes and the corresponding lncRNA was drawn (B). An overview of the aberrant lncRNAs by analyzing the microarray data (C).

Figure 7 Real-time PCR were conducted to verify the differentially expressed lncRNAs and mRNAs. One downregulated and four upregulated lncRNAs and one upregulated mRNA were randomly selected and studied. Summary data showed the relative lncRNA level of ENSRNOG00000022535 (A) and NR_027478 (B) and ENSRNOG00000042027 (C) and NR_030675 (D). Data showed the relative of mRNA levels of Apobec3b (E) in the dorsal horn of the spinal cord. Data were expressed as mean ± SEM (n =6 per strain). *p < 0.05, compared to the vehicle group.

Discussion

Neuropathic pain is a prevalent clinical problem with limited treatment options. In the SDH, systemic RTX exposure can induce severe tactile allodynia and sprouting of myelinated afferent fibers.^3,4 It was revealed that RTX-induced neuropathy was associated with α2δ-1 up-regulation in the dorsal root ganglion (DRG) and increased physical interaction between α2δ-1 and GluN1 in the spinal cord synaptosomes.³¹ Besides, α2δ-1–bound NMDA receptors at presynaptic terminals of sprouting myelinated afferent nerves contribute to RTX-induced potentiation of nociceptive input to the spinal cord and tactile allodynia.³² Microarray technology represents a potentially powerful method to assess the changed expression levels of genes in the nervous system.³³ The genes enriched in the spinal dorsal horn could be grouped into various functional categories such as neuropeptides, ion channels, receptors, calcium/calmodulin-binding proteins, transcription factors, and synaptic proteins.³⁴ In the spinal dorsal horn, changes of gene expression are involved in the development of neuropathic pain,³⁵ however, changes in gene expression relevant to neuropathic pain remain uncertain. The transcription factors zfp and signaling molecules Galpha (olf) are enriched in the dorsal horn of adult rats.³⁴ Our study demonstrated that zfp and olf were profoundly upregulated in the dorsal horn of RTX-treated rats. It was found that Zinc finger E box binding protein-1 (ZEB1) and myeloid zinc finger protein 1 (MZF1) were greatly increased in rats with chronic constriction injury (CCI), and knockdown of ZEB1 and MZF1 could reduce neuropathic pain development.^36,37 MZF1 is the upstream regulator of KCNA2-AS, and microinjection of KCNA2-AS into DRG reduced total Kv current and increased the excitability of DRG neurons, producing neuropathic pain symptoms.¹⁵ As a direct target of miR-128-3p, ZEB1 can be significantly repressed by LV-miR-128-3p. Moreover, miR-128-3p rescued the effects of ZEB1 on neuropathic pain progression via inhibiting neuroinflammation, and it was implied that miR-128-3p can alleviate the progression of neuropathic pain via modulating ZEB1.³⁸ In previous research, the in-vivo models of neuropathic pain were constructed and the mRNA levels in the spinal dorsal horn were determined by using microarray technology or investigating certain transcripts.^34,39,40 LncRNAs were found to regulate pain and excitability of sensory neuronal.^15,31 For instance, MEG3, a novel lncRNA, plays an essential role in the development of neuropathic pain.⁴¹ To explore molecular mechanisms of PHN, we investigated expression profiling of mRNAs and lncRNAs in dorsal horn tissues extracted from rats after RTX treatment to explore new RNAs (eg, mRNA isoforms and lncRNAs) involved. We found differential expression patterns of lncRNAs and mRNAs in RTX-treated rats, in which 745 mRNAs and 139 lncRNAs were upregulated, whereas 590 mRNAs and 140 lncRNAs were downregulated in spinal dorsal horn tissues.

Protein network analysis indicated that the genes with elevated expression levels in SDH were related to injury-induced responses. We identified lncRNAs targeted genes from LNCipedia.org, an online comprehensive lncRNA database. RNAs, such as Apobec3b, Kcnd3, Netrin-1, were new targets in neuropathic pain development.^42–44 The differential expression levels of lncRNA (ENSRNOG00000022535, NR_027478, ENSRNOG00000042027 and NR_030675) and mRNA (Apobec3b) in the dorsal horn after RTX treatment were validated by RT-qPCR. These data were consistent with the microarray data. Additionally, Apobec3b was upregulated by lncRNA (NR_030675), which was obtained by lncRNA target gene prediction in the spinal cord dorsal horn. The possible association between high expression levels of NR_030675 and Apobec3b is notable The transcriptional change in Apobec3b is strongly related to the development of neuropathic pain.³⁸ Thus, upregulation of Apobec3b in the spinal cord may lead to tactile allodynia induced by RTX treatment. These lncRNAs and mRNAs may be novel targets in the treatment of neuropathic pain.

Emerging studies have shown that lncRNAs play crucial roles in the occurrence and development of neuropathic pain by regulating ion channels and neuroinflammation, and the two key features drive the pathogenesis of neuropathic pain. Notably, circulating lncRNAs with rich differences could also be found in the plasma of patients and were positively correlated with pain scores, indicating that lncRNAs could be a new source of biomarkers for identifying and monitoring patients with neuropathic pain.⁴⁵ However, the association between patients’ pain scores and lncRNA or mRNA levels in the plasma remains to be confirmed. Importantly, more translational studies are needed to maximize their use as biomarkers and therapeutic targets in neuropathic pain.

Conclusion

Ours study found that 745 mRNA and 139 lncRNAs were upregulated and 590 mRNA and 140 lncRNAs were downregulated in spinal dorsal horn tissues after RTX exposure. lncRNA (ENSRNOG00000022535, ENSRNOG00000042027, NR_027478, NR_030675) and Apobec3b mRNA were expressed differently in spinal cord tissues, indicating the association between NR_030675 and Apobec3b levels. These specifically and differentially expressed molecule targets are enriched in the spinal cord dorsal horn, which provide new information for further investigation on the mechanisms and therapeutics of neuropathic pain.