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 package35 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 mice14 (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 databases18,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) catalogue7 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.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 set9) and day 28 (Baskozos [B] data set5) 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.
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 set31) and day 21 (Baskozos rat [Br] data set5) 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 algorithm1 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
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-Seq14 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).
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.
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 pain24,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 Grin3a26 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 Drd211,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.
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, top7). 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, middle8). 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, bottom46). 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.
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.
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.
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-analysis17 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/P41). 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 plausible46 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 pain28 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.
Appendix A. Supplemental digital content
Supplemental digital content associated with this article can be found online at http://links.lww.com/PAIN/A972.
The authors thank Prof Caleb Webber and Adam Handel for valuable comments on the manuscript.
The study has been funded by Oxford BRC and Oxford Cannabinoid Technologies Ltd.
Author contributions: A. Pokhilko performed analysis of transcriptional signatures; A. Nash performed drug analysis. All authors contributed to the writing of the manuscript. A. Pokhilko and M.Z. Cader conceived and designed the study. M.Z. Cader supervised the study.
. Alexa A, Rahnenfuhrer J, Lengauer T. Improved scoring of functional groups from gene expression data by decorrelating GO graph structure. Bioinformatics
. Arner S, Meyerson BA. Lack of analgesic effect of opioids on neuropathic and idiopathic forms of pain. PAIN 1988;33:11–23.
. Bangash MA, Alles SRA, Santana-Varela S, Millet Q, Sikandar S, de Clauser L, Ter Heegde F, Habib AM, Pereira V, Sexton JE, Emery EC, Li S, Luiz AP, Erdos J, Gossage SJ, Zhao J, Cox JJ, Wood JN. Distinct transcriptional responses of mouse sensory neurons in models of human chronic pain conditions. Wellcome Open Res 2018;3:78.
. Basbaum AI, Bautista DM, Scherrer G, Julius D. Cellular and molecular mechanisms of pain. Cell 2009;139:267–84.
. Baskozos G, Dawes JM, Austin JS, Antunes-Martins A, McDermott L, Clark AJ, Trendafilova T, Lees JG, McMahon SB, Mogil JS, Orengo C, Bennett DL. Comprehensive analysis of long noncoding RNA expression in dorsal root ganglion reveals cell-type specificity and dysregulation after nerve injury. PAIN 2019;160:463–85.
. Borsani E. Aquaporins in sensory and pain transmission. Curr Neuropharmacol 2010;8:122–7.
. Buniello A, MacArthur JAL, Cerezo M, Harris LW, Hayhurst J, Malangone C, McMahon A, Morales J, Mountjoy E, Sollis E, Suveges D, Vrousgou O, Whetzel PL, Amode R, Guillen JA, Riat HS, Trevanion SJ, Hall P, Junkins H, Flicek P, Burdett T, Hindorff LA, Cunningham F, Parkinson H. The NHGRI-EBI GWAS Catalog of published genome-wide association studies, targeted arrays and summary statistics 2019. Nucleic Acids Res 2019;47:D1005–12.
. Chen EY, Tan CM, Kou Y, Duan Q, Wang Z, Meirelles GV, Clark NR, Ma'ayan A. Enrichr: interactive and collaborative HTML5 gene list enrichment analysis tool. BMC Bioinformatics
. Cobos EJ, Nickerson CA, Gao F, Chandran V, Bravo-Caparros I, Gonzalez-Cano R, Riva P, Andrews NA, Latremoliere A, Seehus CR, Perazzoli G, Nieto FR, Joller N, Painter MW, Ma CHE, Omura T, Chesler EJ, Geschwind DH, Coppola G, Rangachari M, Woolf CJ, Costigan M. Mechanistic differences in neuropathic pain modalities revealed by correlating behavior with global expression profiling. Cell Rep 2018;22:1301–12.
. Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, Batut P, Chaisson M, Gingeras TR. STAR: ultrafast universal RNA-seq aligner. Bioinformatics
. Donello JE, Guan Y, Tian M, Cheevers CV, Alcantara M, Cabrera S, Raja SN, Gil DW. A peripheral adrenoceptor-mediated sympathetic mechanism can transform stress-induced analgesia into hyperalgesia. Anesthesiology 2011;114:1403–16.
. Guo YP, Zhi YR, Liu TT, Wang Y, Zhang Y. Global gene knockout of Kcnip3 enhances pain sensitivity and exacerbates negative emotions in rats. Front Mol Neurosci 2019;12:5.
. Hartmann B, Ahmadi S, Heppenstall PA, Lewin GR, Schott C, Borchardt T, Seeburg PH, Zeilhofer HU, Sprengel R, Kuner R. The AMPA receptor subunits GluR-A and GluR-B reciprocally modulate spinal synaptic plasticity and inflammatory pain. Neuron 2004;44:637–50.
. Hu G, Huang K, Hu Y, Du G, Xue Z, Zhu X, Fan G. Single-cell RNA-seq reveals distinct injury responses in different types of DRG sensory neurons. Sci Rep 2016;6:31851.
. Jaggi AS, Jain V, Singh N. Animal models of neuropathic pain. Fundam Clin Pharmacol 2011;25:1–28.
. Jamieson DG, Moss A, Kennedy M, Jones S, Nenadic G, Robertson DL, Sidders B. The pain interactome: connecting pain-specific protein interactions. PAIN 2014;155:2243–52.
. LaCroix-Fralish ML, Austin JS, Zheng FY, Levitin DJ, Mogil JS. Patterns of pain: meta-analysis of microarray studies of pain. PAIN 2011;152:1888–98.
. Lacroix-Fralish ML, Ledoux JB, Mogil JS. The Pain Genes Database: an interactive web browser of pain-related transgenic knockout studies. PAIN 2007;131:3.e1–4.
. Liao Y, Smyth GK, Shi W. The Subread aligner: fast, accurate and scalable read mapping by seed-and-vote. Nucleic Acids Res 2013;41:e108.
. Liu S, Liu WT, Liu YP, Dong HL, Henkemeyer M, Xiong LZ, Song XJ. Blocking EphB1 receptor forward signaling in spinal cord relieves bone cancer pain and rescues analgesic effect of morphine treatment in rodents. Cancer Res 2011;71:4392–402.
. Liu S, Tang Y, Shu H, Tatum D, Bai Q, Crawford J, Xing Y, Lobo MK, Bellinger L, Kramer P, Tao F. Dopamine receptor D2, but not D1, mediates descending dopaminergic pathway-produced analgesic effect in a trigeminal neuropathic pain mouse model. PAIN 2019;160:334–44.
. Llorian-Salvador M, Gonzalez-Rodriguez S. Painful understanding of VEGF. Front Pharmacol 2018;9:1267.
. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol 2014;15:550.
. Maldonado R, Banos JE, Cabanero D. Usefulness of knockout mice to clarify the role of the opioid system in chronic pain. Br J Pharmacol 2018;175:2791–808.
. Marker CL, Stoffel M, Wickman K. Spinal G-protein-gated K+ channels formed by GIRK1 and GIRK2 subunits modulate thermal nociception and contribute to morphine analgesia. J Neurosci 2004;24:2806–12.
. Mohamad O, Song M, Wei L, Yu SP. Regulatory roles of the NMDA receptor GluN3A subunit in locomotion, pain perception and cognitive functions in adult mice. J Physiol 2013;591:149–68.
. Ng KY, Wong YH, Wise H. The role of glial cells in influencing neurite extension by dorsal root ganglion cells. Neuron Glia Biol 2010;6:19–29.
. North RY, Li Y, Ray P, Rhines LD, Tatsui CE, Rao G, Johansson CA, Zhang H, Kim YH, Zhang B, Dussor G, Kim TH, Price TJ, Dougherty PM. Electrophysiological and transcriptomic correlates of neuropathic pain in human dorsal root ganglion neurons. Brain 2019;142:1215–26.
. Palazzo E, Romano R, Luongo L, Boccella S, De Gregorio D, Giordano ME, Rossi F, Marabese I, Scafuro MA, de Novellis V, Maione S. MMPIP, an mGluR7-selective negative allosteric modulator, alleviates pain and normalizes affective and cognitive behavior in neuropathic mice. PAIN 2015;156:1060–73.
. Penas C, Navarro X. Epigenetic modifications associated to neuroinflammation and neuropathic pain after neural trauma. Front Cell Neurosci 2018;12:158.
. Perkins JR, Antunes-Martins A, Calvo M, Grist J, Rust W, Schmid R, Hildebrandt T, Kohl M, Orengo C, McMahon SB, Bennett DL. A comparison of RNA-seq and exon arrays for whole genome transcription profiling of the L5 spinal nerve transection model of neuropathic pain in the rat. Mol Pain 2014;10:7.
. Perkins JR, Lees J, Antunes-Martins A, Diboun I, McMahon SB, Bennett DL, Orengo C. PainNetworks: a web-based resource for the visualisation of pain-related genes in the context of their network associations. PAIN 2013;154:2586 e2581–2512.
. Pettingill P, Weir GA, Wei T, Wu Y, Flower G, Lalic T, Handel A, Duggal G, Chintawar S, Cheung J, Arunasalam K, Couper E, Haupt LM, Griffiths LR, Bassett A, Cowley SA, Cader MZ. A causal role for TRESK loss of function in migraine mechanisms. Brain 2019;142:3852–67.
. Reinhold AK, Batti L, Bilbao D, Buness A, Rittner HL, Heppenstall PA. Differential transcriptional profiling of damaged and intact adjacent dorsal root ganglia neurons in neuropathic pain. PLoS One 2015;10:e0123342.
. Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, Smyth GK. Limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res 2015;43:e47.
. Rodriguez Parkitna J, Korostynski M, Kaminska-Chowaniec D, Obara I, Mika J, Przewlocka B, Przewlocki R. Comparison of gene expression profiles in neuropathic and inflammatory pain. J Physiol Pharmacol 2006;57:401–14.
. Subramanian A, Narayan R, Corsello SM, Peck DD, Natoli TE, Lu X, Gould J, Davis JF, Tubelli AA, Asiedu JK, Lahr DL, Hirschman JE, Liu Z, Donahue M, Julian B, Khan M, Wadden D, Smith IC, Lam D, Liberzon A, Toder C, Bagul M, Orzechowski M, Enache OM, Piccioni F, Johnson SA, Lyons NJ, Berger AH, Shamji AF, Brooks AN, Vrcic A, Flynn C, Rosains J, Takeda DY, Hu R, Davison D, Lamb J, Ardlie K, Hogstrom L, Greenside P, Gray NS, Clemons PA, Silver S, Wu X, Zhao WN, Read-Button W, Wu X, Haggarty SJ, Ronco LV, Boehm JS, Schreiber SL, Doench JG, Bittker JA, Root DE, Wong B, Golub TR. A next generation connectivity map: L1000 platform and the first 1,000,000 profiles. Cell 2017;171:1437–52.e1417.
. Tandrup T, Woolf CJ, Coggeshall RE. Delayed loss of small dorsal root ganglion cells after transection of the rat sciatic nerve. J Comp Neurol 2000;422:172–80.
. Thakur M, Crow M, Richards N, Davey GI, Levine E, Kelleher JH, Agley CC, Denk F, Harridge SD, McMahon SB. Defining the nociceptor transcriptome. Front Mol Neurosci 2014;7:87.
. Xiao HS, Huang QH, Zhang FX, Bao L, Lu YJ, Guo C, Yang L, Huang WJ, Fu G, Xu SH, Cheng XP, Yan Q, Zhu ZD, Zhang X, Chen Z, Han ZG, Zhang X. Identification of gene expression profile of dorsal root ganglion in the rat peripheral axotomy model of neuropathic pain. Proc Natl Acad Sci U S A 2002;99:8360–5.
. Xu B, Cao J, Zhang J, Jia S, Wu S, Mo K, Wei G, Liang L, Miao X, Bekker A, Tao YX. Role of MicroRNA-143 in nerve injury-induced upregulation of Dnmt3a expression in primary sensory neurons. Front Mol Neurosci 2017;10:350.
. Yu G, Wang LG, Han Y, He QY. clusterProfiler: an R package for comparing biological themes among gene clusters. Omics 2012;16:284–7.
. Zhang J, Liang L, Miao X, Wu S, Cao J, Tao B, Mao Q, Mo K, Xiong M, Lutz BM, Bekker A, Tao YX. Contribution of the suppressor of variegation 3-9 homolog 1 in dorsal root ganglia and spinal cord dorsal horn to nerve injury-induced nociceptive hypersensitivity. Anesthesiology 2016;125:765–78.
. Zhou XL, Yu LN, Wang Y, Tang LH, Peng YN, Cao JL, Yan M. Increased methylation of the MOR gene proximal promoter in primary sensory neurons plays a crucial role in the decreased analgesic effect of opioids in neuropathic pain. Mol Pain 2014;10:51.
. Zimmer A, Zimmer AM, Baffi J, Usdin T, Reynolds K, Konig M, Palkovits M, Mezey E. Hypoalgesia in mice with a targeted deletion of the tachykinin 1 gene. Proc Natl Acad Sci U S A 1998;95:2630–5.
. Zorina-Lichtenwalter K, Parisien M, Diatchenko L. Genetic studies of human neuropathic pain conditions: a review. PAIN 2018;159:583–94.