# Analysis of Crucial Genes and Pathways Associated with Spared Nerve Injury-Induced Neuropathic Pain

Authors: Dong Mao; Huang Zhai; Gang Zhao; Jingyi Mi; Yongjun Rui
Journal: Neural Plasticity
 (2020)

---

## Abstract

Purpose. The study was aimed at elucidating the molecular mechanism underlying neuropathic pain induced by spared nerve injury (SNI). Methods. The microarray data of GSE30691 were downloaded from the Gene Expression Omnibus database, including sciatic nerve lesion samples at 3, 7, 21, and 40 days after SNI and sham control samples at 3, 7, and 21 days. Differential analysis along with Mfuzz clustering analysis was performed to screen crucial clusters and cluster genes. Subsequently, comprehensive bioinformatic analyses were performed, including functional enrichment analysis, protein-protein interaction (PPI) network and module analysis, and transcription factor- (TF-) gene and miRNA-target interaction predictions. Moreover, the screened differentially expressed genes (DEGs) were corroborated using two other microarray datasets. Results. Three clusters with different change trends over time after SNI were obtained. Protein kinase CAMP-activated catalytic subunit beta (Prkacb), complement C3 (C3), and activating transcription factor 3 (Atf3) were hub nodes in the PPI network, and fibroblast growth factor 9 (Fgf9) was found to interact with more TFs. Prkacb and Fgf9 were significantly enriched in the MAPK signaling pathway. Moreover, rno-miR-3583-5p was targeted by Fgf9, and rno-miR-1912-3p was targeted by neuregulin 1 (Nrg1). Key genes like Nrg1 and Fgf9 in cluster 1, Timp1 in cluster 2, and Atf3 and C3 in cluster 3 were screened out after corroborating microarray data with other microarray data. Conclusions. Key pathways like the MAPK signaling pathway and crucial genes like Prkacb, Nrg1, Fgf9, Timp1, C3, and Atf3 may contribute to SNI-induced neuropathic pain development in rats.


---

## Body

## 1. Introduction

Neuropathic pain refers to chronic pain originating from neurological pathology, and it affects approximately 7–10% of the global population [1, 2]. It is characterized by spontaneous hyperalgesia, dysesthesia, and allodynia [3, 4]. Neuropathic pain can negatively affect quality of life, and most neuropathic pain patients may suffer from negative moods like depression and anxiety disorders [3, 5, 6]. Despite great progress in understanding the pathogenesis of neuropathic pain and patient prognosis, little is known about the genetic basis and mechanism underlying this disease, and many patients respond poorly to current therapies. A better understanding of the molecular mechanism of neuropathic pain will be important for further effective therapies.The pathogenesis of neuropathic pain is complex, and elucidation of specific molecular alterations helps understand the mechanisms involved in neuropathic pain development. Several biological alterations like ion channel or inflammatory mediator expression, extracellular proteins, and epigenetic influences have been implicated in neuropathic pain [7]. In recent years, microarray data have been widely used to globally assess gene expression signatures that provide new insights into disease pathophysiology [8]. Spared nerve injury (SNI) is a robust neuropathic pain model [9]. Since the development of bioinformatics, massive microarray data have been used to extensively investigate the candidate molecules associated with SNI-induced neuropathic pain and help identify potential targets for disease diagnosis and treatment. For instance, crucial genes like C-X-C motif chemokine receptor 2 protein coding (CXCR2) and G protein-coupled receptor kinase 1 (GRK1) and miRNAs like miR-208a-5p and miR-135a-2-3p are differentially expressed in SNI based on microarray data, suggesting a possible role in neuropathic pain [10]. In addition to genes and miRNAs, several pathways like focal adhesion are also revealed to be involved in neuropathic pain development [11]. However, the possible mechanism underlying neuropathic pain remains incompletely understood, and reliable biomarkers for diagnosis and treatment are lacking.In previous studies, the microarray data GSE30691 developed by Costigan et al. [12] has been used to identify novel pain-related genes by mining expression profiling data in three rodent neuropathic pain models using an iteratively reweighted least squares outlier-resistant regression method and weighted gene coexpression network analysis, followed by analysis of the associations between polymorphisms in the gene and pain phenotypes in human cohorts by a combination of bioinformatic analysis of transcriptional changes in rodent models and human gene polymorphism association studies. As a result, a neuropathic potassium channel modulatory subunit (also called Kv9.1) was downregulated in all three neuropathic pain models, and a common amino acid-altering KCNS1 polymorphism is associated with the pain phenotype in five of six independent cohorts [12]. Moreover, the microarray data GSE30691 are used to screen candidate genes associated with neuropathic pain using functional and weighted coexpression modular analysis [13] or differential analysis with a random walk with restart [14]. In contrast to these previous studies, we also downloaded the microarray data GSE30691 from the NCBI Gene Expression Omnibus (GEO) [15] and reanalyzed by other bioinformatic methods, aiming to find more neuropathic pain-related genes and pathways. In detail, differential analysis along with Mfuzz clustering analysis was utilized to screen crucial clusters and cluster genes. Subsequently, comprehensive bioinformatic analyses were performed, including functional enrichment analysis, protein-protein interaction (PPI) network and module analysis, and transcription factor- (TF-) gene and miRNA-target interaction predictions. Moreover, the screened DEGs were corroborated using two other microarray datasets. Our findings will help to elucidate key molecular mechanisms associated with neuropathic pain and discover new potential targets for disease therapies.


## 2. Materials and Methods

### 2.1. Data Sources

The microarray data GSE30691 deposited in the NCBI GEO (http://www.ncbi.nlm.nih.gov/geo/) database by Costigan et al. [12], which were generated on the GPL85 [RG_U34A] Affymetrix Rat Genome U34 Array platform, was downloaded. This dataset contained adult rat L4 and L5 dorsal root ganglion (DRG) ipsilateral samples from different sciatic nerve lesions, including SNI, spinal nerve ligation, and chronic constriction injury. Expression profiling of sciatic nerve lesion samples at 3, 7, 21, and 40 days after SNI and 3, 7, and 21 days of sham control samples (n=3 per time point; total=24) were extracted. The data were downloaded in May 2020.


### 2.2. Differential Expression Analysis

The differential expression of probes between SNI samples and sham controls was analyzed at time points of 3, 7, and 21 days post-SNI using the online tool GEO2R (http://www.ncbi.nlm.nih.gov/geo/geo2r/). The probes that did not map to gene symbols, or probes mapping to different genes were removed. The threshold value was set to p value < 0.05 and ∣logfoldchange∣>0.585. The number of differentially expressed probes and genes was counted. The unions of DEGs that had the same expression change tendency at three time points were collected for further analysis.


### 2.3. Mfuzz Clustering Analysis

To study the expression pattern of DEGs with the time-course injury, the expression values of DEGs in 0, 7, 21, and 40 days of SNI samples were extracted before using clustering analysis with Mfuzz (version 2.42.0,http://bioconductor.org/packages/release/bioc/html/Mfuzz.html) [16] in the R environment. In the clustering analysis, the optimal cluster number was calculated using default parameters. Meanwhile, the min score (membership) threshold was set to 0.6. For a given cluster, if the membership of two genes was higher, their expression trend was more similar, representing the importance of the gene in the cluster to a certain extent. Clusters with different change trends of genes over time were used for subsequent analysis. Cluster genes in each cluster were obtained using the online tool Metascape (http://metascape.org) [17].


### 2.4. Functional Enrichment Analysis

To investigate the function of cluster genes in each cluster, Gene Ontology (GO) [18] biological process (BP) and Kyoto Encyclopedia of Genes and Genomes (KEGG) [19] pathway enrichment analyses were performed, which are widely used for the functional annotation of large-scale data. The species was set as Rattus norvegicus, and the parameters were default, including min overlap: 3, p value cutoff: 0.01, and min enrichment: 1.5. After obtaining the terms conforming to the above parameters, further term clustering was performed based on the similarity of genes enriched in each term>0.3, and the term with the greatest statistical significance (p value minimum) was selected to represent the term cluster. The top 20 term clusters with the most significance were displayed. In addition, the interaction network of terms was constructed to further capture the relationship between terms in term clusters. The selection criteria for terms in this network were as follows: the terms with the best p value in each cluster of the top 20 clusters were selected, there were ≤15 terms in each cluster, ≤250 terms in total, and similarity of >0.3. After obtaining the interaction relations of terms, the interactive network of terms was established and visualized using Cytoscape (version 3.4.0, http://chianti.ucsd.edu/cytoscape-3.4.0/) [20].


### 2.5. PPI Network Construction and Module Analysis

The Search Tool for the Retrieval of Interacting Genes (STRING) database [21] was used to analyze interactions between the gene coding proteins. The PPI pairs between cluster genes were analyzed using a medium confidence of 0.4. The species was rat. The PPI network was constructed using Cytoscape software. The degree centrality of nodes was excavated using the CytoNCA plugin (version 2.1.6, http://apps.cytoscape.org/apps/cytonca) [22] without weighting, and the hub nodes were identified by ranking the node degree. In addition, the modules of the PPI network with nodes>3 were screened out using the Molecular Complex Detection (MCODE) (version 1.5.1, http://apps.cytoscape.org/apps/mcode) [23] plugin. The parameters were set to default, including DegreeCutoff=2, NodeScoreCutoff=0.2, k‐core=2, and maxdepth=100. KEGG pathway enrichment analysis for DEGs in modules were conducted using clusterProfiler (version:3.8.1,http://bioconductor.org/packages/release/bioc/html/clusterProfiler.html) [24]. The p value was adjusted using the Benjamini and Hochberg (BH) approach [25], and p adjust < 0.05 was defined as significant.


### 2.6. TF-Gene Interaction Prediction

TFs that could interact with cluster genes were predicted using Overrepresentation Enrichment Analysis (ORA) with the online tool WebGestalt (version WebGestalt 2019,http://www.webgestalt.org/) [26]. The species was Rattus norvegicus. The threshold values were set as follows: enrichedgenes>5 and falsediscoveryrateFDR<0.05. The TF-target network was established by Cytoscape.


### 2.7. miRNA-Target Interaction Prediction

We used miRWalk2.0 (http://zmf.umm.uniheidelberg.de/apps/zmf/mirwalk2/generetsys-self.html) [27] to predict the upstream miRNAs that could target cluster genes. The species was rat. Meanwhile, the miRNA-target pairs should be simultaneously presented in miRWalk, miRanda, miRmap, RNAhybrid, and TargetScan databases. The miRNA-target regulatory network was visualized using Cytoscape.


### 2.8. DEG Corroboration with More Microarray Data

The microarray data GSE15041 deposited by Vega-Avelaira et al. [28] and GSE18803 deposited by Costigan et al. [29], which were developed on the [Rat230_2] Affymetrix Rat Genome 230 2.0 Array and GPL341 [RAE230A] Affymetrix Rat Expression 230A Array platforms, respectively, were also downloaded from NCBI GEO. The GSE15041 dataset contained 16 neonates (P10) or adult (8-12 wk) rat L4 and L5 DRG ipsilateral or contralateral seven days post-SNI/sham control samples. GSE18803 contained 24 neonates (P10) or adult (8-12 wk) rat DRG ipsilateral seven days post-SNI/sham control samples. The expression profiling data of three adult rat L4/L5 DRG ipsilateral SNI samples and three sham controls in the GSE15041 dataset and the expression profiling data of six adult rat DRG ipsilateral SNI samples and six sham controls in the GSE18803 dataset were selected to ensure sample source and age consistency. The data were downloaded in May 2020. The expression data of three adult rat L4/L5 DRG ipsilateral SNI samples and three sham controls in the GSE15041 dataset have been used to analyze the key genes associated with neuropathic pain in several studies [28, 30]. Similarly, the expression profiling data of six adult rat DRG ipsilateral SNI samples and six sham controls in the GSE18803 dataset have also been utilized to identify crucial genes involved in neuropathic pain development [29, 31–33]. We thus selected the two microarray data for DEG corroboration. In this study, the differential expression of probes between SNI samples and sham controls in the GSE15041 and GSE18803 databases were also analyzed using this method. The intersection of cluster genes obtained from the GSE30691 dataset and DEGs obtained from the GSE15041 and GSE18803 datasets were acquired to corroborate the differential expression of key cluster genes.


## 2.1. Data Sources

The microarray data GSE30691 deposited in the NCBI GEO (http://www.ncbi.nlm.nih.gov/geo/) database by Costigan et al. [12], which were generated on the GPL85 [RG_U34A] Affymetrix Rat Genome U34 Array platform, was downloaded. This dataset contained adult rat L4 and L5 dorsal root ganglion (DRG) ipsilateral samples from different sciatic nerve lesions, including SNI, spinal nerve ligation, and chronic constriction injury. Expression profiling of sciatic nerve lesion samples at 3, 7, 21, and 40 days after SNI and 3, 7, and 21 days of sham control samples (n=3 per time point; total=24) were extracted. The data were downloaded in May 2020.


## 2.2. Differential Expression Analysis

The differential expression of probes between SNI samples and sham controls was analyzed at time points of 3, 7, and 21 days post-SNI using the online tool GEO2R (http://www.ncbi.nlm.nih.gov/geo/geo2r/). The probes that did not map to gene symbols, or probes mapping to different genes were removed. The threshold value was set to p value < 0.05 and ∣logfoldchange∣>0.585. The number of differentially expressed probes and genes was counted. The unions of DEGs that had the same expression change tendency at three time points were collected for further analysis.


## 2.3. Mfuzz Clustering Analysis

To study the expression pattern of DEGs with the time-course injury, the expression values of DEGs in 0, 7, 21, and 40 days of SNI samples were extracted before using clustering analysis with Mfuzz (version 2.42.0,http://bioconductor.org/packages/release/bioc/html/Mfuzz.html) [16] in the R environment. In the clustering analysis, the optimal cluster number was calculated using default parameters. Meanwhile, the min score (membership) threshold was set to 0.6. For a given cluster, if the membership of two genes was higher, their expression trend was more similar, representing the importance of the gene in the cluster to a certain extent. Clusters with different change trends of genes over time were used for subsequent analysis. Cluster genes in each cluster were obtained using the online tool Metascape (http://metascape.org) [17].


## 2.4. Functional Enrichment Analysis

To investigate the function of cluster genes in each cluster, Gene Ontology (GO) [18] biological process (BP) and Kyoto Encyclopedia of Genes and Genomes (KEGG) [19] pathway enrichment analyses were performed, which are widely used for the functional annotation of large-scale data. The species was set as Rattus norvegicus, and the parameters were default, including min overlap: 3, p value cutoff: 0.01, and min enrichment: 1.5. After obtaining the terms conforming to the above parameters, further term clustering was performed based on the similarity of genes enriched in each term>0.3, and the term with the greatest statistical significance (p value minimum) was selected to represent the term cluster. The top 20 term clusters with the most significance were displayed. In addition, the interaction network of terms was constructed to further capture the relationship between terms in term clusters. The selection criteria for terms in this network were as follows: the terms with the best p value in each cluster of the top 20 clusters were selected, there were ≤15 terms in each cluster, ≤250 terms in total, and similarity of >0.3. After obtaining the interaction relations of terms, the interactive network of terms was established and visualized using Cytoscape (version 3.4.0, http://chianti.ucsd.edu/cytoscape-3.4.0/) [20].


## 2.5. PPI Network Construction and Module Analysis

The Search Tool for the Retrieval of Interacting Genes (STRING) database [21] was used to analyze interactions between the gene coding proteins. The PPI pairs between cluster genes were analyzed using a medium confidence of 0.4. The species was rat. The PPI network was constructed using Cytoscape software. The degree centrality of nodes was excavated using the CytoNCA plugin (version 2.1.6, http://apps.cytoscape.org/apps/cytonca) [22] without weighting, and the hub nodes were identified by ranking the node degree. In addition, the modules of the PPI network with nodes>3 were screened out using the Molecular Complex Detection (MCODE) (version 1.5.1, http://apps.cytoscape.org/apps/mcode) [23] plugin. The parameters were set to default, including DegreeCutoff=2, NodeScoreCutoff=0.2, k‐core=2, and maxdepth=100. KEGG pathway enrichment analysis for DEGs in modules were conducted using clusterProfiler (version:3.8.1,http://bioconductor.org/packages/release/bioc/html/clusterProfiler.html) [24]. The p value was adjusted using the Benjamini and Hochberg (BH) approach [25], and p adjust < 0.05 was defined as significant.


## 2.6. TF-Gene Interaction Prediction

TFs that could interact with cluster genes were predicted using Overrepresentation Enrichment Analysis (ORA) with the online tool WebGestalt (version WebGestalt 2019,http://www.webgestalt.org/) [26]. The species was Rattus norvegicus. The threshold values were set as follows: enrichedgenes>5 and falsediscoveryrateFDR<0.05. The TF-target network was established by Cytoscape.


## 2.7. miRNA-Target Interaction Prediction

We used miRWalk2.0 (http://zmf.umm.uniheidelberg.de/apps/zmf/mirwalk2/generetsys-self.html) [27] to predict the upstream miRNAs that could target cluster genes. The species was rat. Meanwhile, the miRNA-target pairs should be simultaneously presented in miRWalk, miRanda, miRmap, RNAhybrid, and TargetScan databases. The miRNA-target regulatory network was visualized using Cytoscape.


## 2.8. DEG Corroboration with More Microarray Data

The microarray data GSE15041 deposited by Vega-Avelaira et al. [28] and GSE18803 deposited by Costigan et al. [29], which were developed on the [Rat230_2] Affymetrix Rat Genome 230 2.0 Array and GPL341 [RAE230A] Affymetrix Rat Expression 230A Array platforms, respectively, were also downloaded from NCBI GEO. The GSE15041 dataset contained 16 neonates (P10) or adult (8-12 wk) rat L4 and L5 DRG ipsilateral or contralateral seven days post-SNI/sham control samples. GSE18803 contained 24 neonates (P10) or adult (8-12 wk) rat DRG ipsilateral seven days post-SNI/sham control samples. The expression profiling data of three adult rat L4/L5 DRG ipsilateral SNI samples and three sham controls in the GSE15041 dataset and the expression profiling data of six adult rat DRG ipsilateral SNI samples and six sham controls in the GSE18803 dataset were selected to ensure sample source and age consistency. The data were downloaded in May 2020. The expression data of three adult rat L4/L5 DRG ipsilateral SNI samples and three sham controls in the GSE15041 dataset have been used to analyze the key genes associated with neuropathic pain in several studies [28, 30]. Similarly, the expression profiling data of six adult rat DRG ipsilateral SNI samples and six sham controls in the GSE18803 dataset have also been utilized to identify crucial genes involved in neuropathic pain development [29, 31–33]. We thus selected the two microarray data for DEG corroboration. In this study, the differential expression of probes between SNI samples and sham controls in the GSE15041 and GSE18803 databases were also analyzed using this method. The intersection of cluster genes obtained from the GSE30691 dataset and DEGs obtained from the GSE15041 and GSE18803 datasets were acquired to corroborate the differential expression of key cluster genes.


## 3. Results

### 3.1. DEG Analysis

GSE30691 dataset analysis using the GEO2R online tool screened out 180 (109 up- and 71 downregulated genes), 230 (133 up- and 97 downregulated genes), and 153 (124 up- and 29 downregulated genes) DEGs between SNI samples and sham controls at the time points of 3, 7, and 21 days post-SNI. The unions of DEGs that had the same expression change tendency at three time points were analyzed, and a total of 306 DEGs were obtained for subsequent analysis.


### 3.2. Mfuzz Clustering Analysis

The expression trends of 306 DEGs at 0, 3, 7, 21, and 40 days post-SNI were analyzed using Mfuzz clustering analysis, and three clusters with different change trends over time were screened out (Figure1). Cluster 1 contained 77 DEGs and presented a trend of decreasing before increasing. The expression value reached the lowest level from three to seven days, and the overall expression of the SNI group was decreased compared to the control group (0 days after injury). Cluster 2 contained 13 DEGs and presented a trend of increasing before decreasing. The expression level rose to the highest level on the third day before gradually declining, and the overall SNI group expression was higher than that of the control group. Cluster 3 contained 49 DEGs and presented a trend of increasing before leveling off. The expression level rose to the highest level on day three before gradually declining, and the overall SNI group expression was higher than that of the control group.Figure 1
The results of Mfuzz clustering analysis. Three clusters with different change trends over time were screened out. The horizontal axis represents the days after spared nerve injury (SNI), and the vertical axis represents the relative expression value of cluster genes. The color from blue to red indicates that the membership value of this gene increases, which represents the gene’s importance in the cluster to some extent.


### 3.3. Functional Enrichment Analysis

We performed GO BP and KEGG pathway analyses for cluster genes in three clusters. In total, 344 GO BP and 33 KEGG pathway terms were significantly enriched by cluster genes in cluster 1, one GO BP and no KEGG pathway term was enriched by cluster genes in cluster 2, and 187 GO BP and 20 KEGG pathway terms were enriched by cluster genes in cluster 3. The top 20 enriched cluster terms are shown in Figure2(a). Moreover, several terms were simultaneously enriched by cluster genes in different clusters like GO:0009611: response to wounding, GO:0043408: regulation of MAPK cascade, and rno04010:MAPK signaling pathway. Subsequently, the term interaction network was constructed using important terms in clusters (Figure 2(b)).Figure 2
The significantly enriched cluster terms. (a) The top 20 enriched cluster terms by cluster genes in three clusters. The darker the color, the smaller thep value. (b) The term interaction network constructed by important terms in clusters. The color of different nodes represents different clusters, and the connection line represents the genetic similarity between terms.
(a)(b)


### 3.4. PPI Network Construction

The PPI network of cluster genes in three clusters included 117 nodes and 243 interactions (Figure3(a)). Notably, the nodes coded by cluster 1 genes were all downregulated, and the nodes coded by cluster 2 and 3 genes were all upregulated in this PPI network. Based on connectivity degree analysis, the nodes with greater connectivity degrees were protein kinase CAMP-activated catalytic subunit beta (Prkacb), complement C3 (C3), synaptosome-associated protein 25 (Snap25), activating transcription factor 3 (Atf3), and protein phosphatase 3 regulatory subunit B, alpha (Ppp3r1), which were considered as PPI network hub nodes. Notably, Prkacb, Ppp3r1, and Mapk9 were revealed to be remarkably enriched in the rno04010:MAPK signaling pathway. Further module analysis identified five modules with node>3 (Figure 3(b)). The KEGG pathway showed that genes in modules 1-5 were dramatically enriched in the synaptic vesicle cycle, complement and coagulation cascades, focal adhesion, wnt signaling pathway, and phagosome (Figure 3(c)).Figure 3
Protein-protein interaction (PPI) network constructed by cluster genes in three clusters, five modules, and KEGG pathway analysis for module genes. (a) The PPI network and (b) five modules. Red nodes represent cluster 1 genes, blue nodes represent cluster 2 genes, and green nodes represent cluster 3 genes. Circle nodes represent upregulated genes, and rhombus nodes represent downregulated genes. Node size represents the degree of connectivity. (c) KEGG pathway analysis for module genes. Bubble size indicates the number of enriched genes. The redder the bubble color, the more significant this pathway.
(a)(b)(c)


### 3.5. TF-Gene Regulatory Network Analysis

TFs that could interact with cluster genes were predicted using WebGestalt. The results showed that 15 TFs that could interact with cluster 1 genes were obtained, whereas TFs that could interact with cluster 2 and 3 genes were not predicted. The TF-gene regulatory network containing 15 TFs and 54 target cluster 1 genes was constructed (Figure4). Most of cluster 1 genes in this network were downregulated. According to the number of target TFs, the important cluster 1 genes that could interact with more TFs were protein phosphatase 2 regulatory subunit B gamma (Ppp2r2c), visinin-like 1 (Vsnl1), and fibroblast growth factor 9 (Fgf9). Fgf9 was also found to be enriched in the rno04010:MAPK signaling pathway.Figure 4
The transcription factor- (TF-) gene regulatory network. Red nodes indicate the cluster 1 genes, and yellow octagon nodes indicate predicted TFs. Rhomboid nodes represent downregulated genes and circular nodes represent upregulated genes. Arrow lines indicate TF-target relationship.


### 3.6. miRNA-Target Regulatory Network Analysis

In total, 4546 miRNA-target interactions, including 681 miRNAs and 103 target genes, were obtained using miRWalk2.0. Since there were so many miRNAs, the miRNAs withtargetclustergenes≥15 were screened out again. The miRNA-target regulatory network was established, containing 14 miRNAs and 70 target cluster genes (Figure 5). Among them, rno-miR-3583-5p was found to be targeted by Fgf9 and rno-miR-1912-3p was targeted by neuregulin 1 (Nrg1).Figure 5
The miRNA-target regulatory network. Red nodes represent cluster 1 genes, blue nodes represent cluster 2 genes, and green nodes represent cluster 3 genes. Circle nodes represent upregulated genes, and rhombus nodes represent downregulated genes. The pink triangle nodes represent the predicted miRNAs with the number oftargetgenes≥15. The gray lines represent target gene regulation by miRNA.


### 3.7. DEG Corroboration with More Microarray Data

To verify the important cluster genes, microarray data GSE15041 and GSE18803 were also downloaded from NCBI GEO and used to identify DEGs between SNI samples and sham controls at seven days post-SNI. The results showed that 830 (498 up- and 332 downregulated genes) and 188 (184 up- and four downregulated genes) DEGs were screened out between SNI samples and sham controls based on the GSE15041 and GSE18803 datasets, respectively. Since most of cluster 1 genes were downregulated according to the aforementioned analysis, we investigated the intersection of cluster 1 genes and downregulated genes obtained from GSE15041 and GSE18803 databases. The results showed that 14 intersected genes between cluster 1 genes and downregulated genes obtained from the GSE15041 database were obtained, includingNrg1 and Fgf9 (Figure 6(a)). Since the overall SNI group expression in cluster 2 and cluster 3 was higher than that in the control group, the intersection of cluster 2 genes or cluster 3 genes and upregulated genes obtained from GSE15041 and GSE18803 databases were analyzed. The results showed that only one intersected gene, named TIMP Metallopeptidase Inhibitor 1 (Timp1), was obtained between cluster 2 genes and upregulated genes obtained from GSE15041 and GSE18803 databases (Figure 6(b)). Nine intersected genes, including Atf3 and C3, were screened between cluster 3 genes and upregulated genes obtained from GSE15041 and GSE18803 databases (Figure 6(c)).Figure 6
The intersection of cluster genes with differentially expressed genes obtained from GSE15041 and GSE18803 databases. (a) The intersection of cluster 1 genes and downregulated genes obtained from GSE15041 and GSE18803 databases. (b) The intersection of cluster 2 genes and upregulated genes obtained from GSE15041 and GSE18803 databases. (c) The intersection of cluster 3 genes and upregulated genes obtained from GSE15041 and GSE18803 databases.
(a)(b)(c)


## 3.1. DEG Analysis

GSE30691 dataset analysis using the GEO2R online tool screened out 180 (109 up- and 71 downregulated genes), 230 (133 up- and 97 downregulated genes), and 153 (124 up- and 29 downregulated genes) DEGs between SNI samples and sham controls at the time points of 3, 7, and 21 days post-SNI. The unions of DEGs that had the same expression change tendency at three time points were analyzed, and a total of 306 DEGs were obtained for subsequent analysis.


## 3.2. Mfuzz Clustering Analysis

The expression trends of 306 DEGs at 0, 3, 7, 21, and 40 days post-SNI were analyzed using Mfuzz clustering analysis, and three clusters with different change trends over time were screened out (Figure1). Cluster 1 contained 77 DEGs and presented a trend of decreasing before increasing. The expression value reached the lowest level from three to seven days, and the overall expression of the SNI group was decreased compared to the control group (0 days after injury). Cluster 2 contained 13 DEGs and presented a trend of increasing before decreasing. The expression level rose to the highest level on the third day before gradually declining, and the overall SNI group expression was higher than that of the control group. Cluster 3 contained 49 DEGs and presented a trend of increasing before leveling off. The expression level rose to the highest level on day three before gradually declining, and the overall SNI group expression was higher than that of the control group.Figure 1
The results of Mfuzz clustering analysis. Three clusters with different change trends over time were screened out. The horizontal axis represents the days after spared nerve injury (SNI), and the vertical axis represents the relative expression value of cluster genes. The color from blue to red indicates that the membership value of this gene increases, which represents the gene’s importance in the cluster to some extent.


## 3.3. Functional Enrichment Analysis

We performed GO BP and KEGG pathway analyses for cluster genes in three clusters. In total, 344 GO BP and 33 KEGG pathway terms were significantly enriched by cluster genes in cluster 1, one GO BP and no KEGG pathway term was enriched by cluster genes in cluster 2, and 187 GO BP and 20 KEGG pathway terms were enriched by cluster genes in cluster 3. The top 20 enriched cluster terms are shown in Figure2(a). Moreover, several terms were simultaneously enriched by cluster genes in different clusters like GO:0009611: response to wounding, GO:0043408: regulation of MAPK cascade, and rno04010:MAPK signaling pathway. Subsequently, the term interaction network was constructed using important terms in clusters (Figure 2(b)).Figure 2
The significantly enriched cluster terms. (a) The top 20 enriched cluster terms by cluster genes in three clusters. The darker the color, the smaller thep value. (b) The term interaction network constructed by important terms in clusters. The color of different nodes represents different clusters, and the connection line represents the genetic similarity between terms.
(a)(b)


## 3.4. PPI Network Construction

The PPI network of cluster genes in three clusters included 117 nodes and 243 interactions (Figure3(a)). Notably, the nodes coded by cluster 1 genes were all downregulated, and the nodes coded by cluster 2 and 3 genes were all upregulated in this PPI network. Based on connectivity degree analysis, the nodes with greater connectivity degrees were protein kinase CAMP-activated catalytic subunit beta (Prkacb), complement C3 (C3), synaptosome-associated protein 25 (Snap25), activating transcription factor 3 (Atf3), and protein phosphatase 3 regulatory subunit B, alpha (Ppp3r1), which were considered as PPI network hub nodes. Notably, Prkacb, Ppp3r1, and Mapk9 were revealed to be remarkably enriched in the rno04010:MAPK signaling pathway. Further module analysis identified five modules with node>3 (Figure 3(b)). The KEGG pathway showed that genes in modules 1-5 were dramatically enriched in the synaptic vesicle cycle, complement and coagulation cascades, focal adhesion, wnt signaling pathway, and phagosome (Figure 3(c)).Figure 3
Protein-protein interaction (PPI) network constructed by cluster genes in three clusters, five modules, and KEGG pathway analysis for module genes. (a) The PPI network and (b) five modules. Red nodes represent cluster 1 genes, blue nodes represent cluster 2 genes, and green nodes represent cluster 3 genes. Circle nodes represent upregulated genes, and rhombus nodes represent downregulated genes. Node size represents the degree of connectivity. (c) KEGG pathway analysis for module genes. Bubble size indicates the number of enriched genes. The redder the bubble color, the more significant this pathway.
(a)(b)(c)


## 3.5. TF-Gene Regulatory Network Analysis

TFs that could interact with cluster genes were predicted using WebGestalt. The results showed that 15 TFs that could interact with cluster 1 genes were obtained, whereas TFs that could interact with cluster 2 and 3 genes were not predicted. The TF-gene regulatory network containing 15 TFs and 54 target cluster 1 genes was constructed (Figure4). Most of cluster 1 genes in this network were downregulated. According to the number of target TFs, the important cluster 1 genes that could interact with more TFs were protein phosphatase 2 regulatory subunit B gamma (Ppp2r2c), visinin-like 1 (Vsnl1), and fibroblast growth factor 9 (Fgf9). Fgf9 was also found to be enriched in the rno04010:MAPK signaling pathway.Figure 4
The transcription factor- (TF-) gene regulatory network. Red nodes indicate the cluster 1 genes, and yellow octagon nodes indicate predicted TFs. Rhomboid nodes represent downregulated genes and circular nodes represent upregulated genes. Arrow lines indicate TF-target relationship.


## 3.6. miRNA-Target Regulatory Network Analysis

In total, 4546 miRNA-target interactions, including 681 miRNAs and 103 target genes, were obtained using miRWalk2.0. Since there were so many miRNAs, the miRNAs withtargetclustergenes≥15 were screened out again. The miRNA-target regulatory network was established, containing 14 miRNAs and 70 target cluster genes (Figure 5). Among them, rno-miR-3583-5p was found to be targeted by Fgf9 and rno-miR-1912-3p was targeted by neuregulin 1 (Nrg1).Figure 5
The miRNA-target regulatory network. Red nodes represent cluster 1 genes, blue nodes represent cluster 2 genes, and green nodes represent cluster 3 genes. Circle nodes represent upregulated genes, and rhombus nodes represent downregulated genes. The pink triangle nodes represent the predicted miRNAs with the number oftargetgenes≥15. The gray lines represent target gene regulation by miRNA.


## 3.7. DEG Corroboration with More Microarray Data

To verify the important cluster genes, microarray data GSE15041 and GSE18803 were also downloaded from NCBI GEO and used to identify DEGs between SNI samples and sham controls at seven days post-SNI. The results showed that 830 (498 up- and 332 downregulated genes) and 188 (184 up- and four downregulated genes) DEGs were screened out between SNI samples and sham controls based on the GSE15041 and GSE18803 datasets, respectively. Since most of cluster 1 genes were downregulated according to the aforementioned analysis, we investigated the intersection of cluster 1 genes and downregulated genes obtained from GSE15041 and GSE18803 databases. The results showed that 14 intersected genes between cluster 1 genes and downregulated genes obtained from the GSE15041 database were obtained, includingNrg1 and Fgf9 (Figure 6(a)). Since the overall SNI group expression in cluster 2 and cluster 3 was higher than that in the control group, the intersection of cluster 2 genes or cluster 3 genes and upregulated genes obtained from GSE15041 and GSE18803 databases were analyzed. The results showed that only one intersected gene, named TIMP Metallopeptidase Inhibitor 1 (Timp1), was obtained between cluster 2 genes and upregulated genes obtained from GSE15041 and GSE18803 databases (Figure 6(b)). Nine intersected genes, including Atf3 and C3, were screened between cluster 3 genes and upregulated genes obtained from GSE15041 and GSE18803 databases (Figure 6(c)).Figure 6
The intersection of cluster genes with differentially expressed genes obtained from GSE15041 and GSE18803 databases. (a) The intersection of cluster 1 genes and downregulated genes obtained from GSE15041 and GSE18803 databases. (b) The intersection of cluster 2 genes and upregulated genes obtained from GSE15041 and GSE18803 databases. (c) The intersection of cluster 3 genes and upregulated genes obtained from GSE15041 and GSE18803 databases.
(a)(b)(c)


## 4. Discussion

Neuropathic pain is chronic pain with an elusive mechanism. To discover the possible mechanism, this study utilized a comprehensive bioinformatic approach to screen candidate genes associated with SNI-induced neuropathic pain. Consistent with previous findings by Costigan et al. [12], we also found that KCNS1 was downregulated in the SNI-induced neuropathic pain model, suggesting that our results were reliable. In addition to this, our analysis further revealed three clusters with different change trends over time after SNI. Prkacb, C3, and Atf3 were hub nodes in the PPI network, and Fgf9 was found to interact with more TFs. Prkacb and Fgf9 were significantly enriched in the MAPK signaling pathway. Moreover, rno-miR-3583-5p was found to be targeted by Fgf9, and rno-miR-1912-3p was targeted by neuregulin 1 (Nrg1). After corroborating microarray data with two other microarray datasets, key genes, like Nrg1 and Fgf9 in cluster 1, Timp1 in cluster 2, and Atf3 and C3 in cluster 3, were screened out to be implicated in SNI-induced neuropathic pain.Increasing evidence has revealed that nerve injury leads to p38 MAPK pathway activation in the spinal cord, consequently resulting in neuropathic pain development by regulating proinflammatory cytokine production [34–36]. Moreover, the p38 MAPK pathway is considered a promising therapeutic target for neuropathic pain [37]. Here, the MAPK signaling pathway was simultaneously enriched by cluster genes in different clusters, confirming the key role of this pathway in SNI-induced neuropathic pain. Moreover, Prkacb and Fgf9 were significantly enriched in the MAPK signaling pathway. Prkacb is a subunit of the cAMP-dependent protein kinase, which has a catalytic role in numerous cellular processes like cell proliferation, gene transcription, and differentiation [38]. Fgf9 is a member of the highly conserved FGF family, which is reportedly important for glial cell development in the nervous system [39]. Fgf9 silencing recapitulates the inhibitory effect of miR-182 overexpression on Schwann cell proliferation at an early stage following SNI [40], suggesting the potential role of Fgf9 after SNI. Although the key role of Prkacb and Fgf9 in neuropathic pain has not been fully disclosed, we speculate that Prkacb and Fgf9 may be implicated in SNI-induced neuropathic pain in rats by being involved in the MAPK signaling pathway.In addition toFgf9, other important genes, like Nrg1 in cluster 1, Timp1 in cluster 2, and Atf3 and C3 in cluster 3, were also found to be differentially expressed in SNI samples after corroborating microarray data with two other microarray datasets. Nrg1 reportedly plays a pivotal role in neural development and plasticity [41]. Wang et al. demonstrated that Nrg1 upregulation reversed the signs of SNI-induced neuropathic pain in rats, and modulating Nrg1 might exhibit therapeutic value for neuropathic pain treatment [42]. Timp1 is well known as an inhibitor of matrix metalloproteinases (MMPs) that are widely involved in pain development following a variety of injury and inflammatory conditions [43, 44]. Knight et al. revealed that Timp1 played a crucial role in pathological pain states associated with inflammation [45]. Moreover, Timp1 is upregulated in DRG and spinal cord tissues after tibial nerve transection, and it may be related to neuropathic pain following peripheral nerve injury [46]. A recent microarray analysis also showed that hub genes like C3 and Timp1 are closely related to SNI-induced neuropathic pain development, providing the theoretical basis for treatment of this disease [14]. In addition, ATF3 has been shown to be implicated in the Bortezomib-induced painful peripheral neuropathy [47]. Atf3 expression was found to be significantly increased after SNI operation, and it may be a key regulator in neuropathic pain development [48], which agrees with our bioinformatic results. Based on our results, we speculate that these genes may be involved in neuropathic pain development and could serve as the therapeutic targets of this disease. Furthermore, rno-miR-3583-5p was found to be targeted by Fgf9, and rno-miR-1912-3p was targeted by Nrg1. However, there are few reports about the roles of rno-miR-3583-5p and rno-miR-1912-3p in neuropathic pain. Given the potential role of Fgf9 and Nrg1, we speculate that rno-miR-3583-5p and rno-miR-1912-3p may also contribute to neuropathic pain development in rats.However, there are some limitations in this study. First, the sample size was small which influences the sample power stability. Second, we only used other microarray data to corroborate the identified key genes; no experimental validations were conducted to confirm our findings. In the future, more studies with large samples and experimental validations are needed to verify the role of critical genes and pathways in neuropathic pain development.


## 5. Conclusions

Our results reveal that key pathways like the MAPK signaling pathway and crucial genes likePrkacb, Nrg1, Fgf9, Timp1, C3, and Atf3 may be implicated in SNI-induced neuropathic pain development in rats. Our findings will provide new insight for designing effective neuropathic pain therapies.



---

*Source: 8822001-2020-11-21.xml*