key: cord-1048111-qkdq1h4k
authors: Zhang, Huai-Gen; Liu, Li; Song, Zhi-Ping; Zhang, Da-Ying
title: Bioinformatics Analysis of the MicroRNA-Metabolic Gene Regulatory Network in Neuropathic Pain and Prediction of Corresponding Potential Therapeutics
date: 2021-09-27
journal: J Mol Neurosci
DOI: 10.1007/s12031-021-01911-w
sha: b9b2a4725b255b250677e850bc690938e35a83d3
doc_id: 1048111
cord_uid: qkdq1h4k

Neuropathic pain (NP) involves metabolic processes that are regulated by metabolic genes and their non-coding regulator genes such as microRNAs (miRNAs). Here, we aimed at exploring the key miRNA signatures regulating metabolic genes involved in NP pathogenesis. We downloaded NP-related data from public databases and identified differentially expressed microRNAs (miRNAs) and mRNAs through differential gene expression analysis. The miRNA target prediction was performed, and integration with the differentially expressed metabolic genes (DEMGs) was used for constructing the miRNA-DEMG network. Subsequently, functional enrichment analysis and protein–protein interaction (PPI) analysis were performed to explore the role of DEMGs in the regulatory network. The drug prediction was performed based on the DEMGs in the miRNA-DEMG network. A total of 8251 differentially expressed mRNAs (4193 upregulated and 4058 downregulated), and 959 differentially expressed miRNAs (455 upregulated and 504 downregulated) were identified. Moreover, after target gene prediction, a miRNA-DEMG network composed of 22 miRNAs and 113 mRNAs was constructed. The network was constituted of 135 nodes and 236 edges. We found that DEMGs in the network were mainly enriched in metabolic pathways and metabolic processes. A total of 1200 drugs were predicted as potential therapeutics for NP based on the differentially expressed genes, while 170 drugs were predicted for the DEMGs in the miRNA-DEMG network. Conclusively, our study predicted drugs that may be effective against the metabolic changes induced by miRNA dysregulation in NP. This information will help further reveal the pathological mechanism of NP and provide more treatment options for NP patients. SUPPLEMENTARY INFORMATION: The online version contains supplementary material available at 10.1007/s12031-021-01911-w.

Neuropathic pain, the main form of chronic pain, is caused by damage to the nervous system and dysfunction (Finnerup et al. 2016; Gaskin and Richard 2012) . Currently, the treatment of NP is limited to symptomatic treatment, and its prognosis is poor. To better understand the pathogenesis of NP, it is essential to develop effective prevention strategies and improve the efficacy of NP treatment. This necessitates the deepening in the comprehension of the molecular pathogenesis of NP.

miRNAs are transcripts of about 20-25 nucleotide (nt)long; they regulate a variety of processes, including DNA methylation, transcription, and post-transcriptional RNA processing. The miRNAs participate in various gene transcriptional processes via interacting with transcription factors, coactivators, and/or inhibitors. For example, miRNAs can bind mRNAs and affect the expression of downstream pathways. Accumulating evidences show that miRNAs play an important role in the pathogenesis of nervous system diseases (Liu et al. 2016; Wei et al. 2018) . For instance, Zhou et al. (2017) found that a total of 134 lncRNAs, 12 miR-NAs, 188 circRNAs, and 1066 mRNAs were significantly regulated after spared nerve injury (SNI) surgery. Wei et al. (2018) found that miR-154-5p inhibition promotes NP progression in rats via upregulation of TLR5S and proposed this pathway as a novel therapeutic target for NP treatment. It is known that miRNAs regulate gene expression and can be used as potential biomarkers. Dayer and colleagues (Dayer et al. 2019) reported that chronic pain may induce the abnormal and specific dysregulation of miRNA expression and that hsa-miR-320a and hsa-miR-98-5p (two circulating miRNA signatures) can serve as biomarkers for pain-type classification. The epigenetic intervention of miRNAs may also be a new therapeutic approach for complications such as injurious hypersensitivity caused by peripheral nerve injury. Pan and colleagues (Pan et al. 2018) reported that miR-23a may be involved in NP through the TXNIP/NLRP3 inflammasome axis in spinal glial cells by directly targeting CXCR4.

Metabolism plays a major role in multiple biological processes. Studies state that metabolic dysregulation plays a central role in painful peripheral neuropathies, and there is cross-communication between metabolism, inflammation, and immune response in NP (Navia-Pelaez et al. 2021) . Previous studies have shown that CCI-induced NP causes metabolic changes in serum and spinal cord (Chen et al. 2021) ; changes in metabolism of thalamic neurotransmitters have also been reported ). However, a systemic analysis showing the overall profile of genes involved in metabolism has not yet been reported in PN. In addition, although many studies indicate the importance of microR-NAs in gene regulation, the interaction between metabolic genes and miRNAs has not yet been elucidated, which constitutes an important field to explore.

With the rapid development of high-throughput sequencing, data mining, and the wide application of precision medicine, it becomes more feasible to extract miRNA and mRNA information of NP from the microarray datasets. Previous studies have reported the key genes and pathways associated with NP based on expression patterns Zhu et al. 2019 ). The contributors of GSE24982, Von and colleagues, conducted a miRNA and mRNA expression profiling study of dorsal root ganglion (DRG) tissue from rats with spinal nerve ligation (SNL), and they found that the expression level of 63 miRNAs was significant by t-test (P-value < 0.001) (von Schack et al. 2011). GSE24982 was analyzed in other studies as well: Gao and colleagues obtained 123 upregulated differentially expressed genes and identified p53 as a candidate prognostic biomarker for NP (Gao et al. 2018 ); a total of 206 differentially expressed genes was obtained in GSE24982 by unpaired t-test with a fold change ≥ 2, and the genes were enriched in DNA binding, cell cycle ); Zhu and colleagues analyzed the differentially expressed genes of high and low pain compared with the normal control group and obtained 1243 upregulated and 1533 downregulated genes in GSE24982 by R "limma" package with a cutoff of |log2FC|> 1 and P-value < 0.05 (Zhu et al. 2019) . A study reported a miRNA expression profiling study in DRG after peripheral nerve injury and found that a total of 220 miRNAs were upregulated in DRGs following three types of peripheral nerve injury (Chang et al. 2017) . Although previous studies have investigated the role of the differentially expressed miRNAs and mRNAs in the NP, the interaction and regulatory network of miRNAs and metabolic genes in NP have not been elucidated.

In this study, two datasets (GSE145199 and GSE24982) were downloaded to identify differentially expressed mRNAs and miRNAs. Then, we performed GO and KEGG pathway analysis of differentially expressed mRNAs and constructed a PPI network to investigate gene interactions. Moreover, the metabolic process related genes were downloaded from the GSEA database and used to identify the differentially expressed metabolic genes (DEMGs). Then, the miRNA-DEMG network was constructed to explore the interaction of miRNAs and DEMGs in NP. Finally, drugs were predicted based on the DEMGs in the regulatory network. Our study will contribute to the understanding of NP pathogenesis and provide new strategies for targeted drugs and metabolic therapies.

Raw data were obtained from Gene Expression Omnibus (GEO, https:// www. ncbi. nlm. nih. gov/ geo/) with the following search strategy: "Neuropathic Pain" AND ("miRNA" OR "mRNA"). Eligible datasets met the following inclusion criteria: (1) the dataset contains NP samples and control samples; (2) each sample has a group label; (3) the microarray data type has been specified; (4) gene symbol or GenBank ID in the dataset file for each probe has been provided; (5) the original data is available. Finally, the included data accessions were GSE145199 and GSE24982. GSE145199 is a miRNA expression profiling of prelimbic cortex (PL), including 3 SNI samples and 3 control samples (Rattus norvegicus). GSE24982 is an mRNA expression profiling in DRG that underwent spinal nerve ligation (SNL) or a sham operation, including 20 SNL samples and 20 SHAM samples (R. norvegicus). In this study, all the samples from the nerve injury were compared with the SHAM samples or control samples in the abovementioned datasets.

Before differential expression analysis, average expression of a gene was kept in the expression matrix when its gene symbol mapped multiple probes, and genes containing the missing value (or zero value) were removed. Data preprocessing and differentially expressed analyses were performed with the R limma package (Ritchie et al. 2015) . Data normalization was implemented by the quantile method with the R limma package and data scaling was performed by logarithmic conversion in R. Differentially expressed miRNAs and mRNAs were identified with the significance thresholds of |log2FC|> 0.5 and P-value < 0.05. The R ggplot2 package (Ginestet C 2011) was used for data visualization. The ten most significant genes (sorted by P-values) of the upregulated cluster and downregulated cluster were extracted to generate a heatmap using the R pheatmap package (https:// cran.rproje ct. org/ web/ packa ges/ pheat map/ index. html). To identify the DEMGs, we first downloaded metabolic process-related gene sets in the GSEA database and combined them in a gene list. This gene list containing metabolic genes and the list of differentially expressed mRNAs were used for Venn diagram analysis. The intersection between both gene lists was considered the list of differentially expressed metabolic genes (DEMGs). The Venn diagram analysis was performed online at the Bioinformatics & Evolutionary Genomics domain (http:// bioin forma tics. psb. ugent. be/ webto ols/ Venn/).

The R multiMiR package (Ru et al. 2014 ) was used to predict target mRNAs regulated by differentially expressed miRNAs. Then, the Pearson correlation analysis was performed to evaluate the correlation between the expression of miRNAs and mRNAs. Target mRNAs with negative correlation (a negative correlation coefficient lower than − 0.8) with miRNA and overlapping with the DEMGs were used for network construction. The miRNA-DEMG network was visualized by Cytoscape software (Shannon et al. 2003) , and the MCODE plugin in Cytoscape was used to identify hub miRNAs and DEMGs.

Search Tool for the Retrieval of Interacting Genes (STRING, version 11.0, https:// string-db. org/) is a database of known and predicted protein-protein interactions (Szklarczyk et al. 2019) . Currently, the STRING database contains 24,584,628 proteins from 5090 organisms, and interactions in the database are derived from genomic context predictions, high-throughput lab experiments, automated text mining, co-expression, and previous knowledge in databases. The identified miRNA target DEMGs and differentially expressed DEMGs were input into the STRING database to construct the PPI network, in which individual non-connected proteins were eliminated.

Functional enrichment analysis of sets of genes was performed using the online platform g:Profiler (https:// biit. cs. ut. ee/ gprofi ler/ gost). R. norvegicus was chosen as organism, and g:SCS threshold of 0.05 was chosen as significance cutoff threshold for the enrichment terms. The results of the top ten terms were used to draw the bubble chart using a customized script in R software with the ggplot2 library.

The Drug Gene Interaction Database (DGIdb, www. dgidb. org) has been used to integrate, organize and display drugs, gene interactions, and gene-drug information from published articles and web resources (Cotto et al. 2018) . It is easily accessed through an intuitive Web user interface, an application programming interface (API), and publicly cloud-based server images. Here, we uploaded the miRNA target DEMGs in the miRNA-DEMG network, differentially expressed mRNAs, and hub genes in the protein-protein network of the differentially expressed mRNAs onto the DGIdb to predict their potential targeting drugs effective for NP. Moreover, the L1000CDS 2 pharmacogenetic search tool (https:// maaya nlab. cloud/ l1000 cds2/#/ index) was applied to DEMGs to further retrieve potential perturbant drugs targeting these genes.

After data preprocessing, the raw data of each dataset was normalized. The boxplots in Fig. 1 A and B show the difference of samples before and after data normalization. The volcano plots of miRNA and mRNA expression profiles in GSE145199 and GSE24982 were as depicted in Fig. 2 A and B. A total of 959 differentially expressed miRNAs were identified in GSE145199, of which 455 miRNAs were upregulated, and 504 were downregulated (Supplementary File S1). In GSE24982, we identified 8251 differentially expressed mRNAs (4193 upregulated and 4058 downregulated mRNAs) (Supplementary File S2). The heatmaps in Fig. 2 C and D displayed the top ten differentially expressed mRNAs sorted by P-values in the upregulated and downregulated clusters in GSE24982.

Functional analysis has a wide range of applications in bioinformatics and can explain biological mechanisms and functional pathways in genomics and transcriptomics. In order to identify the functional importance of differentially expressed mRNAs, functional enrichment analysis of the 8251 differentially expressed mRNAs was performed. We found that the downregulated differentially expressed mRNAs were significantly enriched in 661 GO terms in the category of biological process (GO: BP) (Supplementary File S3). Based on adjusted P-values, the most significantly enriched biological processes were localization, regulation of biological quality, and positive regulation of biological process while on the basis of gene ratio, the largest sets of genes were involved in cellular process, biological regulation, and metabolism-related biological processes such as metabolic process, organic substance metabolic process, cellular metabolic process, and primary metabolic process ( Fig. 3 ; Supplementary File S3). Protein binding, binding, ion binding, identical protein binding, small molecule binding, and catalytic activity were the most enriched in the category of molecular function ( Fig. 3 ; Supplementary File S3), whereas cytoplasm, intracellular anatomical structure, intracellular membrane-bounded organelle, membrane-bounded organelle, organelle, intracellular organelle, endomembrane system, cell junction, and synapse were the terms mostly enriched in the category of cellular component ( Fig. 3 ; Supplementary File S3). Additionally, we found that the downregulated differentially expressed mRNAs mainly participated in metabolic pathways, glutamatergic synapse, MAPK signaling pathway, calcium signaling pathway, cocaine addiction, neuroactive ligand-receptor interaction, PI3K-Akt signaling pathway, and cAMP signaling pathway in the KEGG pathway analysis ( Fig. 3 ; Supplementary File S3). More interestingly, in the WikiPathways (WP) analysis, spinal cord injury was found as the most significantly enriched pathway (Supplementary File S3). For the upregulated differentially expressed mRNAs, we found that the biological process terms of localization, positive regulation of biological process, transport, developmental process, and multicellular organism development were the most enriched ( Fig. 4 ; Supplementary File S4). For molecular function category, the terms mostly enriched for the upregulated genes were protein binding, binding, identical protein binding, and enzyme binding ( Fig. 4 ; Supplementary File S4), whereas cytoplasm and intracellular membrane-bounded organelle were the most enriched terms in the category of cellular component ( Fig. 4; Supplementary File S4) . The neuroactive ligand-receptor interaction, "growth hormone synthesis, secretion and action", calcium signaling pathway, proteoglycans in cancer, MAPK signaling pathway, prolactin signaling pathway, cholinergic synapse, and GABAergic synapse were the most enriched pathways ( Fig. 4 ; Supplementary File S4) obtained from the KEGG database, while VEGF-receptor signal transduction and IL-2 signaling pathway were the most significantly enriched in the WP database (Supplementary File S4). 

The PPI network based on the top 2000 differentially expressed mRNAs indicated strong interactions between these genes. The PPI network contained 1122 nodes and 6012 edges, and the average number of neighbors was 11.399 (Supplementary Figure S1 ). Using MCODE, we identified a significant cluster with the highest score of 45.43 containing 79 genes as hub genes for the PPI; this representative subnetwork was as depicted in Fig. 5 and contained 72 nodes and 1772 edges. Functional enrichment analysis indicated that the hub genes were primarily implicated in the biological processes of G protein-coupled receptor signaling pathway, signal transduction, signaling, cell communication, regulation of cytosolic calcium ion concentration, and cellular calcium ion homeostasis ( Fig. 6 ; Supplementary File S5). The predominant molecular functions of the hub genes were G protein-coupled receptor activity, transmembrane signaling receptor activity, and molecular transducer activity ( Fig. 6 ; Supplementary File S5), while the most enriched cellular component terms were integral component of membrane, intrinsic component of membrane, plasma membrane, and cell periphery ( Fig. 6; Supplementary File S5) .

Furthermore, in the KEGG database, we identified neuroactive ligand-receptor interaction, calcium signaling pathway, and chemokine signaling pathway as the most significantly enriched pathways for the hub genes ( Fig. 6 ; Supplementary File S5).

To identify the miRNA-mRNA interactor pairs, the Pearson correlation between the differentially expressed miR-NAs and the differentially expressed mRNAs was analyzed, and miRNA-mRNA pairs (24 miRNAs and 305 mRNAs) showing negative correlation coefficient lower than − 0.8 (Supplementary File S6) were chosen for further analysis. Next, the 24 miRNAs were used for miRNA-target prediction with the multimir package in R based on R norvegicus, Mus musculus and Homo sapiens as species. The prediction results were summarized in Supplementary File S7, and subsequent intersection analysis indicated that the 305 differentially expressed mRNAs were targets for the 24 miRNAs. As shown in Supplementary Figure S2 , the intersection between the 304 mRNAs negatively correlated with the 24 miR-NAs, the metabolic process gene list and predicted miRNA Figure S3) .

The miRNA-DEMG pairs were used for miRNA-DEMG network visualization in Cytoscape. The miRNA-DEMG network showing miRNA-DEMG interactions and PPI interactions of the DEMGs in the network was depicted in Fig. 7 , whereas the network parameters were summarized in Supplementary File S8. The network was constituted of 135 nodes and 236 edges (Fig. 7) . The average number of neighbors was 3.496, while the network diameter and radius were nine and five, respectively (Supplementary File S8). MCODE analysis indicated miR-3613-3p, Polr2d, Nudt21, Sf3b3, Ddx6, Xrn1, Cks2, Mcm4, Mcm5, Mcm6, Zfp36, Oxa1l, Mrpl17, and ENSRNOG00000047563 as hub genes in the network (colored in orange) (Fig. 7) . The top miRNAs with the highest number of interactors were miR-940, miR-5192, miR-2277-3p, miR-6882-3p, miR-5698, and miR-6133. To explore the biological function of DEMGs in the miRNA-DEMG network, we performed functional enrichment analysis ( Fig. 8 ; Supplementary File S9). We found that the DEMGs in the network were enriched in the biological processes of cellular metabolic process, metabolic process, and regulation of cellular protein metabolic process (Fig. 8, Supplementary File S9) . For molecular function category, organic cyclic compound binding, RNA binding, and carboxylic acid binding were the most enriched terms ( Fig. 8; Supplementary File S9) , 

The DGIdb database provides the association between genes and their matching known or potential drugs. All the differentially expressed mRNAs were used for drug prediction via the DGIdb database, and the results were presented in Supplementary File S10. A total of 1200 drugs were predicted. Among these drugs, fostamatinib had the highest number (66) of target genes. Other predicted drugs included cisplatin (37 target genes), copper (28 target genes), paclitaxel (26 target genes), dasatinib (24 target genes), erlotinib (23 target genes), artenimol (22 target genes), sorafenib (22 target genes), fluorouracil (21 target genes), gemcitabine (21 target genes), alcohol (18 target genes), carboplatin (18 target genes), gefitinib (18 target genes), docetaxel (17 target genes), imatinib (17 target genes), methotrexate (15 target genes), dexamethasone (14 target genes), everolimus (14 target genes), methyldopa (14 target genes), palbociclib (14 target genes), progesterone(14 target genes), and tamoxifen (14 target genes) (Supplementary File S10). For the 79 hub genes in the PPI network, drugs could be predicted for 60 of them, and For the DEMGs targeted by miRNAs, 170 drugs were predicted from the DGIdb database, among which each of Adriamycin, bevacizumab, carboplatin, cisplatin genistein, imatinib, irinotecan, pioglitazone, and triamcinolone were found to have two targets in the set of DEMGs in the miRNA-DEMG network (Supplementary File S12). All of the other drugs had only one target as DEMG (Supplementary File S12). Furthermore, the characteristic direction signature search engine L1000CDS 2 was applied to the 110 DEMGs for drug repurposing. Depending on the expression trend (upregulated or downregulated) of the 110 DEMGs in the network, compounds or drugs that potentially reversed the gene signatures were retrieved. The heatmap of potential drugs and their modulated signatures were shown in Fig. 9 . Upregulated consensus DEMGs were downregulated by these compounds, while downregulated consensus DEMGs were upregulated after stimulation by the compounds. Rottlerin, BRD-K31912990, and BRD-K08448573 were the top-3 drugs obtained based on the overlap score values. Rottlerin downregulated nine upregulated DEMGs (CKS2, FBLN5, MCM4, MCM5, MCM6, NUP210, RPS6KA3, SFPQ, and SHMT1) among the input DEMG list. It is worth noting that Cks2, Mcm4, Mcm5, and Mcm6 are hub genes in the miRNA-DEMGs network (Fig. 5) ; these genes were targeted by the vast majority of the predicted drugs (Fig. 9) . Among the downregulated DEMGs, DUSP1 and GADD45A were both activated by 23 of the predicted drugs, followed by DUSP6 (16 drugs) (Fig. 9 ). 

Neuropathic pain is a debilitating condition caused by damage to the nervous system and chronic disease. As a kind of non-coding RNA, miRNA-related research has become a hot spot in genetic research. As short non-coding single-stranded RNA molecules, miRNAs participate in the regulation of post-transcriptional gene expression in animals and plants (Singh et al. 2020; Yanai et al. 2020) and are involved in many cellular activities such as epigenetic regulation, cell cycle regulation, cell differentiation, and metabolic processes. At present, many miRNAs have been reported in neurodegenerative diseases (Eyileten et al. 2020) , cancers (Kapoor et al. 2020; Wilczynski et al. 2020) , bowel disease (Ambrozkiewicz et al. 2020) , and diabetic heart disease (DHD) (Lew et al. 2020) . Although some miRNAs (Kalpachidou et al. 2020; Tang et al. 2020) have been reported in NP, reports on how miRNAs and mRNAs jointly regulate the pathogenesis of NP are still limited and have only captured partial insights of the complex pathogenetic mechanisms. In addition, studies have indicated the shift in metabolic processes in NP and nervous system-related diseases, but, to the best of our knowledge, no study has explored the regulatory relationship between miRNAs and metabolic genes, which would be important for understanding the pathogenesis of NP and developing essential metabolic therapies for this condition. Thus, in the present study, we identified the miRNA signature regulating metabolic genes in NP and constructed the corresponding miRNA-DEMG network. Furthermore, in order, to uncover drugs that could revert the deleterious effects induced by the miRNAs, drug-gene interactions and drug repositioning analysis were performed and a list of potential drugs were predicted. This study is the first of its kind to establish the miRNA-DEMG network and predict the potential drugs for this network in NP.

At present, the pathological mechanism of NP is still unknown and might be caused by differentially expressed mRNAs. Therefore, we screened differentially expressed mRNAs in the GSE24982 dataset and functional enrichment analysis showed that the main pathways for these differentially expressed mRNAs in NP were neuroactive ligand-receptor interaction, PI3K-Akt signaling pathway, MAPK signaling pathway, cAMP signaling pathway, and calcium signaling pathway. A previous study reported that jct-801 relieves paclitaxel-induced neuropathic pain through the PI3K/Akt pathway . Our study suggested that differentially expressed mRNAs are involved in NP possibly through regulating the MAPK signaling pathway, which is consistent with a previous study (Galan-Arriero et al. 2015) . The cyclic adenosine monophosphate (cAMP) signaling pathway is a key contributor to the development of chronic pain, and an existing study indicated that knockdown of the cAMP effector can relieve the pain-like responses in chronic pain models (Berkey et al. 2020 ). In addition, the PPI network analysis of differentially expressed mRNAs indicated that proteins encoded by these genes were involved in strong interactions, which indicated that these genes work in concert to occasion the deleterious effects associated with NP. A robust cluster of 79 hub genes were identified and was found to be mainly involved in the biological process of G protein-coupled receptor signaling pathway. The finding corroborated previous studies indicating that the G protein-coupled receptors are involved in pain transmission (Pan et al. 2008) and that targeting these molecules may be vital for the relief of pain . The hub genes were also involved in other processes and pathways associated with signal transmission such as neuroactive ligand-receptor interaction, calcium signaling pathway, and chemokine signaling pathway, suggesting that targeting these hub genes may be vital for treating NP. Therefore, our study provided a supplementary contribution for the molecular understanding of the NP pathological mechanism.

miRNAs and mRNAs often coordinate the occurrence and prognosis of diseases in organisms. Up to date, only few studies have analyzed the mRNA-miRNA network in NP (Zhou et al. 2017) . Various studies have provided new mechanisms for the molecular roles of miRNAs in the pathogenesis of NP (Chang et al. 2020; Phạm et al. 2020; Wilkerson et al. 2020 ). However, whether there is a regulatory miRNA signature controlling the metabolic genes in NP is unknown. Herein, we identified a 24-miRNA signature regulating DEMGs in NP and constructed the corresponding miRNA-DEMG network. These miRNAs have not been reported in metabolic processes associated with NP in the past. For the first time, we reported that miR-940 and miR-2277-3p may play an important role in metabolic changes associated with NP pathology based on the miRNA-DEMG regulatory network of NP; these miRNAs were previously reported mainly in cancer research (Fan et al. 2018; Gao et al. 2020; Zhou et al. 2019) . Additionally, miR-5689, identified in the miRNA-DEMG regulatory network, was also reported as a biomarker for predicting the development of new distant metastasis (Satomi-Tsushita et al. 2019) . MCODE analysis pinpointed miR-3613-3p as a hub miRNA in the miRNA-DEMG network; this finding corroborated the findings of previous researchers putatively showing that miR-3613-3p is a regulator of pain-associated genes such as GABRB3, NPY1R, GRIN3A, TRPV1, and SCN9A (Linnstaedt et al. 2015) . The hub DEMGs in the miRNA-DEMG regulatory network were Polr2d, Nudt21, Sf3b3, Ddx6, Xrn1, Cks2, Mcm4, Mcm5, Mcm6, Zfp36, Oxa1l, Mrpl17, and ENSRNOG00000047563 . However, the functions of these genes have not been systematically reported in NP. Thus, we anticipate that these genes might be new potential therapeutic targets for NP. Further experimental studies are required to confirm their role in NP and explore whether their modulation might be beneficial for the treatment of NP. These genes might also be essential for the development of metabolic therapies for NP. In short, we have discovered for the first time some miRNAs and mRNAs that may be involved in NP, but their specific function in NP needs more validation studies.

To further explore the biological function of the miRNA-DEMG network, we performed the GO and KEGG pathway analysis. We found that the DEMGs in the network were involved in the pathways of RNA degradation and cell cycle. These results indicated that the miRNA-DEMG network may participate in the pathogenesis of NP by regulating these pathways. Studies have reported that the cell cycle is activated in spinal cord injury (SCI), which may regulate chronic pain after SCI (Wu et al. 2013 (Wu et al. , 2016 . RNA degradation plays significant roles in neurodegenerative diseases as previously reviewed (Weskamp and Barmada 2018) . However, the implication of RNA degradation in NP is not-well documented; the enrichment of the DEMGs in the miRNA-mRNA network in this pathway indicated that RNA degradation is highly active in NP and this pathway could be a potential therapeutic target for NP.

Up to date, there is still no effective treatment drug for NP. As indicated above, significant dysregulation of miRNAs and mRNAs is associated with NP. Thus, finding drugs targeting these genes may help correct the deleterious effect encountered in NP. For this reason, we planned to uncover drugs targeting the differentially expressed mRNAs and DEMGs in the miRNA-DEMG network using the DGIdb. DGIdb database can infer the targeted drugs of genes in diseases based on existing resources. The drugs included in this database contain the US Food and Drug Administration (FDA) certified drugs, which will provide more professional and reliable assurance for target drug prediction. For example, Nambou and Anakpa (2020) recently discovered nicotinamide adenine dinucleotide (NAD) and CHEMBL1161866 with potential therapeutic value for coronavirus disease 2019 (COVID-19) through the DGIdb database, which provided novel insight for the treatment of COVID-19. In addition, Yang and colleagues predicted candidate drugs for lung adenocarcinoma (LUAD) via the DGIdb database as well ). Here, we predicted a series of drugs that could be potentially used for NP treatment based on differentially expressed mRNAs and DEMGs in the regulatory networks of NP. It is reported that inhibition of spleen tyrosine kinase (Syk) can alleviate mechanical allodynia (Choi and Yang 2019) . Fostamatinib is an oral Syk inhibitor that has been approved by the FDA for the adult treatment of chronic immune thrombocytopenia (ITP) (Bussel et al. 2018) . Although there are differences between the pathological processes of mechanical pain and neuropathic pain, the number of genes targeted by this drug is high (66 target genes), which deserves more clinical studies. Dexamethasone is a synthetic glucocorticoid with anti-inflammatory activity and minimal glucocorticoid effect, which is widely used in the treatment of various inflammation disorders (Coutinho and Chapman 2011) . It has been reported that low doses of ibuprofen and dexamethasone have a synergistic therapeutic effect on trigeminal NP in rats and can significantly inhibit mechanical allodynia (Park et al. 2019) . Dexamethasone is effective in the treatment of NP caused by tumor-related spinal cord compression, and pregabalin is effective for malignant painful radiculopathy (Tagami et al. 2020) . Nevertheless, the long-term use of dexamethasone also has certain side effects (Lieberthal et al. 2015) . In the past, opioids such as morphine were often ineffective in the treatment of neuropathic pain, but a study showed that the combination of imatinib and morphine can completely relieve pain (Donica et al. 2014 ), which provides more treatment options for NP patients. A recent study reported that in rats of infraorbital nerve ligation trigeminal neuralgia, c-Abl expression was significantly increased, and the downstream activation product p38 was also abnormally activated. Interestingly, imatinib mesylate (STI571), a specific c-Abl family kinase inhibitor, can reduce the expression of p38 and reduce the loss of dopaminergic neurons, suggesting that imatinib may alleviate the symptoms of NP through the c-abl-p38 signaling pathway (Fu et al. 2020 ). Additionally, our study predicted the treatment ability of progesterone (PROG) in NP and was supported by a previous study, demonstrating that PROG may provide new strategies for the treatment of NP (Jarahi et al. 2014) . NP caused by chemotherapy can reduce the prognosis of patients with the quality of life. Tamoxifen can alleviate NP induced by paclitaxel, vincristine, and bortezomib in chemotherapy through inhibiting protein kinase C/extracellular signal-regulated kinase pathway (Tsubaki et al. 2018 ). Although we have only discussed the therapeutic effects of the drugs above, we cannot ignore the potential effects of other drugs that have not been mentioned. These predicted drugs can be verified in future studies, either individually or in combination. The number of genes corresponding to the targeting drugs was small, indicating that the therapeutic effect of these drugs may still be relatively weak. The drug repurposing for DEMGs in the miRNA-DEMG network allowed the identification of drugs that could reverse the expression levels of these genes; these drugs are potentially those that can possibly correct the abnormal expression induced by miRNAs in the pathogenesis of NP. However, all the abovementioned therapeutic effects of drugs in NP were based on bioinformatic analysis and literature reports. The specific therapeutic effects of drugs still need more experimental and clinical studies for validation.

Although this study found several miRNAs and mRNAs that may be related to NP through the public datasets, and analyzed the pathways that may participate in the pathological mechanism of NP through GO and KEGG pathway analysis, our research still has some limitations. Circular RNAs (circR-NAs) and lncRNAs are non-coding RNA molecules that have recently been reported to regulate the function of miRNAs and affect a variety of complex human diseases (Fan et al. 2020; Harrison et al. 2020; Li et al. 2020b) . Although there have been some studies on the circRNA regulation in NP (Li et al. 2020a; Wei et al. 2020; Zhang et al. 2020) , the present public databases still do not include the circRNA expression data of NP. In addition, the annotation of the available lncRNA data is not complete, and the prediction of lncRNA target miR-NAs and mRNAs is critical. Therefore, we used microarray sequencing data, including the sequencing results of miRNA and mRNA in NP, which may not fully reflect the pathological mechanism of NP. This study found some miRNAs and mRNAs that may be related to NP, but their relationship with NP is still unclear. We predicted several drugs that may have therapeutic effects on NP, but the effects of these drugs in clinical treatment have not yet been reported. It should be noted that the information of miRNA and mRNA is taken from the datasets obtained by different samples, which may affect the reliability of the results. In the future, we will collect samples from human NP patients, and obtain lncRNA, miRNA, mRNA, and circRNA data of NP patients through high-throughput sequencing technology. We plan to verify the functions of the genes and drugs identified here through more in vitro and in vivo experiments and strive to further reveal the pathological mechanism of NP.

Our study predicted drugs to treat NP based on the NP regulatory miRNA-DEMG network. Some of the predicted drugs have been reported to alleviate NP in previous studies, but the role of other drugs in NP therapy remains unknown. Future research will require more studies to validate these drugs. In vitro and in vivo experiments based on cell and animal models of NP, and clinical trials are necessary to validate these drugs.

The online version contains supplementary material available at https:// doi. org/ 10. 1007/ s12031-021-01911-w.

Author Contribution HZ analyzed the data, prepared figures, and wrote the manuscript. LL and ZS analyzed the data and performed the experiments. DZ conceived and designed the study. All authors have read and approved the final version of the manuscript.

All data generated or analyzed during this study are included in this published article. Data is available at NCBI GEO: GSE145199 and GSE24982.

Ethics Approval and Consent to Participate Not applicable.

Competing Interests The authors declare no competing interests.

In search for interplay between stool microRNAs, microbiota and short chain fatty acids in Crohn's disease -a preliminary study

EPAC1 and EPAC2 promote nociceptor hyperactivity associated with chronic pain after spinal cord injury

Fostamatinib for the treatment of adult persistent and chronic immune thrombocytopenia: results of two phase 3, randomized, placebo-controlled trials

) miRNA expression change in dorsal root ganglia after peripheral nerve injury

Upregulation of miR-133a-3p in the sciatic nerve contributes to neuropathic pain development

Identification of key genes and pathways associated with neuropathic pain in uninjured dorsal root ganglion by using bioinformatic analysis

Alterations in the gut microbiota and metabolite profiles in the context of neuropathic pain

Role of macrophage inducible C-type lectin/ spleen tyrosine kinase signaling in the development of pain hypersensitivity state in rat

DGIdb 3.0: a redesign and expansion of the drug-gene interaction database

The anti-inflammatory and immunosuppressive effects of glucocorticoids, recent developments and mechanistic insights

Differences in the miRNA signatures of chronic musculoskeletal pain patients from neuropathic or nociceptive origins

Platelet-derived growth factor receptor-β antagonism restores morphine analgesic potency against neuropathic pain

The relation of the brain-derived neurotrophic factor with MicroRNAs in neurodegenerative diseases and ischemic stroke

CircRNAs: a new chapter in oral squamous cell carcinoma biology

MiR-940 promotes the proliferation and migration of gastric cancer cells through up-regulation of programmed death ligand-1 expression

Neuropathic pain: an updated grading system for research and clinical practice

2020) c-Abl-p38α signaling pathway mediates dopamine neuron loss in trigeminal neuralgia

Early treatment with UR13870, a novel inhibitor of p38α mitogenous activated protein kinase, prevents hyperreflexia and anxiety behaviors, in the spared nerve injury model of neuropathic pain

Functional passenger-strand miRNAs in exosomes derived from human colon cancer cells and their heterogeneous paracrine effects

Bioinformatics analysis identifies p53 as a candidate prognostic biomarker for neuropathic pain

The economic costs of pain in the United States

A circle RNA regulatory axis promotes lung squamous metastasis via CDR1-mediated regulation of Golgi trafficking

JTC-801 alleviates mechanical allodynia in paclitaxel-induced neuropathic pain through the PI3K/Akt pathway

Effects of progesterone on neuropathic pain responses in an experimental animal model for peripheral neuropathy in the rat: a behavioral and electrophysiological study

Non-coding RNAs in neuropathic pain

Evaluating the use of microRNA blood tests for gastric cancer screening in a stratified population-level screening program: an early modelbased cost-effectiveness analysis

Exercise tegulates microRNAs to preserve coronary and cardiac function in the diabetic heart

Targeting G protein-coupled receptor for pain management

Circ-ZNF609 aggravates neuropathic pain via miR-22-3p/ENO1 axis in CCI rat models

Prediction of circRNA-disease associations based on inductive matrix completion

Inflammation in joint injury and post-traumatic osteoarthritis

MicroRNA circulating in the early aftermath of motor vehicle collision predict persistent pain development and suggest a role for microRNA in sex-specific pain differences

LncRNA NONRATT021972 siRNA regulates neuropathic pain behaviors in type 2 diabetic rats through the P2X7 receptor in dorsal root ganglia

Deciphering the co-adaptation of codon usage between respiratory coronaviruses and their human host uncovers candidate therapeutics for COVID-19

Normalization of cholesterol metabolism in spinal microglia alleviates neuropathic pain

Modulation of pain transmission by G-protein-coupled receptors

miRNA-23a/CXCR4 regulates neuropathic pain via directly targeting TXNIP/NLRP3 inflammasome axis

Co-administered low doses of ibuprofen and dexamethasone produce synergistic antinociceptive effects on neuropathic mechanical allodynia in rats

2020) miRNA 146a-5p-loaded poly(d, l-lactic-co-glycolic acid) nanoparticles impair pain behaviors by inhibiting multiple inflammatory pathways in microglia

Limma powers differential expression analyses for RNAsequencing and microarray studies

The multiMiR R package and database: integration of microRNA-target interactions along with their disease and drug associations

Serum microRNA-based prediction of responsiveness to eribulin in metastatic breast cancer

Cytoscape: a software environment for integrated models of biomolecular interaction networks

In silico identification and validation of miRNA and their DIR specific targets in Oryza sativa Indica under abiotic stress

STRING v11: proteinprotein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets

The current clinical use of adjuvant analgesics for refractory cancer pain in Japan: a nationwide crosssectional survey

Micro-RNAs in the spinal microglia serve critical roles in neuropathic pain

Tamoxifen suppresses paclitaxel-, vincristine-, and bortezomib-induced neuropathy via inhibition of the protein kinase C/extracellular signal-regulated kinase pathway

Dynamic changes in the microRNA expression profile reveal multiple regulatory mechanisms in the spinal nerve ligation model of neuropathic pain

Altered thalamic neurotransmitters metabolism and functional connectivity during the development of chronic constriction injury induced neuropathic pain

Downregulated circular RNA zRANB1 mediates Wnt5a/β-catenin signaling to promote neuropathic pain via miR-24-3p/LPAR3 axis in CCI rat models

LncRNA X inactive specific transcript contributes to neuropathic pain development by sponging miR-154-5p via inducing toll-like receptor 5 in CCI rat models

RNA degradation in neurodegenerative disease

ggplot2: elegant graphics for data analysis

MiRNA-103/107 in primary high-grade serous ovarian cancer and its clinical significance

Alterations in mouse spinal cord and sciatic nerve microRNAs after the chronic constriction injury (CCI) model of neuropathic pain

TrkB.T1 contributes to neuropathic pain after spinal cord injury through regulation of cell cycle pathways

Cell cycle inhibition limits development and maintenance of neuropathic pain following spinal cord injury

Quantitative real-time PCR evaluation of microRNA expressions in mouse kidney with unilateral ureteral obstruction

Predictions of the dysregulated competing endogenous RNA signature involved in the progression of human lung adenocarcinoma

Circ_0005075 targeting miR-151a-3p promotes neuropathic pain in CCI rats via inducing NOTCH2 expression

Identification of the spinal expression profile of non-coding RNAs involved in neuropathic pain following spared nerve injury by sequence analysis

miR-940 potentially promotes proliferation and metastasis of endometrial carcinoma through regulation of MRVI1

Identification of novel therapeutic targets for neuropathic pain based on gene expression patterns