Screening and Identification of Key Genes, Pathways, and Drugs Associated with Neuropathic Pain in Dorsal Horn: Evidence from Bioinformatic Analysis

Introduction

Neuropathic pain represents a complex chronic condition occurring post-nervous damage and affecting the somatosensory nervous system.¹ The pain is usually described as a kind of persistent, unbearable feeling and seriously compromises the patient’s quality of life.^2,3 The prevalence of neuropathic pain in the general population may report up to 10%, and around 20% to 25% of them will evolve into lifetime chronic pain, which becomes a tremendous social‐economic burden.^4–6 Unfortunately, current therapies only work for less than 40% of neuropathic pain patients, and treatment outcomes vary from person to person with high instability.^7–,9 Therefore, identifications of key genes and potential effective therapeutics toward neuropathic pain are urgently needed.

In recent years, plenty of studies focused on the potential mechanism of neuropathic pain,^10–12 but most of them only restricted to dorsal root ganglia (DRG) without the involvement of spinal cord like dorsal horns.¹³ Meanwhile, it was verified that hundreds of pain-related genes existed substantial expression changes in the spinal cord tissue, and most of these genes were tested in neuronal, instead of microglial. Interestingly, neuropathic pain was considered as immune reactions in spinal cord regulated by resident immune cells such as microglia or astrocytes.¹⁴ The activation of microglia in spinal dorsal horns contributed to the central nerve’s sensitization, which was thought to be the key mechanism for neuropathic pain development.^15–17 Although progresses like the caspase-6 signal pathway, activation of the NMDA receptor or ERK phosphorylation were revealed in neuropathic pain mediated by microglia,^15,18,19 the precise and comprehensive molecular mechanisms have not yet been fully understood.

In recent decades, with the development of high-throughput technology and bioinformatic analysis database, millions of gene expression profiles for neuropathic pain have been gradually performed. However, few comprehensive comparative analyses have been carried out on gene expression profiling of uninjured spinal dorsal horns, and potential interactions among DEGs and pathways in microglia have not been taken into consideration. Notably, Jeong et al originally submitted the GSE60670 dataset to the Gene Expression Omnibus (GEO) database.¹⁴ In this study, we will reanalyze it using a bioinformatic method to identify the DEGs, finding the crucial genes related to biological processes (BPs) and molecular functions (MFs) at the microglia level and exploring potential pathways in the uninjured spinal dorsal horn.

Materials and Methods

Microarray Data

The microarray data of GSE60670 were downloaded from the GEO database (https://www.ncbi.nlm.nih.gov/geo/), and its microarray platform was GPL6246 (Affymetrix Mouse Gene 1.0 ST Array). GSE60670 was a group of microglia-enriched microarrays from male C57 mice between nerve-injury and sham surgery. The nerve-injured model was generated by L4 spinal nerve transection (SNT), and the experiment contained 23 samples, including 12 samples in the SNT group and 11 sham samples.

Identification of DEGs

Firstly, the downloaded raw matrix file was imported into R Studio software (v 4.3). All data were separated into two parts, namely, the SNT group and the sham group. Bayes test in the limma package was used to identify DEGs between SNT and Sham groups.³ Genes with conditions of |LogFC| >1.5 and P-values <0.05 were considered to be DEGs, which were visualized by generating heat map and volcano plot via R heat map and ggplot2 packages. DEGs with condition of |LogFC| >1.5 were considered as upregulated genes in the SNT group compared with the Sham group. After that, all gene probes were transformed into the corresponding gene symbols using database mogene10sttranscriptcluster.db and package org.Mm.eg.DBDB in R Studio as what was done previously.

Gene Ontology and Pathway Enrichment Analyses

Significant DEGs were combined with the Gene Ontology (GO) database or The Kyoto Encyclopedia of Genes and Genomes (KEGG) database for further comprehensive analysis. GO provided a specific way to illustrate the function of the selected genes in three main domains, namely, cellular component (CC), molecular function (MF), and biological process (BP). KEGG database contained abundant pathway information, which could help detect relevant potential signal pathways of DEGs. Both GO analysis and KEGG analysis were performed in R software with clusterProfiler package and visualized by ggplot2 package in R Studio.

Protein to Protein Interaction Network and Module Analysis

Protein to protein interaction (PPI) network analysis was aimed to evaluate relationships among DEGs from protein aspect, and the network was constructed using the Search Tool for the Retrieval of Interacting Genes (STRING; http://STRING‐db.org), an online database. The STRING results were visualized by Cytoscape software, and the MCODE plugin was utilized to screen the most significant modules in the whole PPI network. The criteria for selection were as follows: MCODE scores >5, node score cut‐off = 0.2, degree cut‐off = 2, Max depth = 100, and k‐score = 2.⁶ At last, the genes with counts over 8 were chosen as the hub gene,³ which were screened using cytoHubba, another plugin in Cytoscape.

Transcription Factors Matching Prediction

The DEGs were uploaded to NetworkAnalyst (HTTPS:// www.networkanalyst.ca/faces /home.XHTM) to obtain the transcription factors (TFs) paired with DEGs. These TFs were predicted from the ENCODE database and the generated list of Transcription factors-Differential expression genes (TFs‐DEGs) pairs outcomes. The results were visualized in Cytoscape.

Related Drugs Prediction

The hub genes analyzed from the PPI network served as candidate targets to find therapeutic drugs in the DGIdb database (http://dgidb.genome.wustl.edu/). The discovery of gene–drug interactions was also visualized in Cytoscape, and the filter of DGIdb was set as “FDA approved.”

Expression Confirmation of Pain-Related DEGs in Neuropathic Pain Model

Rodent Pain Model Establishment

Health C57 female mice (weighing around 0.018–0.021 kg) were provided by the Nanfang Hospital, Southern Medical University. All mice were housed at separate cages in a constant temperature room, with ad libitum feeding and 12 hours diurnal rhythms. The study was approved by the Animal Ethics Committee of Southern Medical University. The ethics number related to this study is as follows: NFYY-2020-98 Model establishment referred to IASP (International Association for Study of Pain) guideline, every possible effort was made to minimize unnecessary suffering for animals.

All mice were randomly assigned into the SNT group and the sham group (15 mice in each group). The SNT model was constructed according to the study by DeLeo.²⁰ Specifically, following the pentobarbital anesthesia, a midline incision of 1 cm longitudinal covering the L4-S2 level was made to expose L4 and L5 spinal nerves. The left paraspinal muscles and L4 transverse processes were separated carefully. The L4 spinal nerve was transected without damaging the L3 spinal nerve. The wound was then closed by suturing the muscle and skin layers. For the sham group, the same procedures were performed except that the L4 spinal nerve was only isolated but not transected.¹⁴

Behavior Test

To ascertain the effect of the neuropathic pain model, 50% mechanical withdrawal palm thresholds (50% MWTs) were evaluated based on the up-down method.⁸ For all mouse, 50% MWTs were tested before the operation as the basal line (day 0), and the same processes were undertaken respectively on the 1st, 7th and 14th days after operation. To reduce the impact of environment, mouse were placed on the test platform at least 30 minutes before the behavioral test, and all tests were carried out between 9:00 and 13:00 under silent atmosphere. Same von-Frey monofilaments were used to stab the third or fourth digit of mouse paw and observe whether palms withdraw or not. The final 50% MWTs were transformed into pain thresholds by an experimenter who was blinded to the test procedure.

Real-Time Polymerase Chain Reaction (RT-PCR)

Whole L4-L6 lumbar enlargement of spinal cord were collected in 14 days after SNT or sham operation was immediately frozen on dry ice and stored at −80°C.³ Total RNA from snap-frozen spinal tissues was isolated using Trizol reagent (Thermo Fisher Scientific, Waltham, MA, USA). The cDNA reverse transcription was performed by reverse transcript enzyme (Toyobo, Shanghai, China) and quantitative analysis was conducted with quantitative RT-PCR (qRT-PCR) using SYBR Green Realtime PCR (Toyobo, Shanghai, China). 18s work as the housekeeping gene. Data of RT-PCR were calculated using the 2-ΔΔCt method,²¹ and the forward and reverse primer sequences were exhibited in Table S1.

Statistical Analysis

All data were analyzed using R software and GraphPad Prism 8. Data regarding the MWTs and RT-PCR were analyzed statistically using the Student Test. P < 0.05 was considered statistically significant.

Results

Data Normalization

A total of 34,831 genes were detected in 23 samples. The original microarray has undergone Robust Multichip Average (RMA) normalization in the production, and authors removed batch effects using a two-stage least-squares method for estimation of location and scale parameters due to nonbiological variability. Researchers performed a series of quality controls to ensure that the chip quality met the study criteria. In our study, data normalization was assessed through biological variability and cross-comparability. Principal component analysis (PCA) was applied to verify biological variability between different samples. SNT samples and Sham samples were separated into two groups with the density of plots and fitted circles indicating globally distinct expression profiles (Figure 1A). A box plot (Figure 1B) illustrated the gene expression of each mouse microglia sample, where the black lines were almost on the same level, indicating the standardization and accuracy of subsequent data processing.

Figure 1 The distribution of samples after normalization. (A) Principal component analysis (PCA) of each sample. The result presented that SNT and nontreatment sham samples were grouped separately, indicating globally distinct expression profiles. (B) Box plot of each sample, the horizontal axis represents the sample name, while the vertical axis represents the expression value after normalization. The black line in the box represents the median of value, which can stand for the degree of normalization. They were all in the same line, indicating the normalization was effective.

DEGs Identification

According to the criteria that P < 0.05 and |logFC| >1.5, 100 DEGs were detected, including 78 downregulated and 12 upregulated DEGs. For all DEGs, volcano plot and heat map were generated using ggplot2 package (Figure 2A and B). The top 20 ranked genes are summarized in Table 1.

Table 1 Results of the Top 20 Ranked DEGs

Figure 2 The volcano plot (A) and heat map of DEGs. (A) The volcano plot of differentially expressed genes(DEGs), the horizontal axis represents the fold change (SNT versus Sham), and the vertical axis represents the p-value of DEGs, genes with absolute value of logFC >1.5 and p-value <0.05 were confirmed as DEGs. (B) The heat map of DEGs, the horizontal axis represents the name of each sample, while the left vertical axis represents the degree of gene clustering. Red represented the upregulated genes, while blue stands for the downregulated genes.

KEGG Pathway Analysis and GO Annotation

A total of 22 KEGG-enriched pathways were detected for all DEGs, which revealed that DEGs were significantly enriched in the signaling pathways shown in Figure 3A, such as the malaria signal and toll-like receptor signal. Top KEGG-enriched pathways in 28 upregulated and 72 downregulated DEGs were analyzed separately, and results are presented in Figure 3B.

Figure 3 GO and KEGG enrichment analysis of the DEGs. (A) the comprehensive analysis of DEGs in the KEGG database. (B) the separative KEGG analysis for upregulated DEGs and downregulated DEGs, grey bars stand for upregulated and black bars represented downregulated genes. (C–E) represented the biological processes (BPs), cellular components (CCs) and molecular functions (MFs) analysis for whole DEGs.

GO enrichment analysis of these 100 DEGs was performed to identify relevant biological processes (BPs), molecular functions (MFs), and cellular components (CCs), and the result found 86 GO terms with P-value <0.05 (69 BP items, 10 CC items, and 7 MF items). Specifically, the DEGs were mostly enriched in the regulation of inflammatory response for BP terms, the secretory granule for CC terms, and the cytokine activity for MF terms (Figure 3C–E). Dominated relationships between DEGs and KEGG pathways or GO annotation were exhibited in Table 2.

Table 2 Significant Enrichment Results of KEGG and GO Analysis

PPI Network Analysis and Hub-Gene Extraction

Base on the database of STRING, the PPI network was established to determine the relationship between proteins expressed by DEGs. The current PPI network contained 60 nodes and 143 protein-pairs (Figure 4A). The plugin of cytoHubba was employed to define hub genes, and 9 hub genes, including Il1b, Tnf-a, Myd88, Lyz2, Clec7a, Cxcl10, Tlr4, Tlr7 and Tlr13, were defined as counts over 8 in 12 mathematical calculation methods. Based on the module analysis with Cytoscape MCODE, three modules were acquired (Figure 4B–D). The KEGG pathway enrichment and GO annotation of DEGs in these three clusters were presented in Figure S1A–D.

Figure 4 PPI network of screened DEGs from string database. (A) The comprehensive analysis of whole DEGs. (B–D) The important modules after filtered by MCODE. Brown nodes stand for upregulated genes, while grey nodes stand for downregulated. The size of each nodes is determined by p-value. The width of each edge is determined by the combined score between the related two nodes.

Establishment of TF-DEGs Network and DEGs-Drugs Network

To discover more potential functions of DEGs, TF–DEGs and DEGs–Drugs relationships were further analyzed. According to the ENCODE database, 277 DEG–TF pairs were found among 54 DEGs and 32 TF. As shown in Table 3, MYB and ZMZ1 regulated 16 DEGs separately. Cytoscape was utilized to visualize TF-DEGs and the vital modules were identified, which were displayed in Figure S2A and B. To limit the scope of drugs prediction, only nine hub genes were considered. And the finding indicated that more than 100 candidate drugs may have effect on hub genes, and 22 drugs may play essential roles in neuropathic pain alleviation (Figure 5, Table S2) after eliminating drugs without confirmed effective ways (eg, antagonist or inhibitor).

Table 3 Result of TF-DEGs Connection

Figure 5 The drug–hub gene interaction network. The grey circle nodes represent the hub genes, while the blue nodes represent the relevant drugs. Grey diamond nodes highlight drug may take effect on least two hub genes.

Behavior Test and DEGs Expression Confirmation for SNT Model

Following the transection of L4 spine, the nociceptive response to mechanical stimulation increased significantly. Figure 6 illustrated that the 50%MWTs for two groups were similar from the statistical perspective before the operation on day 0, but decreased significantly in SNT groups compared with Sham groups after the operation. And the difference between the two groups was statistically significant on day 7 and day 14.

Figure 6 Time-course curves of 50% MWT induced by SNT. 50% MWT was measured before (day 0) and 7, 14 days after surgery. **P<0.01, ***P< 0.001, SNT versus sham (at each time point, n=6 per group). Abbreviations: MWT, mechanical withdrawal threshold; SNT, spinal nerve transection.

Given the fact that Ctdnep1, Rom1, Susd2, Gtpbp8 and Fggy were the top five most significantly upregulated DEGs, and Tmem176a, Cdc86, Scn9a, Ctsc and Sucl2 were top five most significantly downregulated DEGs, the mRNA expressions of these 10 genes were then detected in both SNT groups and Sham groups within 14 days. Results confirmed that SNT led to remarkably enhanced expressions of Ctdnep1, Rom1, Susd2 and Fggy, but Tmem176a, Cdc86, Scn9a, Ctsc and Sucl2 had a significant reduction in SNT groups in comparison with sham groups (Figure 7). However, only little increase was observed in Gtpbp8 expression in the SNT group compared with the sham group, and the difference was not statistically significant (Figure 7).

Figure 7 mRNA expression of top five upregulated and downregulated DEGs in the L4 lumbar enlargement of spinal cord. *P<0.05, **P<0.01, SNT versus sham (at each time point, at least n=6 per group). Abbreviations: MWT, mechanical withdrawal threshold; SNT, spinal nerve transection.

Discussion

Complex processes have been demonstrated in the mechanism of neuropathic pain development, but the exact mechanism has not been completely elucidated. Although some gene expression changes in dorsal horns have been found in the neuropathic pain model, and accumulating evidence indicates that microglia, as a kind of immune cell in the nervous system, contributes a lot to the occurrence or development of neuropathic pain, no genome-level analysis has been carried out in this aspect to date.^13,22 For this reason, our study is the first bioinformatic analysis that focuses on the gene variation in dorsal horns’ microglia.

After the integrated interpretation of differential gene expression, the top 20 DEGs were deemed as the most valuable candidate genes that may affect neuropathic pain development. Notably, several reports in recent years indicated that some of them played a vital role in mediating neuropathic pain. For example, accumulating evidence reported that Tmem176a was involved in the remodeling of the dendritic cell after spinal cord injury.²³ Dendritic cells play a central role not only in the development of innate and adaptive immunity but also in the induction and maintenance of immune tolerance.²⁴ Previous studies showed that Tmem176a and Tmem176b could inhibit dendritic cell maturation and activation, together with abnormal functioning of cell‐mediated immunity.^23,24 As such, the decreased expression of Tmem176a downregulated the inhibition of dendritic cell and accelerated the remodeling process of dendritic cells.

Further studies indicated that in neuropathic models, disruption of Rac1, a kind of GTP-binding protein and known as an antagonism to regulate actin filament reorganization, could reduce dendritic remodeling and attenuate neuropathic pain following peripheral nerve injury.^25–27 We speculated that it might explain why Ctpbp8, another GTP-binding protein, has higher expression. Additionally, it was also known that Scn9a positively encoded the expression of voltage-dependent sodium 1.7 channels (NaV 1.7 channels), and the increased expression of NaV1.7 channels contributed to the exacerbation of pain behavior.^28,29 Besides, both Suclg2 and Fggy were demonstrated to have a tight relationship with energy metabolisms like the tricarboxylic acid cycle or carbohydrate metabolic process.^30,31 The question whether there is a potential link between energy metabolism and neuropathic pain remains to be explored.

KEGG pathway analysis and GO annotation were then utilized to obtain a comprehensive understanding of these DEGs. GO terms of DEGs for biological process analysis indicated more than 10 DEGs enriched in inflammation response or inflammatory factors production like IL1 and IL6. IL-1β is an essential mediator in microglia after activating the p38 pathway, which can further regulate pre-synaptic and post-synaptic function by affecting the N-methyl-D-aspartic acid receptor (NMDAR).^10,29 Furthermore, the KEGG pathway analysis result revealed that DEGs were mainly enriched in malaria inflammatory-related pathway and Toll-like receptor signal pathway, which were consistent with the BP term in GO annotation. Interestingly, previous studies found that some clinically used drugs in the treatment of malaria also alleviated neuropathic pain, like tetracyclines and mefloquine,^33,34 which suggested that there may exist overlapped mechanisms between malaria and neuropathic pain. Some studies verified that drugs for malaria might inhibit the microglial activation, attenuate the expression of IL-1β, PGE2 or BDNF, suppress the activation of the p38 pathway, and make an effect on the synaptic plasticity mediated by microglia.^34,35 All these processes have the potential to alleviate neuropathic pain.

PPI network and key modules identified several hub genes, including inflammatory factors like IL-1β and TNF, Toll-like receptors like TLR7, TLR4, TLR13, and other genes such as Myd88 and Cxcl10. Similar to IL-1β introduced above, TNF was also a vital mediator produced by microglia after the activation of p38/ERK/NF-kappaB pathways.^13,32,36,37 The toll-like receptors (TLRs) has played an emerging role in neuropathic pain. Microglial activation via p38 phosphorylation was common after the activation of cell-surface receptors on microglia like TLR4.^13,38,39 By the way, TLR4 was a specific pattern recognition receptor on microglia and induced proinflammatory mediator production via TLR4/AP-1 pathway.^3,40 Besides TLR4, TLR7 also contributed to neuropathic pain by activating the NF-kappaB pathway in spinal level.^41,42 Additionally, a previous study also explained that depressed expression of Myd88 could alleviate neuropathic pain.⁴³

The DGIdb database makes it possible to screen the potential drugs that might relieve neuropathic pain. Some malaria drugs such as chloroquine and hydroxychloroquine were also included in the list of potential therapeutics, which was consistent with the KEGG pathway enrichment analysis Former studies confirmed that Adalimumab and Infliximab, the inhibitors of TNF or IL-1β as monoclonal antibody, had been clinically used for neuropathic pain therapy.^3,44,45 However, other potential antibodies in our analysis, like Golimumab and Certolizumab, were still under research. Besides, accumulated evidence showed that tissue plasminogen activator (t-PA) had effects on neuropathic pain, and the underlying mechanism was considered related to Annexin A2 (ANX2), a calcium-binding protein.^46–48 Whether alteplase, as a kind of t-PA, can alleviate neuropathic pain in clinical was unclear right now. Interestingly, several studies and case reports have demonstrated that atorvastatin, a common drug for lipid disorders, also has an effect on neuropathic pain, and the NF-kappaB pathway may contribute to this effect.^49,50 The question whether another potential drug, pravastatin, has similar effect is also worth discovering.

Several limitations should be claimed in our study. Firstly, the sample size is relatively small, and more fitted microarrays are waited to add in. Meanwhile, the original microarray was submitted in 2016, due to objective reasons, raw data of the microarray is not available at the present, it’s not available for us to conduct more detailed assessment of the chip’s quality, more updated microarrays are waited to refresh it. Secondly, most results of our study are only based on the bioinformatic analysis, and even if experiments have been conducted to verify top DEGs, more confirmative experiments for hub-gene, viral pathway, and potential drugs are needed. Thirdly, due to the sample limitation, we merged the samples of 1 day post-operation and that of 7 days, so our study cannot separate genes work for acute neuropathic pain or chronic pain specifically. At last, the original research employed near single-cell microarray of microglia from spinal cord slices, but our study used L4 lumbar enlargement of spinal cord for DEGs expression confirmation, explaining why we cannot get statistically increased expression in gene Gtpbp 8. The expression variation may be only limited in microglia instead of the whole spinal cord, Considering the existing conclusion cannot completely rule out the interference of neurons and other glial cells (such as astrocytes). However, from the perspective of activation timespan, the sequence of activity time of different central nerve cells in the development of neuropathic pain is relatively distinct, so we supposed that the mixing of other types of cells would not cause the external validation to lose its full credibility, but more precise empirical experiments on microglia level are wanted in future researches.

Conclusion

Multiple studies have suggested that activated microglia is a critical contributor to the initiation and development of neuropathic pain. Our current study has focused on the microglia in the dorsal horn, and demonstrated that Tmem176a and Scn9a may play essential roles in neuropathic pain development. More potential mechanisms involved in toll-like receptors and inflammatory reaction are worth to be studied in the future, and several drugs like alteplase and pravastatin may take a part in alleviating neuropathic pain. Thus, further experiments are required to probe the function of these genes and drugs.