Bioinformatics Analysis of the MicroRNA-Metabolic Gene Regulatory Network in Neuropathic Pain and Prediction of Corresponding Potential Therapeutics.

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.

Introduction

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 miRNAs, 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 (Wang et al. 2020). 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 microRNAs 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 (Chen et al. 2017; 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 (Chen et al. 2017); 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.

Methods

Raw Data Collection

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.

Raw Data Preprocessing and Screening for Differentially Expressed miRNAs and mRNAs and DEMGs

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.r-project.org/web/packages/pheatmap/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://bioinformatics.psb.ugent.be/webtools/Venn/).

Construction of miRNA-DEMG Network

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.

Protein–Protein Interaction Network

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

Functional enrichment analysis of sets of genes was performed using the online platform g:Profiler (https://biit.cs.ut.ee/gprofiler/gost). R. norvegicus was chosen as organism, and g:SCS threshold of 0.05 was chosen as significance cut-off 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.

Drug Prediction and Repurposing for the Set of Genes

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 L1000CDS2 pharmacogenetic search tool (https://maayanlab.cloud/l1000cds2/#/index) was applied to DEMGs to further retrieve potential perturbant drugs targeting these genes.

Results

Detection of Differentially Expressed miRNAs and mRNAs

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.

Fig. 1

Box plot of each sample distribution in GEO data before and after normalization. The blue boxplots represent the distribution of the original data, while the red boxplots represent the distribution of the normalized data. A GSE145199 and B GSE24982

Fig. 2

Differential analysis of miRNAs and mRNAs in NP. Differentially expressed genes were screened by the criteria of |log2FC|> 1 and P-value < 0.05, the green spots represented the downregulated molecules, while the red spots represented the upregulated molecules. Volcano plot (A) and heatmap of top 10 differentially expressed miRNAs (B). Volcano plot (C) and heatmap of top 10 differentially expressed mRNAs (D)

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).

Fig. 3

Bubble chart of GO and KEGG pathway analysis of the downregulated differentially expressed mRNAs. The bubble size represented the count of differential expressed mRNAs, and the different colors of the bubble represented the indicated functional category

Fig. 4

Bubble chart of GO and KEGG pathway analysis of the upregulated differentially expressed mRNAs. The bubble size represented the count of differential expressed mRNAs, and the different colors of the bubble represented the indicated functional category

Protein–Protein Interaction Network of the DEGs

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).

Fig. 5

Protein–protein interaction network of 79 hub genes identified in the protein–protein interaction network of differentially expressed mRNAs. Different colors represent different clusters, and the clustering result was obtained by the k-mean clustering method. The green nodes represent downregulated genes, while the red nodes represent upregulated genes

Fig. 6

Bubble chart of GO and KEGG pathway analysis of the 79 hub genes identified in the protein–protein interaction network of differentially expressed mRNAs. The bubble size represented the count of mRNAs, and the depth of color of the bubble represented the significance of the mRNAs based on the P-values

Target Genes of Differentially Expressed miRNAs and Construction of miRNA-DEMG Networks

To identify the miRNA-mRNA interactor pairs, the Pearson correlation between the differentially expressed miRNAs 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 miRNAs, the metabolic process gene list and predicted miRNA targets allowed the identification of 110 metabolic process-related differentially expressed genes (DEMGs) as targets of 22 differentially expressed miRNAs. These 22 miRNAs were considered as the metabolism-related signature miRNAs in NP. The PPI network of the 110 DEMGs targeted by miRNAs was constructed in the STRING database and revealed strong interactions between proteins encoded by the 110 DEMGs. We eliminated the unconnected nodes in the network and finally got a network with 83 nodes and 126 edges (Supplementary Figure S3). The average number of neighbors was 3.477, while the network diameter and network radius were 10 and 5, respectively. The MCODE identification of hub genes indicated Cks2, Mcm4, Mcm5, Mcm6, Oxa1l, Mrpl17, and ENSRNOG00000047563 as the most significant hub genes in the PPI network of the 110 DEMGs (Supplementary 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), whereas for cellular component, intracellular organelle lumen, membrane-enclosed lumen, and organelle lumen were the most enriched GO terms (Fig. 8; Supplementary File S9). In the KEGG pathway enrichment analysis, RNA degradation, cell cycle, and MAPK signaling pathway were the most enriched pathways (Fig. 8; Supplementary File S9).

Fig. 7

Construction of the miRNA-DEMG network. The ellipses represent the DEMGs and the triangles represent miRNAs. The orange color indicate hub miRNAs and DEMGs; the green color indicates non-hub miRNAs; the violet color indicates non-hub DEMGs

Fig. 8

Bubble chart of GO and KEGG pathway analysis of the DEMGs in the miRNA-DEMG network. The bubble size represented the count of mRNAs, and the depth of color of the bubble represented the significance of the mRNAs based on the P-values

Drug Prediction for NP

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 the genes targeted by most of the drugs were DRD2 (177 drugs), OPRM1 (124 drugs), HTR2 (146 drugs), CHRM3 (100 drugs), OPRK1 (89 drugs), AGTR1 (49 drugs), and ADRA1D (47 drugs) (Supplementary File S11).

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 L1000CDS2 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).

Fig. 9

Drug repositioning analysis of the DEMGs in the miRNA-DEMG network

Discussion

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 (Huang et al. 2020). 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 (Li et al. 2017). 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, 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 (Yang et al. 2020). 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 (circRNAs) 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 miRNAs 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.

Conclusions

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.

Below is the link to the electronic supplementary material.

Supplementary file1 (TIF 21061 KB)

Supplementary file2 (TIF 12338 KB)

Supplementary file3 (TIF 16251 KB)

Supplementary file4 (XLSX 76 KB)

Supplementary file5 (XLSX 788 KB)

Supplementary file6 (CSV 2950 KB)

Supplementary file7 (CSV 7081 KB)

Supplementary file8 (CSV 38 KB)

Supplementary file9 (XLSX 133 KB)

Supplementary file10 (XLSX 1648 KB)

Supplementary file11 (XLSX 23 KB)

Supplementary file12 (CSV 75 KB)

Supplementary file13 (XLSX 196 KB)

Supplementary file14 (XLSX 73 KB)

Supplementary file15 (XLSX 17 KB)