Identification of the key genes associated with neuropathic pain.

Neuropathic pain is a chronic pain state associated with multiple etiologies that results in considerable social and economic burden. The identification of key genes associated with neuropathic pain is important for the development of novel therapies. Therefore, the present study downloaded the gene expression profile GSE15041 from the Gene Expression Omnibus database. The unverified gene chip was removed and the microarray data was normalized following quality control. The limma package in R was used to screen the differentially expressed genes (DEGs), followed by Gene Ontology and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis. Furthermore, a protein‑protein interaction (PPI) network based on the identified DEGs was constructed to select hub proteins, and reverse transcription‑quantitative polymerase chain reaction was performed to detect the expression of these proteins in a mouse model of neuropathic pain. In total, 86 common DEGs were identified. DEGs were significantly enriched in ̔extracellular space̓ and KEGG pathway enrichment analysis demonstrated that the DEGs were significantly enriched in inflammatory diseases and the mitogen‑activated protein kinase signaling pathway. The PPI network consisted of 27 nodes (proteins) and 47 PPI edges (interactions). Interleukin (IL)‑6, transcription factor AP‑1 (c‑Jun) and urikinase‑type plasminogen activator (Plau) were identified as hub proteins and key genes in neuropathic pain. The mRNA expression of these hub proteins was significantly increased in the neuropathic pain model, compared with the sham group. IL‑6, c‑Jun, and Plau may be involved in development of neuropathic pain and further research investigating the exact role of these key genes is required.

Introduction

Neuropathic pain is a chronic pain state associated with multiple etiologies, including trauma, diabetes, cancer, infection and autoimmune pathology (1). The incidence of neuropathic pain is approximately 7–10% worldwide (2,3). Patients with neuropathic pain induced by nerve injury exhibit varying degrees of mechanical and heat hyperalgesia. Additionally, sleep and anxiety disorders are more prevalent in individuals with neuropathic pain, and quality of life is more seriously impaired, compared with individuals with chronic non-neuropathic pain (4,5). Chronic neuropathic pain also results in considerable social and economic burden, leading to high medical costs and reduced productivity (6,7).

The mechanism underlying neuropathic pain remains unclear and it is still difficult to treat, despite the development of various novel therapies (8). The rapid development of gene chip technology has become a technical basis for biological research; thousands of genes are able to be detected simultaneously with gene chip analysis. At present, large amounts of data have not been effectively explored, despite the availability of large public gene chip databases. Bioinformatics has revealed a considerable number of complex biological interactions through the comprehensive utilization of biology, computer science and statistics (9).

Vega-Avelaira et al (10) used a microarray expression profile to identify differentially expressed genes (DEGs) associated with neuropathic pain. However, the quality of one gene chip studied was not verified. Therefore, for quality control, it should be removed for reanalysis of the normalized unscaled standard errors (NUSE) and RNA degradation curves in the microarray expression profile, as performed by the affy package. The present study reanalyzed the gene expression profile from the Vega-Avelaira et al (10) study to further identify the DEGs associated with neuropathic pain. Function and pathway enrichment analyses of the DEGs were performed, followed by protein-protein interaction (PPI) network construction. Furthermore, reverse transcription-quantitative polymerase chain reaction was performed in a mouse model of neuropathic pain to verify the results obtained from the bioinformatics analysis.

Materials and methods

Microarray data

The gene expression profile GSE15041 was downloaded from the Gene Expression Omnibus (GEO) of the National Center for Biotechnology Information (http://www.ncbi.nlm.nih.gov/gds/) on the Affymetrix Rat Genome 230 2.0 Array (Thermo Fisher Scientific, Inc., Waltham, MA, USA). A total of 9 adult rat samples were analyzed in the present study, including 3 samples from the ipsilateral spinal cord of the sham group (Sham group), 3 samples from the contralateral spinal cord of the surgery group (Contra group) and 3 samples from the ipsilateral spinal cord of the surgery group (Ipsi group). Additionally, the probe annotation package (rat2302.db) of GSE15041 was downloaded from Bioconductor to convert the probe names to gene names (11).

Data preprocessing

Quality control of the gene chip was performed through weighting and residual plot, relative log expression (RLE) box plot, NUSE box plot, principal components plot (PCA) and RNA degradation curve to remove unqualified samples (12,13). The Robust Multi-array Average integrative algorithm was chosen for preprocessing the microarray data, including background correction, quartile data normalization and probe summarization. Following this, the boxplot of the raw and normalized data was drawn (14). The annotate package in R (15) was subsequently used to convert probe names to gene names. For each sample, a maximum expression value was used for the genes with multiple probe names.

Identification of DEGs

The limma package in R (16) was used to screen out the DEGs (17). Genes with adjusted P value of P<0.05 and |fold change| >2 were selected as DEGs. Then, cluster analysis and heat map construction was performed using the Pheatmap package (https://cran.r-project.org/package=pheatmap) in R to further analyze the DEGs.

Enrichment analysis of DEGs

Gene function annotation is usually performed by Gene ontology (GO) analysis (18). The GO database consists of three categories, including biological process, molecular function and cellular component. The Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways database is a comprehensive and recognized database encompassing numerous biochemical pathways. In order to gain further insight into the function and pathways of the identified DEGs, the Database for Annotation, Visualization and Integrated Discovery (DAVID; version 6.8) tool (19) was used to perform the GO and KEGG pathway enrichment analysis. The threshold was set as P<0.05 and enrichment counts at >2.

PPI network analysis of DEGs

The Search Tool for the Retrieval of Interacting Genes (STRING; version 10.0; http://www.string-db.org/) provides comprehensive and easily accessible interaction information derived from known and predicted protein data (20). Interactions with a combined score of >0.4 were selected and a network was built using the STRING online tool. Networks were visualized using Cytoscape software (21) to identify hub proteins and key genes associated with neuropathic pain. Nodes with higher degrees (connected nodes) were considered hub proteins in the PPI network and may have an essential role in the network.

Animals

Male 8-week-old C57BL/6J (weight, 20–23 g; n=10) purchased from Model Animal Research Center of Nanjing University (Nanjing, China) were used in the present study and were randomly divided into the sham group and surgery group. All animals were housed in a climate-controlled room (temperature: 22±2°C; humidity: 50–70%) with a 12 h light/dark cycle for at least 2 weeks, with ad libitum access to food and water. All animal experiments were approved by the Experimental Animals Welfare and Ethical Inspection of Nanjing Drum Tower Hospital (Nanjing, China).

Neuropathic pain model

The most widely used model of neuropathic pain is the unilateral sciatic nerve chronic constriction injury (CCI), developed by Bennet et al (22). Briefly, animals were anesthetized with an intraperitoneal injection of 80 mg/kg sodium pentobarbital (Sigma-Aldrich; Merck KGaA, Darmstadt, Germany). The skin of the right hind limb was prepared at the mid-thigh level in a lateral position prior to surgery. The sciatic nerve was carefully isolated following skin incision and blunt dissection, two ligatures (5-0 chromic gut suture) were knotted loosely around the sciatic nerve. The distance between ligatures was approximately 1 mm. The surgery was the same in the sham group, however the nerve was only exposed without ligation.

RNA extraction

Mice were deeply anesthetized with sodium pentobarbital and were rapidly decapitated 7 days post-surgery. The lumbar (L4-L5) segments of spinal cords were transversely sectioned and hemi-dissected along the midline. The ipsilateral (CCI-Ipsi) and contralateral (CCI-Contra) lumbar spinal cord of the surgery group and the ipsilateral (Sham-Ipsi) lumbar spinal cord of the sham group was collected. Total RNA was extracted using an RNA extraction kit (Bioteke Corporation, Beijing, China) according to the manufacturer's protocol. RNA concentration was determined using NanoDrop ND-1000 Spectrometer (NanoDrop Technologies; Thermo Fisher Scientific, Inc., Waltham, MA, USA).

RT-qPCR

RNA was reverse transcribed into cDNA using Takara PrimeScript RT master mix (RR036A; Takara Biotechnology Co., Ltd., Dalian China). RT-qPCR was performed using an ABI StepOne Plus Real-Time PCR system (Thermo Fisher Scientific, Inc.) and SYBR Premix Ex Taq II master mix (Takara Biotechnology Co., Ltd.) according to the manufacturer's protocol. The reaction system (10 µl) consisted of: cDNA (1 µl), forward primers (10 µM; 0.2 µl), reverse primers (10 µM; 0.2 µl), ROX reference dye (0.2 µl), RNase-free water (3.4 µl) and SYBR-Green mixture (5 µl). The thermocycling conditions were as follows: Initial denaturation, 95°C for 30 sec, followed by 40 cycles of 95°C for 5 sec and 60°C for 30 sec. GAPDH was used as a housekeeping gene. The relative expression of genes was calculated using the 2−∆∆Ct method (23). The primer sequences used for RT-qPCR were as follows: GAPDH forward, 5′-TGTGTCCGTCGTGGATCTGA-3′ and reverse, 5′-TTGCTGTTGAAGTCGCAGGAG-3′; interleukin (IL)-6 forward, 5′-AGACAAAGCCAGAGTCCTTCA-3′ and reverse, 5′-GTGACTCCAGCTTATCTCTTGGT-3′; transcription factor AP-1 (c-Jun) forward, 5′-CCTTCTACGACGATGCCCTCAA-3′ and reverse, 5′-GGGGTCGGTGTAGTGGTGATGT-3′; urikinase-type plasminogen activator (Plau) forward, 5′-AGTGTGGCCAGAAGGCTCTA-3′ and reverse, 5′-GCTGCTCCACCTCAAACTTC-3′.

Statistical analysis

Data are presented as the mean ± standard deviation. One-way analysis of variance followed by the Bonferroni post hoc test was used to determine the statistical significance, with SPSS software, version 19.0 (IBM Corp., Armonk, NY, USA). P<0.05 was considered to indicate a statistically significant difference.

Results

One disqualified gene chip sample (GSM375686) was removed following quality control of the microarray expression profile. A total of eight samples based on the GPL1355 platform were included for research; this consisted of 2 samples from the Sham group, 3 samples from the Contra group and 3 samples from the Ipsi group. As a result, 13,686 genes were obtained from GSE15041. The box plot of the log expression values for all genes in each sample before and after normalization were drawn (Fig. 1A and B). The median values of each sample were extremely similar, indicating that the data should be further analyzed.

Figure 1.

Data preprocessing and DEG screening. Boxplots of log expressed values (A) before and (B) after normalization. (C) Volcanic plot of the DEGs for the Sham group vs. the Ipsi group (D) and the Contra group vs. Ipsi group. Significantly upregulated genes are indicated in red, downregulated in green, and black indicates no significant change. DEG, differentially expressed gene; Sham, ipsilateral spinal cord sham group; Ipsi, ipsilateral spinal cord surgery group; Contra, contralateral spinal cord surgery group.

Identification of the DEGs

The limma package in R was used to perform the DEG analysis. No DEGs were identified between the Sham and Contra group, 90 DEGs were identified between the Sham and Ipsi group and 191 DEGs were identified between the Contra and Ipsi group. The volcanic plot of the Sham vs. Ipsi groups (Fig. 1C) and the Contra vs. Ipsi groups were drawn (Fig. 1D). In total, 86 common DEGs were selected by Venn diagram (Fig. 2A). A heat map of the identified DEGs of the Sham, Contra and Ipsi groups was subsequently constructed to perform cluster analysis (Fig. 2B). The DEG clustering was predominantly separated into two clusters; one cluster was the Ipsi group and the other was the Sham and Contra group.

Figure 2.

Venn diagram and cluster analysis of identified DEGs. (A) A total of 86 common DEGs were identified from the Venn diagram of the Sham vs. Ipsi (90 DEGs) and Contra vs. Ipsi groups (191 DEGs). (B) Heat map and the clustering pattern of the common DEGs. Red indicates the upregulated genes and green indicates the downregulated genes. Samples were separated into the Ipsi cluster and the Sham and Contra cluster. DEG, differentially expressed gene; Sham, ipsilateral spinal cord sham group; Ipsi, ipsilateral spinal cord surgery group; Contra, contralateral spinal cord surgery group.

Function and pathway enrichment analysis

A total of 86 common DEGs were significantly enriched in 77 GO terms (Table I) and 17 KEGG pathways (Table II). Notably, function enrichment analysis revealed that in cellular component terms, DEGs were significantly enriched in ‘extracellular space’; in biological processes, DEGs were enriched in response to lipopolysaccharides and cytokines and ‘wound healing’; and in molecular function terms, DEGs were enriched in neuropeptide hormone activity. In KEGG pathway enrichment analysis, DEGs were significantly enriched in tuberculosis, rheumatoid arthritis, human T-cell lymphotropic virus type 1 (HTLV-1) and the mitogen-activated protein kinase (MAPK) signaling pathway.

Table I.

Functional enrichment analysis for the common differentially expressed genes.

GO category	GO number	GO term	Number of differentially expressed genes	P-value
CC	GO:0005615	Extracellular space	23	1.77×10‒8
BP	GO:0071222	Cellular response to lipopolysaccharide	8	1.02×10‒5
MF	GO:0005184	Neuropeptide hormone activity	5	1.15×10‒5
BP	GO:0034097	Response to cytokine	6	8.39×10‒5
BP	GO:0043066	Negative regulation of apoptotic process	11	1.06×10‒4
BP	GO:0007565	Female pregnancy	6	1.44×10‒4
BP	GO:0042060	Wound healing	6	2.61×10‒4
BP	GO:0032496	Response to lipopolysaccharide	8	2.67×10‒4
BP	GO:0032355	Response to estradiol	7	2.92×10‒4
BP	GO:0006954	Inflammatory response	8	3.81×10‒4

CC, cellular component; BP, biological process; MF, molecular function; GO, Gene Ontology.

Table II.

KEGG pathway enrichment analysis for the common differentially expressed genes.

KEGG number	KEGG term	Number of differentially expressed genes	P-value
rno05152	Tuberculosis	8	2.79×10‒5
rno05323	Rheumatoid arthritis	5	1.01×10‒3
rno04145	Phagosome	6	2.77×10‒3
rno05166	Human T-cell lymphotropic virus type 1 infection	7	2.99×10‒3
rno05321	Inflammatory bowel disease	4	4.10×10‒3
rno05202	Transcriptional misregulation in cancer	5	8.45×10‒3
rno04010	Mitogen-activated protein kinase signaling pathway	6	8.98×10‒3
rno04514	Cell adhesion molecules	5	1.04×10‒2
rno04612	Antigen processing and presentation	4	1.24×10‒2
rno05146	Amoebiasis	4	1.82×10‒2

KEGG, Kyoto Encyclopedia of Genes and Genomes.

Protein-protein interaction network analysis of DEGs

In the present study, 27 nodes (proteins) and 47 PPI edges (interactions) were obtained following deletion of the disconnected nodes in the network (Fig. 3). IL-6, c-Jun and Plau were considered as hub proteins and key genes (Table III) that may have an important role in the development of neuropathic pain.

Figure 3.

Protein-protein interaction network of the identified differentially expressed genes. The degree of color indicates the degree of fold change (lighter grey indicate a lower fold change and darker grey a higher fold change) and the size of circle represents the node degree. IL-6, c-Jun and Plau were identified as key genes. IL-6, interleukin-6; c-Jun, transcription factor AP-; Plau, urikinase-type plasminogen activator.

Table III.

Hub proteins in the protein-protein interaction network.

Proteins	Degree
IL-6	20
Jun	9
Plau	8
Gadd45α	6
Atf3	5
Cdkn1a	5

IL-6, interleukin 6; c-Jun, transcription factor AP-1; Plau, urikinase-type plasminogen activator; Gadd45α, growth arrest and DNA damage-inducible protein GADD45 α; Atf3, cyclic AMP-dependent transcription factor Atf3; Cdkn1a, cyclin-dependent kinase inhibitor 1.

RT-qPCR validation

RT-qPCR was used to detect the expression of the identified key genes to confirm the reliability of bioinformatics analysis. Compared with the CCI-Ipsi group, the expression levels of IL-6, c-Jun and Plau were significantly decreased (P<0.05) in the Sham-Ipsi and CCI-Contra group 7 days post-surgery; no significant difference (P>0.05) was detected between the Sham-Ipsi group and CCI-Contra group (Fig. 4).

Figure 4.

Expression of key genes in a mouse model of neuropathic pain. Reverse transcription-quantitative polymerase chain reaction was performed to detect the mRNA expression of IL-6, c-Jun and Plau. *P<0.05 vs. Sham-Ipsi group, #P<0.05 vs. CCI-Contra group. IL-6, interleukin 6; c-Jun, transcription factor AP-1; Plau, urikinase-type plasminogen activator; Sham, sham group; Ipsi, ipsilateral lumbar spinal cord; Contra, contralateral lumbar spinal cord; CCI, chronic constriction injury.

Discussion

At present, the mechanisms underlying neuropathic pain remain unclear. The development of chronic pain is thought to be due to peripheral and/or central nerve injury. Several studies have reported that nerve injury induces alteration in the peripheral immune system and glial cells (24–26). Additionally, it has been demonstrated that the activation of microglia is associated with the early stages of neuropathic pain in the spinal dorsal horn where the nerve is injured (27–29). Therefore, identifying key genes involved may be helpful in further understanding these mechanisms underlying the development of neuropathic pain.

Vega-Avelaira et al (10) used the microarray expression profile GSE15041 to identify DEGs associated with neuropathic pain. However, the expression profile required re-analysis, as the gene chip quality was not verified. The present study obtained different results following correct quality control. Microarray quality control has a significant impact on the results of the analysis and unverified microarrays may lead to incorrect results. In the present study, mRNA expression was analyzed in a mouse model of neuropathic pain 7 days post-surgery to confirm the bioinformatics analysis results. RT-qPCR was performed following 7 days as mice have the lowest paw withdrawal threshold and this is accompanied by significant changes in molecular biology (30).

A total of 86 common DEGs among the Sham, Contra and Ipsi groups were identified following quality control of the microarray and removal of the disqualified sample. The pivotal roles of IL-6, Jun and Plau in neuropathic pain were subsequently identified.

IL-6 is an inflammatory cytokine with a wide range of biological effects. It has been reported that unilateral sciatic nerve injury increases the expression of IL-6 mRNA (31), and mechanical allodynia induced by spinal nerve lesion is attenuated and delayed in IL-6 knockout mice (32). Furthermore, intrathecal injection of anti-IL-6 antibody reduces mechanical allodynia induced by spinal nerve transection (33). These data indicate the involvement of IL-6 in the initiation and development of neuropathic pain. Immediate early genes (IEGs) may be transiently and rapidly activated in response to a wide variety of cellular stimuli. Well-characterized IEGs include c-Jun, and the proto-oncogenes c-Fos and c-Myc. Double labeling with the tracers fast blue and horseradish peroxidase-labeled gold underwent retrograde transport from the site of injury and revealed that c-Jun protein expression was confined to injured neurons (34). Activated c-Jun in the dorsal root ganglion (DRG) contributes to the pathogenesis of neuropathic pain induced by CCI (35), and neuropathic pain is alleviated by inhibiting the mitogen-activated protein kinase 8 (JNK)/c-Jun signal pathway (36,37). Furthermore, Plau mRNA is abundantly expressed by many neurons in the peripheral and central nervous system (38). Plau expression is upregulated in rat DRGs following nerve injury (39). The aforementioned data suggests the involvement of IL-6, c-Jun and Plau in neuropathic pain.

Common DEGs were revealed to be enriched in the injury and inflammation-associated biological processes through GO enrichment analysis. Pathway analysis revealed that DEGs were typically enriched in tuberculosis, rheumatoid arthritis, HTLV-1 and MAPK signaling pathways. Spinal tuberculosis may induce varying degrees of damage to the spinal cord nerve roots (40), resulting in neuropathic pain. A clinical study revealed that a third of patients with rheumatoid arthritis have symptoms of neuropathic pain (41), and a review reported that approximately 53.1% of HTLV-1 carriers have symptoms of neuropathic pain (42). The conventional MAPKs include the p38 isoforms, extracellular signal-regulated kinases 1/2 (ERK1/2) and JNKs (43). MAPKs regulate numerous physiological activities, including inflammation, apoptosis, cancer, tumor cell invasion and metastasis. Active MAPK signaling also has an important role in neuropathic pain. It has been reported that IL-6 is involved in the maintenance of neuropathic pain by activating ERK; and the JNK/c-Jun pathway also participates in the maintenance of neuropathic pain (43). It has been demonstrated that p38 MAPK activation in the ipsilateral spinal dorsal horn is significantly increased following CCI (44). Additionally, MAPK activity is significantly decreased following an intrathecal injection of anti-IL-6 antibody (45). These results indicate that IL-6 and c-Jun have an important role in neuropathic pain that is mediated through the MAPK signaling pathway. This is consistent with the bioinformatics findings of the present study and reflect the reliability of the bioinformatics analysis performed.

In the present study, the results of the bioinformatics analysis revealed that there was a significant difference between the Sham and Ipsi group (P<0.05), however no difference was observed between the Sham and Contra group (P>0.05). RT-qPCR analysis was consistent with this. The sciatic nerve was only injured in the CCI-Ipsi group, not the CCI-Contra or Sham-Ipsi group, and the nerve was only exposed in the Sham-Ipsi group. Sciatic nerve injury induces the inflammatory response in the ipsilateral spinal cord, whereas the contralateral spinal cord was not, compared with the sham group, as previously reported (46). Nerve injury induces a set of progressive alterations in ipsilateral spinal cord microglia with increased secretion of cytokines. This does not occur in the contralateral spinal cord following nerve injury, or the uninjured ipsilateral spinal cord. Therefore, significant changes in neurophysiology in spinal cord of the Sham and Contra group may not be induced. DEG analysis revealed that there were numerous DEGs detected in the Ipsi group compared with the Contra and Sham groups, and not between the Contra and Sham groups. This may be due to only the Ipsi group exhibiting nerve injury, whereas the other two groups did not. However, DEGs may be detected between the Contra and Sham group only in the local tissue of the surgical incision, rather than the spinal cord. Thus, many DEGs were identified in the Ipsi group compared with Sham and Contra group by bioinformatics analysis.

In summary, 86 DEGs were identified and IL-6, c-Jun, Plau were identified as key genes in neuropathic pain through pathway and PPI network analysis, which may aid in further understanding of the underlying pain mechanisms. Following the bioinformatics analysis, the expression of key genes was confirmed with RT-qPCR in a mouse model of neuropathic pain. However, the current study has certain limitations. Although the results revealed that IL-6, c-Jun and Plau have an important role in neuropathic pain, further studies are required to confirm the results at the molecular level.