Identification of candidate genes and miRNAs associated with neuropathic pain induced by spared nerve injury.

Neuropathic pain (NP) is a complex, chronic pain condition caused by injury or dysfunction affecting the somatosensory nervous system. This study aimed to identify crucial genes and miRNAs involved in NP. Microarray data (access number GSE91396) were downloaded from the Gene Expression Omnibus (GEO). Murine RNA‑seq samples from three brain regions [nucleus accumbens, (NAc); medial prefrontal cortex, (mPFC) and periaqueductal gray, (PAG)]were compared between the spared nerve injury (SNI) model and a sham surgery. After data normalization, differentially expressed RNAs were screened using the limma package and functional enrichment analysis was performed with Database for Annotation, Visualization and Integrated Discovery. The microRNA (miRNA/miR)‑mRNA regulatory network and miRNA‑target gene‑pathway regulatory network were constructed using Cytoscape software. A total of 2,776 differentially expressed RNAs (219 miRNAs and 2,557 mRNAs) were identified in the SNI model compared with the sham surgery group. A total of two important modules (red and turquoise module) were found to be related to NP using weighed gene co‑expression network analysis (WGCNA) for the 2,325 common differentially expressed RNAs in three brain regions. The differentially expressed genes (DEGs) in the miRNA‑mRNA regulatory network were significantly enriched in 21 Gene Ontology terms and five pathways. A total of four important DEGs (CXCR2, IL12B, TNFSF8 and GRK1) and five miRNAs (miR‑208a‑5p, miR‑7688‑3p, miR‑344f‑3p, miR‑135b‑3p and miR‑135a‑2‑3p) were revealed according to the miRNA‑target gene‑pathway regulatory network to be related to NP. Four important DEGs (CXCR2, IL12B, TNFSF8 and GRK1) and five miRNAs (miR‑208a‑5p, miR‑7688‑3p, miR‑344f‑3p, miR‑135b‑3p and miR‑135a‑2‑3p) were differentially expressed in SNI, indicating their plausible roles in NP pathogenesis.

Introduction

Neuropathic pain (NP) is a severe chronic condition caused by injury or dysfunction affecting the somatosensory nervous system (1). A total of >100 million Americans are thought to be affected by chronic pain (2). NP is characterized by a wide range of sensory, cognitive and affective symptoms, such as tactile allodynia (burning pain resulting from noxious stimuli), hyperalgesia, and spontaneous pain (3). Most NP patients may suffer from depression, anxiety disorders, or other negative moods and various therapies offer only partial relief to 40-60% of NP patients (4). Despite great progress in understanding NP's prognosis, a number of patients respond poorly to current therapies. Therefore, it is necessary and urgent to explore NP's molecular mechanisms, especially the association between NP and anxiety, depression, and other mood disorders.

Chronic pain is associated with adaptations in several brain networks involved in mood, motivation and reward. Evidence indicates that the nucleus accumbens (NAc) is necessary for expressing NP-like behavior (5), whereas the medial prefrontal cortex (mPFC) is linked to a wide variety of cognitive functions critical for social behavior, including personal traits (6). NP can lead to mPFC remodeling, which associates it with emotional regulation of chronic pain (7). Periaqueductal gray (PAG) is the primary center for descending pain modulation (8) that receive schronic pain and temperature signals from the spinomesencephalic tract, hinting at its key role in NP progression (9). Microarray profiling was previously used to explore the altered gene expression of the murine nervous system following nerve injury (10,11). This led to the identification of IL-6, c-Jun and Plau as crucial genes (12) and DNA binding, cell cycle, and forkhead box protein as a major signaling pathway involved in NP (13). Despite these investigations profiling post nerve injury gene expression using microarray analysis, NP's pathological mechanisms remain poorly understood.

Recently, Descalzi et al (11) used RNA-sequencing technology to explore NP; they found NP can affect the expression of multiple genes in three distinct brain regions. However, they did not explore differentially expressed mRNAs (DEMs) and microRNAs (miRNAs/miR), especially the interactions between them. Based on the microarray data deposited by Descalzi et al, DEMs and differentially expressed miRNAs were screened from three distinct brain regions (NAc, mPFC, and PAG) of the spared nerve injury (SNI) and sham surgery murine model.

miRNA is a subset of non-coding, small RNAs, about 22nt long, that can combine with the 3′ untranslated region of messenger RNA (mRNA) to regulate post-transcription gene expression (14). miRNA binds to mRNA, forming a complex regulatory network that plays a vital role in a number of biological processes, such as cellular proliferation, apoptosis, differentiation and metabolism (15). Biological networks have provided a systems biology approach using data from DNA microarray, RNA-seq, miRNA and signaling pathways (16). Gene Ontology (GO) annotation and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were used to analyze biological networks to reveal potential synergism between gene and the organism (17). To achieve this, a miRNA-mRNA regulatory network was constructed by integrating DEMs, related miRNAs and major signaling pathways.

The present study results might provide further understanding of molecular mechanisms involved in NP progression and how differentially expressed mRNAs and miRNAs may serve as potential targets for NP therapy.

Materials and methods

Microarray database

The microarray data under access number GSE91396 (11) were downloaded from the National Center of Biotechnology Information (NCBI) Gene Expression Omnibus (GEO; www.ncbi.nlm.nih.gov/) (18). In GSE91396, RNA-seq samples were derived from three brain regions (NAc, mPFC and PAG) of animals two and a half months after sham surgery or SNI on the sciatic nerve. Bilateral punches were collected from the NAc, mPFC and PAG regions from 12 adult male C57BL/6 mice in the sham or SNI surgery group. The punches were pooled from two mice per sample. Thus, 36 samples in GSE91396 were processed, including 18 (6 mPFC, 6 NAc and 6 PAG) from sham surgery and 18 (6 mPFC, 6NAc and 6 PAG) from the SNI surgery group. The SNI model of neuropathic pain was performed under Avertine general anesthesia according to the study of Descalzi et al (11). Briefly, skin and muscle incisions were made on the left hind leg at mid-thigh level, revealing the sciatic nerve and its three branches. The common peroneal and sural nerves were carefully ligated with 6.0 silk sutures, transected, and 1-2 mm sections of each of these nerves were removed. The tibial nerve was left intact. Skin was then sutured with silk 4.0 sutures. Sham surgery mice underwent the same procedure, but all nerves were left intact.

Gene expression profiles in TXT format were downloaded from the NCBI database. The R3.4.1 preprocessCore (19) version 1.40.0 (bioconductor.org/packages/2.4/bioc/html/preprocessCore.html) was used to process the fragments per kilobase per million mapped reads (FPKM) values. The limma package (20) version 3.32.5 (bioconductor.org/packages/release/bioc/html/limma.html) was used to identify differentially expressed RNAs (mRNAs and miRNAs) in the three brain regions that are significantly different (P<0.05) and |log2 fold change (FC)| >0.585 were considered as thresholds. The numbers of miRNAs and genes were calculated. Subsequently, pheatmap (21) version 1.0.8 package (cran.r-project.org/web/packages/pheatmap/index.html) in R3.4.1 software were used for unsupervised hierarchical clustering analysis based on a correlation algorithm. Finally, the Database for Annotation, Visualization and Integrated Discovery (DAVID) (22,23) version 6.8 (david.ncifcrf. gov/) was used to perform GO (Biology Process, Cellular Component, and Molecular Function) and KEGG pathway enrichment analysis for the identified differentially expressed genes (DEGs).

Analysis of the common differentially expressed RNAs for all the three brain regions

Venn diagram (24) (cran.r-project. org/web/packages/VennDiagram/index.html) version 1.6.17 in R3.4.1 software was used to visualize differentially expressed RNAs. The common differentially expressed RNAs from all the three brain regions were selected for further analysis.

The weighed gene co-expression network analysis (WGCNA) algorithm provided topology properties for co-expression network analysis (25). This allows construction of scale-free networks, defines the co-expression matrix and adjacency function, calculates different node coefficients, and identifies functional modules associated with disease from high through put data (26). The WGCNA package (27) version 1.61 in R3.4.1 software (cran.r-project.org/web/packages/WGCNA/) was used to analyze the correlation with disease status for common differentially expressed RNAs. Genes in major modules were screened by GO functional and KEGG pathway enrichment analysis; P<0.05 was considered to indicate a statistically significant difference.

Construction and analysis of miRNA-mRNA regulatory network

To construct an miRNA-mRNA regulatory network and explore the correlations between miRNA and mRNA for mining for the potential roles of these RNAs in NP, the TargetScan Release7.1 (28) (www.targetscan.org/vert_71/) database was used to predict the target mRNAs of miRNAs. The miRNA-mRNA interactions involving common DEMs were selected. According to co-expression WGCNA results, the negative miRNA-mRNA interactions were used to construct the miRNA-mRNA regulatory network that was visualized using Cytoscape3.5.1 (29) (www.cytoscape.org/). GO annotation and KEGG pathway analysis was performed for target genes in the regulatory network.

miRNA-mRNA-pathway network construction

The Comparative Toxicogenomics Database (30) (CTD; ctd.mdibl. org/), a public website and research tool launched by Mount Desert Island Biological Laboratory, elucidates relationships between genesproteins, diseases, phenotypes, GO annotations, pathways, and interaction modules (31). CTD's primary objective is to advance understanding of how gene-environment interactions affect human health (32). 'Neuropathic pain' was used as the keyword to search the KEGG pathways related to NP disease in CTD. Pathways thus retrieved were compared with enriched pathways for target genes in the regulatory network. The overlapping pathways and related genes were hypothesized to play important roles in NP pathogenesis.

Results

Data preprocessing and screening of differentially expressed RNA

The transcriptional profiles downloaded from the GEO database were preprocessed and, after data normalization, 21,193 RNAs were identified, including 1,123 miRNAs and 20,070 mRNAs (Fig. 1A and B). RNAs with a zero median value were excluded, yielding 383 miRNAs and 17,654 mRNAs. A total of 2,776 differentially expressed RNAs were identified (219 miRNAs and 2,557 mRNAs) in SNI samples compared with sham surgery samples for further analysis (Table I). The clustering heatmaps showed significant differences in expression levels of differentially expressed RNAs between SNI and sham surgery groups pertaining to the three distinct brain regions (NAc, mPFC and PAG; Fig. 1C-E). These results show that clustering heatmaps of differentially expressed RNAs, identified from microarray datasets, may accurately distinguish samples from SNI and sham surgeries.

Figure 1

Data normalization and cluster analysis for differentially expressed RNAs. (A) Data before normalization. (B) Data after normalization. Volcano and clustering heatmaps for the differential expressed RNAs in three distinct regions of brain samples (C) NAc and (D) mPFC. Volcano and clustering heatmaps for the differential expressed RNAs in three distinct regions of brain samples (E) PAG. The horizontal axis in Volcano represents log2FC; vertical axis refers to −1*log10 (P-value). Red triangles and blue triangles represent the significantly upregulated and downregulated RNAs respectively in SNI samples. Black dots represent the RNAs that are not differentially expressed. Crosswise red dashed lines represent threshold of P=0.05 and vertical red dashed lines refer to the threshold of |log2FC|=0.585. FPKM, fragments per kilobase per million mapped reads; NAc, nucleus accumbens; mPFC, medial prefrontal cortex; PAG, periaqueductal gray; FC, fold-change; SNI, spared nerve injury.

Table I

Differentially expressed mRNAs and miRNAs in three distinct brain regions (NAc, mPFC, PAG) related to neuropathic pain.

Type	NAc		mPFC		PAG
	Down	Up	Down	Up	Down	Up
miRNA	50	50	49	56	7	7
mRNA	489	491	437	422	95	623
Total	539	541	486	478	102	630
	1,080		964		732

miRNA, microRNA; NAc, nucleus accumbens; PAG, periaqueductal grey; medial prefrontal cortex.

GO functional and KEGG pathway enrichment analysis for the DEGs

GO functional and KEGG pathway enrichment analyses were performed for the DEGs identified above (Fig. 2). The DEGs in NAc were mainly involved in DNA binding (GO:0003677, n=99), extracellular region part (GO:0044421, n=82), transcription regulator activity (GO:0030528, n=75), transcription factor activity (GO:0003700, n=64) and extracellular space (GO:0005615, n=57; Fig. 2A). The DEGs in mPFC were mainly involved in immune response (GO:0006955, n=41), transcription factor activity (GO:0003700, n=41), extracellular region part (GO:0044421, n=35), sequence-specific DNA binding (GO:0043565, n=32) and positive regulation of the macromolecule metabolic process (GO:0010604, n=32; Fig. 2B). The DEGs in PAG were mainly involved in DNA binding (GO:0003677, n=83), transcription regulator activity (GO:0030528, n=59), regulation of transcription from RNA polymerase II promoter (GO:0006357, n=34), protein kinase activity (GO:0004672, n=31) and structural molecule activity (GO:0005198, n=30; Fig. 2C).

Figure 2

GO functional and KEGG pathway enrichment analysis for the differentially expressed genes related to neuropathic pain. The pie chart slice sizes correspond to the number of annotated differentially expressed genes. GO functional enrichment analysis results for the differentially expressed genes in three distinct regions of brain (A) NAc, (B) mPFC and periaqueductal gray, (C) PAG. The yellow, pink and green colors represent Biological Process, Molecular Function and Cellular Component, respectively. (D) KEGG pathway enrichment analysis results for the differentially expressed genes in three distinct regions of brain (NAc, mPFC and PAG). The pink, yellow and green represent NAc, mPFC and PAG groups, respectively, and black dots represent -log2 (P-value). GO, gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; NAc, nucleus accumbens; mPFC, medial prefrontal cortex; PAG, periaqueductal gray.

In addition, nine, six and four significant pathways were identified for the DEGs identified in NAc, mPFC, and PAG, respectively (Fig. 2D). There were five overlapping pathways for NAc and mPFC, including cytokine-cytokine receptor interaction (mmu04060), chemokine signaling pathway (mmu04062), neuroactive ligand-receptor interaction (mmu04080), JAK-signal transducer and activator of transcription signaling pathway (mmu04630), and metabolism of xenobiotics by cytochrome P450 (mmu00980). The DEGs in PAG were mainly enriched in vascular endothelial growth factor (mmu04370) and Notch (mmu04330) signaling pathways.

Analysis of critical gene modules related to NP

A total of 2,325 common differentially expressed RNAs were identified in all the three brain regions (NAc, mPFC and PAG), including 173 miRNAs and 2,152 mRNAs. These common differentially expressed RNAs were analyzed using the WGCNA algorithm to identify disease-related RNAs. The topological matrix's scale-free distribution was calculated based on the GSE91396 dataset. Once the square value of the correlation coefficient reached 0.9 for the first time, the corresponding power value (power=7) was selected to calculate community dissimilarity of RNAs. After constructing the clustering dendrogram, the minimum gene number was set as 50 and the cut Height= 0.99. This process resulted in the identification of six modules-module-red, module-turquoise, module-blue, module-green, module-brown and module-yellow (Fig. 3).

Figure 3

Differentially expressed RNAs-identification and critical functional modules analysis. (A) Venn diagram of differentially expressed RNAs in three brain regions (NAc, mPFC and PAG). (B) The adjacency function definition for the RNAs: Left chart represents the power selection diagram of the adjacency matrix weight parameter. The horizontal axis represents the weight parameters power, while the vertical axis represents the square values of correlation coefficient between log (k) and log [p (k)]. The red line represents the standard as the square value reached 0.9. The right chart represents the mean connectivity of RNAs under a different adjacency matrix weight parameter. The red line shows that the average connection degree of the node is under differential values of power parameter. (C) Identification of major functional modules. The left chart refers to a cluster dendrogram based on the dynamic tree. Each dendrogram color represents a unique module. The right chart represents the heatmap for the correlation between each module and clinical factors. The number in the grid represents the correlation coefficient while the number in parentheses represents P-value. NAc, nucleus accumbens; mPFC, medial prefrontal cortex; PAG, periaqueductal gray.

Correlations between the RNA module and disease characteristics (different regions and disease states) were calculated and are shown in Table II. As shown, genes in module-red (correlation=0.33, P=4.0×10-66) and module-turquoise (correlation=0.82, P<0.001) were positively correlated with disease characteristics. Furthermore, the correlation values were significantly increased compared with the module-gray (correlation with disease=−0.10, P=2.0×10−07). Thus, the differentially expressed RNAs in the red and turquoise modules were used for further analysis.

Table II

Differentially expressed RNAs in gene modules related to neuropathic pain.

Module color	Correlation with disease	P-value	#RNA	#miRNA	#mRNA
Blue	−0.32	2.0×10−63	412	11	401
Brown	−0.79	<0.0001	209	11	198
Green	−0.05	6.0×10−3	87	0	87
Grey	−0.10	2.0×10−7	665	118	547
Red	0.33	4.0×10−66	71	2	69
Turquoise	0.82	<0.0001	781	31	750
Yellow	−0.31	7.0×10−59	100	0	100

Correlations between the RNA module and disease were calculated using weighed gene co-expression network analysis package version 1.61 in R3.4.1 software (cran.r-project.org/web/packages/WGCNA/).

Pathway enrichment analyses by GO functional and KEGG for DEGs in red and turquoise modules (Table III) showed 18 GO terms and seven pathways (Fig. 4). The DEGs in the two modules were mainly enriched in immune response (GO:0006955), defense response (GO:0006952) and response to wounding (GO:0009611). The major pathway categories for DEGs in the two modules were cytokine-cytokine receptor interaction (mmu04060) and neuroactive ligand-receptor interaction (mmu04080).

Table III

GO terms and KEGG pathways for the differential expressed mRNAs in module-red and module-turquoise.

A, Biology process
Term	Count	P-value
GO:0006955 immune response	33	6.36×10−5
GO:0007218 neuropeptide signaling pathway	11	2.01×10−4
GO:0009611 response to wounding	25	4.02×10−4
GO:0006952 defense response	29	7.03×10−4
GO:0006814 sodium ion transport	13	7.38×10−4
GO:0042742 defense response to bacterium	12	7.42×10−4
GO:0050954 sensory perception of mechanical stimulus	11	7.70×10−4
GO:0007605 sensory perception of sound	10	1.55×10−3
GO:0009617 response to bacterium	14	1.78×10−3
GO:0006954 inflammatory response	17	2.72×10−3
GO:0006811 ion transport	38	3.05×10−3
GO:0007601 visual perception	10	5.42×10−3
GO:0050953 sensory perception of light stimulus	10	5.79×10−3
GO:0002684 positive regulation of immune system process	15	7.24×10−3
GO:0007155 cell adhesion	30	9.03×10−3
GO:0022610 biological adhesion	30	9.42×10−3
GO:0007423 sensory organ development	17	9.56×10−3
GO:0030182 neuron differentiation	23	1.06×10−2

GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; ECM, extracellular matrix; TGF, transforming growth factor.

Figure 4

Functional and pathway enrichment analysis for differentially expressed mRNAs in turquoise and red modules. (A) GO functional and (B) KEGG pathway analyses of the differentially expressed mRNAs in turquoise and red modules. The horizontal axis refers to the number of mRNAs in functional analysis, while the vertical axis represents GO terms and KEGG pathway categories. The colors range from red to green, representing the change in P-values from low to high. GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes.

The mRNA-miRNA regulatory network construction

Target gene prediction of 33 miRNAs in red (2 miRNAs) and turquoise (31 miRNAs) modules using TargetScan Release 7.1 database resulted in overlapping target genes, which were matched with 819 mRNAs in red (69 mRNAs) and turquoise (750 mRNAs) modules, involving 558 pairs of negative miRNA-mRNA regulatory interactions. The correlation coefficients >0.6 for these relationships, based on the WGCNA algorithm, yielded 59 miRNA-mRNA interactions and 172 mRNA-mRNA interactions that facilitated construction of a miRNA-mRNA regulatory network (Fig. 5), which consisted of 58 nodes, including 10 miRNAs and 48 mRNAs.

Figure 5

Construction and analyses of miRNA-mRNA regulatory network. (A) The miRNA-mRNA regulatory network: The squares and circles represent miRNAs and mRNAs respectively, and the node colors from green to red indicate change from downregulated to upregulated in gene expression. The red 'T' shape represents a negative connection between miRNA and mRNA. The gray solid line refers to a positive connection. (B) GO analysis for the differential expressed mRNAs in miRNA-mRNA regulatory network. The horizontal axis refers to the number of mRNAs in functional analysis, while the vertical axis represents GO terms. The colors range from red to green, representing the change in P-values from low to high. (C) The pie chart of KEGG pathway analysis results: Each part represents one pathway category, while the number refers to the number of mRNAs. miRNA/miR, microRNA; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes.

The GO functional and KEGG pathway enrichment analysis showed that DEGs in the miRNA-mRNA regulatory network were enriched in 21 GO terms and five pathways. The GO terms included regulating T cell proliferation (GO:0042129) and T cell differentiation (GO:0030217). The major pathway was cytokine-cytokine receptor interaction (mmu04060; Table IV).

Table IV

GO terms and KEGG pathways for the critical genes in mRNA-microRNA regulatory network.

A, Biology process
Term	Count	P-value	Genes
GO:0042129 regulation of T cell proliferation	3	0.011	FOXJ1, IL12B, SPN
GO:0042981 regulation of apoptosis	6	0.016	COL18A1, PLEKHF1, GFRAL, MGMT,
			COL2A1, SPN
GO:0043067 regulation of programmed cell death	6	0.017	COL18A1, PLEKHF1, GFRAL, MGMT, COL2A1, SPN
GO:0010941 regulation of cell death	6	0.017	COL18A1, PLEKHF1, GFRAL, MGMT, COL2A1, SPN
GO:0030217 T cell differentiation	3	0.018	CD3D, IL12B, SPN
GO:0032944 regulation of mononuclear cell proliferation	3	0.020	FOXJ1, IL12B, SPN
GO:0050670 regulation of lymphocyte proliferation	3	0.020	FOXJ1, IL12B, SPN
GO:0070663 regulation of leukocyte proliferation	3	0.021	FOXJ1, IL12B, SPN
GO:0007601 visual perception	3	0.031	OPN5, PDE6G, CNGB3
GO:0050953 sensory perception of light stimulus	3	0.031	OPN5, PDE6G, CNGB3
GO:0050863 regulation of T cell activation	3	0.035	FOXJ1, IL12B, SPN
GO:0002820 negative regulation of adaptive immune response	2	0.035	FOXJ1, SPN
GO:0002704 negative regulation of leukocyte mediated immunity	2	0.037	FOXJ1, SPN
GO:0002707 negative regulation of lymphocyte mediated immunity	2	0.037	FOXJ1, SPN
GO:0030098 lymphocyte differentiation	3	0.038	CD3D, IL12B, SPN
GO:0042110 T cell activation	3	0.040	CD3D, IL12B, SPN
GO:0045061 thymic T cell selection	2	0.040	CD3D, SPN
GO:0042325 regulation of phosphorylation	4	0.044	DOK7, RAPGEF3, PDE6G, GRK1
GO:0002520 immune system development	4	0.046	CD3D, FOXJ1, IL12B, SPN
GO:0051174 regulation of phosphorus metabolic process	4	0.048	DOK7, RAPGEF3, PDE6G, GRK1
GO:0019220-regulation of phosphate metabolic process	4	0.048	DOK7, RAPGEF3, PDE6G, GRK1

Overlapping pathways identified by comparing pathways selected from Comparative Toxicogenomics Database with enriched KEGG pathways. GO, Gene Ontology; KEGG, Kyoto Encyclopedia of genes and genomes.

Crucial miRNAs, mRNAs and pathways related to NP

Keyword ('neuropathic pain') search of the CTD resulted in 88 pathways. Three pathways overlapped when the 88 pathways were compared to previously-identified enriched pathways for the DEGs in the mRNA-miRNA regulatory network. These pathways, cytokine-cytokine receptor interaction, chemokine signaling pathway and endocytosis were involved in four DEGs, including interleukin (IL)-8 receptor β (CXCR2), subunit β of interleukin (IL)-12 (IL12B), CD153 (also known as TNFSF8), and rhodopsin kinase (GRK1). The CXCR2 can be targeted by miR-7688-3p, IL12B by miR-208a-5p, GRK1 by miR-135b-3p and TNFSF8 by miR-344f-3p, miR-135b-3p, and miR-135a-2-3p, respectively (Fig. 6).

Figure 6

mRNA-miRNA-pathway regulatory network. The square, circle and rhombus shapes represent miRNA, mRNA, and pathway respectively. The variation of node colors from green to red indicates change from down-regulated to up-regulated in gene expression. The red 'T' shape represents a negative connection between miRNA and mRNA. The black arrow represents a connection between mRNA and the participating pathway. miRNA/miR, microRNA.

Discussion

In this study, 2,776 differentially expressed RNAs (219 miRNAs and 2,557 mRNAs) were identified in SNI compared with the sham surgery samples. In the three brain regions (NAc, mPFC and PAG), there were 2,325 common differentially expressed RNAs (173 miRNAs and 2,152 mRNAs). Two important modules (red and turquoise module) were identified as related to NP using WGCNA for the common differentially expressed RNAs. The miRNA-mRNA regulatory network was constructed based on the differentially expressed RNAs in the red and turquoise modules. The DEGs in miRNA-mRNA regulatory network were enriched in 21 GO terms and 5 pathways. Three pathways in the miRNA-target gene-pathway regulatory network, including cytokine-cytokine receptor interaction, chemokine signaling pathway and endocytosis, comprised four important DEGs (CXCR2, IL12B, TNFSF8, and GRK1) that were related to NP.

Cytokines are a group of small proteins (5-20 kDa) produced by a range of cells, including immune cells-macrophages, B lymphocytes, T lymphocytes, mast cells, as well as endothelial cells and various stromal cells (33). Cytokines have important roles in the immune system (34). Cytokines, along with chemokines, interferons (IFN), ILs, colony stimulation factor and tumor necrosis factor (TNF) are involved in responsiveness to trauma, pain, and infection (35). In the present study, three DEGs (CXCR2, IL12B and TNFSF8) participated in cytokine-cytokine receptor interaction. IL12B located on chromosome 5q31-33 encodes the p40 subunit of IL-12, an immunomodulatory cytokine (36). IL12B polymorphisms are associated with asthma and psoriasis (37). A recent study showed that significant decreases in systemic concentrations of chemokines, IL-12 and IFNγ, were observed in nerve-injured Foxp3+ regulatory T cell-depleted transgenic mice; decrease in IL-12 promoted pain hypersensitivity in this model (38). The immune system, particularly T cells, plays a key role in mediating NP. IL-12p70, also referred to as IL-12, is a pro-inflammatory cytokine secreted by activated hematopoietic phagocytic cells; Chen et al (39) investigated pain response following systemic administration of IL-12p70 and IL-12p40 homodimer and found that IL-12p40 exhibited significant anti-nociceptive effects in a rat model of chronic NP. The present study found that IL12B was dysregulated in brain tissues of NP mice, which is consistent with previous studies (37-39). This dysregulation may be effectuated by miRNA, one of which, miR-208a-5p, was found in the miRNA-target gene-pathway regulatory network, to target IL12B. Findings in the literature show miR-208a to regulate cardiac hypertrophy and conduction in mice (40), whereas the expression level of miR-208b progressively declined after spinal cord injury in humans (41). The results of the present study suggest miR-208a-5p may play a vital role in NP pathogenesis by regulating IL12B expression levels.

TNFSF8, a ligand of cluster of differentiation (CD)30 (or CD153), exhibits polymorphisms that showed significant associations with spondylarthritis in a French cohort (42). In addition, TNFSF8 is a susceptibility gene in excessive inflammatory responses (43). Since inflammation is a key pathophysiological process in NP (44), TNFSF8 might be a major regulator in NP inflammatory responses. The miRNA-target gene-pathway regulatory network revealed differentially expressed TNFSF8 could be targeted by three miRNAs (miR-344f-3p, miR-135b-3p and miR-135a-2-3p). It has been reported that miR-135a modulates inflammatory molecules, IL-6, IL-1β and TNF-α, which enhances inflammatory responses of vascular smooth muscle cells involved in vascular disease complications (45,46). Another study showed that overexpressing the miR-344b-1-3p inhibitor in alveolar macrophages significantly increased the expression of TNF-α, IL-1β and macrophage inflammatory protein (MIP)-2 (47). Together, the results of the present study indicate that differentially expressed IL12B targeted by miR-208a-5p and TNFSF8 targeted by miR-344f-3p, miR-135b-3p and miR-135a-2-3p, might play critical roles in NP through cytokine-cytokine receptor interaction.

CXCR2 and GRK1 are chemokines that were identified in this study to be involved in NP progression. CXCR2 (or IL 8 receptor β) is a chemokine receptor whose interaction with MIP elicits chronic neuroinflammation through neutrophil accumulation and hyperacetylation of histone H3 leading to NP (48).

CXCL1/CXCR2 signaling plays a vital role in pathological pain, including peripheral and central sensitization (49). CXCL1 sensitizes primary peripheral neurons by directly acting on CXCR2 (50). In the central nervous system, CXCL1/CXCR2 signaling increases N-methyl-D-aspartate receptor currents in neurons and promotes expression of genes related to neuroplasticity, contributing to prolonged chronic pain (51).

The enzyme phosphorylating rhodopsin receptor, rhodopsin kinase (GRK1), was identified in the late 1970s as controlling vision (52). Nerve injury, high norepinephrine concentration and abnormal β2-adrenoreceptor functions indicate high sympathetic nerve activity and dysfunction (53). Phosphorylation by GRK1/2, under high norepinephrine concentrations and high sympathetic nerve activity, induced β2-adrenoreceptor internalization by recruiting β-arrestin-1 to the receptor, leading adrenoreceptors to either recycle to the membrane or traffic to lysosomes for degradation (54). Taken together, it was speculated that CXCR2 targeted by miR-7688-3p and GRK1 targeted by miR-135b-3p might be involved in NP progression through chemokine signaling pathway and endocytosis.

The present study has limitations. The small sample size available for analysis may have obscured rare interactions. Additional validation of the roles of DEGs and miRNAs in NP progression is warranted.

In conclusion, 2,776 differentially expressed RNAs (219 miRNAs and 2,557 mRNAs) were identified in the NP SNI model. Differentially expressed IL12B, targeted by miR-208a-5p, as well as TNFSF8 dysregulated by miR-344f-3p, miR-135b-3p and miR-135a-2-3p, might play critical roles in NP through cytokine-cytokine receptor interaction. CXCR2 targeted by miR-7688-3p and GRK1 targeted by miR-135b-3p might be involved in NP progression through chemokine signaling pathway and endocytosis. The present study has provided new insights into the regulatory mechanisms of NP, which merits their potential as candidates for therapies for this debilitating condition.