Common transcriptional signatures of neuropathic pain.

ABSTARCT: The dorsal root ganglia (DRG) are key structures in nociception and chronic pain disorders. Several gene expression studies of DRG in preclinical pain models have been performed, but it is unclear if consistent gene changes are identifiable. We, therefore, compared several recent RNA-Seq data sets on the whole DRG in rodent models of nerve injury. Contrary to previous findings, we show hundreds of common differentially expressed genes and high positive correlation between studies, despite model and species differences. We also find, in contrast to previous studies, that 60% of the common rodent gene response after injury is likely to occur in nociceptors of the DRG. Substantial expression changes are observed at a 1-week time-point, with smaller changes in the same genes at a later 3- to 4-week time-point. However, a subset of genes shows a similar magnitude of changes at both early and late time-points, suggesting their potential involvement in the maintenance of chronic pain. These genes are centred around suppression of endogenous opioid signalling. Reversal of this suppression could allow endogenous and exogenous opioids to exert their analgesic functions and may be an important strategy for treating chronic pain disorders. Currently used drugs, such as amitriptyline and duloxetine, do not seem to appropriately modulate many of the critical pain genes and indeed may transcriptionally suppress endogenous opioid signalling further.

1. Introduction

The mechanisms of neuropathic pain are complex and may vary considerably depending on the nature of the original nerve injury. These mechanisms have been investigated through different preclinical pain models (eg, injury or immune pain) and in different species (eg, mouse and rat).[3,17,36] Multiple groups of genes, including neuropeptides, ion channels, and genes involved in neuroimmunity, stress, and regeneration, are expressed in the peripheral nervous system (dorsal root ganglia [DRG]) and central nervous system (spinal cord and brain) and have been implicated in neuropathic pain. However, effective therapies for neuropathic pain remain elusive, perhaps due to differences in pain mechanisms between human and animals. In addition, tissue heterogeneity increases the complexity of interpreting neuropathic pain gene signatures. Thus, only 10% to 15% of DRG cells are sensory neurons, with the rest being glia and other cell types,[27,39] such that pain-specific expression change in nociceptors may be masked. Furthermore, pain-related genes, including ion channels, may be expressed at low levels, requiring sensitive methods for their detection, such as deep RNA-Seq.[31]

Changes in expression of DRG genes related to neuropathic pain have been extensively studied in rodent models using microarrays and subject to meta-analysis, providing early data on the transcriptional response to nerve injury.[17] However, microarray analysis has limited coverage in comparison with RNA-Seq, and to explore further the transcriptional signatures of neuropathic pain, we compared several recently published RNA-Seq data sets. We used rodent data on nerve injuries that were shown to have robust and reproducible chronic pain behaviours.[15] We chose the data sets where raw data were available and the sequencing had high coverage, aiming to assess whether a common gene response in neuropathic pain was evident. We found a highly consistent differential expression of DRG genes after nerve injury and explored functional classification of the common genes and their network. Also, we used the previously developed pain interactome approach to explore the protein–protein interaction (PPI) networks of pain-related genes.[16] We then focused on a subset of genes with similar magnitude of differential expression at early (1 week) and later (3 to 4 weeks) time-points after the nerve injury. This revealed a highly connected network centred around opioid signalling. Furthermore, we observed a substantial overlap between rodent pain network genes and human genes involved in chronic pain, suggesting potential usefulness of the large collection of rodent signatures in exploratory analysis and clinical applications to human neuropathic pain. Finally, we were able to use the rodent signatures to identify potential repurposing candidate drugs for the treatment of neuropathic pain.

2. Materials and methods

2.1. Analysis of differential gene expression

We analysed 4 rodent DRG RNA-Seq data sets on nerve injury: 2 mouse data set—C (from Cobos et al.[9]) and B (from Baskozos et al.[5]) and 2 rat data sets—P (from Perkins et al., 2014) and Br (from Baskozos et al.[5]). In these data sets, 2 injury models were used: spared nerve injury (SNI) for mouse data sets and spinal nerve transection (SNT) for rat data sets. The data were analysed from raw SRA data (for C (GSE102937) and P (GSE53861) data sets) or from reported counts (for B, Br data sets). The raw RNA-Seq data consisted of paired-end 69 bp reads for C and single-end 34 bp reads for P. After trimming adapters with cutadapt (M. Martin, https://doi.org/10.14806/ej.17.1.200), C and P samples had mean sequencing quality (Phred) scores per base between 31.4 and 40.3. The reads were mapped to the references genomes (mm10 for mouse; rn5 for rat) using STAR aligner.[10] Seventy-nine percent and 86% of reads in all samples were uniquely mapped for C and P. The mapped reads were counted to genes using the FeatureCounts function of the Subread package.[19] Differential expression was analysed using DESeq2 package in R.[23] Only genes with common names were retained for the further comparison between the data sets. In all cases, the false discovery rate (FDR) adjusted P values for multiple-hypothesis testing with the Benjamini–Hochberg method was ≤0.05 and absolute value of log2 fold changes, log2fc ≥ 0.5. C and P data sets had 3 samples per each condition (injured and control animals; C—one animal per sample; P—4 pooled animals per sample). Cobos samples had 4 runs each (17 millions reads), and these technical replicates were collapsed into the total sums of counts after the mapping; the samples from other rodent DRG data sets (P, B, and Br) had deep coverage (∼50 millions reads). B data set (bulb mouse strain) had 6 samples for control and 4 for injured animals (each sample—one animal); Br data set had 4 samples (3 pooled animal each) for both conditions. All data sets were based on male animals, except B, having both male and female mice. Two outlier samples (one injured and one control) were removed from the analysis of the Br data set based on the principle component analysis (PCA) plot of the Br samples. The retained samples (3 injured and 3 controls in both rat data sets) showed clear separation between injured and control groups (Fig. S1A, available at http://links.lww.com/PAIN/A972). After the outlier removal, all rat samples were analysed together using Deseq2 (similarly, mouse samples were analysed together, Fig. S1B, available at http://links.lww.com/PAIN/A972).

The technical variations between rat or mouse data sets/batches were minimized by Deseq2 package, using the design formula ∼ batch + condition. This correction of technical variations was visualized (Figs. S1 and S2, available at http://links.lww.com/PAIN/A972) by using the removeBatchEffect function from the limma package[35] on variance stabilizing transformed counts (using vst function). The expression values correlated between batches for both control and injured samples, which was improved after the correction of technical variations (Figs. S2C and D, http://links.lww.com/PAIN/A972). The correction had little effect on the correlation between control and injured conditions. Unsurprisingly, given the greater number of differentially expressed genes (DEGs) in the rat chronic pain model, the correlation in gene expression between control vs injured samples is slightly lower for rat than for mouse data sets (DEGs are red dots on Figs. S2C and D, available at http://links.lww.com/PAIN/A972). The genes, which are differentially expressed upon injury had high correlation between log2fc for both rat and mouse data sets (red dots on Figs. S2E and F, available at http://links.lww.com/PAIN/A972). However, genes not identified as differentially expressed showed a greater dispersion (black dots), which was especially pronounced in mouse data sets due to the lower proportion of DEGs in the whole transcriptome and the prevalence of non-DEGs with highly variable fold changes between data sets (Figs. S2E and F, http://links.lww.com/PAIN/A972).

To demonstrate the significance of the overlaps between gene lists, we performed permutation tests without replacement, using all genes with nonzero expression in at least one sample in any of the tested conditions. Thus, 19,312 expressed genes were used in the permutation test for mouse data sets and 13,940 genes for rat data sets. The differences in numbers of expressed genes between the data sets might be related to library composition and differences between species/models. The P value of the overlap was calculated using a hypergeometric distribution.

Nociceptor-related genes were identified in the whole DRG data sets using reported results of single-cell RNA-Seq analysis of DEGs in DRG sensory neurons during SNT injury in mice[14] (log2fc ≤1). The comparison of log2fc in expression of nociceptor-related DEGs between single-cell and bulk data sets demonstrated nearly-proportional changes between data sets, supporting nociceptor origin of the expression of these genes (Fig. S3, http://links.lww.com/PAIN/A972). To mitigate the differences in gene expression between bulk and single-cell data sets, we ranked the genes by the expression changes in single-cell data set multiplied to the differences in expression changes between single-cell and the whole DRG data set (Fig. S3A, http://links.lww.com/PAIN/A972). This ranking resulted in nearly-monotonous curves for the expression changes in both data sets, allowing visual comparison. We also used the above pipeline to analyse bulk RNA-Seq data on pooled DRG neurons during SNI injury in mice.[14] The pooled neurons data set had 6 samples per condition. After trimming of adapters, the mean quality scores were between 30 and 40.3, with ∼85% of the reads being uniquely mapped.

2.2. Enrichment and protein–protein interaction analysis

The analysis of significantly enriched GO (Gene Onthology) terms over-represented among DEG lists was performed with the ClusterProfiler package in R.[42] The significance of the overlap between DEG sets from the top GO terms was calculated using a hypergeometric distribution. This analysis demonstrated significant differences between most of the top DEG sets. The functional annotation clustering was performed with the TopGO package in R using conservative Elim algorithm, which minimizes false-positive results by removing annotated genes from ancestor GO terms.[1] The genefilter package was used to find background genes with expression similar to DEGs. The list of pain genes was constructed based on PainNetworks and PainGene databases[18,32] and literature curation. A PPI network was drawn using R package igraph based on the list of reported 1002 PPIs involved in pain.[16] Single edges not connected to the main network were removed. The list of pain-associated human genes was constructed based on the Genome-Wide Association Study (GWAS) catalogue[7] and recent review on genetic association studies.[46] The enrichment of the orthologous rat genes with pain-related PPI human networks was analysed using the Online Mendelian Inheritance in Man (OMIM)-extended library of human genetic diseases from the Enrichr web server.[8]

2.3. Analysis of differential genes against cellular responses to perturbation using connectivity map

The gene expression signatures of amitriptyline and duloxetine across the cell lines available on CMAP (connectivity map) (kidney, liver, epithelial, and others), for different doses and exposure times, were exported as gct files from the Touchstone database using the clue.io online service.[37] In-house R scripts were then used to import the saved gct files and remove all Z-scores that were not statistically significant (an absolute Z-score approximately greater than 1.65). The 37 rat persistent genes were extracted from cell lines treated with 10 μM of each drug of interest. Each gene was scored with either 1 for an anti-correlated gene expression direction, −1 for a correlated gene expression direction, or 0 if the gene was missing from that combination of dose and cell line or the Z-score was not significant. The magnitude of gene expression was not taken into account. The correlation results were filtered further for an exposure-to-drug time of 6 hours.

3. Results

3.1. Highly correlated gene signatures in 4 different rodent RNA-Seq dorsal root ganglion data sets of nerve injury models

We first examined 2 recent RNA-Seq data sets on DRG tissues derived from mouse SNI model of neuropathic pain. In these data sets, DRGs were collected at day 7 (Cobos [C] data set[9]) and day 28 (Baskozos [B] data set[5]) after the onset of the injury. The variations related to the technical differences between B and C data sets were corrected during the differential expression analysis using the Deseq2 package in R. To visualize the effect of technical variation, we used the removeBatchEffect function from the limma package.[35] Plotting the gene loadings before the correction showed clear separation between top PC1 and PC2 genes (Fig. S1C and Table S1, available at http://links.lww.com/PAIN/A972). The genes with highest PC1 loadings are mainly related to technical variations (Fig. S1A, available at http://links.lww.com/PAIN/A972), and 497 out of 500 top PC1 genes were not present among the top 500 PC1 genes after the correction (Fig. S1D, left, available at http://links.lww.com/PAIN/A972). Noticeably, the top PC2 genes related to the effect of injury before the correction become the top 500 PC1 genes after the correction (Fig. S1D, right; Fig. S1A, available at http://links.lww.com/PAIN/A972). The PCA plot showed clear separation between control and injured groups of samples (Fig. 1A and Fig. S1A, available at http://links.lww.com/PAIN/A972). Differential expression analysis was performed with DESeq2, using FDR-adjusted P ≤ 0.05 and log2fc ≥ 0.5, identifying 966 DEGs for C and 474 for B. The Venn diagram of Figure 1B demonstrates a substantial overlap between the data sets, with 233 common genes (Table S2, http://links.lww.com/PAIN/A972), significantly larger than expected by chance (Fig. S2A, http://links.lww.com/PAIN/A972). We found that changes in differential expression (log2fc of the mean) of the common DEGs were highly positively correlated between the data sets (Fig. 1C, R ∼0.9, P < 2.2·10−16). Similarly, there was a high correlation (R ∼0.89, P < 2.2·10−16) between log2fc calculated using the spread of log2fc between replicates (Fig. S2G, http://links.lww.com/PAIN/A972). The correlation of log2fc between the genes identified as DEGs in at least one data set was lower (R ∼0.6, Fig. S2I, http://links.lww.com/PAIN/A972); therefore, we used only overlapping DEGs for downstream analysis. The correlation across all genes was even lower (R ∼0.4, Fig. S2E, http://links.lww.com/PAIN/A972), with many genes having high fold changes in early (C), but not in a late (B) data set.

Figure 1.

Comparison on the common DEGs in rodent DRG data sets on the injury pain. (A) Principle component analysis plot of variance-stabilized transformed counts for mouse: Cobos (C), day 7, and Baskozos (B), day 28. (D) Principle component analysis plot of variance-stabilized transformed counts for rat: Perkins (P), day 7, and Baskozos (Br), day 21. To illustrate the correction of technical variations between the data sets, which was performed internally by Deseq2, the PCAs were generated after applying the removeBatchEffect function from the limma package. (B and E) Venn diagrams of DEG numbers in mouse (B) and rat (E) data sets, showing the numbers of overlapping genes. (C and F) Linear correlation between log2fc in expression levels relative to control noninjured DRG for the common genes (adjusted P ≤ 0.05) in mouse (C) and rat (F) data sets. (G) Linear correlation of the expression changes in common genes between at least one mouse and rat data set (|log2fc| ≥ 0.5; adjusted P ≤ 0.05). The sets of the common genes were split into correlated and anticorrelated groups as indicated in the legend. DEG, differentially expressed gene; DRG, dorsal root ganglia; PCA, principle component analysis.

We next examined RNA-Seq data from 2 independent SNT nerve injuries in rats. In these data sets, DRG tissues were collected on day 7 (Perkins [P] data set[31]) and day 21 (Baskozos rat [Br] data set[5]) after injury. For differential expression analysis, we used the same criteria as for the mouse data sets (adjusted P ≤ 0.05; log2fc ≥ 0.5). After the removal of technical variations (Figs. S1B, C, and E and Table S1, available at http://links.lww.com/PAIN/A972), PCA plots again showed a clear separation between control and injured group of samples (Fig. 1D and Fig. S1B, http://links.lww.com/PAIN/A972). The total numbers of DEGs were 2950 for Br and 5154 for P, with 2421 common genes between the data sets (Fig. 1E and Table S3, available at http://links.lww.com/PAIN/A972). The size of the overlap between the rat data sets was significantly larger than would be expected by chance (Fig. S2B, http://links.lww.com/PAIN/A972). There was again a very high positive correlation (R ∼0.9) between expression changes of the common rat DEGs (Fig. 1F and Figs. S2 H and J, available at http://links.lww.com/PAIN/A972). The correlation between expression changes of the whole transcriptome was lower (Figs. S2E and F, available at http://links.lww.com/PAIN/A972), due to the effect of non-DEGs. The ∼10 times higher number of DEGs in rat compared with mouse data sets leaded to higher proportion of DEGs in the whole genome and, hence, higher correlation for all genes in rat compared with mouse (Figs. S2E and F, http://links.lww.com/PAIN/A972). The higher number of DEGs in rat data sets was possibly related to differences in the technical variations between laboratories, the injury models, or less likely species differences.

The comparison of mouse and rat DEGs showed that 89 out of 233 common mouse genes were also differentially expressed in both rat data sets (Table S4, http://links.lww.com/PAIN/A972), with 655 genes being present in at least one mouse and one rat data sets (Table S5, http://links.lww.com/PAIN/A972). The expression changes of the common 655 mouse/rat genes were linearly correlated (R ∼0.7, P < 2.2·10−16; Fig. 1G). Interestingly, the changes in gene expression were more drastic in rat compared with mouse, possibly due to more severe SNT trauma. Between rat and mouse models, 68 of the 655 genes (10%) had opposite changes in gene expression (Fig. 1G and Table S5, http://links.lww.com/PAIN/A972). The STRING analysis of these negatively correlated genes identified several connections between them, involving Calb1, Gabrg1, and Slc1a1 and, interestingly, hormonal regulatory genes (Crh, Mc4r, and Tshr) (Fig. S4, http://links.lww.com/PAIN/A972).

3.2. Nerve injury dorsal root ganglion differentially expressed genes are enriched for ion channels and neuronal membrane function

To further characterize the functional properties of the common rodent DEGs (655 genes from Table S5, http://links.lww.com/PAIN/A972), we performed an enrichment analysis of ontology terms. The analysis of the functional clusters using ClusterProfiler showed over-representation of neuron- and membrane-related GO terms in the top GO category of the cell component (CC) (Fig. 2A and Fig. S5E, http://links.lww.com/PAIN/A972). The inflammatory response, response to wounding, blood circulation, sensory, and membrane processes were among the highly represented GO terms in the category of the biological processes (Fig. 2B and Fig. S5 F, http://links.lww.com/PAIN/A972). Also, functional annotation clustering using the TopGO package with the more conservative Elim algorithm[1] showed that the top GO categories include the integral component of membrane, extracellular space, neuron cell body, and plasma membrane (Fig. 2C and Table S7, http://links.lww.com/PAIN/A972). The common mouse (Figs. S5A, B, G, and H, http://links.lww.com/PAIN/A972) or common rat (Figs. S5C, D, I, and J, http://links.lww.com/PAIN/A972) DEGs were also enriched in neuronal- and membrane-related GO terms in the CC category and processes related to membrane transport, immunity, cell motility, and secretion in the biological process category. We concluded that rodent DRG signatures are enriched in neuron- and membrane-related genes, in agreement with previous observations.[9,40]

Figure 2.

The GO enrichment analysis of the common rat/mouse DEGs in nerve injury models (Table S5, http://links.lww.com/PAIN/A972). (A and B) The ClusterProfiler analysis of significantly enriched GO terms within the cell component (A) and biological processes (B) categories; FDR-adjusted P ≤ 0.05. (C) TopGO analysis of the common genes with Elim method. The top categories with FDR-adjusted P < 0.01 are shown in bright colours. DEG, differentially expressed gene.

3.3. Mouse dorsal root ganglion genes, which are differentially expressed after nerve injury, are largely expressed in nociceptors

The DRG is a complex tissue constituted by several cell types including sensory neurons, fibroblasts, immune cells, satellite glial cells, and Schwann cells, with neurons only representing ∼10% to 15% of DRG cells.[27,39] This poses a question about the relative contribution of nociceptor cells in the whole DRG transcriptome. The transcriptional responses of sensory neurons to SNT injury (day 3) were recently characterised in mice using single-cell RNA-Seq[14] on individually isolated DRG neurons, where neurons of different sizes were manually isolated and further classified based on the expression of the known markers. In this study, 2386 DEGs were reported (with log2fc ≥1) for 3 main subtypes of DRG sensory neurons: small nonpeptidergic (NP) nociceptors, medium peptidergic (PEP) nociceptors, and large myelinated (LM) sensory neurons. We used these data to estimate the numbers of genes related to nociceptive and sensory neurons among the common DEGs we identified above. We found that 59% of the common 233 mice DEGs from Table S2 (http://links.lww.com/PAIN/A972) were nociceptor (NP or PEP)-related. Nonpeptidergic genes represented the largest group (118 genes), followed by PEP (55 genes), with many genes belonging to several types of neurons (Fig. 3A). Only a small number of DEGs (29) belong to LM sensory neurons (Fig. 3A), consistent with predominant role of NP and PEP nociceptors in neuropathic pain.[14,38] Furthermore, the expression changes of nociceptor-related genes in the whole DRG data sets correlated with the reported expression changes in individually isolated nociceptor neurons (Figs. 3B),[14] supporting a significant contribution of nociceptors to the whole DRG transcriptomic response to injury. The STRING network analysis of the genes differentially expressed in each type of sensory neurons showed an interconnected network of genes implicated in nerve injury and regeneration (eg, Atf3, Sox11, Jun, Npy, and Gal) (Fig. 3C and Table S1, available at http://links.lww.com/PAIN/A972). The interactions in PEP- and LM-related genes involved less genes (Figs. 3D and E) but were centred around the subsets of the same genes (Gal, Atf3, Nts, and Csf1).

Figure 3.

Sensory neurons genes are differentially expressed in mouse DRG during nerve injury. The set of 233 common mouse DEGs from Table S2 (http://links.lww.com/PAIN/A972) was analysed for nociception-related genes using reported lists of mouse DEGs in DRG sensory neurons.[14] (A) Venn diagram of the DEG numbers related to NP, PEP, and LM sensory neurons, shown in blue, magenta, and green, respectively. (B) Linear correlation of the whole DRG expression changes in the nociception genes with those in individual sensory neurons (NP, PEP, and LM) during nerve injury.[14] (C–E) Protein–protein interaction analysis of biological interactions within common mouse DEGs related to NP (C), PEP (D), or LM (E) neurons. The analysis was performed with STRING, using high confidence of interactions (score 0.7) and not more than 5 interactors in the first shell. The line colour indicates the type of interaction evidence. The nodes were clustered using the MCL algorithm. Those genes that did not have a direct partner in the network are not shown. DEG, differentially expressed gene; DRG, dorsal root ganglia; NP, nonpeptidergic nociceptor; PEP, peptidergic nociceptor.

For the 2421 common rat DEGs, the proportion of genes that were nociceptor-related (by mouse classification) seems to be substantially smaller, at 25%. However, neuron-related GO terms were enriched among the remaining 74% of the rat DEGs that were not present in the results of the single-cell analysis of mouse neuronal DEGs after nerve injury (Fig. S6 and Table S8, http://links.lww.com/PAIN/A972). Unsurprisingly, 62% of the core 89 DEGs (common for all 4 rodent data sets) were related to nociceptors (Table S4, http://links.lww.com/PAIN/A972). Hence, overall, our analysis demonstrates that nociceptors or neurons significantly contribute to the transcriptomic response to injury, although neurons represent only 10% to 15% of DRG cells.[27,39]

3.4. Most of the differentially expressed DRG genes show consistent changes at early and late time-points after nerve injury

On the PCA plots, the separation between pain and control conditions was more pronounced for early time points (day 7) compared with the later time-points (day 21 or 28; Figs. 1A and D). Examining the magnitude of expression changes between early and late time-points, we found that in most cases, the absolute fold change of the gene expression was higher at early time-points (Figs. 1C and F), in agreement with the PCA plots. Interestingly, a subset of DEGs showed a similar magnitude of expression changes at the early and later time-points (| log2fc | at later time-point ≥ |log2fc| at early time-point). We found that these “persistent” gene sets were different for the mouse and rat data sets—40 out of 233 common mouse genes and 949 out of 2421 common rat genes, which might be attributed to the more severe SNT nerve injury. In most cases, genes were downregulated (22 for mouse and 662 for rat). A total of 50% of the mouse and 22% of the rat persistent genes were classified as neuronal by the mouse classification.[14] The enrichment analysis of the remaining 78% of rat persistent genes showed that the top DEGs (|log2fc| ≥ 1) are enriched in neuron-related GO terms (Fig. 4), while the DEGs with more modest expression changes (0.5 ≤ |log2fc| < 1) are enriched in mitochondrial terms (Fig. S7 and Table S8, http://links.lww.com/PAIN/A972). Only 6 genes showed persistent expression in all 4 rodent data sets: 2 upregulated genes—Lipn (a lipase) and Ms4a7 (potentially involved in signal transduction)—and 4 downregulated genes—Slca5a1 (an organic anion transporter), Frmpd4 (a positive regulator of dendritic spine morphogenesis), Vsnl1 (a neuronal calcium sensor), and Glb1l2 (a galactosidase). Overall, however, it seems that the mouse SNI and rat SNT persistent gene changes are divergent.

Figure 4.

Enrichment analysis of the persistent rat DEGs with pronounced expression changes, which are not present in mouse data set of sensory neurons. ClusterProfiler analysis of significantly enriched GO terms among DEGs with |log2fc| at 28 days (data set Br) ≥ than at 7 days (data set P). Only highly DEGs (|log2fc| ≥ 1; 304 genes total) were used. GO terms are shown within the CC (A) and BP (B) categories; FDR-adjusted P < 0.05. BP, biological process; CC, cell component; DEG, differentially expressed gene; FDR, false discovery rate.

3.5. Downregulation of the endogenous opioid system is a central feature in the rat pain model

To explore the pain-related PPI network, we used the previously published pain interactome—a PPI network of 611 interconnected pain-related proteins.[16] Examining all common rat DEGs revealed an interconnected network of 257 genes (Fig. S8, http://links.lww.com/PAIN/A972), including several hub genes related to immunity (Il6, IIl1a, Cx3cr1) that were mainly upregulated at day 7, but lower at day 21, which may be consistent with a time-dependent decrease of inflammation. Several regenerative hub genes of this network also show an initial increase followed by a decrease with time. This includes galanin (Gal), Sprr1a, Adcyap1; transcription factors Jun, Fos, and Sox11; and Avp (vasopressin growth factors) (Tables S1 and S2, available at http://links.lww.com/PAIN/A972). Temporally upregulated DEGs also include Pdyn—a precursor of κ-opioids, having antinociceptive functions in inflammatory, but not in neuropathic pain.[24] Likewise, transcription factor Atf3 is upregulated only at the early, but not late time-point (Table S4, http://links.lww.com/PAIN/A972). By contrast, we found that persistently downregulated genes include Oprm1 (µ-opioid receptor) and Penk (proenkephalin, encoding δ-opioids), both having antinociceptive effects in neuropathic pain[24,46], as well as upregulated Cck (cholecystokinin, encoding the opioid-antagonizing peptide). Importantly, Oprm1, Penk, and Cck were the key hub genes (Fig. S8, http://links.lww.com/PAIN/A972) in this network.

We next investigated the 949 persistent rat genes and found that this core network consists of 148 genes, including 37 persistent DEGs (13 upregulated and 24 downregulated) (Fig. 5), linking to a further 111 interaction partners (19 of the partners were nonpersistent DEGs, Table S9, http://links.lww.com/PAIN/A972). The key hub genes of this network were again Oprm1, Penk, and Cck. In addition to the global downregulation of the opioid system, the persistent network included several other pronociceptive changes, such as downregulation of antinociceptive glutamate receptor subunit Grin3a[26] and potassium channel regulators Kcnip3 and Kcnj3.[12,25] However, many changes in the persistent genes from the pain network involved antinociceptive mechanisms, such as upregulation of antinociceptive adrenergic and dopamine receptors Adra2a and Drd2[11,21] and downregulation of pronociceptive sodium channels Scn9a, Scn10a, Scn11a, Scn8a,[4] glutamate receptors Gria2 and Grm7,[13,29] aquaporin Aqp1,[6] ephrin receptor Ephb1,[20] tachykinin Tac1,[45] and vascular growth factor Vegfa.[22] Therefore, both antinociceptive and pronociceptive changes were involved in the rat pain network based on the persistent/chronic genes.

Figure 5.

The pain-related PPI network of the persistent rat DEGs. The network was built based on the pain interactome,[16] using rat DEGs with persistent expression changes, together with neighbours from the pain interactome network (DEGs and non-DEGs). The coloured arrows mark the type and direction of interactions; coloured nodes mark the expression changes in rat data sets (up/down/persistent/nonpersistent/non-DEG partner), as explained in the legend. DEG, differentially expressed gene; PPI, protein–protein interaction.

3.6. Rat genes with persistent changes are involved in human neuropathic pain

To explore the potential relevance of rat genes to human pain, we analysed several different sources of information about human pain genes. First, we found that 16 orthologues of the rat persistent genes and their neighbours on the pain network are associated with human pain in the GWAS catalogue (Table 1, top[7]). These include 9 persistent DEGs and 7 neighbour genes on the pain network (3 DEGs and 4 non-DEGs). Second, we found that 20 rat orthologues (7 persistent, 3 nonpersistent, and 10 neighbours) are involved in the human PPI pain network, based on the Online Mendelian Inheritance in Man-extended library in EnrichR (Table 1, middle[8]). Third, 5 orthologous genes (4 persistent and 1 nonpersistent) were found to be associated with monogenic disorders with mutations implicated in human neuropathic pain (Table 1, bottom[46]). And finally, 2 rat genes DRD2 and KCNJ3 were present among the 6 pain-associated genes, which were recently reconfirmed using the UKBiobank cohort.[46] In total, 42 orthologous rat genes (20 persistent) are implicated in human neuropathic pain.

Table 1

Persistent rat genes and their neighbours on the pain network implicated in chronic pain conditions in human.

3.7. Drugs used in human neuropathic pain modulate rodent chronic pain genes

We next explored the effect upon our identified persistent pain genes to drugs currently used for treatment of neuropathic pain in humans, such as amitriptyline and duloxetine. Because nociceptor-specific transcriptomic drug responses are not available, for this analysis, we used the gene expression reference data set Touchstone and the online clui.io toolset to examine drug effect in 9 available non-neuronal cell lines. The 37 rat persistent genes were isolated from these Touchstone expression signatures and partitioned by whether the drug response was correlated or anticorrelated (Figs. S9–S12, available at http://links.lww.com/PAIN/A972) to the rat-persistent gene expression (Fig. 6). Both amitriptyline and duloxetine demonstrated consistent correlated expression with many of the downregulated rat genes including those genes that may have an antinociceptive effect when downregulated (Scn10a, Scn11a, Scn9a, Tac1, Vegfa, Gria2, Grm7, Aqp1, and Ephb1). These drugs on the other hand were anticorrelated to the upregulated rat-persistent genes. The direction of the drug effects on the expression of the persistent genes was similar for most of the cell lines (Fig. 6), suggesting consistent responses, that may also occur in nociceptors.

Figure 6.

Cell-line responses to 10-μm amitriptyline and duloxetine for 6 hours relative to the differential gene expression of rat. For each rat gene, the number of anticorrelated and correlated cell-line responses from the drugs amitriptyline (A) and duloxetine (B) were collated. Genes in the rodent model whose changes are expected to have antinociceptive effects are indicated with an asterisk.

We next attempted to identify other drugs that had an anticorrelated signature to the rat pain gene expression changes. We first removed the 12 genes that we considered as having pain ameliorating effects (marked by * on Fig. 6) because we wanted to identify drugs that ideally did not modify these genes. Several drugs were identified including a number that increased endogenous opioid signalling, such as daunorubicin, doxorubicin, and others (Fig. S13, available at http://links.lww.com/PAIN/A972). However, all those drugs that increased opioid signalling genes also increased pronociceptive genes, particularly sodium channels. Future experiments on drug effects on nociceptor cell lines would allow for verification of our computational predictions and screen for drugs upregulating opioid-related genes, but not pronociceptive genes.

4. Discussion

Neuropathic pain is a highly disabling disorder, for which therapies are inadequate, and puts a considerable economic burden on society and the health service. Drug discovery efforts are costly and high risk because of a lack of clear understanding of the mechanisms of the development and maintenance of pain.[4] Transcriptomic studies have become a common approach to understand diseases and drug mechanisms, but direct applications of transcriptomics to study human pain are hindered by the inaccessibility of relevant cells and tissues (eg, DRGs). Preclinical animal models, therefore, have become an important tool for studying neuropathic pain, but the reproducibility of findings and the relevance to human mechanism have been called into question.

Previous comparison of transcriptional changes in neuropathic pain models has been performed based on microarray studies using different pain models (eg, nerve injury and inflammatory).[3,17] A meta-analysis[17] suggested only a narrow set of common DEGs, dominated by the immune genes.[3,17] Furthermore, there was a minimal correlation between gene expression in different models and species (rat/mouse).[17] Although inflammation may represent an essential component of neuropathic pain, many other studies indicate an important role for changes in neuronal function after nerve injury.[34,40] Moreover, genome-wide association and related genetic studies implicate genes involved in neurotransmission, metabolism, and other nonimmune functions in human neuropathic pain.[46] In agreement with these studies, we found that transcriptome changes involve a major neuronal component underlying neuropathic pain after nerve injury.

As far as we are aware, we present the first study that compares multiple RNA-Seq data sets on DRG gene expression after nerve injury. Contrary to the analysis of microarray data, we find a substantial similarity in transcriptional signatures of injury pain in rodents, suggesting common mechanisms. Thus, 233 genes were common between 2 mice modes of the SNI injury and 2421 genes were common between 2 rat models of the SNT injury. The larger overlap in rat data sets was related to larger number of DEGs and large fold changes, which is likely related to species differences and also differences between models, because SNT injury is more severe compared with SNI. To verify the cell-type origin of the identified DEGs, we used previously published data on gene expression in individual DRG neurons after nerve injury.[14] We should note that this classification is based on relative enrichment of particular genes and as such is dependent on the analytical method and the fold-change cutoff. Nevertheless, the Hu et al. (2016) assignments do seem to be consistent with other published studies.

Genes identified as differentially expressed may well be shared between neurons and other cell types. However, postinjury changes in the expression of neuronal genes show strong correlation between whole DRG and individually isolated nociceptor neurons (Fig. 3B), with a relatively stable ratio of expression changes between these data sets (Fig. S3, http://links.lww.com/PAIN/A972), supporting a largely nociceptor origin of the gene expression changes in neuropathic pain. The ∼ 5-fold reduction of the magnitude of expression changes in the whole DRG relative to single neurons is probably related to the “diluting” effect of other cell types. Overall, our results demonstrate that most of the overlapping DEGs across data sets seem to be driven by neuronal responses, despite neurons representing less than 15% of DRG cells.[27,39] This suggests that targeting of DRG peripheral neurons is a relevant strategy for developing protective therapies after nerve injury and for neuropathic pain treatments.[33]

We find that early changes at 1 week after injury are pronounced and affect genes involved in neuronal function, nerve regeneration, and immunity, presumably as a protective transcriptional response to injury, in agreement with previous observation.[9,14] Interestingly, however, we found that expression changes in most of these genes have much lesser magnitude after 3 to 4 weeks of injury. It seems likely that many of these transiently altered genes are not involved in chronic pain but are the appropriate response to nerve injury. Since each of the studies we examined had only one time-point, it is possible that the differences between the early and later time-points are simply technical differences between studies. However, it is notable that the greater magnitude of gene changes at the early time-points were evident in both mouse and rat studies, and it is biologically plausible that nerve injury induces dramatic gene expression changes early, which then declines after the initial response.

Genes that show persistent changes at later weeks after injury may be more relevant to chronic pain mechanisms. Further studies with multiple time-points after nerve injury are now required to explore this further. Noticeably, the immune genes, Il1 receptor antagonist (Il1rn), and Il16, blocking cell cycle progression in resting T-cells, are present among the upregulated persistent genes as well as Nfkbia and Rela, suggesting that at least these immune genes may have relevance to later stages after injury. These genes might be involved in neuroimmune interactions, which are known to be important in neuropathic pain.[4] Interestingly, we find many of these persisting genes, eg, sodium channels Snc9a, Scn10a, Scn11a, and glutamate receptors Gria2 and Grm7, have previously been implicated in pronociceptive pain mechanism. Downregulation of these genes suggest an antinociceptive, protective transcriptional response. This highlights a challenge when interpreting these data because in some cases, it may be difficult to distinguish protective from deleterious changes. Antinociceptive genes of opioid signalling (µ- and δ-opioid receptors Oprm1 and Oprd1; and Penk, coding δ-opioids) are persistently downregulated, while Cck, encoding opioid-antagonising peptides, is persistently upregulated. The downregulation of the opioid analgesic system agrees with multiple observations on lack of efficacy of exogenous opioids in the treatment of human neuropathic pain.[2] Our analysis of the rat pain network suggests that suppression of opioid signalling is central in the network, with Oprm1 being a top hub gene (Fig. 5). The mechanism of downregulation of Oprm1 expression has been shown to be dependent on methylation of its promoter by Dnmt3a methyltransferase, which in turn is dependent on decreased expression of miRNA-134.[41,44] Consistently, we find Dnmt3a expression is upregulated after nerve injury in the 2 rat data sets (log2fc = 0.74/0.63 for Br/P[41]). In addition, histone methylation by Suv39h1 methyltransferase was shown to be involved in Oprm1 silencing,[43] and again, we find increased expression in 2 rat data sets (log2fc = 0.54/0.9 for Br/P). Therefore, several epigenetic mechanisms likely participate in the downregulation of Oprm1 in nerve injury models. The degree of network connections of opioid signalling system might be partially exaggerated due to the opioid system being well-studied. However, the downregulation of opioid signalling in neuropathic pain is biologically plausible[46] and aligns with the clinical experience that opioids are not effective in neuropathic pain.

We identified a number of genes from the rat SNT pain network that are implicated in human pain and genetic studies. Several of these genes have strategic locations on the rat pain network, representing hub (SCN10A, OPRD1, FOS, IL6, TNF, IL1A, NGF, and DLG2) or bottleneck (SCN11, OPRD1, and DRD2) genes (Figs. 5 and Fig. S8, http://links.lww.com/PAIN/A972). Notably, a recent analysis of DRG samples from operated cancer patients with neuropathic pain[28] identified downregulation of κ-opioid receptor Oprk1. Therefore, drugs upregulating expression of opioid-related genes and their upstream regulators might represent strong candidates for the treatment of neuropathic pain.

Finally, we explored transcriptional drug signatures of 2 commonly prescribed therapeutics for neuropathic pain, duloxetine and amitriptyline. Both drugs were shown to drive changes in pain-related genes, but mostly, they downregulate pronociceptive genes (eg, Scn11a) that are already driven in the appropriate direction to reduce pain responses after nerve injury in the rodent models. Interestingly, these drugs are able to revert the persistently upregulated neuropathic pain genes identified in our preceding analysis, including the immune-related genes, in most of the cell lines evaluated. However, they do not seem to upregulate opioid signalling. Further studies are now required to assess whether these drugs have similar effects on DRG cells. Collectively, the analysis suggests that at a transcriptional level, these drugs do not seem ideal to counter the peripheral nerve gene expression changes, and this may explain why they have limited benefit in neuropathic pain disorders.

Our attempt to identify other drugs that might reverse the rat-persistent pain signature proved to be challenging. Although we found that drugs, such as daunorubicin and doxorubicin, increased endogenous opioid genes, they did so at the cost of increasing sodium channels, which may, therefore, attenuate any pain-relieving effects. This suggests that sodium channel and opioid receptors may be coregulated. Indeed, multiple epigenetic mechanisms (eg, histone acetylation) have been found to be involved in silencing of opioid receptors (eg, Oprm1) and sodium channels (eg, Scn10a) in neuropathic pain.[30] Further efforts are required to investigate potential therapeutic approaches that can upregulate the endogenous opioid pathway in the DRG without increasing pronociceptive molecules such as voltage-gated sodium channel expression. We conclude that information derived from animal preclinical models of neuropathic pain can provide useful insights and relevance to human pain disorders and suggest therapeutic strategies to help alleviate this disabling disorder.

Conflict of interest statement

M.Z. Cader has received honoraria from Novartis, TEVA, Orion, and Eli Lilly. The remaining authors have no conflicts of interest to declare.

Supplemental digital content associated with this article can be found online at http://links.lww.com/PAIN/A972.