Analysis of Potential Hub Genes for Neuropathic Pain Based on Differential Expression in Rat Models

Bai, Jie; Geng, Bin; Wang, Xingwen; Wang, Shenghong; Xia, Yayi

doi:https://doi.org/10.1155/2022/6571987

Pain Research and Management

On this page

Abstract Introduction Materials and Methods Results Discussion Conclusions Data Availability Ethical Approval Conflicts of Interest Authors’ Contributions Acknowledgments Supplementary Materials References Copyright Related Articles

Research Article | Open Access

Volume 2022 | Article ID 6571987 | https://doi.org/10.1155/2022/6571987

Analysis of Potential Hub Genes for Neuropathic Pain Based on Differential Expression in Rat Models

Jie Bai,^1,2Bin Geng,^2,3Xingwen Wang,^2,3Shenghong Wang,^2,3and Yayi Xia^2,3

Academic Editor: Vahid Rakhshan

Received02 Nov 2021

Revised03 Jan 2022

Accepted21 Jan 2022

Published03 Mar 2022

Abstract

Objective. Neuropathic pain (NP) is a type of intractable chronic pain with complicated etiology. The exact molecular mechanism underlying NP remains unclear. In this study, we searched for molecular biomarkers of NP. Methods. Differentially expressed genes (DEGs) were predicted by analyzing three NP-related microarray datasets in Gene Expression Omnibus with robust rank aggregation. A weighted gene coexpression network analysis was conducted to construct a network of differentially expressed genes, followed by the evaluation of correlations between gene sets and the determination of hub genes. The candidate genes from the key module were identified using a gene set enrichment analysis. Results. In total, 353 upregulated and 383 downregulated genes were obtained, among which five hub genes were determined to be related to pain phenotypes. Reverse transcription-quantitative polymerase chain reaction was performed to verify the expression of these hub genes in the dorsal root ganglia of rats with spared nerve injury, which revealed the decreased expression of EMC4. Hence, EMC4 was defined as a biomarker for NP development. Conclusions. The results of this study form a basis for further research into the mechanism of NP development and are expected to aid in the development of novel therapeutic strategies.

1. Introduction

As a type of chronic pain with complex etiology, neuropathic pain (NP) is characterized by hyperalgesia, numbness, and allodynia [1]. More than 6% of patients experience debilitating NP-related physical and emotional trauma. Damage to the sensory system afflicts the transmission of sensory signals, thus resulting in hyperalgesia symptoms [2]. Owing to the complex pathogenesis of NP, there is still no effective treatment. With the development of bioinformatics, gene expression datasets such as those in the Gene Expression Omnibus (GEO) have been widely utilized to construct genetic networks and identify the potential roles and functions of differentially expressed genes (DEGs) in the development of NP [3].

However, the study of molecular mechanisms has significant implications in the treatment of NP. Dorsal root ganglia (DRGs) have previously attracted the attention of researchers. DRGs consist of pseudounipolar neurons that transmit signals from the peripheral nerves to the dorsal horns via the neuronal cell bodies [4]. Evidence shows that genetic variations in DRGs are related to pain phenotypes [5]. Therefore, it is of importance to analyze gene expression changes in the DRGs after peripheral nerve injury for the understanding of the molecular mechanism underlying NP, which may contribute to the development of an effective therapeutic regimen.

A large number of gene studies on the molecular basis of NP have been performed, but most are based on the analysis of peripheral blood samples from patients to identify DEGs. Therefore, the non-spinal cord in situ genetic information obtained does not truly reflect the lesion site. Additionally, it is difficult to obtain a stable and reproducible phenotype owing to the complex background of patient samples [6].

Some investigators have utilized rat models with more stable phenotypes to map genetic changes at the site of NP lesions. For example, Yu H et al. screened for hub genes and measured the expression of these genes in a rat model of NP. As a result, the authors identified DEGs substantially enriched in “extracellular space,” and their Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis showed that DEGs were enriched in inflammatory disease and mitogen-activated protein kinase signaling pathways [7]. However, these findings were not assessed in rat models for further validation, and thus, there is still a lack of reliable analyses of molecular mechanisms and hub genes in NP.

In this context, we proposed a research hypothesis that potential core genes for NP could be identified by DEGs and replicated in animal models. This work was conducted to detect the expression of DEGs and determine their reproducibility in an animal model of NP. Because in situ information could not be obtained from human samples, a phenotypically stable and reproducible animal model was employed to ascertain the pathogenesis of NP. Our results provide novel insights into central genes and biological pathways in the pathogenesis of NP using bioinformatics analyses.

2. Materials and Methods

2.1. Selection of Gene Microarray Datasets

Microarray data were downloaded from the GEO database (http://www.ncbi.nlm.nih.gov/geo/) for expression profiling of neuralgia. Three independent microarray datasets, including GSE24982, GSE30691, and GSE63442, were selected. GSE24982 comprised data from 20 spinal nerve ligation-induced NP model animals and 20 sham cohort samples generated using the GPL1355 microarray platform (Affymetrix Rat Genome 230 2.0 Array). GSE30691 included data from 56 samples of three independent NP models, which were generated using the GPL85 microarray platform (Affymetrix Rat Genome U34 Array). GSE63442 consisted of data from six spinal nerve ligation model animals and six sham control samples produced using the GPL341 microarray platform (Affymetrix Rat Expression 230A Array). Raw data were preprocessed using R software (version 4.0.2; http://www.R-project.org/).

2.2. Identification of DEGs

Analysis of DEGs was carried out to retrieve genes with differential expression in the treatment group compared to the control group. Microarray data were normalized to predict DEGs between the control and treatment groups using the Limma package [8]. Robust rank aggregation (RRA) was applied to identify reliable DEGs, minimize inconsistencies, and integrate the results of the three datasets. Firstly, the genes in each dataset were ranked as per their fold-changes among the groups. Then, the “RRA” package was utilized to integrate the list of all ranked genes. The adjusted P value of all genes was obtained to indicate the possibility of their ranking in the resultant gene list. DEGs were screened out with the threshold of log2 (fold-change) >0.3 and adjusted < 0.05 as screening criteria, followed by the construction of a new data frame. Finally, a “pheatmap” package in R was used to visualize and plot the top 40 DEGs (the top 20 upregulated and downregulated genes were selected based on the adjusted value).

2.3. Functional Enrichment Analyses

Methods and tools have been developed to analyze and filter these datasets to produce smaller, more meaningful, and biologically relevant gene/protein lists. The functional enrichment analysis is an approach that can identify genes enriched in the datasets of molecular functions, biological processes, and pathways of interest. Functional enrichment analysis contributes to the focus of researchers on a specific gene of interest or a specific biological problem. Considering a threshold of < 0.05, the clusterProfiler package [9] was applied for Gene Ontology (GO) and KEGG pathway analyses. Then, the enrichment results were visualized using the GOplot package [10] to further analyze the biological functions of genes related to NP.

2.4. Identification and Verification of Hub Genes

To identify the modules most closely correlated with NP phenotypic traits, a weighted gene coexpression network analysis (WGCNA) was performed to construct a scale-free network, define a coexpression matrix and adjacency functions, and calculate the coefficients of different nodes. The GSE24982 dataset was selected for the identification of coexpression modules. Based on the RRA analysis results, the top 4000 genes (according to the P value) from GSE24982 were extracted, followed by WGCNA of expression data to obtain the key modules with the highest correlation with NP. Modules and candidate genes were acquired using the WGCNA package [11]. An unsigned topological overlap matrix was employed to construct the WGCNA network and determine coexpressed gene modules. The soft-thresholding power was set to 12 with the threshold for cutting height as 0.25, and the minimum number of gene modules as 30. The degree of correlation between genes and modules was determined as module membership (MM), and the degree of correlation between genes and clinical information was regarded as gene significance (GS). Genes with high connectivity tended to have crucial functions. Therefore, genes with high correlation were defined as hub genes in the candidate modules. The hub genes in the modules met the criteria of MM > 0.80 and GS > 0.70. The top 50 genes with the highest connectivity were selected to confirm the expression data in GSE63442 and GSE30691 using an independent t-test.

2.5. Gene Set Enrichment Analysis (GSEA) and Gene Set Variation Analysis (GSVA)

GSEA and GSVA refer to the analyses based on gene sets. As the name implies, a gene set is a collection of genes, and any number of genes together can be called a gene set. However, the gene set applied for analyses must have a certain biological significance. The most widely used gene set databases are GO and KEGG, one of which classifies genes according to GO, and the other integrates related genes as per metabolic pathways. In addition, genes can be assembled into biologically significant gene sets through transcription factor regulatory networks, coexpression networks, and lists of marker genes that define biological states. The R packages “clusterProfiler” and “GSVA” were employed to perform GSEA on candidate hub genes to yield the biological pathways related to these genes. The six samples in GSE63442 were grouped according to the median expression of the candidate hub genes: high and low expression groups. < 0.05 was considered to be statistically significant.

2.6. Animal Experiments

Our study was ratified and reviewed by the Animal Care and Use Committee of the Second Hospital of Lanzhou University (D2019-003). Twelve adult male Sprague Dawley rats (weighing 200–220 g) were purchased from the Lanzhou Institute of Veterinary Medicine, Chinese Academy of Agricultural Sciences. There were 2 groups in this study, spared nerve injury group (SNI, n = 6) and sham-operated group (Sham, n = 6). Animals were housed in cages with a 12-h light/dark cycle and food and water ad libitum.

2.7. Establishment of a Rat Model of NP

To establish an NP model, spared nerve injury (SNI) was induced in the left sciatic nerve of rats following a previously published protocol [12]. Briefly, after intraperitoneal injection with sodium pentobarbital (40 mg/kg), the skin of the left lateral thigh was incised, and the biceps femoris muscle was bluntly dissected to expose the terminal branches of the sciatic nerve. At the bifurcation point, the common peroneal nerve and tibial nerve were tightly ligated with 4–0 silk that was then severed at the distal end of the ligature. Approximately 3–5 mm of the distal end of the nerve was removed, and contact with the sural nerve was avoided during the surgery. Following the confirmation of complete hemostasis, the wound was sutured. In sham-operated rats, the sciatic nerve was exposed without ligation and incision. Six rats each were utilized as biological replicates in the Sham and SNI groups. The sample size was not predetermined based on a priori power calculations but was estimated based on previous literature and the authors’ experience.

2.8. Measurement of Mechanical Hyperalgesia

Von Frey filaments (Aesthesio, Danmic Global, USA) were adopted to assess the withdrawal threshold for mechanical hyperalgesia in rats as previously described [13]. Rats were placed in transparent boxes with a wire mesh platform to acclimatize for at least 30 min. Then, the lateral edge in the left hind paw of the rats was stimulated with von Frey filaments, and the stimulation intensity started from 2 g. The paw withdrawal threshold (PWT) was determined using the up-and-down method described by Chaplan et al. [13]. The measurement was repeated three times for each rat, and the response frequency of each filament force was recorded and expressed as a percentage. The mechanical withdrawal threshold of all rats was tested 3 days before surgery and 3, 7, 10, and 14 days after surgery. On day 14 after SNI surgery, all animals were euthanized by administering an overdose of sodium pentobarbital after the assessment of mechanical hyperalgesia. Then, the histological samples were collected.

2.9. Quantitative Real-Time Polymerase Chain Reaction (RT-qPCR) Analysis

L_3-5 DRG neurons were collected and frozen in liquid nitrogen after SNI or sham surgery (n = 6 in each group). Total RNA was extracted using a TRIzol reagent (Invitrogen, Carlsbad, USA) following the manufacturer’s protocols. SYBR Premix Ex Taq^TM was applied for qPCR. Gene expression was calculated using the 2^–△△CT method with glyceraldehyde-3-phosphate dehydrogenase (GAPDH) as a normalizer. The primers used for qPCR are listed in Table 1. The ggstatsplot package was used to perform t-test and created graphics.

2.10. Statistical Analysis

GraphPad Prism 8 was employed for statistical analysis. The data are summarized as the mean ± standard deviation, and the Shapiro–Wilk test was used to examine the normal distribution of data. The sample size was not predetermined based on a priori power calculations [14]. It was estimated based on previous literature and the authors’ experience. Repeated measures analysis of variance (ANOVA) (factors: group and time) was utilized to compare PWT evoked by mechanical stimulus between ipsilateral and contralateral paw of rats in the Sham and SNI groups. Student’s t-test was carried out to compare the expression of five hub genes between the Sham and SNI groups. A value < 0.05 was considered statistically significant. P value was adjusted based on Benjamini–Hochberg FDR using the R package fdrtool.

3. Results

3.1. Identification of DEGs

Figure 1 depicts the process flow of identification and functional analyses of DEGs. Three datasets in the GEO database were used for RRA analysis, which integrated DEGs from three GSE datasets and retrieved intersected DEGs. According to the RRA analysis results, 353 upregulated and 383 downregulated genes were filtered out based on a threshold of < 0.05 (Supplementary 1). The most significantly upregulated gene was REG3B ( = 3.85E − 07 and adjusted = 2.96E − 03), followed by ATF3 ( = 4.33E-07 and adjusted = 2.96E-03). The most significantly downregulated gene was DDYSL4 ( = 5.94E − 08 and adjusted = 8.12E-04), followed by KCNS3 ( = 1.92E-07 and adjusted = 1.31E-03).

The top 20 upregulated and downregulated genes are shown in a heatmap in Figure 2, where the specific FDR-adjusted P value for each gene is shown in Supplementary 4. Following RRA analysis, a PPI network was constructed to define the relationship between the proteins expressed by the DEGs (Supplementary Figure 1), which increases our understanding of protein functions and relationships. Research has extensively explored the significant roles of some genes identified, such as REG3B, ATF3, SPRR1A, GAL, and JUN, in NP.

3.2. Functional Enrichment Analyses of DEGs

The DEGs identified by RRA analysis were subjected to GO and KEGG functional enrichment analyses. “Regulation of ion transmembrane transport,” “sensory perception of pain,” “response to axon injury,” “regulation of membrane potential,” and “potassium ion transport” were the most substantially activated biological functions of the DEGs (Figures 3(a)–3(c)). KEGG analysis showed that the DEGs were remarkably enriched in 27 signaling pathways, including “complement and coagulation cascade,” “neuroactive ligand-receptor interaction,” and “ECM-receptor interaction” (Figure 3(d)). The pathways identified are closely related to the development of NP, suggesting that the obtained DEGs are indeed related to NP.

(a)

(b)

(c)

(d)

3.3. Identification and Validation of Hub Genes

A hub gene is a gene that plays a critical part in a biological proc gene. WGCNA was adopted to investigate functional modules and genes associated with clinical traits. The key modules were identified by setting the soft-thresholding power as 12, the cutting height as 0.25, and the minimum module size as 30 followed by the determination of six modules with sizes from 77 to 2209 genes (labeled with different colors in Figure 4). The correlations between module eigengenes and clinical traits are visualized in Figures 4(c)–4(e). From the heatmap of module–trait relationships, it was found that the turquoise module, which contained 2, 209 genes, shared the strongest correlation with clinical traits (Figures 4(d) and 4(f)). Based on a heatmap of module–feature correlations, the turquoise module was considered the key module (correlation coefficient = 0.71 and = 2E − 06; Figure 4(d)). The enrichment of DEGs of the turquoise module in GO and KEGG pathways is shown in Supplementary 3. As shown in Figure 4(d), genes in the blue module (correlation coefficient = 0.53 and = 9E − 04) were also correlated with NP traits. The top 50 genes with the highest connectivity (Supplementary 2) were extracted based on the screening criteria for candidate hub genes in the key modules (MM > 0.80 and GS > 0.70). Seven genes (MKI67, VOM2R75, TJP1, EXT1, FOXP1, RNASEH2C, and EMC4) were retrieved as candidate hub genes. Among these, five genes were upregulated, whereas only RNASEH2C and EMC4 were downregulated in the NP group (Figure 5).

(a)

(b)

(c)

(d)

(e)

(f)

Figure 4

WGCNA. (a) The Cluster dendrogram of the module eigengenes. (b) The selection of soft-thresholding powers through scale-free topology fit index and mean connectivity among genes. (c) The cluster dendrogram of all DEGs based on the dissimilarity measure and the assignment modules. (d) The heatmap of module-trait correlations. The number in each small box represented the corresponding correlation and value. (e) The identification of modules related to NP. The horizontal axis indicated the name of modules, and the vertical axis indicated the gene significance value. (f) The scatter plot of GS and MM for genes in the turquoise module. DEG, differentially expressed gene; WGCNA, weighted gene coexpression network analysis.

Figure 5

The selection and verification of the candidate hub genes. The expression of MKI67, VOM2R75, TJP1, EXT1, FOXP1, EMC4, and RNASEH2C differed between the two groups. The genes with the highest connectivity were screened out by WGCNA. To validate the expression data of these genes, the GSE63442 and GSE30691 datasets were selected for pairwise validation by independent t-test. The ggstatsplot package was used to perform t-test and plot graphs. GEO, Gene Expression Omnibus. Statistical analysis was conducted using the independent t-test. Plots represented mean ± 95% confidence interval (CI).

3.4. Screening of the Five Hub Genes by GSEA and GSVA

To ascertain the biological functions of the seven candidate hub genes involved in NP development, GSEA and GSVA were conducted on the GSE63442 dataset. As exhibited in Figures 6 and 7, and Supplementary Figure 2, the genes in the groups with high expression of MKI67, TJP1, EXT1, RNASEH2C, or EMC4 were enriched in “sodium channel activity” in accordance with the GO term enrichment analysis. The genes of the MKI67, TJP1, RNASEH2C, and EMC4 high-expression groups were enriched in “cation channel activity” and “cation channel complex.” Further, TJP1, RNASEH2C, and EMC4 were all enriched in the pathways “cytosolic DNA sensing pathway,” “ECM-receptor interaction,” and “focal adhesion.” Therefore, these five hub genes (MKI67, TJP1, EXT1, RNASEH2C, and EMC4) were selected for further verification.

(a)

(b)

(c)

(d)

(e)

(f)

(a)

(b)

(c)

(d)

(e)

(f)

3.5. Behavioral Changes of Rats following SNI

Peripheral nerve injury can cause hyperalgesia and allodynia in affected rats [15]. Following SNI, the rats developed mechanical allodynia-like behavior, which was evidenced by the reduction of the ipsilateral withdrawal threshold from days 3 to 14 after injury (n = 6, male) (Figure 8(a), Table 2 and Supplementary 5). The ipsilateral withdrawal threshold of rats in the SNI group exhibited conspicuously lower values than that in the sham group (3 days, = 0.045; 10 days, = 0.008; 7, 14 days, < 0.001). From preoperative stage to 14 days after injury, there were no apparent differences between the Sham and SNI groups in terms of the withdrawal threshold of the contralateral paw (3 days, = 0.965; 7 days, = 0.992; 10 days, = 0.922; 14 days, = 0.783). The results illustrate the successful establishment of the NP rat model.

(a)

(b)

Figure 8

PWT changes after SNI surgery in rats and experimental validation of the expression of hub genes. (a) Compared with the sham-operated rats, the PWT of SNI rats decreased from day 3 after surgery. The withdrawal threshold was evaluated by von Frey filaments as a response evoked by a mechanical stimulus over time. The data in the sham and SNI groups were normally distributed. Repeated measures ANOVA followed by Tukey’s multiple comparisons test was used to evaluate mechanical hyperalgesia, and data were expressed as mean ± standard deviation. , < 0.05, , < 0.01, , < 0.001, compared with the sham-ips group on the corresponding days. (b) The expression of the five hub genes in DRG tissue after SNI surgery was validated by RT-qPCR. The data in the two groups were normally distributed. Unpaired Student’s (t) test was carried out, and the data were expressed as mean ± standard deviation. , < 0.01, , < 0.001 compared with the Sham group. DRG, dorsal root ganglia; SNI, spared nerve injury; PWT, paw withdrawal threshold.

3.6. Validation of Hub Genes Using RT-qPCR

RT-qPCR was performed to measure the mRNA levels of the five hub genes. Compared with the sham-operated rats, MKI67 expression increased ( = 0.001; Figure 8(b)), and EMC4 expression clearly decreased in rats after SNI surgery ( = 0.002; Figure 8(b), Table 3 and Supplementary 6). These results are consistent with those of the microarray analyses. However, the expression of the other genes was not appreciably different between the SNI and sham-operated rats (TJP1, = 0.793; EXT1, = 0.053; RNASEH2C, = 0.151; Figure 8(b), Table 3, and Supplementary 6).

4. Discussion

NP is attributable to either nervous system dysfunction or neuronal damage, thus leading to abnormal pain. As a complex disorder, NP is often accompanied by abnormal gene expression [16].

The management of NP remains a challenge for clinicians. Therefore, there is an urgent need to deepen our understanding of the genetic changes behind the pathogenesis of NP to find novel potential genetic targets for NP treatment.

In recent years, bioinformatics technology has been commonly utilized to dissect the molecular mechanisms underlying NP and has revealed several relevant DEGs [17]. ATF3, JUN, and GPR151 are mainly correlated with cytokine–cytokine receptor interaction and p53 signaling in the SNI model [18]. CCL2, FCER1G, NF-κB1, RAC2, and C1Q [19], p53 [20], CCL3, CTLA2B, ATF3, PLEK, and TGIF1 [21] are associated with the pathogenesis of NP.

In the present work, three gene expression datasets (GSE24982, GSE30691, and GSE63442) were integrated for RRA analysis to screen out DEGs. RRA analysis is a widely utilized tool for the integration of genome-wide gene expression data from different datasets and the identification of the key genes that are most likely to be implicated in the development of the disease under investigation [22].The DEGs identified comprise REG3B [23] and ATF3 [19, 24], which have been reported to assume a major role in the etiology of NP. After RRA analysis, enrichment analyses were conducted on 736 DEGs, which elucidated that they were mainly enriched in ion transmembrane transport [25], membrane potential regulation [26], sensory perception of pain [27], response to axon injury [28], and potassium ion transport [29]. GO and KEGG enrichment analyses uncovered that these functions were associated with the occurrence and development of NP. Additionally, complement and coagulation cascade [30], neuroactive ligand-receptor [31], and ECM-receptor interactions [32] were signaling pathways related to various functions mediated by NP. The genes associated with these pathways were also observed in our analysis. Based on the above findings, we further evaluated the roles of these genes in NP development.

WGCNA and coexpression network analysis were applied to find the hub genes involved in the pathogenesis of NP. The turquoise module that showed the highest connectivity with traits was chosen to further identify the genes in the module. Seven candidate hub genes (MKI67, VOM2R75, TJP1, EXT1, FOXP1, RNASEH2C, and EMC4) were retrieved after filtering for connectivity, GS, and MM values. To study the biological functions of these hub genes, they underwent GSEA and GSVA. The results show that most were associated with cation channel activity. Finally, MKI67, TJP1, EXT1, RNASEH2C, and EMC4 were selected as hub genes for further confirmation.

The expression of the five hub genes was determined in an SNI rat model. As reflected by RT-qPCR results, MKI67 expression was elevated, and EMC4 expression was diminished in SNI model rats compared to the sham-operated rats, which are consistent with the results of microarray analysis. There were differences in animal models and experimental conditions, and only a one-time point was selected for verification in this study, which might cause differences in the final results. Therefore, more rigorous experiments are warranted in the future for the comprehensive exploration of the role of these genes. However, our results suggest that EMC4 is a very sensitive factor and is potentially an important indicator of prognosis or prediction.

NP results from nerve injury, in which glial cell activation is one of the most prominent characteristics. Proliferation, upregulated cell surface markers and receptors, and functional changes are typical features associated with glial cell activation [33].The nuclear protein MKI67 is a well-known proliferation marker that is employed to assess cell proliferation. The number of Ki67-positive astrocytes and microglia is enhanced following peripheral nerve injury in rats, and the proliferation of glial cells may contribute to central sensitization [34]. Thus, MKI67 may participate in the pathological process of NP and may be the basis for the development and maintenance of hyperalgesia [35].

There are very few reports on the role of EMC4 and no reports on its role in NP in PubMed. EMC4 is a subunit of the endoplasmic reticulum (ER) membrane protein complex (EMC). EMC is a highly conserved oligomeric complex located on the ER membrane, which is essential for the folding and lipid transport of transmembrane proteins [36]. In all organisms evaluated, EMC destruction results in a pleiotropic phenotype. The phenotypes associated with EMC4 disruption are stress response activation in organisms and cells [37].

The ER is a pivotal organelle related to maintaining Ca²⁺ homeostasis, cell death signaling, and posttranslational modification [38]. Numerous factors, such as cellular stress (glucose deficiency and depletion of ER Ca²⁺ reserves), can trigger imbalances in the structure and function of ER, thereby leading to ER stress [39]. The ER stress response may be a crucial factor in the formation of a wide range of inflammatory diseases and neuroinflammation [40]. Nerve injury in an NP model induces an ER stress response in the DRG [41]. The repression of ER stress effectively relieves NP, which is regarded as an indicator for ER stress [42]. Therefore, we speculated that ER stress was correlated with the induction of NP. Our data show that EMC4 expression was dramatically reduced on day 14 after SNI surgery. This decrease in EMC4 expression may be associated with post-SNI pain, and ER stress may exert effects. Whether and how classic pain targets respond to ER stress or ER stress, which is alleviated, remains an active research topic. Hence, further work is needed to explore the related mechanisms.

The hub genes in NP were explored in this study. Unlike in prior studies, we did not return to the database for validation. The screened hub genes were validated in an NP rat model, which likely yields more reliable results. Relying on the reproducibility and phenotypic consistency of the animal model, we can continue to map the subsequent downstream signaling molecular pathways. Transgenic animals were also obtained through knockout or knock-in techniques to further verify whether the hub genes we screened facilitate or suppress NP.

However, our study has some limitations. The lack of sample size predetermination based on a priori power calculations and the small sample size available for analyses may obscure rare interactions, and further studies are required to dissect the relevant mechanisms.

5. Conclusions

EMC4 was downregulated in the SNI rat model and was a sign of NP occurrence. This finding provides novel insights into the mechanism of NP development and the associated therapeutic targets. In the future, we will ascertain the specific action mechanism of EMC4 in the development of NP to clarify its function in the pathogenesis of this condition.

Data Availability

The datasets of GSE24982, GSE30691, and GSE63442 can be downloaded from the website https://www.ncbi.nlm.nih.gov/geo/and other datasets analyzed during the current study are available from the corresponding authors upon reasonable request.

Ethical Approval

All institutional and national guidelines were followed.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this work.

Authors’ Contributions

Jie Bai conceived, designed, and performed the experiments, analyzed and assembled the data, prepared the figures and tables, reviewed drafts, and approved the final draft. Bin Geng conceived and designed the experiments. Xingwen Wang and Shenghong Wang prepared the figures and tables.

Acknowledgments

The authors thank Dr. Yi Sun for advice and assistance with the bioinformatics analysis. This work was supported by the National Natural Science Foundation of China under (grant nos. 81874017, 81960403, and 82060405), the Natural Science Foundation of Gansu Province of China (grant nos. 20JR5RA320 and 20JR10RA726), and the Cuiying Scientific and Technological Innovation Program of Lanzhou University Second Hospital (grant nos. CY2017-ZD02 and CY2018-QN18).

Supplementary Materials

Supplementary Figure 1. PPI network of DEGs from RRA analysis. DEGs, differentially expressed genes; PPI, protein-protein interaction; RRA, robust rank aggregation. Supplementary Figure 2. GSEA of candidate hub genes in the GEO dataset. (A) C, (E) and (G) GSEA of the single candidate hub genes in GO terms. (A) VOM2R75, (C) TJP, (E) FOXP1, (G) EXT1, and(B, (D) F, and (H) GSEA of the single candidate hub genes in KEGG pathways. (B) VOM2R75, (D) TJP, (F) FOXP1, and (H) EXT1. GEO, Gene Expression Omnibus; GO, Gene Ontology; GSEA, gene set enrichment analysis. Supplementary 1. Upregulated and downregulated genes obtained by RRA analysis. Supplementary 2. Top 50 genes with the highest connectivity. Supplementary 3. Enrichment of DEGs in the turquoise module in GO and KEGG pathways. Supplementary 4.Heatmap of the top 20 upregulated and downregulated genes identified by RRA analysis. Supplementary 5. Repeated measure ANOVA followed by Tukey’s multiple comparison test was used to evaluate the paw withdrawal threshold (PWT) between the Sham and SNI groups in Figure 8(a) and Table 2. Supplementary 6.Expression of the five hub genes in the Sham and SNI groups measured using RT-qPCR in Figure 8(b) and Table 3. (Supplementary Materials)

References

B. P. Murnion, “Neuropathic pain: current definition and review of drug treatment,” Australian Prescriber, vol. 41, no. 3, pp. 60–63, 2018.
View at: Publisher Site | Google Scholar
D. Burke, B. M. Fullen, D. Stokes, and O. Lennon, “Neuropathic pain prevalence following spinal cord injury: a systematic review and meta-analysis,” European Journal of Pain, vol. 21, no. 1, pp. 29–44, 2017.
View at: Publisher Site | Google Scholar
J. Zhou, Y. Fan, and H. Chen, “Analyses of long non-coding RNA and mRNA profiles in the spinal cord of rats using RNA sequencing during the progression of neuropathic pain in an SNI model,” RNA Biology, vol. 14, no. 12, pp. 1810–1826, 2017.
View at: Publisher Site | Google Scholar
J. M. Jacobs, R. M. Macfarlane, and J. B. Cavanagh, “Vascular leakage in the dorsal root ganglia of the rat, studied with horseradish peroxidase,” Journal of the Neurological Sciences, vol. 29, no. 1, pp. 95–107, 1976.
View at: Publisher Site | Google Scholar
M. Parisien, S. Khoury, A.-J. Chabot-Doré et al., “Effect of human genetic variability on gene expression in dorsal root ganglia and association with pain phenotypes,” Cell Reports, vol. 19, no. 9, pp. 1940–1952, 2017.
View at: Publisher Site | Google Scholar
H. H. Andersen, M. Duroux, and P. Gazerani, “MicroRNAs as modulators and biomarkers of inflammatory and neuropathic pain conditions,” Neurobiology of Disease, vol. 71, pp. 159–168, 2014.
View at: Publisher Site | Google Scholar
H. Yu, P. Zhang, Y.-R. Chen et al., “Temporal changes of spinal transcriptomic profiles in mice with spinal nerve ligation,” Frontiers in Neuroscience, vol. 13, p. 1357, 2019.
View at: Publisher Site | Google Scholar
M. E. Ritchie, B. Phipson, Y. F. Hu, and C. W. Law, “limma powers differential expression analyses for RNA-sequencing and microarray studies,” Nucleic Acids Res, vol. 43, no. 7, p. e47, 2015.
View at: Publisher Site | Google Scholar
G. Yu, L.-G. Wang, Y. Han, and Q.-Y. He, “clusterProfiler: an R package for comparing biological themes among gene clusters,” OMICS: A Journal of Integrative Biology, vol. 16, no. 5, pp. 284–287, 2012.
View at: Publisher Site | Google Scholar
W. Walter, F. Sánchez-Cabo, and M. Ricote, “GOplot: an R package for visually combining expression data with functional analysis: fig. 1,” Bioinformatics, vol. 31, no. 17, pp. 2912–2914, 2015.
View at: Publisher Site | Google Scholar
P. Langfelder and S. Horvath, “WGCNA: an R package for weighted correlation network analysis,” BMC Bioinformatics, vol. 9, no. 1, p. 559, 2008.
View at: Publisher Site | Google Scholar
I. Decosterd and C. J. Woolf, “Spared nerve injury: an animal model of persistent peripheral neuropathic pain,” Pain, vol. 87, no. 2, pp. 149–158, 2000.
View at: Publisher Site | Google Scholar
S. R. Chaplan, F. W. Bach, J. W. Pogrel, J. M. Chung, and T. L. Yaksh, “Quantitative assessment of tactile allodynia in the rat paw,” Journal of Neuroscience Methods, vol. 53, no. 1, pp. 55–63, 1994.
View at: Publisher Site | Google Scholar
F. Faul, E. Erdfelder, A.-G. Lang, and A. Buchner, “GPower 3: a flexible statistical power analysis program for the social, behavioral, and biomedical sciences,” Behavior Research Methods, vol. 39, no. 2, pp. 175–191, 2007.
View at: Publisher Site | Google Scholar
G. Yang, L. Sun, and G. Yang, “Spared nerve injury model of neuropathic pain in mice,” Bio-protocol, vol. 8, no. 6, 2018.
View at: Publisher Site | Google Scholar
C. A. von Hehn, R. Baron, and C. J. Woolf, “Deconstructing the neuropathic pain phenotype to reveal neural mechanisms,” Neuron, vol. 73, no. 4, pp. 638–652, 2012.
View at: Publisher Site | Google Scholar
J. A. Yang, J. M. He, J. M. Lu, and L. J. Jie, “Jun, Gal, Cd74, and C1qb as potential indicator for neuropathic pain,” Journal of Cellular Biochemistry, vol. 119, no. 6, pp. 4792–4798, 2018.
View at: Publisher Site | Google Scholar
S. Tang, H. Jing, Z. Huang et al., “Identification of key candidate genes in neuropathic pain by integrated bioinformatic analysis,” Journal of Cellular Biochemistry, vol. 121, no. 2, pp. 1635–1648, 2019.
View at: Publisher Site | Google Scholar
J. Wang, S.-H. Ma, R. Tao, L.-J. Xia, L. Liu, and Y.-H. Jiang, “Gene expression profile changes in rat dorsal horn after sciatic nerve injury,” Neurological Research, vol. 39, no. 2, pp. 176–182, 2016.
View at: Publisher Site | Google Scholar
Y. Gao, N. Sun, L. Wang et al., “Bioinformatics analysis identifies p53 as a candidate prognostic biomarker for neuropathic pain,” Frontiers in Genetics, vol. 9, 2018.
View at: Publisher Site | Google Scholar
G. Zhang and P. Yang, “Bioinformatics genes and pathway analysis for chronic neuropathic pain after spinal cord injury,” BioMed Research International, vol. 2017, pp. 1–11, 2017.
View at: Publisher Site | Google Scholar
R. Kolde, S. Laur, P. Adler, and J. Vilo, “Robust rank aggregation for gene list integration and meta-analysis,” Bioinformatics, vol. 28, no. 4, pp. 573–580, 2012.
View at: Publisher Site | Google Scholar
J. Li, H. Shi, H. Liu et al., “Nerve injury-induced neuronal PAP-I maintains neuropathic pain by activating spinal microglia,” Journal of Neuroscience, vol. 40, no. 2, pp. 297–310, 2020.
View at: Publisher Site | Google Scholar
I. Bravo‐Caparrós, M. C. Ruiz‐Cantero, G. Perazzoli et al., “Sigma‐1 receptors control neuropathic pain and macrophage infiltration into the dorsal root ganglion after peripheral nerve injury,” The FASEB Journal, vol. 34, no. 4, pp. 5951–5966, 2020.
View at: Publisher Site | Google Scholar
M. Devor, “Sodium channels and mechanisms of neuropathic pain,” The Journal of Pain, vol. 7, no. 1, pp. S3–S12, 2006.
View at: Publisher Site | Google Scholar
S. G. Waxman and G. W. Zamponi, “Regulating excitability of peripheral afferents: emerging ion channel targets,” Nature Neuroscience, vol. 17, no. 2, pp. 153–163, 2014.
View at: Publisher Site | Google Scholar
L. Meyer, O. Taleb, C. Patte-Mensah, and A.-G. Mensah-Nyagan, “Neurosteroids and neuropathic pain management: basic evidence and therapeutic perspectives,” Frontiers in Neuroendocrinology, vol. 55, p. 100795, 2019.
View at: Publisher Site | Google Scholar
A. J. Davies, S. Rinaldi, M. Costigan, and S. B. Oh, “Cytotoxic immunity in peripheral nerve injury and pain,” Frontiers in Neuroscience, vol. 14, p. 142, 2020.
View at: Publisher Site | Google Scholar
C. Gu and J. Barry, “Function and mechanism of axonal targeting of voltage-sensitive potassium channels,” Progress in Neurobiology, vol. 94, no. 2, pp. 115–132, 2011.
View at: Publisher Site | Google Scholar
Y. Gu, Z. Qiu, N. Cheng, C. Chen, Z. Hei, and X. Li, “Identification of potential mechanism and hub genes for neuropathic pain by expression‐based genome‐wide association study,” Journal of Cellular Biochemistry, vol. 120, no. 4, pp. 4912–4923, 2019.
View at: Publisher Site | Google Scholar
Z. Du, S. Yin, X. Song et al., “Identification of differentially expressed genes and key pathways in the dorsal root ganglion after chronic compression,” Frontiers in Molecular Neuroscience, vol. 13, p. 71, 2020.
View at: Publisher Site | Google Scholar
J. Zhou, Q. Xiong, H. Chen, C. Yang, and Y. Fan, “Identification of the spinal expression profile of non-coding RNAs involved in neuropathic pain following spared nerve injury by sequence analysis,” Frontiers in Molecular Neuroscience, vol. 10, p. 91, 2017.
View at: Publisher Site | Google Scholar
S. Echeverry, X. Q. Shi, and J. Zhang, “Characterization of cell proliferation in rat spinal cord following peripheral nerve injury and the relationship with neuropathic pain,” Pain, vol. 135, no. 1, pp. 37–47, 2008.
View at: Publisher Site | Google Scholar
J. L. Watson, T. J. Hala, R. Putatunda, D. Sannie, and A. C. Lepore, “Persistent at-level thermal hyperalgesia and tactile allodynia accompany chronic neuronal and astrocyte activation in superficial dorsal horn following mouse cervical contusion spinal cord injury,” PLos One, vol. 9, no. 9, p. e109099, 2014.
View at: Publisher Site | Google Scholar
R. A. Slivicki, S. S. Mali, and A. G. Hohmann, “Voluntary exercise reduces both chemotherapy-induced neuropathic nociception and deficits in hippocampal cellular proliferation in a mouse model of paclitaxel-induced peripheral neuropathy,” Neurobiology of Pain, vol. 6, p. 100035, 2019.
View at: Publisher Site | Google Scholar
G. Savidis, W. M. McDougall, P. Meraner et al., “Identification of zika virus and dengue virus dependency factors using functional genomics,” Cell Reports, vol. 16, no. 1, pp. 232–246, 2016.
View at: Publisher Site | Google Scholar
P. J. Chitwood and R. S. Hegde, “The role of EMC during membrane protein biogenesis,” Trends in Cell Biology, vol. 29, no. 5, pp. 371–384, 2019.
View at: Publisher Site | Google Scholar
R. J. Kaufman, “Stress signaling from the lumen of the endoplasmic reticulum: coordination of gene transcriptional and translational controls,” Genes & Development, vol. 13, no. 10, pp. 1211–1233, 1999.
View at: Publisher Site | Google Scholar
E. Zhang, M.-H. Yi, N. Shin et al., “Endoplasmic reticulum stress impairment in the spinal dorsal horn of a neuropathic pain model,” Scientific Reports, vol. 5, no. 1, p. 11555, 2015.
View at: Publisher Site | Google Scholar
A. M. Deslauriers, A. Afkhami-Goli, A. M. Paul et al., “Neuroinflammation and endoplasmic reticulum stress are coregulated by crocin to prevent demyelination and neurodegeneration,” The Journal of Immunology, vol. 187, no. 9, pp. 4788–4799, 2011.
View at: Publisher Site | Google Scholar
L. Dong, A. O. Odeleye, K. L. Jordan-Sciutto, and B. A. Winkelstein, “Painful facet joint injury induces neuronal stress activation in the DRG: implications for cellular mechanisms of pain,” Neuroscience Letters, vol. 443, no. 2, pp. 90–94, 2008.
View at: Publisher Site | Google Scholar
C. Hetz, E. Chevet, and H. P. Harding, “Targeting the unfolded protein response in disease,” Nature Reviews Drug Discovery, vol. 12, no. 9, pp. 703–719, 2013.
View at: Publisher Site | Google Scholar

Copyright

Copyright © 2022 Jie Bai et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

PDF Download Citation

Download other formats

Order printed copies

Views

385

Downloads

585

Citations