Transcriptomic analysis of long noncoding RNAs and mRNAs expression profiles in the spinal cord of bone cancer pain rats
Molecular Brain volume 13, Article number: 47 (2020)
Bone cancer pain (BCP) is one of the most common types of chronic cancer pain and its pathogenesis has not been fully understood. Long non-coding RNAs (lncRNAs) are new promising targets in the field of pain research, however, their involvements in BCP have not been reported. In the present study, we established the BCP model by implantation of Walker 256 carcinoma cells into rats’ tibial medullary cavity and performed transcriptome sequencing of the ipsilateral lumbar spinal cord to explore changes in expression profiles of lncRNA and mRNA. We identified 1220 differently expressed mRNAs (DEmRNAs) (1171 up-regulated and 49 down-regulated) and 323 differently expressed lncRNAs (DElncRNAs) (246 up-regulated and 77 down-regulated) in BCP model, among which 10 DEmRNAs (5 up-regulated and 5 down-regulated) and 10 DElncRNAs (5 up-regulated and 5 down-regulated) were validated the expression by RT-qPCR. Then, we performed Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis on the expression of DEmRNAs and DElncRNAs, showing that they were mainly enriched in inflammatory and immunologic processes/pathways. Finally, we constructed a co-expression network and a ceRNA network of DEmRNAs and DElncRNAs to exhibit a potential regulatory mechanism of DElncRNAs, directly regulating protein coding gene expression in cis or in trans and indirectly regulating protein coding gene expression by sponging miRNA. In conclusion, our study provided a landscape of dysregulated lncRNA and mRNA in spinal cord of bone cancer pain and detected novel potential targets for treatment in the future.
Cancer pain, caused by primary cancer itself or metastases , often deteriorates the quality of life of cancer patients, with a prevalence rate of 50.7% in all cancer stages and 66.4% in advanced stage . Bone cancer pain (BCP), in which the majority is caused by metastatic tumors, is the most common type of cancer pain . Clinically, BCP is one of the most intractable pain, with its complex manifestation such as ongoing pain and breakthrough pain . As many as 70% of breast or prostate cancer patients with skeletal metastases would subsequently develop into BCP . Opioid analgesics were considered as the golden standard in the treatment of cancer-related pain, but nonmedical opioid use and potential social disorder in cancer patients complicate the situation , besides the direct adverse reactions of opioids. It is of great clinical demand to clarify the underlying mechanism of BCP, thus determining new potential therapeutic targets.
The development of BCP is accompanied by changes in numerous genes expression in the peripheral and central nervous system, which may account for this dysfunctional nociceptive perception [6,7,8]. Our previous studies have shown that dysregulated HDACs and its potential downstream targets in the lumbar spinal cord contribute to the development of BCP [9, 10]. In fact, endogenous noncoding RNAs are another key regulators to modify the expression of mRNAs or protein targets, thereby contributing to the development of pathologic pain .
Long non-coding RNA (lncRNA) is a class of non-coding RNAs with sequence length greater than 200 nucleotides. LncRNAs are proved to be able to regulate the expression of ion channel related genes , purinergic receptor related genes , etc. in neuropathic pain. Although the pivotal role of lncRNA in neuropathic pain has attracted extensive attention [14, 15], the involvement of lncRNA in BCP pathogenesis has rarely been explored.
Thus, in the present study, we performed transcriptome sequencing of ipsilateral lumbar spinal cord of BCP rats, induced by Walker 256 carcinoma cells inoculation into tibia cavity, to explore changes in expression profiles of lncRNA and mRNA, and functional analysis of differential genes to provide new insights on the mechanism of BCP.
Materials and methods
Female Wistar rats from the experimental animal center of Central South University (Changsha, Hunan, China), were used in all animal experiments. Rats were group housed in clear plastic cages with sawdust bedding with room temperature at 22 ± 2 °C, light/dark cycle of 12/12 h and freely accessible to food and water. All efforts were made to minimize suffering and to reduce the number of animals used to the minimum required for statistical accuracy. All experimental procedures were approved by the Animal Committee of Xiangya Hospital, Central South University and were performed in accordance with the guidelines of the International Association for the Study of Pain (IASP).
Modeling of bone cancer pain
The rat model of BCP was established by implantation of Walker 256 carcinoma cells into the medullary cavity of the tibia as we previously described . Briefly, Walker 256 breast cancer cell suspension (Dingguo Changsheng, China), syngeneic to Wistar rats, was thawed, rinsed, and then centrifuged at 1000 rpm for 5 min, followed by intraperitoneal inoculation in Wistar rats (weighed 80–100 g). 7–10 days later, malignant ascites was extracted, rinsed, centrifuged and diluted into 1 × 107 cells/mL for inoculation.
Adult Wistar rats (weighed 220–250 g) were anesthetized with intraperitoneal injection of 50 mg/kg phenobarbital sodium. After shaved and disinfected, a 1 cm transverse superficial incision was made in the skin overlying ligamentum patellae, then a 22-gauge needle was used to perforate the bone cortex with entry point at intercondylar eminence. After removal of puncture needle, a 50-μL microsyringe was inserted into the needle hole and 10 μL of Walker 256 cell suspension was injected into the medullary cavity of tibia slowly, while sham-operated rats were injected with heat-killed cancer cells instead. The puncture point was sealed with bone wax and the incision was irrigated with normal saline and interruptedly sutured.
Radiological and pathological examination of the tibia
To verify the bone destruction caused by tumor inoculation, X-ray radiography image of the left hind limb was obtained at 14 days after BCP modeling. In addition, the left tibia bone tissue was harvested, preserved in 4% paraformaldehyde and decalcified in 10% EDTA for another 21 days. Finally, the tibia was embedded in paraffin, sectioned into 5-μm-thick slice and stained with hematoxylin and eosin (H&E).
Mechanical hyperalgesia was measured by von Frey filaments as previously described . In brief, each rat was allowed to acclimate to a Plexiglas box (22 × 12 × 12 cm3) with a wire mesh floor for 30 min before test; then, a series of calibrated blunt-pointed von Frey filaments with logarithmically incremental stiffness (ranging from 0.6 g to 26 g) were stabbed perpendicularly against the plantar surface of hind paw with even force. A positive response was defined as rapid paw withdrawal and/or paw flinching accompanied by head turning, biting, or licking upon application of von Frey filament. The incremental von Frey filament would be applied until a positive response aroused. The paw withdrawal mechanical threshold (PWMT) was recorded as the lowest force evoking three positive responses out of five applications. Each rat was tested three times and the data were averaged with the experimenter blind to animal grouping.
Tissue collection and RNA extraction
On the 14th day posterior to BCP modeling, all rats were deeply anesthetized, followed by decapitation and exsanguination. The lumbar enlargement spinal cord was quickly isolated out, removed the meninges, and dissected longitudinally on ice; then the ipsilateral spinal cord was snap-frozen and preserved in liquid nitrogen. Total RNA was isolated using RNeasy mini kit (Qiagen, Germany) and checked for a RIN (RNA Integrity Number) to inspect RNA integrity by an Agilent Bioanalyzer 2100 (Agilent technologies, USA), followed by purified by RNAClean XP Kit (Beckman Coulter, USA) and RNase-Free DNase Set (QIAGEN, Germany) according to manufacturers’ instructions, respectively.
Construction of cDNA libraries and high-throughput sequencing
Strand specific libraries were prepared using the TruSeq® Stranded Total RNA Sample Preparation kit (Illumina, USA). Briefly, ribosomal RNA was removed from total RNA by Ribo-Zero rRNA removal beads. Then fragmented RNA was copied into first strand cDNA using reverse transcriptase, followed by second strand cDNA synthesis using DNA Polymerase I and RNase H. After undergoing an end repair process, the addition of a single ‘A’ base, and then ligation of the adapters, cDNA fragments were purified and enriched with PCR to create the final cDNA library followed by quantified with Qubit® 2.0 Fluorometer (Life Technologies, USA) and validated by Agilent 2100 bioanalyzer (Agilent Technologies, USA) to confirm the insert size and calculate the concentration. Cluster was generated by cBot with the library diluted into 10 pM and sequenced on the Illumina HiSeq X-ten (Illumina, USA). The library construction and sequencing were performed at Shanghai Biotechnology Corporation.
Quantification of gene expression and differential expression analysis
Sequencing raw reads were filtered by Seqtk followed by mapping to the rat Rnor_6.0 reference genome using Hisat2 (version:2.0.4) . Stringtie (version:1.3.0) [17, 18] was used to calculate Fragments Per Kilobase of exon model per Million mapped reads (FPKM) with a reference annotation, then edgeR was used to identify differentially expressed mRNAs (DEmRNAs) and differentially expressed lncRNAs (DElncRNAs) , with filter criteria set as P < 0.05, false discovery rate (FDR) < 0.05  and fold change (FC) > 2. Gffcompare (version: 0.9.8) was used to predict novel lncRNA not included in NONCODE v5  and Ensembl database. The location of total identified mRNAs and lncRNAs were schematized on rat chromosome karyotype, among which DEmRNAs and DEmRNAs were indicated, with the use of Idiographica .
Real-time quantitative PCR (RT-qPCR)
Total RNA was extracted from the spinal cord sample, purified and quality inspected ibidem. Then, cDNA was synthesized with ReverTra Ace qPCR Kit (TOYOBO, FSQ-101, Japan) according to the manufacturer’s instructions. Real-time fluorescent quantitative PCR was performed on QuantStudio 5 Real-Time PCR System (Applied Biosystems, USA) with total volume of 20 μl containing 10 μl 2X SYBR Green PCR buffer (ABI, 4368708, USA), 1 μl PCR forward primer (10 μM), 1 μl PCR reverse primer (10 μM), 2 μl template and 6 μl nuclease-free water in each reaction. Three-step program for PCR reaction by QuantStudio™ Design & Analysis Software was as follows: 1 cycle for initial denaturation at 50 °C for 2 min and 95 °C for 10 min, then 40 cycles for PCR amplification with denaturation at 95 °C for 15 s, annealing and extension at 60 °C for 1 min, followed by a melting curve analysis. All experiments included no-template controls and all samples were analyzed independently in triplicate. Expression data were normalized to the expression of Gapdh with 2-ΔΔCt method. The sequences of the used primers were listed in Table 1.
Target prediction of DElncRNAs
LncRNA could regulate gene expression in cis and in trans manner. The potential target coding genes in trans were predicted by the sequence similarity of lncRNA and mRNA with BLAST and their complementary energy, which is above the threshold, computed by RNAplex ; while the coding genes transcribed within 10 kb upstream or downstream of lncRNA in genomic localization were considered as potential target genes in cis.
Gene ontology (GO) and Kyoto encyclopedia of genes and genomes (KEGG) enrichment analysis
GO analysis for DEmRNAs and targets of DElncRNAs was performed to construct gene annotations from the perspective of biological process, cellular component and molecular function. GO analysis was accomplished by DAVID Bioinformatics Resources 6.8 with P values and FDR used to test the reliability of the analysis  and the top 10 terms (according to P value) in each class were presented in bar graph. KEGG pathway analysis was performed to understand the function and interactions among differentially expressed genes and the top 30 pathways were presented in bubble diagram. The enrichment factor was the value ratio between the sequenced gene and all annotated genes enriched in the pathway.
Cellular expression analysis of DEmRNA and comparative analysis, protein-protein interaction (PPI) network construction and functional enrichment of top 200 DEmRNAs
Based on the expression data reported by Ye Zhang et al. , we constructed expression heatmap of DEmRNA among 7 types of cells (i.e., neuron, microglia, astrocytes, oligodendrocyte precursor cell, newly formed oligodendrocyte, myelinating oligodendrocytes, and endothelial cells). The top 200 DEmRNAs (according to P value) were selected to compare with the results of another two researches focused on gene expression variation in the spinal cord in rats chronic pain model, with one by microarray assay  and the other by RNA-seq assay . The comparative analyses of three datasets were performed by a web-based portal, Metascape . The PPI network between the top 200 DEmRNA-encoded protein was constructed by the online database STRING (v11.0, http://string-db.org/) and visualized by Cytoscape 3.7.1  and was further analyzed by Molecular Complex Detection (MCODE)  clustering algorithm to construct subnetwork. The GO biological process and KEGG pathway enrichment of the top 200 DEmRNA were performed and visualized by Metascape .
Co-expression analysis of DElncRNAs/mRNAs
The Pearson correlation coefficient (PCC) of DElncRNAs and DEmRNAs expression value was calculated, and |PCC| > 0.99 and P < 0.001 was set as filter to construct co-expression network of DElncRNAs and DEmRNAs, however, there were too many nodes to exhibit. So, we took the intersection of DEmRNAs and target mRNA of DElncRNAs, and constructed a co-expression network of DElncRNAs and the intersectional DEmRNAs and Cytoscape 3.7.1 was used to exhibit the network .
Competing endogenous RNA (ceRNA) analysis of DElncRNAs and DEmRNAs
The possible target binding of lncRNA/mRNA and miRNA was predicted by miRanda, with a maximum binding free energy less than − 20 . In consideration of the scale of the figure, we set the threshold of |PCC| of relative expression value of DElncRNAs and DEmRNAs more than 0.98 to exhibit the network and Cytoscape 3.7.1 was used to exhibit the network .
The behavioral test and RT-qPCR data were presented as the mean ± standard deviation (SD). For behavioral test data, repeated measures ANOVA was performed to detect overall differences followed by Bonferroni post hoc analysis for multiple comparisons. For RT-qPCR data, student’s t-test was used to determine the differences between groups. P < 0.05 was regarded as statistically significant. GraphPad Prism 8 Software (Graph-Pad, USA) was used for all statistical analyses and statistical graph drawing.
Rats developed mechanical hypersensitivity and osteolytic lesions after tumor cells inoculation
There was no significant difference in the PWMT baseline between sham-operated rats and BCP rats; the PWMT of BCP rats declined consistently until 10 days after tumor inoculation and maintained at relatively low level, while PWMT of sham rats remained unchanged. The differences in PWMT between the two groups were significant on day 6, 8, 10, 12 and 14 after modeling (P < 0.01 and P < 0.001 respectively), indicating the development of mechanical hypersensitivity in BCP rats (Fig. 1a, Additional file 1 Table S1). Radiological and pathological examination of tibia bone demonstrated consistently severe cortical bone destruction and disarrangement of bone trabecula in BCP rats compared with sham-operated rats indicating osteolytic lesions after modeling (Fig. 1b, c).
Transcriptomic high-throughput sequencing revealed changes in gene expression profiles
A total of 506,546,810 raw reads were obtained using the Illumina HiSeq X-ten platform and 485,363,772 clean reads were generated after removing adaptor, short fragment reads and other low-quality sequences with clean ratio more than 95% in all samples. After ribosome RNA reads trimmed, 440,181,683 reads were mapped to the rat Rnor_6.0 reference genome among 479,470,030 reads in all with an average mapping rate of 91.8% (Additional file 2 Table S2). The raw data and processed data have been deposited in the National Center for Biotechnology Information’s Gene Expression Omnibus (GEO), with the accession number GSE137326.
In each sample, the gross transcripts were found to be almost equally distributed on all chromosomes (Fig. 2a). Pearson’s correlation coefficient of whole transcripts expression levels among all samples were calculated and visualized by heatmap of inter-sample correlation, demonstrating the obvious difference between BCP and sham-operated rats, but slight difference within each group (Fig. 2b). In total, 30,957 mRNAs and 22,783 lncRNAs were identified. With the filtering criteria of P < 0.05, FDR < 0.05 and fold change > 2, 1220 differently expressed mRNAs (DEmRNAs) (1171 up-regulated and 49 down-regulated) and 323 differently expressed lncRNAs (DElncRNAs) (246 up-regulated and 77 down-regulated) were detected after quantification of gene expression and comparison between the two groups (Fig. 2c). Hierarchical clustering of the expression of mRNA and lncRNA showed obvious discrimination between BCP and sham-operated rats (Fig. 2d, e). The count and fold change of DEmRNAs and DElncRNAs were visualized by volcano plot (Fig. 2f, g) and the detailed information of top 20 up-regulated and top 20 down-regulated mRNAs and lncRNAs were listed in Table 2 and Table 3. The total identified mRNAs/DEmRNA and lncRNAs/DElncRNA covered location of all chromosomes (Additional file 3 Figure S1 and Additional file 4 Figure S2). These results indicated that tumor inoculation in tibia shifted expression profiles of lncRNA and mRNA in the ipsilateral lumbar spinal cord, which may account for the development of hyperalgesia of BCP rats.
Validation of expression changes of DElncRNAs and DEmRNAs by RT-qPCR assay
To validate the reliability of sequencing quantification, 10 DEmRNAs (5 up-regulated and 5 down-regulated) and 10 DElncRNAs (5 up-regulated and 5 down-regulated) were selected to undergo RT-qPCR on another batch of samples independent of the ones used for RNA-seq (n = 6 per group). As shown in Fig. 3 and Additional file 5 Table S3, quantification by RT-qPCR of all the ten DEmRNAs were consistent with those of RNA sequencing; Class II, major histocompatibility complex, transactivator (Ciita), Cd74, Cytochrome b-245 beta chain (Cybb), Transient receptor potential cation channel, subfamily M, member 2 (Trpm2), annexin A3 (Anxa3) were all significantly up-regulated and AABR07031184.1, cytokine receptor-like factor 1 (Crlf1), heat shock protein 1B (Hspa1b), troponin C type 2 (fast) (Tnnc2), CDC42 effector protein (Rho GTPase binding) 2 (Cdc42ep2) significantly down-regulated. And all the five DElncRNAs, namely NONRATT007487.2, NONRATT003582.2, NONRATT026544.2, NONRATT004661.2, NONRATT008764.2 were up-regulated significantly and four out five DElncRNAs, namely MSTRG.12616.2, MSTRG.13351.2, MSTRG.16806.2, MSTRG.29385.15 were down-regulated, consistent with that of RNA sequencing.
GO and KEGG analysis of DEmRNAs and DElncRNAs
GO enrichment analysis of the DEmRNAs showed that the DEmRNAs were mainly enriched in the following biological process (BP) terms: inflammatory response, immune response, response to lipopolysaccharide, and in the following cellular component (CC) terms: external side of plasma membrane, cell surface, immunological synapse, and in the following molecular function (MF) terms: cytokine receptor activity, receptor activity, chemokine activity, peptide antigen binding (Fig. 4a).
GO analysis of potential target genes of DElncRNAs showed that they were enriched in the following BP terms: defense response, cellular response to interferon-beta, innate immune response, and in the following CC terms: nucleolus, membrane, integral component of endoplasmic reticulum membrane, and in the following MF terms: ATP binding, metal ion binding, nucleic acid binding (Fig. 4b).
Similarly, KEGG analysis of these dysregulated genes showed that neglecting apparently unrelated pathways, DEmRNAs were primarily enriched in allograft rejection, antigen processing and presentation, primary immunodeficiency, etc. (Fig. 4c); and potential target genes of DElncRNAs were primarily enriched in allograft rejection, antigen processing and presentation, natural killer cell mediated cytotoxicity, etc. (Fig. 4d). These results also indicated the involvement of extra−/intra-cellular pathways in inflammatory and immunologic process. The overlap between the GO enrichment and KEGG analysis of DEmRNAs and DElncRNAs suggested that there might be many interactions between DEmRNAs and DElncRNAs in the development of BCP.
The top 200 DEmRNAs were mainly expressed in microglia and interacted extensively
As shown in the expression heatmap of DEmRNA among 7 main types of the cells in central nervous system (Additional file 6 Figure S3), the majority of DEmRNAs, especially the top 200, were highly expressed in microglia, in spite of a few exceptions, indicating the critical role of microglial dysfunction in BCP development. By comparison with another two similar researches, 64 DEmRNAs identically overlapped with peer results and much more genes were functionally related (Fig. 5a). PPI network constructed by the top 200 DEmRNA-coding protein showed extensive interaction among them and four subnetworks were clustered, which may represent more compact functional blocks (Fig. 5b). Functional analysis of these genes showed that they were mainly concentrated in immunological processes/pathways, indicating they were highly involved in the pathogenesis of bone cancer (Fig. 5c, Additional file 7 Table S4).
Co-expression analysis of DElncRNAs/mRNAs
A delicate co-expression network was constructed based on DElncRNAs and their potential target DEmRNAs, with 131 lncRNAs, 137 mRNAs and 350 edges (Fig. 6). Most of DElncRNAs and their trans-regulated DEmRNAs were interweaved with the hub lncRNAs, such as ENSRNOT00000080979, NONRATT012538.2, MSTRG.12616.2, NONRATT026544.2, NONRATT022653.2, NONRATT001625.2, NONRATT017037.2, and the hub mRNAs, such as Cybb, Tap1, Srgn, Ccr5, RT1-DOa, Msn, Dsg3, CD80, CD247, indicating the involvement of inflammatory and immunologic process, oxidative stress and intercellular communication in BCP. Meanwhile, major DElncRNAs and their cis-regulated DEmRNAs were separated pairs, except CD74, Pou2f2, Cxcl13, Anxa3, etc. regulated by two or three lncRNAs.
Competing endogenous RNA (ceRNA) analysis of DElncRNAs, miRNAs and DEmRNAs
As lncRNA may affect the mRNA translation by sponging miRNA, the ceRNA network was constructed based on sequence similarity, maximum binding free energy and expression consistency, and 34 DElncRNAs, 13 miRNAs and 192 DEmRNAs with 287 edges were finally involved (Fig. 7). The majority of ceRNA pairs may bind competitively the following miRNA, miR-383-3p, miR-146-3p, miR-146-5p, miR483-5p, miR-292-5p, miR-18a-3p, miR-135b-3p, miR-130b-5p, indicating their core roles in the potential ceRNA regulation mechanisms.
By high-throughput sequencing of the transcriptome, we found 1220 differently expressed mRNAs (DEmRNAs) (1171 up-regulated and 49 down-regulated) and 323 differently expressed lncRNAs (DElncRNAs) (246 up-regulated and 77 down-regulated) in the lumbar spinal cord of BCP rats established by implantation of Walker 256 carcinoma cells into tibia; then through bioinformatic analysis of differently expressed genes, we found that a considerable proportion of them were highly expressed in microglia and extensively involved in immunological and inflammatory processes/pathways and partial DElncRNAs and DEmRNAs may have direct or indirect regulatory relationship.
As previously reported [32,33,34], the inoculation of Walker 256 breast cancer cells into medullary cavity of tibia could cause significant bone destruction and BCP response in rats, and could induce microglia and astrocyte activation in the spinal cord. As the site of secondary neuron of pain processing, spinal cord dorsal horn integrates and modulates noxious stimuli signal and extensive molecular, cellular and functional changes occur in cancer bone metastasis conditions. LncRNA could regulate coding genes expression and function [35, 36]; however, researches about their roles in BCP are still rare and limited. Given that the biological function of most lncRNA is still unknown, RNA sequencing of lncRNA and mRNA was performed to obtain a global landscape of gene expression alteration in the present study. Furthermore, bioinformatic analyses were conducted to deduce possible functions and regulatory coding genes of DElncRNAs attempting to finding new interventional targets of BCP.
Functional analysis of DEmRNAs and DElncRNAs indicated that they are primarily involved in inflammatory and immunological responses, represented by Class II major histocompatibility complex transactivator (CIITA) and CD74. It has been reported that robust activation of microglia and astrocytes in the spinal cord in BCP model and inhibition of spinal neuroinflammation could significantly attenuate BCP, probably through the means of AMPK activation , reactive oxygen species (ROS) scavenging , leucine-rich repeat and pyrin domain containing protein 3 (NLRP3) inhibition , etc. The central nervous system (CNS) may be not a completely immunological privileged organ especially in certain pathologic conditions; thus CD8+ T cell could infiltrate into the spinal cord in BCP , maybe by the structural and functional CNS lymphatic vessels . In our differential expression list, CIITA is a key regulator of MHC II expression, which mediates antigen presentation and the reactivation of CNS-infiltrated encephalitogenic T cells in experimental allergic encephalomyelitis  and multiple sclerosis . In the study of BCP, Song, Z. et al. has proved CIITA is upregulated by MEK/ERK/STAT1 axis and facilitates MHC II RT1B expression in microglia of spinal cord . Another immune molecule, CD74, as chaperone that regulates antigen presentation for immune response and cell surface receptor for the cytokine macrophage migration inhibitory factor (MIF), is upregulated in our results. While in the study of formalin-induced inflammatory pain, it was reported that CD74 and MIF were both upregulated in the ipsilateral spinal cord dorsal horn . CD74 also participates in sciatic chronic constriction nerve injury-evoked neuropathic pain by activation of ERK1/2-NMDA receptor and/or -PGE2 cascade .
We also identified significant upregulation of a series of oxidative stress related genes in BCP, represented by Cytochrome B-245 Beta Chain (CYBB) and Transient Receptor Potential Cation Channel Subfamily M Member 2 (TRPM2). Oxidative stress, which is mainly carried out by reactive oxygen species (ROS) and heightens during inflammation, is involved in the pathogenesis of BCP . As the main sources of ROS in CNS, NADPH oxidase 2 (Nox2) has been proved to play a critical role in nerve injury induced microglial activation in spinal cord and subsequent expression of proinflammatory cytokines in neuropathic pain . As the catalytic subunit of NADPH oxidase complex 2, transcripts of CYBB (gp91phox) are usually used to evaluate the Nox2 level. TRPM2 is a Ca2+-permeable nonselective cation channel that is directly activated by free cytosolic ADP-ribose and indirectly activated by ROS such as hydrogen peroxide (H2O2) . TRPM2 is expressed in peripheral macrophages and spinal microglia and it mediates the pathogenesis of inflammatory and neuropathic pain through the aggravation of pronociceptive inflammatory responses . As a monitor of oxidative stress, TRPM2 is considered as a new target of a spectrum of neuroinflammation associated with CNS disorders, including chronic pain .
A few of lncRNAs function have been explored in pain research, such as KCNA2-AS, uc.48+, NONRATT021972, MRAK009713, XIST, CCAT1 ; but the function of most lncRNAs remains unknown. We performed GO enrichment and KEGG analysis of potential target genes of DElncRNAs to detect their general biological function, which showed high degrees of consistency with DEmRNAs. Then we constructed co-expression network and ceRNA network to explore the relationship between DElncRNAs and coding genes.
LncRNA could regulate gene expression in cis, which influence the expression and/or chromatin state of nearby genes via the sequence-dependent manner (recruiting regulatory factors), transcription and/or splicing of the lncRNA dependent regulatory manner, or functional DNA elements within lncRNA loci . The co-expression network and validated DEmRNAs and DElncRNAs indicated that NONRATT007487.2 may regulate hydroxyl carboxylic acid receptor type 2 (Hcar2) in cis, which was recently reported to be involved in neuropathic pain . And CD74 and Trpm2 may be regulated by their adjacent lncRNA shown in the co-expression network. LncRNA could also regulate gene expression in trans, which influence gene expression throughout the cell, by regulating chromatin states and distant gene expression (for example, acting as a scaffold), influencing nuclear architecture, directly binding to protein and RNA . The co-expression network and validated DEmRNAs and DElncRNAs suggested NONRATT026544.2 may regulate the expression of Cybb, and NONRATT004661.2 may regulate Ncf1 or ccl9 in trans.
LncRNA has been acknowledged to sponge miRNA and therefore protect their target mRNAs from repression, called competing endogenous RNA (ceRNA). Those lncRNAs harboring multiple binding sites of identical miRNA are more likely candidates [52, 53]. CeRNA network hints that the upregulated lncRNA NONRATT007487.2 may interact with rno-miR-292-5p or rno-miR-383-3p, NONRATT004661.2 may interact with rno-miR-292-5p, consequently influencing the levels of the downstream mRNA.
However, we have to recognize that the quantification of RNA-seq could be affected by many factors though we proceeded strict quality control in each major step and the results were confirmed by traditional RT-qPCR. Moreover, gene expression fluctuations measured by RNA-seq could be partially reflected by cell population fluctuations (e.g. increase in microglia population in pathological pain [54, 55]), so it is necessary to perform single cell RNA sequencing (scRNA-Seq) or RNA-seq for a certain cell type (sorted by FACS, for example) to investigate gene expression in individual cell types. In addition, these bioinformatic analyses and predictions are theoretical, limited and maybe even controversial , therefore rigorous experimental verification is indispensable to obtain accurate conclusions.
To conclude, the sequencing analysis provides a landscape of dysregulated lncRNA and mRNA in the spinal cord of bone cancer pain; and functional analysis indicates differential expressed genes mainly concentrated in inflammatory and immunologic processes. Our present work may provide some insights and lay a foundation for future research to shed light on the pathogenesis of bone cancer pain.
Availability of data and materials
Please contact the author for data requests.
Bone cancer pain
Competing endogenous RNA
Differently expressed lncRNAs
Differently expressed mRNAs
False discovery rate
Kyoto Encyclopedia of Genes and Genomes
Long non-coding RNAs
Pearson correlation coefficient
Paw withdrawal mechanical threshold
Reactive oxygen species
Real-time quantitative PCR
Bennett MI, Kaasa S, Barke A, Korwisi B, Rief W, Treede RD. The IASP classification of chronic pain for ICD-11: chronic cancer-related pain. Pain. 2019;160(1):38–44.
van den Beuken-van Everdingen MH, Hochstenbach LM, Joosten EA, Tjan-Heijnen VC, Janssen DJ. Update on Prevalence of Pain in Patients With Cancer: Systematic Review and Meta-Analysis. J Pain Symptom Manage. 2016;51(6):1070–1090.e1079.
Kane CM, Hoskin P, Bennett MI. Cancer induced bone pain. Bmj. 2015;350:h315.
Figura N, Smith J, Yu HM. Mechanisms of, and adjuvants for, bone pain. Hematol Oncol Clin North Am. 2018;32(3):447–58.
Arthur J, Bruera E. Balancing opioid analgesia with the risk of nonmedical opioid use in patients with cancer. Nat Rev Clin Oncol. 2019;16(4):213–26.
Hu XF, He XT, Zhou KX, Zhang C, Zhao WJ, Zhang T, et al. The analgesic effects of triptolide in the bone cancer pain rats via inhibiting the upregulation of HDACs in spinal glial cells. J Neuroinflammation. 2017;14(1):213.
Chiou CS, Chen CC, Tsai TC, Huang CC, Chou D, Hsu KS. Alleviating bone Cancer-induced mechanical hypersensitivity by inhibiting neuronal activity in the anterior cingulate cortex. Anesthesiology. 2016;125(4):779–92.
Hua B, Gao Y, Kong X, Yang L, Hou W, Bao Y. New insights of nociceptor sensitization in bone cancer pain. Expert Opin Ther Targets. 2015;19(2):227–43.
Hou X, Weng Y, Ouyang B, Ding Z, Song Z, Zou W, et al. HDAC inhibitor TSA ameliorates mechanical hypersensitivity and potentiates analgesic effect of morphine in a rat model of bone cancer pain by restoring μ-opioid receptor in spinal cord. Brain Res. 1669;2017:97–105.
Hou X, Weng Y, Wang T, Ouyang B, Li Y, Song Z, et al. Suppression of HDAC2 in spinal cord alleviates mechanical Hyperalgesia and restores KCC2 expression in a rat model of bone Cancer pain. Neuroscience. 2018;377:138–49.
Bali KK, Kuner R. Noncoding RNAs: key molecules in understanding and treating pain. Trends Mol Med. 2014;20(8):437–48.
Zhao X, Tang Z, Zhang H, Atianjoh FE, Zhao JY, Liang L, et al. A long noncoding RNA contributes to neuropathic pain by silencing Kcna2 in primary afferent neurons. Nat Neurosci. 2013;16(8):1024–31.
Li G, Jiang H, Zheng C, Zhu G, Xu Y, Sheng X, et al. Long noncoding RNA MRAK009713 is a novel regulator of neuropathic pain in rats. Pain. 2017;158(10):2042–52.
Wu S, Bono J, Tao YX. Long noncoding RNA (lncRNA): a target in neuropathic pain. Expert Opin Ther Targets. 2019;23(1):15–20.
Li Z, Li X, Chen X, Li S, Ho IHT, Liu X, et al. Emerging roles of long non-coding RNAs in neuropathic pain. Cell Prolif. 2019;52(1):e12528.
Kim D, Langmead B, Salzberg SL. HISAT: a fast spliced aligner with low memory requirements. Nat Methods. 2015;12(4):357–60.
Pertea M, Pertea GM, Antonescu CM, Chang TC, Mendell JT, Salzberg SL. StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat Biotechnol. 2015;33(3):290–5.
Pertea M, Kim D, Pertea GM, Leek JT, Salzberg SL. Transcript-level expression analysis of RNA-seq experiments with HISAT, StringTie and Ballgown. Nat Protoc. 2016;11(9):1650–67.
Robinson MD, McCarthy DJ, Smyth GK. edgeR: A bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010;26(1):139–40.
Benjamini Y, Yekutieli D. The control of the false discovery rate in multiple testing under dependency. Ann Stat. 2001;29(4):1165–88.
Fang S, Zhang L, Guo J, Niu Y, Wu Y, Li H, et al. NONCODEV5: a comprehensive annotation database for long non-coding RNAs. Nucleic Acids Res. 2018;46(D1):D308–d314.
Kin T, Ono Y. Idiographica: a general-purpose web application to build idiograms on-demand for human, mouse and rat. Bioinformatics (Oxford, England). 2007;23(21):2945–6.
Tafer H, Hofacker IL. RNAplex: a fast tool for RNA-RNA interaction search. Bioinformatics. 2008;24(22):2657–63.
Huang da W, Sherman BT, Lempicki RA. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protoc. 2009;4(1):44–57.
Zhang Y, Chen K, Sloan SA, Bennett ML, Scholze AR, O'Keeffe S, et al. An RNA-sequencing transcriptome and splicing database of glia, neurons, and vascular cells of the cerebral cortex. J Neurosci. 2014;34(36):11929–47.
Costigan M, Moss A, Latremoliere A, Johnston C, Verma-Gandhu M, Herbert TA, et al. T-cell infiltration and signaling in the adult dorsal spinal cord is a major contributor to neuropathic pain-like hypersensitivity. J Neurosci. 2009;29(46):14415–22.
Du H, Shi J, Wang M, An S, Guo X, Wang Z. Analyses of gene expression profiles in the rat dorsal horn of the spinal cord using RNA sequencing in chronic constriction injury rats. J Neuroinflammation. 2018;15(1):280.
Zhou Y, Zhou B, Pache L, Chang M, Khodabakhshi AH, Tanaseichuk O, et al. Metascape provides a biologist-oriented resource for the analysis of systems-level datasets. Nat Commun. 2019;10(1):1523.
Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003;13(11):2498–504.
Bader GD, Hogue CWV. An automated method for finding molecular complexes in large protein interaction networks. BMC Bioinformatics. 2003;4:2.
Betel D, Wilson M, Gabow A, Marks DS, Sander C. The microRNA.org resource: targets and expression. Nucleic Acids Res. 2008;36(Database issue):D149–53.
Dong C, Wu R, Wu J, Guo J, Wang F, Fu Y, et al. Evaluation of bone Cancer pain induced by different doses of Walker 256 mammary gland carcinoma cells. Pain Physician. 2016;19(7):E1063–77.
Shenoy PA, Kuo A, Leparc G, Hildebrandt T, Rust W, Nicholson JR, Corradini L, Vetter I, Smith MT. Transcriptomic characterisation of the optimised rat model of Walker 256 breast cancer cell-induced bone pain. Clin Exp Pharmacol Physiol 2019;46(12):1201–15.
Dai J, Ding Z, Zhang J, Xu W, Guo Q, Zou W, et al. Minocycline Relieves Depressive-Like Behaviors in Rats With Bone Cancer Pain by Inhibiting Microglia Activation in Hippocampus. Anesth Analg. 2019. https://doi.org/10.1213/ANE.0000000000004063.
Kopp F, Mendell JT. Functional classification and experimental dissection of long noncoding RNAs. Cell. 2018;172(3):393–407.
Quinn JJ, Chang HY. Unique features of long non-coding RNA biogenesis and function. Nat Rev Genet. 2016;17(1):47–62.
Song H, Han Y, Pan C, Deng X, Dai W, Hu L, et al. Activation of adenosine monophosphate-activated protein kinase suppresses Neuroinflammation and ameliorates bone Cancer pain: involvement of inhibition on mitogen-activated protein kinase. Anesthesiology. 2015;123(5):1170–85.
Zhou YQ, Liu DQ, Chen SP, Sun J, Zhou XR, Rittner H, et al. Reactive oxygen species scavengers ameliorate mechanical allodynia in a rat model of cancer-induced bone pain. Redox Biol. 2018;14:391–7.
Chen SP, Zhou YQ, Wang XM, Sun J, Cao F, HaiSam S, et al. Pharmacological inhibition of the NLRP3 inflammasome as a potential target for cancer-induced bone pain. Pharmacol Res. 2019;147:104339.
Fu Q, Shi D, Zhou Y, Zheng H, Xiang H, Tian X, et al. MHC-I promotes apoptosis of GABAergic interneurons in the spinal dorsal horn and contributes to cancer induced bone pain. Exp Neurol. 2016;286:12–20.
Louveau A, Smirnov I, Keyes TJ, Eccles JD, Rouhani SJ, Peske JD, et al. Structural and functional features of central nervous system lymphatic vessels. Nature. 2015;523(7560):337–41.
Nikodemova M, Watters JJ, Jackson SJ, Yang SK, Duncan ID. Minocycline down-regulates MHC II expression in microglia and macrophages through inhibition of IRF-1 and protein kinase C (PKC)alpha/betaII. J Biol Chem. 2007;282(20):15208–16.
Liu X, Deng J, Li R, Tan C, Li H, Yang Z, et al. ERβ-selective agonist alleviates inflammation in a multiple sclerosis model via regulation of MHC II in microglia. Am J Transl Res. 2019;11(7):4411–24.
Song Z, Xiong B, Zheng H, Manyande A, Guan X, Cao F, et al. STAT1 as a downstream mediator of ERK signaling contributes to bone cancer pain by regulating MHC II expression in spinal microglia. Brain Behav Immun. 2017;60:161–73.
Wang F, Shen X, Guo X, Peng Y, Liu Y, Xu S, et al. Spinal macrophage migration inhibitory factor contributes to the pathogenesis of inflammatory hyperalgesia in rats. Pain. 2010;148(2):275–83.
Wang F, Xu S, Shen X, Guo X, Peng Y, Yang J. Spinal macrophage migration inhibitory factor is a major contributor to rodent neuropathic pain-like hypersensitivity. Anesthesiology. 2011;114(3):643–59.
Kim D, You B, Jo EK, Han SK, Simon MI, Lee SJ. NADPH oxidase 2-derived reactive oxygen species in spinal cord microglia contribute to peripheral nerve injury-induced neuropathic pain. Proc Natl Acad Sci U S A. 2010;107(33):14851–6.
Jang Y, Cho PS, Yang YD, Hwang SW. Nociceptive roles of TRPM2 Ion Channel in pathologic pain. Mol Neurobiol. 2018;55(8):6589–600.
Haraguchi K, Kawamoto A, Isami K, Maeda S, Kusano A, Asakura K, et al. TRPM2 contributes to inflammatory and neuropathic pain through the aggravation of pronociceptive inflammatory responses in mice. J Neurosci. 2012;32(11):3931–41.
Malko P, Syed Mortadza SA, McWilliam J, Jiang LH. TRPM2 channel in microglia as a new player in Neuroinflammation associated with a Spectrum of central nervous system pathologies. Front Pharmacol. 2019;10:239.
Boccella S, Guida F, De Logu F, De Gregorio D, Mazzitelli M, Belardo C, et al. Ketones and pain: unexplored role of hydroxyl carboxylic acid receptor type 2 in the pathophysiology of neuropathic pain. FASEB J. 2019;33(1):1062–73.
Rashid F, Shah A, Shan G. Long non-coding RNAs in the cytoplasm. Genomics Proteomics Bioinformatics. 2016;14(2):73–80.
Thomson DW, Dinger ME. Endogenous microRNA sponges: evidence and controversy. Nat Rev Genet. 2016;17(5):272–83.
Ding Z, Xu W, Zhang J, Zou W, Guo Q, Huang C, et al. Normalizing GDNF expression in the spinal cord alleviates cutaneous hyperalgesia but not ongoing pain in a rat model of bone cancer pain. Int J Cancer. 2017;140(2):411–22.
Xu F, Huang J, He Z, Chen J, Tang X, Song Z, et al. Microglial polarization dynamics in dorsal spinal cord in the early stages following chronic sciatic nerve damage. Neurosci Lett. 2016;617:6–13.
The authors thank Dr. Jun-Ming Zhang from University of Cincinnati for the comments of the manuscript.
This study was supported by grants from the Nature Science Foundation of Hunan Province (2018JJ3836 to Dr. Song), National Nature Science Foundation of China (81300958 to Dr. Song), Nature Science Foundation of China (81901143 to Dr. Ding) and Hunan Province Clinical Research Center for Anesthesia and Perioperative Medicine(2016Sk4003 to Dr.s Guo).
Ethics approval and consent to participate
All experiments were approved by the Animal Care and Use Committee of Central South University and performed following the guidelines of the International Association for the Study of Pain.
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Average paw withdrawal mechanical threshold (PWMT) of each animal at each timepoint.
Statistical data of high-throughput sequencing for six samples.
Total identified mRNAs and differential expressed mRNAs (DEmRNAs) mapped to the rat genome. Total identified mRNAs were marked by blue lines on the left side of chromosomes and DEmRNAs were marked by blue circles on the right side of chromosomes. The idiogram and GC content background data of rn6 rat genome were from Idiographica (Version 2.4) as illustrated in Materials and Methods.
Total identified lncRNA and differential expressed lncRNA (DElncRNA) mapped to rat genome. Total identified lncRNAs were marked by red lines on the left side of chromosomes and DEmRNAs were marked by red triangles on the right side of chromosomes. The idiogram and GC content background data of the rn6 rat genome were from Idiographica (Version 2.4) as illustrated in Materials and Methods.
The raw data of RT-qPCR validation assay.
The expression level of differential expressed mRNAs (DEmRNAs) in seven types of cells in central nervous system. Expression heatmap was constructed by our DEmRNA list and relative expression data from research of Ye Zhang et al. as illustrated in Materials and Methods. FPKM, Fragments Per Kilobase per Million.
Top 20 clusters with their representative enriched terms (one per cluster).
About this article
Cite this article
Hou, X., Weng, Y., Guo, Q. et al. Transcriptomic analysis of long noncoding RNAs and mRNAs expression profiles in the spinal cord of bone cancer pain rats. Mol Brain 13, 47 (2020). https://doi.org/10.1186/s13041-020-00589-2
- Bone cancer pain
- Long noncoding RNA
- High-throughput RNA sequencing