Identification of lncRNA expression profiles and ceRNA analysis in the spinal cord of morphine-tolerant rats

Morphine tolerance is a challenging clinical problem that limits the use of morphine in pain treatment, but the mechanisms of morphine tolerance remain unclear. Recent research indicates that long noncoding RNAs (lncRNAs) might be a novel and promising target in the pathogeneses of diseases. Therefore, we hypothesized that lncRNAs might play a role in the development of morphine tolerance. Male Sprague-Dawley rats were intrathecally injected with 10 μg morphine twice daily for 7 consecutive days. The animals were then sacrificed for lncRNA microarray tests, and the results were validated by RT-qPCR. Next, functional predictions for the differentially expressed mRNAs (DEmRNAs) were made with the Gene Ontology/Kyoto Encyclopedia of Genes and Genomes (GO/KEGG), and predictions for the differentially expressed lncRNAs (DElncRNAs) were made based on competitive endogenous RNA (ceRNA) analyses. The rats successfully developed morphine tolerance. LncRNA microarray analysis revealed that, according to the criteria of a log2 (fold change) > 1.5 and a P-value < 0.05, 136 lncRNAs and 278 mRNAs were differentially expressed in the morphine tolerance group (MT) compared with the normal saline group (NS). The functions of the DEmRNAs likely involve in the processes of the ion channel transport, pain transmission and immune response. The ceRNA analysis indicated that several possible interacting networks existed, including (MRAK150340, MRAK161211)/miR-219b/Tollip.Further annotations of the potential target mRNAs of the miRNAs according to the gene database suggested that the possible functions of these mRNAs primarily involved the regulation of ubiquitylation, G protein-linked receptors, and Toll-like receptors, which play roles in the development of morphine tolerance. Our findings revealed the profiles of differentially expressed lncRNAs in morphine tolerance conditions, and among these lncRNAs, some DElncRNAs might be new therapeutic targets for morphine tolerance.


Introduction
Morphine tolerance is defined as the diminished analgesic effect and the need for a higher dose to achieve the desired analgesic effect after chronic exposure to morphine [1][2][3][4]. Over the past decades, people have attempted to elaborate the mechanisms of morphine tolerance with minimal success. The recent focus of this area is the role of non-coding RNAs (ncRNAs) [5,6]. Among ncRNAs, some studies have reported that miR-NAs are involved in the development of morphine tolerance [7,8], including the let-7 family, miR-23b [9], miR-133b, miR-339 [10], miR-365 [11] and miR-219-5p [12]. However, the research regarding long non-coding RNAs (lncRNAs) is still in its infancy.
The sequences of lncRNAs range in size from approximately 200 nt to over 100 kb, and lncRNAs can function in novel mechanisms of the modulation of the expression of genes [13,14]. Over the past decades, lncRNAs have been the highlight of some studies of cancer, osteoarthritis, nervous system function and development [15][16][17]. In pain research, recent studies have claimed that Knav2 AS [18], uc.48+ [19] and some other lncRNAs [20] might play roles in the process of the development of neuropathic pain [21,22]. The function of lncRNAs has been illustrated as follows: lncRNAs can interact with mRNAs, bind to transcription factors, modulate chromatin remodeling, and even directly regulate the functions of target proteins. Among lncRNAs, some may act as competitive endogenous RNAs (ceRNAs) that target miRNAs [14].
The ceRNA hypothesis was proposed in 2011 with the aim of further elaborating the relationships among RNAs. Salmena et al. [23] highlighted that some RNAs act as ceRNAs that participate in mutual competition for common binding sites of target miRNAs and thus modify the functions of the target miRNAs. Recent research has indicated that ceRNA analysis might shed a light on functional predictions of the effects of lncRNAs [24].
Therefore, in the present study, we hypothesized that lncRNAs might be differentially expressed and act in direct or indirect manners during the development of morphine tolerance. To address this hypothesis, we attempted to identify the lncRNA expression profiles in the spinal cords of rats under normal and morphine tolerance conditions and to predict the possible functions of differentially expressed lncRNAs (DElncRNAs).

Intrathecal injection of morphine induces a chronic morphine tolerance model in rats
Adult male Sprague-Dawley rats (weight 240-260 g) were obtained from the Hunan SJA Laboratory Animal Company (Hunan, China). The rats were housed in groups and maintained on a 12/12 light-dark cycle at a room temperature of 22 ± 1°C with food and water freely available. The experimental procedures were approved by the Animal Care and Use Committee of Central South University and conducted in strict accordance with the guidelines of the International Association for the Study of Pain [25]. To establish the rat model of morphine tolerance, rats in the morphine tolerance group (MT, n = 8) were intrathecally injected (i.t.) with 10 μg (1 μg/1 μl) of morphine twice daily at 08: 00-09:00 am and 16:00-17:00 pm for 7 consecutive days [3,26]. The normal saline group (NS, n = 8) was injected with equal volumes of normal saline at the same time points.

Tail flick test
The tail flick test was used to measure thermal sensitivity. Before conducting this test, the rats were placed on the plantar surface for 15 min to adapt to the testing environment. Then, we tested one fixed point of the tail 2-3 cm from the tip using the Hargreaves apparatus (Italy, UGO Basile) [26,27]. The results are expressed as the tail withdrawal latency (TWL), which was ultimately converted to the percent of the maximum possible effect (%MPE). The radiant index was set at 90, and the cut-off was 20 s to avoid tissue damage. The tests were conducted 30 min before and after the morning injection of morphine on days 1, 3, 5, and 7.

Tissue collection and RNA isolation
On day 8, after the morning injection with morphine or saline, the rats were deeply anesthetized with pentobarbital sodium (1%) one hour later. Then, we decapitated the rats, collected the lumbar enlargements and placed the collected tissues into liquid nitrogen as quickly as possible for preservation. Next, we extracted the total RNA from the spinal tissues and tested the RNA quantity and quality with a NanoDrop ND-1000 Spectrophotometer (Thermo,USA) and tested the RNA integrity with agarose gel electrophoresis (2%). Subsequently, we stored the remnant RNA at − 80°C for later use. The RNA isolation was performed by Kangcheng Bio-tech (Shanghai, China).

Microarray assay
We employed a rat lncRNA microarray 4 × 44 k,V1.0 (Arraystar) containing approximately 9000 lncRNAs to screen the differentially expressed lncRNAs and mRNAs. The total RNAs of the MT and NS groups (n = 5) were hybridized with the gene chips. The RNA samples were transcribed into fluorescent cRNAs along the entire lengths of the transcripts without a 3′ bias utilizing random primers. The labeled cRNAs were hybridized to the rat lncRNA microarray. Next, the arrays were scanned with an Agilent DNA Microarray Scanner (part number G2505C). The array images were analyzed with Agilent Feature Extraction software (version 11.0.1.1). Quantile normalization and subsequent data processing were performed using the GeneSpring GX v12.1 software package (Agilent Technologies,USA). The microarray hybridization was performed by Kangcheng Bio-tech (Shanghai, China).

Bioinformatics analysis
We used the criteria of a log2 (fold-change) > 1.5 and a P-value < 0.05 to screen for the deregulated RNAs. Hierarchical clustering was performed with Cluster 3.0, and the heat maps were generated in Java Treeview. The DEmRNAs were analyzed according to the pathway annotations of the Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) functional enrichment using CapitalBion. The -log 10 (P-values) of the GO and pathway results are displayed in the histogram. The DElncRNAs were analyzed by ceRNA analyses, which were conducted with Arraystar's homemade miRNA target prediction software, which is based on TargetScan and miRanda [28]. A lncRNA/ miRNA/mRNA interaction network was generated to visualize the interactions using Cytoscape. The NCBI Database was used to annotate the functions of the potential target genes.

Real-time quantitative polymerase chain reaction (RT-qPCR)
The microarray results were confirmed by RT-qPCR. The total RNAs of the MT and NS groups (n = 5 in each group) were reverse transcribed using random hexamer primers (Arraystar Flash RNA Labeling Kit, Arraystar) according to the manufacturer's description. The expression levels of 18 lncRNAs (MRAK080737, MRAK159688, MRAK046606, DQ266361, XR_005988, uc.48+, uc.310-, XR_009527, XR_008662, S66184, MRAK161211, MRAK 150340, XR_006440, AF196267, MRAK077287, MRAK01 4088, MRAK141001 and MRAK038897), as well as 12 DEmRNAs (S100a8, Batf, Ccl7, RT1-Bb, RatNP-3b, Grm2 and Mmp9, Prrx2, Asb2, Fam111a, Kcnv2 and Tmem119) were examined. GAPDH was used as the house-keeping gene. The sequences of all primers are presented in Table 1. We conducted the RT-qPCR tests with a 10 μL reaction system in a ViiA 7 Real-time qPCR System according to the manufacturer's protocol. Melting-curve analysis was performed to monitor the specificity of the production. All experiments were replicated three times. The gene expression levels in the MT and NS groups were analyzed with the 2 −ΔΔCT method.

Statistical analysis
All data were presented as mean ± s.e.m. The statistical significance of differences between groups was analyzed with two-way repeated-measures of ANOVA followed by Bonferroni test or with Student's t-test. P values less than 0.05 were considered statistically significant.

Construction of the rat morphine tolerance model
The tail flick test data revealed that there were no significant changes in the thermal sensitivities of NS group (n = 8 in each group). However, comparisons between the two groups revealed that, on day 1 post-injection, the %MPE of the MT group (n = 8 in each group) was significantly higher than that of the NS group (P<0.05). On day 3 post-injection, the % MPE of the MT group began to decline but remained higher than that of the NS group (P<0.05). On day 5 post-injection, the %MPEs did not significantly differ between the two groups (P>0. 05), and this state continued to day 7 (Fig. 1a), which suggested that a stable morphine tolerance model had been established.

Overview of the lncRNA and mRNA expression profiles
First, we created an overview of the lncRNA and mRNA expression profiles using a using scatter plot, which revealed that large numbers of lncRNAs and mRNAs were differentially expressed between the MT and NS groups (n = 5; Fig. 1d). Next, hierarchical cluster analyses of all of the lncRNAs and mRNAs was applied and revealed that the 5 NS and 5 MT samples clustered independently, and the results also indicated high degrees of consistency in both the NS and MT groups (Fig. 1b, c). All of the microarray results have been uploaded to the GEO database (GSE110115).

Differentially expressed lncRNAs and mRNAs in morphine tolerance
According to the criteria of a log2 (fold change) > 1.5 and a P-value< 0.05, the microarray data identified 136 lncRNAs, including 84 up-regulated and 52 downregulated lncRNAs, which were significantly altered in the MT group compared with the NS group. The lncRNAs that exhibited the greatest up-regulations were XR_005988, DQ266361, and MRAK046606 with XR_ 005988 exhibiting the largest up-regulation [log2 (fold change) =12.4243]. The lncRNAs that exhibited the greatest down-regulations were AF196267, XR_009493, and MRAK150340 with AF196267 exhibiting the largest down-regulation [log2 (fold change) =2.2025]. Detailed information, including the top 40 up-regulated and top 40 down-regulated lncRNAs, is provided in Table 2.
Regarding the DEmRNAs, there were 278 genes (176 up-regulated and 102 down-regulated) whose changes met the criteria. These DEmRNAs contained many genes that are known to be involved in pain processing, including Ccl7, Batf, S100a8, Kcnv2, Rgs1, Prrx2, Mmp9, etc. Detailed information about the top 30 up-regulated and top 30 down-regulated mRNAs is listed in Table 3.

Validation of the lncRNA and mRNA expressions
To validate the reliability of the microarray results, 18 DElncRNAs and 12 DEmRNAs were selected and validated by RT-qPCR. The data shown that 11 of 18 selected DElncRNAs (XR_005988, DQ266361, MRAK159688, XR_ 008662, XR_009527, S66184, MRAK150340, MRAK16 1211, MRAK038897, AF196267 and MRAK141001) and 7 of 12 selected DEmRNAs (Batf, Ccl7, RatNP-3b, Mmp9, Kcnv2, Tmem119 and Asb2) exhibited the same trends in altered expressions and the same significant differences in the microarray and RT-qPCR analyses. On the other hand, the other altered mRNAs and lncRNAs exhibited the trends in changes, but the differences between the two groups were not significant, possibly due to the small sample size (Fig. 2a, b).

Class distributions of the DElncRNAs
The examined lncRNAs were categorized into five groups according to their associations with coding genes: intergenic lncRNAs, antisense overlap lncRNAs, sense overlap lncRNAs, bidirectional lncRNAs, and other. One of the mechanisms by which lncRNAs act is through the interplay with adjacent coding genes [20,29]; therefore, it was important to classify the locations of the lncRNAs.
Our data revealed that, among these DElncRNAs, intergenic lncRNAs and sense overlap lncRNAs accounted for the majority, and only five lncRNAs belonged to the bidirectional category. The concrete data are presented in Fig. 2c.

Functional predictions for the DEmRNAs in the morphinetolerant rats
To explore the molecular functions of the DEmRNAs in morphine tolerance conditions, we further performed GO and pathway analyses genes that differentially regulated in the MT and NS groups. The pathway analyses indicated that the most significantly enriched pathways of the up-regulated genes included the TNF metabolic pathway and phagocytic processes, and the downregulated genes were involved in synaptic vesicle activity, the arachidonic acid metabolic pathway, etc. (Fig. 3a, b). The GO results revealed that the most significantly enriched molecular functions of the up-regulated genes in the MT group were peptidase activity, G-proteincoupled receptor activity, and biological processes concentrated on the cytokine response, defense and the immune response. The cell components primarily belonged to extracellular domains and intercellular domains. The most significantly enriched molecular functions of the down-regulated genes in the MT group were voltagegated channel activity, ion transmembrane activity, and biological processes concentrated on potassium ion transport. The cell components were associated with the sarcolemma, cytoplasmic membrane, and ion channel complexes ( Fig. 3c-h).

Functional predictions of the DElncRNAs lncRNA/miRNA/ mRNA interactions
We performed coding-noncoding gene co-expression (CNC) analysis, but we found no DEmRNA-associated DElncRNAs, which meant that we were unable to make forecasts about the functions of the DElncRNAs according to the related DEmRNAs. However, ceRNA analysis allowed us to predict the possible functions of the  DElncRNAs regardless of their adjacent coding genes. According to the ceRNA analysis, we obtained an overview of the potential lncRNA/miRNA/mRNA interactions (Fig. 4a), and we then further identified several promising networks of lncRNA/miRNA/mRNA interactions (Fig. 4b), which included (MRAK161211, MRAK150340)/miR-219b/Tollip and XR_006440/(miR-365, let7)/(Usp31, Usp42, Clcn4-2) networks, and we used these networks to create functional annotations of the predicted target mRNAs by searching the Gene database. The results indicated that the predicted target mRNAs mainly functioned in the process of ubiquitinylation, the GRCP and TLR signaling pathways, and the modulations of transcription, translation and posttranslational modification, and these functions might constitute the foundation of morphine tolerance.

Discussion
In the present study, we detected 136 DElncRNAs and 278 DEmRNAs overall and we found that compared with normal rats, DElncRNAs and DEmRNAs were present in the spinal cords of morphine-tolerant rats; GO and KEGG pathway analyses revealed that the potential functions of the DEmRNAs may be concentrated on the processes that are thought be involved in the formation of morphine tolerance; and the ceRNA analysis identified several potential lncRNA/miRNA/mRNA interaction networks that might modulate the development of morphine tolerance. To identify as many DElncRNAs and DEmRNAs candidates as possible, we set the criteria at a log2 (fold change) > 1.5 and a P-value < 0.05 [30]. Next, we considered several criteria to select lncRNAs and mRNAs from our microarray data for validation. Firstly, we attempted to select the relevant candidates that might be related to our previous studies of morphine tolerance and miRNAs [11,12]. Secondly, for validation, we chose the lncRNA candidates that exhibited higher fold changes and greater raw expression intensities and had adjacent mRNA that was related to morphine tolerance.
Among the DElncRNAs that we detected, XR_005988 exhibited the most significant up-regulation and was Fig. 1 (a) Continuous injection of morphine for 7 days induced the formation of morphine analgesic tolerance. The data are expressed as the mean ± s.e.m. (n = 8 in each group). **P < 0.01, ***P < 0.001, Two-way repeated-measures of ANOVA followed by Bonferroni test. Ten samples (n = 5 in each group) were subjected to microarray analysis. b-c The entire and partial hierarchical clusterings of the lncRNAs and mRNAs, respectively; the up-and down-regulated genes are colored in red and green, respectively. d Scatter plot displaying the lncRNAs and mRNAs that exhibited expression differences between the MT and NS groups that exceeded 1.5-fold classified as a long intergenic non-coding RNA (lincRNA). Although we obtained no information about XR_005988 or its associated genes through searches of all types of gene databases, it is well known that lncRNAs account for the main portion of lncRNAs and exhibit the most substantial biological functions [31], which indicates that XR_005988 is still a promising lncRNA molecule for further study. Additionally, XR_    005988 might bind to some miRNAs according to the ceRNA analysis, which indicates its potential action as a ceRNA. Of course, additional in vivo and in vitro functional research is needed. MRAK-159688, which was another up-regulated lncRNA identified in our study, is associated with the Fos gene and was named the Fos downstream transcript (FosDT) in a report from Mehta et al. [32] Mehta et al. reported that, in an ischemia/reperfusion model, FosDT interacts with the chromatin-modifying proteins Sin3a and co-repressor of the transcription factor REST (coR-EST) and subsequently represses REST-downstream genes. Moreover, other researchers have indicated that MOR is one of the downstream targets of REST and is negatively modulated by REST in specific neuronal cells [33,34]. Therefore, the dysregulated MRAK-159688 might be involved in the development of morphine tolerance through interactions with REST.
Since the ceRNA hypothesis was proposed, it has been verified in some tumor diseases. For example, the FER1L4/miRNA106a-5p/PTEN pathway constitutes a novel regulatory pathway that is involved in the occurrence and progression of gastric cancer [35]. Currently, ceRNA analysis is a novel method for predicting the functions of lncRNAs. In our study, we found that several possible interacting pathways among ceRNAs exist, including the following: MRAK161211, MRAK150340/ miR-219b/mRNAs (e.g., Tollip, and Ubqln4); XR_ In our previous study, We have confirmed that miR-219-5p can attenuate morphine tolerance by targeting CaMKIIγ [12]. Furthermore, other researchers have also reported that miR-219 is down-regulated following the continuous application of morphine and can regulate NMDA receptor signaling [36]. Whereas in the present study, we predicted that miR-219b, rather than miR-219-5p, would exert this function. Recent research about miR-219 has mainly concentrated on the functions of miR-219-5p, and the other isoform, i.e., miR-219b, has not been studied in detail. However, on the one hand, we identified miR-219-5p and miR-219b, which have both previously been reported to function in the suppression of the proliferation, migration and invasion of cells [37,38]. On the other hand, when we forecasted the downstream targets of miR-219b, we found that their functions were mainly related to the ubiquitination process, G-protein-coupled receptors, alterations in ion channels, and Toll-like receptor signaling pathway regulation, and these function are all involved in the mechanisms of morphine tolerance. Therefore, given its possible target mRNAs, miR-219b is likely to be a new target related to morphine tolerance.
Toll-interacting protein (Tollip) is a potential downstream target of miR-219b and can modulate the expression of the Toll-like receptor (TLR) via its action as an The red nodes mean down-regulated lncRNAs, the gray nodes mean up-regulated lncRNAs, the blue squares mean down-regulated miRNAs we are interested in, the green squares mean up-regulated miRNA we are interested in, and the pink nodes mean the other miRNAs and mRNAs. b The lncRNA/miRNA/mRNA networks that we are interested in are displayed. The green nodes represent lncRNAs, the yellow pentagons represent miRNAs and the pink nodes represent mRNAs we forecasted endogenous inhibitor and regulate the IL-1β-induced activation of NF-κB; thus, miR-219b plays an inhibitory role in inflammatory signaling [39,40]. Previous studies have demonstrated that the TLR-mediated activation of glial cells and TLR4-mediated NF-κB activation might influence the development of morphine tolerance [41,42]. Therefore, further study is required and we presumed that the MRAK150340, MRAK161211 (i.e., down-regulated lncRNAs)/miRNA-219b/Tollip interaction network potentially functions in morphine tolerance.
On the one hand, previous studies have suggested that let-7 can suppress the expression of MOR [43], and our previous results suggest miR-365 can modulate morphine tolerance by targeting the beta-arrestin 2 protein [11]. On the other hand, Law et al [44] reported that opioid receptor agonists (such as morphine) promote the ubiquitylation of scaffold proteins and thereby change the expression pattern of the receptor signaling pathway. In the present study, the potential downstream genes that we predicted, i.e., Usp31 and Usp42, are both ubiquitylation-related genes. Therefore, the XR_006440/ (let7, miR-365)/(Usp31, Usp42) pathway is likely to be responsible for the formation of morphine tolerance. Moreover, there are other interested candidates other than the validated lncRNAs and mRNAs, and further studies will be needed.
In summary, although our study is preliminary and lacks additional functional experiments, our research has revealed that hundreds of lncRNAs, especially XR_ 005988 and MRAK159688, are differentially expressed in the spinal cords of morphine-tolerant rats. Further ceRNA analysis revealed that several possible lncRNA/ miRNA/mRNA interaction networks exist, and among these networks, the (MRAK150340, MRAK161211)/ miR-219b/Tollip network holds the most potential for further studies. Availability of data and materials Please contact author for data requests.
Authors' contributions WZ designed the experiments. JS, JW, JH, CL, YP and WZ performed experiments and analyzed data; JS, QG and WZ drafted the manuscript and finished the final vision of the manuscript. All authors read and approved the final manuscript.

Ethics approval
All experiments were performed in accordance with the Animal Care and Use Committee of Central South University and the guidelines of the International Association for the Study of Pain.