Bioinformatics Analysis of the Regulatory lncRNA-miRNA-mRNA Network in Patients with Neuropathic Pain

Background: Neuropathic pain (NP) is the main form of chronic pain, caused by damage to the nervous system and dysfunction. Methods: Here, we explore the key molecules involved in the development of NP condition via identication of lncRNA-miRNA-mRNA expression pattern of patients with NP. We identied differentially expressed miRNAs, lncRNA and mRNA through a comprehensive analysis strategy. Subsequently, we used bioinformatics approach to perform pathway enrichment analysis on DEGs and protein-protein interaction analysis. Combined with the three datasets, the lncRNA-miRNA-mRNA network was constructed. It will then be used as targets for drug prediction. Results: The results showed that a total of 8,251 DEGs (4,193 upregulated and 4,058 downregulated) were identied from the three microarray datasets, 959 DEmiRs (455 upregulated and 504 downregulated), 2,848 DElncs (1,324 upregulated and 1,524 downregulated). GO analysis showed that DEGs are mainly enriched in blood circulation, regulation of membrane potential and regulation of ion transmembrane transport. KEGG results showed that DEGs are enriched in neuroactive ligand-receptor interaction, PI3K-Akt signaling pathway and MAPK signaling pathway. When the correlation is set to above 0.8, a total of 31 lncRNAs, 36 miRNAs and 24 mRNAs were screened in the lncRNA-miRNA-mRNAs network. The results of drug prediction indicated the targeted drugs mainly include INDOMETHACIN, GLUTAMIC ACID and PIRACETAM. Conclusion: The lncRNA-miRNA-mRNA network has been carried out a comprehensive biological information analysis and predicted the potential therapeutic application of drugs in patients with NP. The corresponding data has a certain reference for studying the pathological mechanism of NP. lncRNA-microRNAs-mRNA


Background
Neuropathic pain (NP), the main form of chronic pain, is caused by damage to the nervous system and dysfunction (1,2). Currently, the treatment of NP is limited to symptomatic treatment. However, the prognosis of NP is poor. To better understand the pathogenesis of NP, it is essential to develop effective prevention strategies and improve the e cacy of NP. To deepen the genetic and various neurobiological bases of NP, it can provide important diagnostic and therapeutic tools for clinicians.
Long non-protein-coding RNA (lncRNA) is a transcript that is more than 200 nucleotide (nt) long.
LncRNAs can regulate a variety of processes, including DNA methylation, staining structure, transcription and post-transcription RNA treatment. LncRNAs can be involved in various gene transcriptional processes via interacting with transcription factors, coactivators and/or inhibitors. Accumulating evidences showed that lncRNAs play roles in the pathogenesis of nervous system disease (3,4). For instance, Zhou, Xiong (5) found that a total of 134 lncRNAs, 12 miRNAs, 188 circRNAs and 1066 mRNAs were signi cantly regulated after SNI surgery. The function analysis has also indicated that the pathogenesis of SNI was enriched in ribosome and ECM-receptor interaction. All this information will contribute to the study on NP pathogenesis. Wei,Li (3) found that lncRNA X inactive-speci c transcript (XIST) promotes NP progression in rats via down-regulation of miR-154-5p and up-regulation of TLR5s, representing a novel therapeutic target in NP treatment.
MicroRNAs (miRNAs), non-coding RNAs, regulate gene expression and can be often used as potential biomarkers. Dayer,Luthi (6) reported that chronic pain may induce the abnormal, speci c dysregulation of miRNA expression and they found that hsa-miR-320a and hsa-miR-98-5p (two c-miRNA signatures) can serve as a biomarker for pain-type classi cation. The epigenetic intervention of miRNA may also be a new therapeutic approach for complications such as injurious hypersensitivity caused by peripheral nerve injury. Pan, Shan (7) reported that miR-23a regulates the NP through the TXNIP/NLRP3 in ammasome axis in spinal glial cells targeting directly with CXCR4, suggesting a potential target for NP treatment.
With the rapid development of high-throughput sequencing, data mining, and the wide application of precision medicine, it becomes more feasible to extract aberrant data of ncRNAs and mRNAs on neuropathic pain from the microarray datasets. Although previous investigations have been reported the role of some ncRNAs in regulation of NP pathological pathways, the pathogenesis of NP has yet been fully understood. Therefore, in this study, in combination with bioinformatics and transcriptomics, three datasets from GEO (GSE107180, GSE145199 and GSE24982) were downloaded to identify differentially expressed mRNAs, miRNAs and lncRNAs. We then performed GO and KEGG functional enrichment analysis of differentially expressed mRNAs and the protein-protein interaction (PPI) to investigate gene interactions. Moreover, the lncRNA-miRNA-mRNA network has been constructed for the interaction of tree transcripts in the NP pathological pathway. Our study will contribute to the understanding of NP pathogenesis and providing new strategies for targeted drugs and therapies.

Raw data collection
The keywords "Neuropathic Pain" and "lncRNA" or "miRNA" or "mRNA" were retrieved using the Integrated Gene Expression Database (GEO, https://www.ncbi.nlm.nih.gov/geo/). Eligible datasets must be met with the following inclusion criteria: 1) A microarray analysis study belonging to human NP patients; 2) The dataset contains NP patients and control samples; 3) Each sample size has a group label; 4) The corresponding gene symbol annotation or gene bank ID of each probe of the microarray in the platform le; 5) The original data is available. The exclusive criteria are as follow: 1) datasets as animal samples or cell columns; 2) non-NP samples 3) non-microarray data types; 4) samples without controls; 5) probes without gene notation or library ID annotation. Finally, this paper lters and downloads GSE107180, GSE145199, and GSE24982 datasets.

Raw data pretreatment and screening for DElncs, DEmiRs and DEGs
The GEO dataset is the mainstream common functional genomics data repository. However, the data uploaded is often crude, incomplete, or noisy. Therefore, the used data are pre-processed before data mining, including data cleaning, data ltering and data normalization processing. For Illumina expression arrays, microarray quality control methods are commonly used for background correction and quantile normalization. This method can eliminate the non-experimental differences caused by technical differences and ensure the comparability of data between different samples. The original signal strength data of GSE107180, GSE145199 and GSE24982 were converted using Illumina package (Version 4.0.1 of R). The limma package was used to identify DElncs, DEmiRs, and DEGs. The gene symbol was annotated to the corresponding probe and calculate the average value of the probes when the gene symbol mapped multiple probes. We used R's Gplot packages to map the volcanoes of DElncs, DEmiRs and DEGs. R software endpoint "SVa" packages were used to eliminate data heterogeneity and batch effects, and changes in sample types were retained for further analysis of DEGs. The "PheatMap" package of R is used to draw a hierarchical clustering diagram. Finally, we used the mixed variance analysis model with the error discovery rate (Benjamin-Hochberg Test (FDR). The expression levels of RNA molecules with signi cant differences were screened by using the adjusted P < 0.05 and |log2Fc| value > 1.
Protein-protein interaction (PPI) network of DEGs, screening for hub genes and lncRNA-miRNA-mRNA network The Database for Annotation, Visualization and Integrated Discovery (DAVID) integrating annotations, visualization, and integration ndings, was used to explore the main functions and pathways of DEGs, including Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment pathway analysis (8, 9). The DAVID is an integrated biological database containing gene function classi cation, pathway mining and clustering tools. These tools were developed to extract systematically biological information from uploaded gene lists. In addition, the Malacards database was used in the analysis in order to ensure the accuracy of the results. Malacards is an integrated database of human diseases with annotations from several noted data sources. Maldacards contains fourteen annotated topics, including symptoms, related diseases, summaries, medications and treatments. We adjusted Pvalue using the Benjamini method or FDR method in multiple test corrections. In addition, P < 0.05 was selected as the threshold and the enrichment score was de ned as conversion log 10 (P-value).
A search tool for the retrieval of interacting Genes (https://string-db.org/) is an online functional protein association network. The association in STrinGd is speci c and biological, including both direct (physical) and indirect (functional) interactions. We uploaded the DEGs required for analysis to the STrinG and obtained the potential PPI network. The hub genes are key genes that play an important role in a variety of biological processes and often play regulatory roles on the expression of other genes in the related pathways. Therefore, the hub gene is also considered as an important target. Hub genes are screened based on the network topology. CytoHubba is an application for an effective Cytoscape version 3.6.1 (http://www.cyto-scape.org/) that can be used to identify hub genes accurately through twelve topological analysis methods. Cytoscape is an open-source software platform to visualize molecular interaction networks. These networks are mixed with annotations of gene expression pro les and other data. We used Maximal Clique Centrality (Mcc) to identify the top 20 hub genes and input the data of paired lncRNA-miRNA-mRNA into Cytoscape software in order to generate the regulatory network image.

Drug prediction of DEGs
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 (10). It can be standardized content from 30 different sources and allows to search and lter information. It is easily accessed through an intuitive Web user interface, an application programming interface (API), and publicly cloud-based server images. We uploaded the DEG list identi ed in the mRNA dataset onto the drugBank to predict their relevance of potential targeted drugs. The aim of this study is to reveal the molecular mechanism of the drugs currently used in NP and provide potential target genes for the development of future drugs.

Detection of DElncs, DEmiRs and DEGs
After data preprocessing, the raw data of each data set were normalized. The boxplot (Fig. 1) showed the difference of each sample before and after data normalization. The volcano plots ( Fig. 2A-2C

GO and KEGG functional enrichment analysis of DEGs
The "clusterPro ler" package was used for GO and KEGG analysis of 8,251 DEGs detected (Fig. 3). GO and KEGG enrichment analysis has a wide range of applications in bioinformatics and can be explained biological mechanisms and functional pathways in genomics and/or transcriptomics. The GO analysis can be divided into "Biological processes" (BP), "Cellular components" (CC), and "Molecular functions" (MF). The functional analysis indicated that the DEGs were enriched in blood circulation, regulation of membrane potential, regulation of ion transmembrane transport, regulation of metal ion transport, metal ion transmembrane transporter activity, amide binding, gated channel activity, ion gated channel activity, synaptic membrane, postsynaptic specialization, neuron to neuron synapse and asymmetric synapse.
Additionally, KEGG analysis showed that they were mainly participated in neuroactive ligand-receptor interaction, PI3K-Akt signaling pathway, MAPK signaling pathway, Rap1 signaling pathway, cAMP signaling pathway and Calcium signaling pathway.

Construction of PPI, hub gene and lncRNA-miRNA-mRNA networks
The PPI analysis is to study the molecular mechanism of diseases from a systematic perspective and discover novel drug targets. The PPI network is run on STRING to explore the interactions between the rst 50 DEGs encoded proteins. We eliminated the unconnected nodes in the network and nally got28 nodes and 45 edges (Fig. 4A) in the network. Additionally, Cytoscape was used to integrate the biomolecular interactions of lncRNA-miRNA-mRNAs into a single network.
When the miRNA-mRNA correlation was set to be greater than 0.8, a total of 24 miRNAs were found, which jointly regulated 381 mRNA targets, among which the top ve miRNAs that regulated the amount of mRNA were miR-940, miR-200C-3p, miR-5192, miR-2277-3p, and miR-5698, respectively. When the miRNA-lncRNA correlation was set to be greater than 0.8, a total of 31 lncRNAs were found and they jointly regulated 36 miRNAs, among which the top ve lncRNAs that regulated the number of miRNAs were CDCA3, SCAMP1, RPL18, NDUFA4, and MGST1, respectively. The results (Fig. 4B) showed that lncRNA and miRNA were highly clustered, among which miRNA-94, miRNA-2277 and miRNA-5192 had the most lncRNA binding sites.

Drug prediction of DEGs
The DGIdbu database provides the association between genes and their matching known or potential drugs. However, it cannot guarantee that any given drug gene association represents an appropriate therapeutic intervention. After processing by the DGIdbu database, we screened the approved drugs related to neuropathic pain from 8,251 results based on the keyword "neuropathic" (Supplementary Table 1). Most of these drugs are some arthritis, antiepileptic drugs and sedatives. The effects of the drug include treatment of osteoarthritis, rheumatoid arthritis, treatment of tonic clonic seizures, regulation of heart rate, body temperature, blood pressure, appetite, attention, mood and reactions related to alertness or alarm conditions.

Discussion
Neuropathic pain is the main form of chronic pain. Neuropathic pain is mainly caused by damage and disorder of the nervous system. Therefore, it is important to investigate the pathological mechanism of NP and identify effective therapeutic targets.
Recently, several studies have been highlighted the function and clinical signi cance of non-coding RNA in NP (6, 7). However, RNA crosstalk on NP has yet been resolved. This work combined with lncRNA, miRNA, and mRNA expression pro les and constructed a lncRNA-miRNA-mRNA regulatory network is to provide new ideas for post-transcriptional gene regulation of NP. In addition, with the help of big data mining and public database annotations, we explored the key molecules of multiple physiological and pathological processes of NP and identi ed several potential NP drug targets by reviewing the previous work.
The identi ed DEmiR and DEGs have not been reported in the previous study related to the pathogenesis of neuropathic pain. MiR-940 is mainly related to proliferation and metastasis of cancer, such as gastric cancer (11), endometrial carcinoma (12), thyroid cancer (13), and cervical cancer (14). Interestingly, miR-940 also regulates the in ammatory response in osteoarthritis (OA) (15) and spinal cord injury (16). lncRNA XIST enhances ANLN by up-regulating miR-200c-3p, promotes the chemical resistance of breast cancer cells to doxorubicin and provides a potential new therapeutic target for the breast (17). The exosomes rich in miR-2277-30p tend to have a high abundance of β4 and lead to a cancer promoting effect (18). The serum levels of miR-5689 were signi cantly correlated with overall survival after initiation of eribulin treatment, providing an evidence that serum miRNA can serve as biomarkers for the predicting the development of new distant metastasis (19). These studies to a certain extent illustrate that the DEmiRs we identi ed have potential research value in diseases, especially cancer. In addition, these results also need to be further veri ed by later work.
Similarly, the lncRNA found in this article has not been fully reported in previous studies. It is reported that knockdown of lncRNA SCAMP1 suppressed malignant biological behaviors of glioma cells through regulating the miR-499a-5p/LMX1A/NLRC5 pathway (20). Additionally, SCAMP1 could be bound with miR-137 and accelerate the progression of ovarian cancer (OC) (21). Zhang, Xu (22) constructed a completing endogenous RNA (ceRNA) to explore the molecular mechanism of cisplatin (DDP) resistance in non-small cell lung cancer (NSCLC) cells. They found that MSTRG51053.2 lncRNA may be used as ceRNA for miR-432-5p to regulate DDP resistance in non-small cell lung cancer.
According to the recommendations from the international association (23), SNRI-duloxetine, venlafaxine, tricyclic antidepressants, and Gabapentin, pregabalin are recommended for rst-line; while Capsaicin 8% patches, Lidocaine (lignocaine) patches, and Tramadol are recommended for second-line; and the strong opioids are recommended for the third line. Most of the drugs that were predicted in our study are arthritis, antiepileptic drugs, and sedatives, indicating that the identi ed genes were associated with the pathogenic mechanism of NP. However, it does not mean that we can use these predicated drugs in the clinic since there needs more research about these genes.
However, there are still some limitations of this study. First of all, we found that lncRNAs and microRNAs were related to neuropathic pain based on the GEO database, but it still needs to be veri ed by PCR analysis or immunoblotting. Secondly, the lncRNA-microRNAs-mRNA network constructed in our work could not indicate whether these RNAs co-express altogether in the same tissue. Finally, the speci city of the identi ed DElnc, DEmiRs and DEGs of neuropathic pain needs to be veri ed in further study.

Conclusions
This study was constructed a lncRNA-microRNAs-mRNA network correlated to neuropathic pain and identi ed potential disease targets through the integration of the microarray data and rigorous analysis.
Our methods can be provided a new idea and references for further study of the pathological mechanism of NP. In addition, this study will contribute to the discovery of new biomarkers or treatments that can improve the prognosis of neuropathic pain and reduce the risk of disease.

Data availability
All data generated or analyzed during this study are included in this published article. Data is available at NCBI GEO: GSE107180, GSE145199, and GSE24982.
HZ analyzed the data, prepared gures 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 nal version of the manuscript.  Differential analysis of lncRNAs, miRNAs, and mRNAs, respectively. Volcano plots of lncRNAs (A), miRNAs (B) and mRNAs (C). 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. The top ten lncRNAs (D), miRNAs (E) and mRNAs (F) with the smallest P-values in the upregulated and downregulated cluster were taken out respectively for the drawing of a heatmap.