Comprehensive Analysis of Differentially Expressed Genes and Key Pathways in Neuropathic Pain by Spinal Nerve Ligation

Background: Neuropathic pain can cause signicant physical and economic burden to people, and there are no effective long-term treatment methods for this condition. We conducted a bioinformatics analysis of microarray data to identify related mechanisms to determine strategies for more effective treatments of neuropathic pain. Methods: GSE24982 and GSE63442 microarray datasets were extracted from the Gene Expression Omnibus (GEO) database to analyze transcriptome differences of neuropathic pain in the dorsal root ganglions caused by spinal nerve ligation. We ltered the differentially expressed genes (DEGs) in the two datasets and Webgestalt was applied to conduct Gene Ontology (GO) functional analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis of the shared DEGs. String Database and Cytoscape software were used to construct the Protein-Protein Interaction (PPI) network to determine the hub genes, which were subsequently veried in the GSE30691 dataset. Finally, miRDB and miRWalk Databases were used to predict potential miRNA of the selected DEGs. Results: A total of 182 overlapped DEGs were found between GSE24982 and GSE63442 datasets. The GO functional analysis and KEGG enrichment analysis showed that the selected DEGs were mainly enriched in infection, transmembrane transport of ion channels, and synaptic transmission. Combining the results of PPI analysis and the verication of the GSE30691 dataset, we identied seven hub genes related to neuropathic pain (Atf3, Aif1, Ctss, Gfap, Scg2, Jun, and Vgf). Predicted miRNA targeting each selected hub genes were identied. Conclusion: Seven hub genes related to the pathogenesis of neuropathic pain and potential targeting miRNA were identied, expanding understanding of the mechanism of neuropathic pain and facilitating treatment development.


Introduction
The International Association for the Study of Pain (IASP) recently de ned pain as "an unpleasant sensory and emotional experience associated with, or resembling that associated with, --actual or potential tissue damage" [1]. Pain is a kind of proprioception, and acute pain is a necessary sensation for humans to avoid risk [2,3]. However, neuropathic pain is a type of chronic pain with a pathological state. The clinical manifestations of neuropathic pain are spontaneous pain, allodynia, and hyperalgesia, and such condition causes serious physical discomfort and potential economic consequences for 11-19% of the adult population in the United States [4][5][6]. There are several potential causes of neuralgia, including mechanical injuries, metabolic or nutritional neurological changes, and neurotransmitter dysfunction [7][8][9]. However, we still lack effective measures to deal with neuropathic pain. Therefore, it is necessary to clarify the mechanisms of neuropathic pain and give a more deeper theoretical basis for the long-term effective treatments of this painful condition.
Many basic and clinical medical studies have speculated that peripheral sensitization and central sensitization are key factors contributing to neuropathic pain [10,11]. Peripheral noxious stimulation could increase the amount of the excitatory amino acid, glutamate, in the spinal cord and enhance the activity of the ionotropic glutamate receptors, leading to central sensitization [12]. Nerve injuries could decrease levels of inhibitory neurotransmitters such as glycine and γ-aminobutyric acid in the central nervous system (CNS) and promote allodynia and hyperalgesia [13]. Injury can also activate microglia cells and astrocytes in CNS, which can also foster the progression of neuralgia [14]. Peripheral sensitization also can be part of the pathogenesis of neuropathic pain. Several cation ion channels in the peripheral nociceptors, such as Nav, Cav and TRP, can be activated due to peripheral stimulation, which increases the generation of neuronal action potentials and overall neuronal excitability [15]. Additionally, injury may result in increased immune mediators such as TNF-α, CCL2, and PGE2, in the Dorsal Root Ganglions (DRGs) to initiate and maintain neuropathic pain [16]. With multiple potential contributions, a deeper understanding of neuropathic pain requires more extensive research into the underlying mechanisms.
High throughput sequencing technology and bioinformatics analysis have been widely used in basic medical research, and can be a powerful tool to identify target molecules and related mechanisms [17]. We obtained several microarray datasets from the Gene Expression Omnibus (GEO) database related to neuropathic pain, GSE24982, GSE63442, and GSE30691 [18,19]. These datasets are differential expression data of the transcriptome in the DRGs in a spinal nerve ligation (SNL) model of rats. We conducted a comprehensive bioinformatics analysis of these datasets to nd genes and functional pathways related to neuropathic pain in the DRGs. The goal was to gain a further understanding of neuropathic pain and provide a foundation for the investigation of mechanisms and development of treatment methods.

Data Retrieval and Identi cation of DEGs
We identi ed two related expression pro les (GSE24982 and GSE63442) via searching the GEO Datasets (Gene Expression Omnibus database, https://www.ncbi.nlm.nih.gov/) [18]. The GEO Datasets is a subdataset of the National Center of Biotechnology Information (NCBI), which includes biological expression data of many species obtained by array, SNP array, high throughput sequencing, and other methods [20]. GSE24982 and GSE63442 datasets both included rat neuropathic pain model caused by SNL surgery and total extracted RNA in L4, L5 DRGs was used for further investigation by array. We used the tool of GEO2R, an online program that comes with the GEO Datasets, to analyze the differentially expressed genes (DEGs) in the sham group and the SNL group. The GEO2R tool provides statistical P values and FDR between the two groups. Compared with the sham group, we de ned DEGs as genes with a P value < 0.05 and a | log (fold change (FC)) | ≥ 1 in the SNL group. LogFC ≥ 1 or logFC ≤ -1 are separately considered as up-regulated or down-regulated genes. After obtaining the DEGs, we constructed two volcano plots of the two GSE datasets based on the original gene sets. We then constructed a Venn diagram to show the intersection of the DEGs in GSE24982 and GSE63442. Gene Enrichment and Pathway Analysis The common up-regulated and down-regulated DEGs in GSE24982 and GSE63442 were subjected to Gene Ontology (GO) analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis using the online tool of Webgestalt(http://www.webgestalt.org/). The set minimum and maximum number of genes for a category in the Webgestalt were 5 and 2000. We used the BH method for multiple test adjustment. Dplyr and ggplots packages in R software were used to further visualize the GO and KEGG results from Webgestalt.

Protein-Protein Interaction Network Construction and Hub Genes Analysis
An online String Database was used to build Protein-Protein Interaction (PPI) networks of all the selected DEGs (https://string-db.org). String Database can predict the possible association networks of a single protein and the functional network between multiple proteins [21]. The minimum required interaction score was set to medium con dence 0.400 and line thickness was simultaneously used to indicate the strength of data support. We used the software Cytoscape v3.7.2 to further elucidate the results obtained from String Database. In Cytoscape, each gene was considered as a node and experimentally determined interactions were used as edge attributes. The MCODE module was used to calculate the core interactions in PPI network. The MCC method in the module of Cytohubba was used to found the top 10 hub genes. After we determined the hub genes, we compared them with the GSE30691 dataset, another GEO array dataset of rat DRGs in the SNL model [19]. Hub genes that met the standards of P value < 0.05 and | log (fold change) | ≥ 1 in GSE30691 were maintained for subsequent analysis.

MiRNA Prediction of Hub Genes
The miRDB online miRNA prediction database (http://mirdb.org/) was used to predict potential miRNA that could target the selected hub genes. This is a commonly used online database that can be used for analysis of human, mouse, and rat data, and can predict miRNA that bind upstream of genes and also miRNA that bind downstream target genes [22]. The miRNA prediction database miRWalk was also used to predict miRNA that could target upstream regions of hub genes. Predicted miRNAs included in the miRNA prediction data from both databases were subjected to further analysis. The bindingp values of these co-predicted miRNA were determined. Bindingp values represent the degree of connection between mRNA and predicted miRNA, and these values were used in the Cytoscape v3.7.2 software to visualize connections.

Identi cation of DEGs
To nd DEGs associated with neuralgia, GSE24982 and GSE63442 datasets were extracted from the GEO database to analyze the transcriptome differences of neuropathic pain in DRGs caused by SNL. Data from GSE24982 were generated using the GPL1335 platform and data from GSE63442 were generated using the GPL341 platform. DEGs of each dataset were found by using the GEO2R tool with cutoff criteria of P < 0.05 and | log (fold change) | ≥ 1. We found 3319 DEGs among the 31100 genes in the GSE24982 (Fig. 1A) dataset and 338 DEGs among the 15924 genes in the GSE15924 (Fig. 1B). A Venn diagram was constructed to show the DEGs common to the two GSE datasets. A total of 182 common DEGs between GSE24982 and GSE63442 were found, of which 95 DEGs were up-regulated and 87 DEGs were down-regulated (Fig. 1C).

GO and KEGG Analysis
To determine possible related gene functions and enrichment pathways of DEGs related to neuralgia, we conducted Gene Ontology (GO) analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis of the common DEGs between GSE24982 and GSE63442 datasets. The overlapped up-regulated and down-regulated DEGs are presented in Additional le 1: Table S1.
The GO functional analysis demonstrated that the overlapped up-regulated DEGs were enriched in response to wounding, response to lipid, response to organic cyclic compounds, response to oxygencontaining compounds, and response to stimulus (Fig. 2A). The GO functional analysis of the overlapped down-regulated DEGs were mainly enriched in the transmembrane transport of ion channels (Fig. 2B).
The KEGG enrichment analysis revealed biological pathways related to the sorted DEGs. Our results demonstrated that the KEGG pathways of the overlapped up-regulated DEGs were enriched in pertussis, complement and coagulation cascades, infection, and MAPK pathway (Fig. 2C). The KEGG enrichment analysis of the down-regulated DEGs were enriched in cAMP signaling pathway, taste transduction, neuroactive ligand-receptor interaction, glutamatergic synapse, serotonergic synapse, and synaptic vesicle cycle (Fig. 2D).

PPI Analysis and Hub Genes Selection
To nd possible interactions between the common DEGs, we used the String Database to construct the PPI network of the DEGs. Cytoscape software was applied to visualize the results from the String Database. Using MCODE module in Cytoscape, we constructed a PPI network with nodes and edges (Fig. 3). The top 10 hub genes were nally determined by using the module of Cytohubba in Cytoscape based on the PPI scores. The top 10 hub genes were: Atf3, Casp3, Aif1, Csfr, Ctss, Ptpn6, Gfap, Scg2, Jun, and Vgf (Table 1). Scg2 was down-regulated, but the other nine hub genes were up-regulated. A network analysis was also performed to show the connection between hub genes (Fig. 4). Table 1 Top ten hub genes from PPI network and the veri cation of these genes in GSE30691. To further evaluate the hub genes, we used another GSE data, GSE30691 dataset, and veri ed the changes of expression of these genes. The GSE30691 dataset is another microarray dataset of DRG expression data for the rat SNL model. We looked at the expression of the identi ed top ten hub genes in the GSE30691 dataset, using cutoff criteria of P value < 0.05 and | log (fold change) | ≥ 1. The results showed that the absolute logFC values of Casp3, Csf1r, and Ptpn6 were less than 1. However, the other seven genes of Atf3, Aif1, Ctss, Gfap, Scg2, Jun, and Vgf showed differential expression in the GSE30691 dataset and were kept for further analysis (Table 1).

MiRNA Prediction of Hub Genes
The above analysis revealed several hub genes related to neuropathic pain. In order to explore the possible regulation by miRNA binding upstream of the selected hub genes, we used the online miRDB Database for prediction. Table 2 lists the top miRNA predicted by miRDB and the full list of predicted miRNA of hub genes are presented in Additional le2: Table S2. We used another miRNA prediction database, miRWalk Database, to check these miRNA results. The prediction for miRNA targeting Ctss was not found by this second analysis, but the predicted miRNA results targeting all the other selected hub genes were also found in the miRWalk Database. The connection score, bindingp values, between mRNA and predicted miRNAs were obtained from miRWalk Database. Cytoscape software was used to depict the interaction network of hub genes and predicted miRNAs based on bindingp values (Fig. 5) .

Discussion
Neuropathic pain is not a single disease but a clinical compound manifestation induced by various diseases or injuries, resulting in diverse and complex pathogenesis [23]. There have been many treatments recently developed to target the central and peripheral sensitization of neuropathic pain, including non-steroidal anti-in ammatory drugs, ion channel blockers, and minimally invasive interventional therapies [24,25]. However, these treatments cannot signi cantly provide long-term relief of pain, so other contributing mechanisms of neuropathic pain should be explored. Here, DEGs were respectively extracted from the GSE24982 and GSE63442 datasets and we identi ed 182 DEGs shared in the two datasets. GO functional analysis and KEGG enrichment analysis demonstrated that most overlapped DEGs were enriched in infection, transmembrane transport of ion channels, and synaptic transmission. We identi ed the seven most related hub genes in the PPI network (Atf3, Aif1, Ctss, Gfap, Scg2, Jun, and Vgf) and predicted potential upstream miRNA for these genes.
Of the seven hub genes, the Atf3 gene showed the greatest interaction with other hub genes. Atf3 encodes a member of the mammalian activation transcription factor/cAMP responsive element-binding (CREB) protein family of transcription factors [26]. Previous work showed that neuropathic pain caused by medial plantar nerve ligation can increase the expression of ATF3 in DRGs [27], suggesting that ATF3 is a signi cant contributor to neuropathic pain. Aif1, allograft in ammatory factor 1, acts in macrophage activation and lymphocyte migration [28], so the effect of Aif1 on neuropathic pain may be due to its involvement in pain-related in ammation. Ctss, Cathepsin S, is the key protease responsible for the removal of the invariant chain from MHC class II molecules [29]. CTSS released from glial cells was shown to relieve neuropathic pain by reversing its expression trend after the partial ligation of the sciatic nerve model [30].
Gfap, glial brillary acidic protein, is a speci c marker of astrocytes [31]. Our study found increased Gfap expression related to neuropathic pain. Tian's research also showed that astrocytes contributed to complex regional pain via the MMP2-JNK1/2 pathway [32]. Scg2, Secretogranin II, was the only downregulated gene of the selected seven hub genes. Scg2 encodes neuroendocrine secretory proteins and packages these neuropeptides into secretory vesicles. Scg2 can also regulate the release of inhibitory neurotransmitters such as r-aminobutyric acid and maintain a high level of neuronal excitability to support the pain state [33]. Jun, AP-1 transcription factor subunit, is an oncogene that can activate the cAMP pathway. Sanna's study showed that HDAC1/C-Jun complex promoted the progression of neuropathic pain via JNK pathway [34]. Vgf, a neuro-endocrine speci c protein, is expressed in neuroendocrine cells and plays many roles in the neuroplasticity of memory, depression, and chronic pain. Previous work showed that central and peripheral Vgf play key role in the chronicity of pain [35].
In addition to the single analysis of the selected hub genes, we also conducted GO functional analysis and KEGG enrichment analysis of all the common DEGs. The DEGs were enriched for infection, transmembrane transport of ion channels, and synaptic transmission. The tissue damage signal carried by the neuropathic pain model may be transmitted from the peripheral nociceptors to the DEGs. DEGs not identi ed as hub genes may still play important roles in the pathogenesis of neuropathic pain. For example, Kcnip3, a DEG not identi ed as a hub gene, is a voltage-gated potassium channel gene. Its expression level decreased in the SNL group, which would reduce the potassium current causing neurons to remain in a high potential state for a longer time, thus promoting neuropathic pain [36].
MicroRNA (miRNA) is a type of non-coding RNA that can bind to the 3'UTR region of mRNA to regulate its translation [37]. In this study, two miRNA prediction databases were used to predict the miRNA of the selected hub genes and we used Cytoscape software to predict the network relationships between the hub genes and the predicted miRNA. The miRNA predicted based on the hub genes might contribute to the regulation of neuropathic pain and this potential regulation should be further investigated.

Conclusions
In conclusion, we identi ed 182 genes exhibiting differential expression in the DRGs in the SNL model. We performed GO functional analysis and KEGG enrichment analysis of these DEGs, and the results indicate potential mechanisms related to neurological pain. We identi ed seven hub genes related to neurological pain and predicted their respective miRNAs. Our research provides a foundation for further investigation of the mechanisms and potential treatment strategies for neuropathic pain. Availability of data and materials The datasets generated and/or analyzed during the current study are available from the corresponding author on reasonable request.  The PPI network of the common DEGs contructed using the Cytoscape software. Spheres represent nodes and lines represent edges.