A Bioinformatics Study of Differentially Expressed Genes in Microarrays of Dorsal Root Ganglia from Rat Models of Neuropathic Pain

Background

Neuropathic pain is a complex clinical condition caused by lesions or diseases in the somatosensory nervous system [1]. Peripheral nerve injury is frequently complicated by neuropathic pain, which causes distress and significantly decreases the quality of life in patients [2]. There are limited treatments available for neuropathic pain. This type of chronic pain is usually poorly controlled, and the treatment is unsatisfactory in patients with peripheral neuropathy [1,2].

The precise molecular mechanisms for the development of neuropathic pain are still unclear. Cell bodies of the primary sensory nerves are located in the dorsal root ganglia (DRG) that transmit the noxious stimuli to the brain. It is well documented that the DRG are implicated in the development of neuropathic pain [3,4]. Following peripheral nerve injury, the DRG undergo dramatic cellular and molecular changes [1]. The spinal nerve ligation (SNL) rat model has been widely used to investigate the molecular mechanisms and therapeutic targets of neuropathic pain [5]. SNL results in significant behavioral signs, including mechanical, heat, and cold hypersensitivity. In addition to behavioral changes, SNL induces pro-inflammatory responses and ion channel alterations in DRG [1]. In recent years, high-throughput next-generation sequencing and transcriptome microarrays have been used to screen differentially expressed genes (DEGs) in the DRG of the SNL rat model [6–8]. These studies were performed in different laboratories, and the DRG samples were collected at various time points after SNL. Therefore, it is challenging to compare and summarize the results of these studies and to make a clear conclusion.

This study aimed to undertake bioinformatics analysis of DEGs in microarrays of DRG from the SNL and sham rat models of neuropathic pain, based on 4 Gene Expression Omnibus (GEO) datasets, GSE15041, GSE38038, GSE2884, and GSE24982. There were 2 recently published bioinformatics analyses of DEGs in microarrays of DRG using other types of animal models of neuropathic pain, including spared nerve injury and spinal nerve transection [9,10]. The present study may improve the understanding of molecular mechanisms in the development of neuropathic pain after peripheral nerve injury.

Material and Methods

DATA COLLECTION:

The microarray datasets GSE15041, GSE38038, GSE2884, and GSE24982 were downloaded from the GEO database (http://www.ncbi.nlm.nih.gov/geo/) [6–8]. The investigators collected the DRG samples on day 7 for GSE15041 and GSE38038 and on day 28 for GSE2884 and GSE24982 datasets. These datasets were generated by performing transcriptome microarray assays in the DRG collected from adult rats that underwent a L4/L5 SNL or a sham operation without ligation. In total, 7 biological repeats in the sham group and 7 in the SNL group were performed on day 7 following the procedure, while 7 biological repeats in the sham group and 9 in the SNL group were performed on day 28. The platforms that were used for gene expression arrays were Affymetrix rat genome 230 2.0 array, Affymetrix rat genome U34 array, and Rosetta/Merck rat 25k v1.2 microarray. Ethics approval was waived for bioinformatic analyses of published datasets.

DATA PROCESSING:

DEGs between SNL and sham DRG specimens were identified via GEO2R online tools (https://www.ncbi.nlm.nih.gov/geo/geo2r/) with log fold change greater than 1.2 or log fold change less than −1.2 and an adjusted P value less than 0.05. Then, we calculated the common DEGs (co-DEGs) among the GSE15041, GSE38038, GSE2884, and GSE24982 datasets. The results were presented in a Venn diagram by using a Venn diagram maker from the Van de Peer Lab.

FUNCTIONAL AND PATHWAY ENRICHMENT ANALYSES:

We performed Gene Ontology (GO) functional enrichment analyses, which includes the biological process (BP), molecular function (MF), and cellular component (CC). The enrichment analyses were performed for the co-DEGs between GSE15041 and GSE38038 for day 7 and the co-DEGs between GSE2884 and GSE24982 for day 28 by using the online bioinformatic tool, the Database for Annotation, Visualization, and Integrated Discovery (DAVID, https://david.ncifcrf.gov/). Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses of the co-DEGs were also performed. P values <0.05 were considered statistically significant. The bubble charts were created by using R software with ggplot2.

PROTEIN–PROTEIN INTERACTION ANALYSES:

We performed Protein–protein interaction (PPI) network analyses of co-DEGs for either day 7 or day 28 by using STRING database of functional protein association networks (http://string-db.org/). Interaction network analysis result charts (TSV files) with a combined score greater than 0.4 were downloaded. The downloaded TSV files were imported into Cytoscape software, an open-source bioinformatics software platform for visualizing molecular interaction networks, to visualize the PPI network. Molecular Complex Detection (MCODE), a Cytoscape plugin that detects network modules, was used to identify the first 3 most significant functional modules for the PPI networks. CytoHubba, a Java plugin for Cytoscape, was used to identify hub genes and sub-networks of the constructed PPI networks. The PPI network charts were created using Cytoscape.

TARGETED DRUG PREDICTION:

The hub genes identified by CytoHubba were submitted to the drug gene interaction database (DGIdb database, http://dgidb.org/) to predict potential drug-gene interactions. The results downloaded from the DGIdb database were visualized in Cytoscape.

Results

IDENTIFICATION OF DEGS:

For the DRG collected on day 7 after SNL, we identified 1207 DEGs from GSE15041, 1959 DEGs from GSE38038, and 75 co-DEGs between GSE15041 and GSE38038. For the DRG collected on day 28 after SNL, we identified 151 DEGs from GSE2884, 1098 DEGs from GSE24982, and 29 co-DEGs between GSE2884 and GSE24982. There was no co-DEG between the DRG collected on day 7 and day 28, suggesting that the molecular mechanisms of neuropathic pain are different between the early and late stages. The results are displayed in a Venn diagram (Figure 1).

FUNCTIONAL ENRICHMENT:

GO enrichment and KEGG analyses of the co-DEGs for days 7 and 28 were performed using DAVID. For the DRG collected on day 7, the enriched biological processes in GO analyses were ubiquitin-dependent protein catabolic process and chromosome segregation (Figure 2A); the enriched cellular components were cytosol and extracellular exosome (Figure 2B); the enriched molecular functions were ATP binding and GTPase activity (Figure 2C); and the enriched KEGG pathways were adherens junction and metabolic pathways (Figure 2D). For the DRG collected on day 28, the enriched biological processes in GO analyses were protein tetramerization and response to drug (Figure 3A); the enriched cellular components were membrane raft and extracellular region (Figure 3B); the molecular functions were protein homodimerization activity and protein binding (Figure 3C); and the enriched KEGG pathways were calcium signaling pathway, proteoglycans, and focal adhesion (Figure 3D).

PPI NETWORK ENRICHMENT:

PPI networks for the DRG collected on day 7 and day 28 after SNL are displayed in Figures 4A and 5A, respectively. The constructed PPI networks were further analyzed by using the MCODE plugin in Cytoscape. For the DRG collected on day 7, the first cluster was enriched in synaptic vesicle trafficking and exocytosis (Figure 4B) and the second cluster was enriched in interleukin-1- and tumor necrosis factor-involved inflammation (Figure 4C). For the DRG collected on day 28, the first cluster was enriched in mRNA splicing (Figure 5B) and the second cluster was enriched in intracellular calcium handling (Figure 5C). Hub genes in the constructed PPI networks were identified by using the CytoHubba plugin in Cytoscape. The hub genes for the DRG collected on day 7 included Snap25, Vamp2, Cplx1, Stx5, Stxbp1, Gosr2, Bet1, Gosr1, Napa, and Stx1a (Figure 4D), which primarily regulate synaptic vesicle trafficking and exocytosis. For the DRG collected on day 28, the hub genes included Sf3b1, Cd2bp2, Lsm4, Sart1, Prpf8, Prpf31, Prpf6, Lsm7, Eftud2, and Sf3b2 (Figure 5D), which primarily regulate mRNA splicing.

TARGETED DRUG PREDICTION:

Potential drugs targeting the hub genes were predicted by using the DGIdb database. Regarding the hub genes for the DRG collected on day 7, six drugs targeting SNAP25 (synaptosome-associated protein of 25 kDa) and VAMP2 (vesicle associated membrane protein 2) were predicted, namely diazoxide, abobotulinumtoxinA, incobotulinumtoxinA, botulinum toxin type A purified neurotoxin complex, botulinum toxin type B, and phorbol myristate acetate (Figure 6A). For the hub genes for the DRG collected on day 28, only Chembl1221944 targeting SF3B1 (splicing factor 3b subunit 1) was predicted (Figure 6B). AbobotulinumtoxinA and incobotulinumtoxinA have been used in clinical practice to treat neck pain due to neck muscle spasm. Diazoxide is a well-known oral medication that is used to treat hyperinsulinism and hypoglycemia but has never been used for pain management.

Discussion

In the present study, we performed comprehensive bioinformatic analyses of 4 transcriptome array datasets that were generated by 4 independent research groups. The datasets compared mRNA expression profiles in the DRG from rats with and without SNL in the early and late stages. We found that there were no common DEGs between the early and late stages. The co-DEGs in the early stage were mainly enriched in synaptic vesicle trafficking and exocytosis and pro-inflammatory responses, while the co-DEGs in the late stage were primarily enriched in intracellular calcium handling and insulin-like growth factor signaling pathways.

Mechanisms for the development and progression of neuropathic pain are not fully understood, which can partly lead to the poor management of patients with neuropathic pain. High-throughput transcriptome array is a promising strategy for screening key genes and pathways involved in the development of neuropathic pain. To reduce the between-laboratory biases, we combined 2 distinct microarray datasets at each time point for bioinformatic analyses. An unpredicted finding was that there were no common DEGs in the DRG between the early (day 7) and late (day 28) stages. This finding suggests that neuropathic pain in the early and late stages has distinct molecular mechanisms. Therefore, further mechanism and functional studies for neuropathic pain in the early and late stages should be performed separately.

In the early stage, co-DEGs of the DRG were majorly enriched in chromosome segregation and protein catabolic processes, cytosol and extracellular exosome components, ATP binding function, and metabolic pathways. These findings suggest that the DRG in the early stage after peripheral nerve injury mainly underwent a metabolically active cellular process. These co-DEGs formed 3 major clusters, including synaptic vesicle trafficking and exocytosis, interleukin-1- and tumor necrosis factor-involved inflammation, and cell autophagy. Snap25, encoding for synaptosome-associated protein 25, was the central hub gene in the cluster of synaptic vesicle trafficking and exocytosis. SNAP25 is a component of the trans-SNARE complex and is involved in the regulation of synaptic vesicle exocytosis during synaptic transmission [11]. This finding is significant as previous studies demonstrated that SNAP25 contributed to neuropathic pain [12]. It was reported that overexpression of SNAP25 in rats induced pain-responsive behaviors and that breaking down SNAP25 by using botulinum toxin A relieved mechanical allodynia and neuropathic pain after peripheral nerve injury in rat models [12–14]. Through drug-gene interaction analyses, the present study found that botulinum toxin and its derivatives are potential drugs for treating neuropathic pain by targeting SNAP25. We also identified Vamp2, encoding for vesicle-associated membrane protein 2, as one of the major hub genes in the cluster of synaptic vesicle trafficking and exocytosis, suggesting it plays a role in neuropathic pain. VAMP2 is also a main component of the SNARE protein complex and is involved in the fusion of synaptic vesicles with the presynaptic membrane [15–17]. A recent microRNA array study revealed that microRNA-34a is a significantly downregulated microRNA in the DRG of rats that underwent chronic constriction injury, and microRNA-34a can target VAMP2 to decrease its expression [18]. These findings suggest that SNAP25 and VAMP2 may play key roles in the development of neuropathic pain.

In the late stage, co-DEGs of the DRG were majorly enriched in protein tetramerization and drug response processes, extracellular and membrane raft components, protein homodimerization and binding functions, and calcium signaling and focal adhesion pathways. These co-DEGs formed 3 major clusters, including mRNA splicing, intracellular calcium handling, and insulin-like growth factor signaling pathways. Sf3b1, encoding for splicing factor 3B subunit 1, was one of the central hub genes in the cluster of mRNA splicing [19]. The role of SF3B1 in the development of neuropathic pain has not been reported to date and warrants further investigation. Tnnc2 and Tnni2, encoding for troponin subtypes, were central hub genes in the cluster of intracellular calcium handling. Tnnc2 and Tnni2 have also been identified in other models of peripheral nerve injury in rat models [20,21], but their role in the development of neuropathic pain has not been defined. The present study found that the co-DEGs were also enriched in insulin-like growth factors pathways, which are responsible for nervous system development and function [22]. Dysregulated IGFBP5 expression may contribute to diabetic axon degeneration [23]. These findings suggest that targeting insulin-like growth factors pathways may attenuate neuropathic pain likely through promoting recovery of the injured nervous system.

All 4 datasets we analyzed in the present study were generated from the SNL rat model of neuropathic pain. There are several different animal models used to explore the molecular mechanisms of neuropathic pain [24,25]. Two recent bioinformatics analyses were performed in microarrays of DRG from neuropathic pain animal models of spared nerve injury and spinal nerve transection, and their identified DEGs are different from those in the present study [9,10]. This inconsistency suggests that different animal models of neuropathic pain may have distinct molecular mechanisms.

One of the limitations of the present study is that only 2 datasets are available for each time point after SNL. The numbers of common DEGs among the microarrays are relatively small, indicating that variations in these datasets might be significant. Future studies involving more datasets may help solve this issue. The microarrays we analyzed in the present study were generated using whole DRG tissue, which contains not only primary sensory neurons but also other types of cells, such as satellite glial cells, vascular cells, and other mesenchymal cells. Therefore, it is not clear which change in which type of cells contributed to the identified DEGs. Future microarrays may need to be performed in purified sensory neurons. Another limitation is that the identified targets were not validated by experiments in the present study because of a funding issue. Validation should be carried out in vitro and in vivo.

Conclusions

In conclusion, the present study deepens our understanding of the molecular mechanism of neuropathic pain. The present study identifies several key genes, including Snap25, Vamp2, and Sf3b1, in neuropathic pain. This study also suggests potential pathways, including pro-inflammatory, intracellular calcium handling, and insulin-like growth factors, involved in the development of neuropathic pain.