Distinct Gene Expression Patterns of Ion Channels and Cytokines in Rat Primary Sensory Neurons During Development of Bone Cancer and Cancer Pain

Abstract

Cancer and cancer pain processes a major clinical challenge and the underlined mechanisms of pathogenesis remain elusive. We examined the specific changes in the transcriptomic profiles in the dorsal root ganglion (DRG) neurons of rats with bone cancer and bone cancer pain (BCP) using RNA sequencing technology. The bone cancer and BCP was induced by tumor cells implantation (TCI) into the tibia bone cavity in adult female rats. One week after treatment, TCI caused up- and down-regulation of thousands of genes in DRG. These genes were mainly involved in the immune process, inflammatory response, and intracellular signaling transduction of carbohydrate and cytokine. The cAMP and calcium signaling pathways were the major processes in the initial responses. Differentially expressed gene (DEG) analysis further showed that the genes for ion channels increased during day 1-7, while the genes for cytokine signaling pathways sustainedly increased during day 7-14 after TCI. The time courses of gene expression for ion channels and cytokines support their distinct roles in the early induction and late maintenance of BCP development. In addition, among the top 500 up- and down-regulated genes, 80-90% were unique for bone cancer pain as well as neuropathic and inflammatory pain, while less than 2% were shared among the three different forms of pain. This study reveals the uniqueness of mechanisms underlying bone cancer with pain, which is, to a large extent, differently from pain after acute inflammatory and nerve injury and provides novel potential targets of DEGs for bone cancer with pain.

Keywords: RNA-sequence; allodynia; cancer pain; cytokines; dorsal root ganglion; hyperalgesia; ion channels; neuropathic pain.

FIGURE 1

TCI treatment-induced painful behaviors and changes in morphology and density of tibia. (A,B) Spontaneous painful behaviors manifested as flinching and movement evoked pain manifested as reduced limb use. Averaged measurements of day 7 and 14. Eight rats included in each group. Student’s t-tests, ** P < 0.01, ***P < 0.001 vs. sham group. (C,D) Mechanical (C) and thermal (D) hypersensitivity manifested as lowered mechanical threshold and shortened latency of thermal withdrawal. Eight rats included in each group. Two-way repeated measure ANOVA, ***P < 0.001 vs. sham group. (E) Morphological changes presented by micro-CT images of 3-dimension reconstruction graph (blue background) and 2-dimension photos (black background) by coronal and sagittal scanning of the ipsilateral and contralateral tibia in TCI rats on day 14. Arrowheads: white, the zone of less dense bone; red, newly formed bone tissue; red arrow: site of cancer cell injection; red line: cross-section. (F-K) Statistical analysis of bone parameters measured by microCT on day 7 (n = 3) and 14 (n = 6), separately. Student’s t-tests between the two groups measured on the same day, * P < 0.05, ** P < 0.01, ***P < 0.001 ipsilateral vs. contralateral to TCI. BV/TV = bone volume/total volume, bone volume fraction (F). Tb.N = trabecular number (G). Tb.Sp = trabecular separation, the thickness of cavities (H). Tb.Th = trabecular thickness, the average thickness of all bone voxels (I). BMD = Bone mineral density (J). Ct.Th = Cortical thickness, the average thickness of cortical bone (K).

FIGURE 2

The quality of RNA-seq data from each of the DRG samples. (A) Summary of raw RNA sequencing data set from16 samples, including raw reads number, clean reads number, clean reads ratio, mapping ratio, and percentage of clean reads as well as Q20 (Phred quality scores Q) and Q30. (B) Distribution of gene expression levels in each of the samples. RPKM: Reads Per Kilobase Million. (C) Heatmap of the correlations between every two samples from each of the groups with the Pearson test. (D) Principal component analysis (PCA) of all the samples.

FIGURE 3

Dynamic alterations of DEGs in DRG ipsilateral to TCI treatment. (A-F) Volcano plot of all the differentially expressed genes from pairs of comparable groups. Log2 (fold change) is plotted as the abscissa and –log10 (Corrected P Value) is plotted as the ordinate. Each dot represents a single gene (red dots = up-regulated genes; green dots = down-regulated genes; gray dots = genes with no significant difference). (G) Comparison of up- and down-regulated DEGs between groups. (H) A Venn diagram presents DEGs that are unique or shared in each of the paired groups.

FIGURE 4

Functional analyses of DEGs in DRG from different stages. (A–D) Gene ontology (GO) analysis showing the enrichment of DEGs in biological process, cellular component, and molecular function from pairs of comparable groups. The top 20 significantly enriched GO terms determined by –log10 (corrected P value) are plotted as the ordinate and the enriched gene number is plotted as the abscissa.

FIGURE 5

KEGG classifications of DEGs in DRG. (A–D) The comparison of pathway enrichment from pairs of comparable groups. The top 20 KEGG terms determined by the rich factor is plotted as the ordinate. The rich factor is plotted as the abscissa. The gene number involved in each pathway is represented by the size of the dots. ARVC: Arrhythmogenic right ventricular cardiomyopathy.

FIGURE 6

Dynamic gene expression profiles of ion channels during the development of BCP. (A) Heatmaps depicting expression patterns for genes encoding ion channels from pairs of comparable groups. Up- or down-regulated genes are presented as the indicated color bars (red to blue). (B) Histogram showing the statistics of DEGs for ion channels. The number and percentage of the altered gene, up- and down-regulated gene in total, with | log2 FC| ≥ 1 and | log2 FC| ≥ 2. (C,D) Venn diagrams shows numbers of up- (C) and down-regulated (D) DEGs that are unique or shared in each paired group. (E-G) Plot data showing the top 10 up- or down-regulated genes for ion channels, ranked by multiple of difference (only genes with | log2 FC| ≥ 1 were plotted). (H) Function analysis of top 20 DEGs for ion channels that are considered associated with pain perception. DEGs with different symbols from the encoded proteins include Kcnu1 (Kca5.1), Cngb3 (cyclic nucleotide-gated cation channel beta-3), Cacng6 (voltage-dependent calcium channel gamma-6 subunit), Kcna7 (Kv1.7), Cacng1 (voltage-dependent calcium channel gamma-1 subunit), Cacna1s (Cav1.1), Kcnk15 (potassium channel subfamily K member 15), Scn4a (Nav1.4), Cnga3 (cyclic nucleotide-gated cation channel alpha-3), Scnn1b (SCNEB), Kcnj16 (Kir5.1), Mcoln3 (TRPML3), Kcnj15 (Kir4.2), Kcnj2 (Kir2.1), Kcnj13 (Kir7.1), Kcnk6 (TWIK-1), Kcnk13 (THIK-1), Kcns2 (Kv9.2).

FIGURE 7

Dynamic gene expression profiles of cytokine-cytokine receptor interaction pathway during the development of BCP. (A) Heatmaps depicting expression patterns for genes encoding cytokine-cytokine receptor interaction pathway from pairs of comparable groups. Up- or down-regulated genes are presented as the indicated color bars (red to blue). (B) Histogram showing the statistics of DEGs for cytokine-cytokine receptor interaction pathway. The number and percentage of the altered gene, up- and down-regulated gene in total, with | log2 FC| ≥ 1 and | log2 FC| ≥ 2. (C,D) Venn diagrams present the number of up- (C) and down-regulated (D) DEGs that are unique or shared in each paired group. (E-G) Plot data showing the top 10 up- or down-regulated genes in cytokine-cytokine receptor interaction pathways, ranked by multiple of difference (only genes with | log2 FC| ≥ 2 were plotted). (H) Function analysis of top 20 DEGs in the cytokine-cytokine receptor interaction pathway that are considered associated with bone cancer pain, pain perception, as well as neuroinflammation or nerve injury. DEGs with different symbols from the encoded proteins include Acvr1c (activin A receptor type-1C), Il18rap (Interleukin 18 receptor accessory protein), Inhbe (inhibin βE subunit), Csf2rb (cytokine receptor common subunit beta).

FIGURE 8

Transcriptomic comparison of the up- and down-regulated genes in rats with different forms of pain each produced by TCI, CFA, and SNL. (A) Venn diagrams showing numbers of the total up-regulated genes that are unique or shared by the three different pain models. (B) Left, Venn diagrams showing the top 500 up-regulated genes that are unique or shared by the three different pain models. Right, details of the shared genes with each comparison are listed in the table. Blue and orange fonts highlight the gene symbols for ion channels and cytokine-cytokine receptor interaction pathways, respectively. (C) Venn diagrams showing numbers of the total down-regulated genes that are unique or shared by the three different pain models. (D) Left, Venn diagrams showing the top 500 down-regulated genes that are unique or shared by the three different pain models. Right, details of the shared genes with each comparison are listed in the table. Orange fonts highlight the gene symbols for cytokine-cytokine receptor interaction pathways. Genes with different symbols from the encoded proteins include Il1r2 (interleukin-1 receptor type 2), Il2ra (interleukin-2 receptor subunit alpha), Ccr1l1 (chemokine C-C motif receptor 1-like 1), Il2rb (interleukin 2 receptor subunit beta), Cacna1s (Cav1,1), Cacng1 (voltage-dependent calcium channel gamma-1 subunit), Tnf (tumor necrosis factor).