Decreased Expression of ENPP6 Predicts the Occurrence of Pain in Malignant Spinal Tumour


 Background: Neuropathic pain (NeP) characterized by neuroplasticity and neuroinflammatory change is a common complication associated with spinal metastasis. However, there are no reliable candidates for diagnosis and treatment. Recently, cancer research has incorporated molecules into the treatment of patients with NeP, therefore, it is necessary to find key molecules of NeP to provide new targets for diagnosis and treatment.Methods: We analyzed RNA-seq data around the expression of ENPP6 based on bioinformatic methods, including differentially expressed genes (DEGs), Gene Ontology (GO), protein-protein interaction (PPI) network, Kyoto Encyclopedia of Genes and Genomes (KEGG) and GSEA analyses, receiver operating characteristic (ROC) curve, immune cell infiltration and mutation analysis.Results: We divided with pain samples into the High and Low ENPP6 expression groups. A total of 231 DEGs were identified. GO and KEGG analysis showed that DEGs were mainly enriched in Inflammation and cancer associated pathways. GSEA analysis showed that DEGs was significantly enriched in ARF3 and P38/MK2, RHO and RAS, and BRAFT and AKT1/E17K pathway. Pearson’s correlation analysis showed that the expression of ENPP6 was significantly correlated with autophagy phenotype and immunophenotype. Immune infiltrating analysis showed that activated NK cells were significantly highly expressed in Low group. ROC analysis of ENPP6 suggested that the area under the ROC curve was 0.925. Mutation sites analysis showed that most of the mutations in ENPP4-7 were phosphorylation sites.Conclusion: This study provides novel insights into molecular mechanisms underlying NeP, and identifying ENPP6 may serve as potentially diagnostic biomarkers and/or therapeutic targets for NeP.


Introduction
Chronic pain has been estimated to affect one-sixth of the population [1]. Cancer-related neuropathic pain (NeP) is a neuroplasticity and neuroin ammatory change caused by damage to the somatosensory system caused by treatment, cancer or paraneoplastic reactions to cancer [2,3]. Spine is a common organ for metastatic cancer, and in patients with advanced cancer, spinal metastasis can lead to signi cant pain or neurological dysfunction, which impacting the quality of life [4,5]. Although advanced radiotherapy and surgical techniques are available for patients with advanced spinal metastases [6,7], some of them have a persistent pain during or shortly after radiotherapy and surgery treatment [8]. Current targeted drugs have many limitations, such as unsatisfactory e cacy and uncontrollable dosage. And few studies have focused on the diagnosis and treatment of patients with cancer-related NeP [2]. Notably, cancer research nowadays has incorporated a variety of molecules into the treatment of cancer patients with NeP, therefore, it is necessary to nd key molecules of NeP to provide new targets for diagnosis and treatment.
Ectonucleotide Pyrophosphatase/Phosphodiesterase6 (ENPP6, ENSP00000296741), which belongs to ENPP family, comprising seven members with structurally related catalytic domains [9]. ENPP1-3 are type II transmembrane proteins, which consist of catalytic domains and nuclease-like domain. ENPP4-7 are type I membrane proteins, which consist of catalytic domains only [10]. ENPP6 is a plasma membrane associated or secreted ectoenzyme that can hydrolyze glycerophosphocholine (GPC) and lysophosphatidylcholine (LPC) and contributes to supplying choline to the cells [11,12]. In mice, ENPP6 mRNA was primarily detected in kidney and brain with a lesser expression in heart [13], and in human it was detected in kidney, ovary and brain [11,14]. Recent research shows that ENPP6 may play a critical role in diseases of the bone, nervous system and cancer [15][16][17][18].
The management of NeP in malignant spinal tumour is extremely complex, so its treatment still faces great challenges. Many promising treatments found in animal models have failed in clinical trials, possibly because of basic cellular and molecular differences between animals and humans [19,20]. It may be possible to identify new drug targets by screening and analyzing the human gene networks associated with pain formation and progression. There has been no report about the correlation between ENPP6 and NeP in patients with cancer. Therefore, this study aims to analyze human RNA-seq data from PAIN Neurobiology Research Group using bioinformatics methods to explore the potential value of ENPP6 in the diagnosis and treatment of patients with spinal tumour.

Data of Patients and Treatment
Datasets and corresponding clinical information were collected of 15 patients undergoing treatment at MD Anderson Cancer Center for malignant tumour involving the spine. The data collection and this study were conducted in compliance with all applicable laws, regulations and policies for the protection of human subjects, and any necessary approvals, authorizations, human subject assurances, informed consent documents [21]. The data type was RNA-seq, and the tissue sample was human dorsal root ganglion (DRG) neuron. This study included 21 samples which from 15 patients, including n = 16 with pain samples (P) (n = 7, 1-6 months; n = 4, 6-12 months; n = 5, > 12 months) and n = 5 no pain samples (NP). Randomly divided n = 16 with pain samples into two groups (P1, P2), and made sure that the samples in each group were from different patients. Each with pain group was compared separately with the no pain group, then a cross validation for each group were performed.
Then, the transcripts per million (TPM) data for 21 samples were used for the following analyses.

Screening Differential and Co-expression Gene.
These raw data already been background corrected and normalized. According to the median value of ENPP6 (TPM = 6.82), the 16 with pain samples were divided into the high (> 6.82; n = 8) and low (≤ 6.82; n = 8) ENPP6 expression groups. Differentially expressed genes (DEGs) and co-expression gene were screened using the R package DEseq [22] and limma [23], according to the fold change (FC) and signi cant difference (p value) between groups. The threshold was set at delta = 1, log 2 FC > 0.5) and p value < 0.05, then we obtained 231 DEGs. Next, volcano plots and heatmaps of differential gene were performed using R package, ggplot2 and pheatmap.

Human Protein Atlas
The Human Protein Atlas (HPA, https://www.proteinatlas.org/) is a database with the aim to map all the human proteins in cells, tissues and organs using an integration of various omics technologies, including antibody-based imaging, mass spectrometry-based proteomics, transcriptomics and systems biology. HPA consists of six separate parts (Tissue, Single Cell Type, Pathology, Blood, Brain, Cell). We obtained the expression analysis and Immunohistochemical (IHC) images of ENPP6 from the Tissue, Pathology, Blood and Brain Atlas.

GO Enrichment and KEGG Pathway Analysis.
DAVID [24,25] was used for GO analysis of DEGs, which threshold was set at FDR < 0.1, then the visualization of results was performed using the ggplot2 R package. For Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis, we used the clusterPro ler package[26] and KEGG Orthology Based Annotation System 3.0 (KOBAS 3.0, http://kobas.cbi.pku.edu.cn) with p < 0.05 were considered to be signi cantly enriched, then produced the relevant bubble and bar charts.

Gene set enrichment analysis
Gene Set Enrichment Analysis (GSEA, version 4.0.3) [27] was used to analyze the DEGs matrix between ENPP6 high and low expression group, with the reference gene set (c2.cp.v7.0.symbols.gmt) and signi cance threshold of p < 0.05. Then, we use the RColorBrewer package for the visualization.

protein-protein interaction (PPI) network
STRING database[28] was used to construct the PPI network of DEGs, with the interaction score was set to con dence of 0.4, and the organism was set to 'Homo sapiens'. Currently, this database contains 18,838 proteins and 25,914,693 core interactions network. After construction of the PPI network, we use MCODE (MCODE scores > 5, node score cut-off = 0.2, degree cut-off = 2, Max depth = 100, and k-score = 2) and cytoHubba in Cytoscape to perform key sub-module analysis to obtain Hub genes and sub-modules with biological signi cance.

Phenotypic correlation analysis
The Human Autophagy Database (HADb) is the rst dedicated human autophagy database, a common repository containing information about human genes related to autophagy described to date. We have acquired autophagy related genes list from the Hadb (http://www.autophagy.lu/index.html) and immune related genes list from the Immunology Database and Analysis Portal Database (ImmPort, https://immport.niaid.nih.gov/home). Subsequently, 12 genes were obtained by intersected with autophagy/immune-related genes respectively for the 231 DEGs. Then, the 'plot' function in R was used to analyze the expression correlation between ENPP6 and the 12 genes, and visualized the results with signi cantly correlated.

CIBERSORT and LM22
CIBERSORT [29] is an analytical tool widely used to characterize the immune cell composition by gene expression values in tumors. LM22 signature matrix is a special genetic marker containing 547 genes that distinguish the 22 immune cell subtypes downloaded from the CIBERSORT (http://cibersort.stanford.edu/). In this study, CIBERSORT package and LM22 algorithm was used to analyze the abundance and expression divergence of 22 kinds of immune cells in ltrates between the groups with high and low ENPP6 expression in 16 pain samples.

Area Under the ROC Curve
Receiver operating characteristic (ROC) curve, which takes true (sensitivity) and false (1-speci city) positive rate as ordinate and abscissa respectively, R package "pROC" was applied to the analysis whether the expression of ENPP6 can distinguish with or without pain in 21 cases of samples, and determine the highest likelihood ratio of the optimal cut-off value to judge ENPP6 recognition threshold of pain. Area Under the ROC Curve (AUC) re ects the diagnostic value of biomarkers for diseases. In the case of AUC > 0.5, the closer the AUC to 1, the better the authenticity of diagnostic is.

Gene expression and mutations
The gene list of ENPP family members was obtained by using HUGO Gene Nomenclature Committee (HGNC) database [30]. Then we analyzed the mutations of ENPP family members in 32 TCGA pan-cancer data sets using the cBioPortal database. cBioPortal database is an open-source and open-data resource for interactive exploration of multiple cancer genomics database.

Statistical analysis
Data analysis was performed using R statistical software (version 4.0.2), and p < 0.05 was considered statistically signi cant.

Downregulation of ENPP6 in malignant tumors of the spine patients with pain
To explore the expression of ENPP6 in pain patients with malignant tumors of the spine in dorsal root ganglion tissues, we rst analyzed previously published datasets [21], including 21 (n = 16 with pain samples (P) and n = 5 no pain samples (NP)) samples of DRG from 15 patients which clinical features were summarized in Supplementary Table 1. Results showed that the expression of ENPP6 in pain group was signi cantly lower than that in no pain group ( Fig. 1. A, p = 0.003). However, there was no signi cant difference in other members of the ENPP family when compared the pain and no pain groups (Supplementary material Fig. 1).

Identi cation of DEGs in differential ENPP6 expression groups
After screened co-expressed genes by R package DEseq and limma (adjusted log 2 FC > 0.5 and p value < 0.05), 231 differentially expressed genes were obtained, Including 216 upregulated and 15 downregulated genes. The volcano plot was used to represent the DEGs ( Fig. 1. B), Subsequently, according to the median value of ENPP6 (TPM = 6.82), we divided 16 with pain samples into the high (H, > 6.82; n = 8) and low (L, ≤ 6.82; n = 8) ENPP6 expression groups. The heatmap was used to represent the distribution of the differential genes in the two groups with different ENPP6 expressions ( Fig. 1. C).

Expression of ENPP6 in human brain
In the HPA database, we found that ENPP6 was highly expressed in the human brain, kidney and reproductive system (

GO and KEGG pathway analyses correlated with DEGs
Using DAVID database, a total of 231 DEGs were subsequently applied to conduct GO and KEGG pathway enrichment analysis. In the "biological process" category, 47 enriched GO terms were observed. Categorization by "cellular component" revealed 21 enriched GO terms, such as neuronal cell body, endosome membrane and postsynaptic density. Moreover, the "molecular function" category revealed 19 signi cant enrichment, such as DNA − binding transcription factor binding, transcription corepressor activity, protein serine/threonine kinase activity. The details of the GO entries were described in Table 1. Results from GO enrichment analysis indicated that the DEGs were existed in endosome membrane of neuronal cell and were primarily associated with the immunoin ammation-related and regulation of trans-synaptic signaling GO terms, such as response to lipopolysaccharide, mononuclear cell migration, Ras protein signal transduction, modulation of chemical synaptic transmission and positive regulation of neurogenesis ( Fig. 3. A-B). To re ect the internal interactions between these GO items, we constructed the GO interaction network as shown in Fig. 3. C-D. showed that the DEGs were mainly enriched in In ammation and cancer associated pathways (Fig. 4 & Table 2), including chemokine, IL-17, MAPK, TNF, NF-κB and Toll-like receptor signaling pathway. The bubble and bar diagram were plotted by the clusterPro ler package, and the length of the column and size of the circle represents the number of genes enriched in the pathway and their color stand for p value, which increases gradually from red to blue (Fig. 4. A-B). Also, the pathways with the most signi cant number of genes represented were showed as Fig. 4. C.  We used the DEGs of two groups between high and low ENPP6 expression groups as the expression matrix, and 'c2.cp.v7.0.symbols.gmt' selected as the reference gene set for GSEA analysis with the standard thresholds p value < 0.05. RColorBrewer package were used for visualization of the GSEA results which showed that DEGs was enriched in ARF3 and P38/MK2 pathway in metabolism category of Pathway Interaction Database (PID), RHO and RAS pathway in BIOCARTA, and BRAFT and AKT1/E17K pathway in Reactome. And these pathways were signi cantly enriched in the ENPP6 low expression group (Fig. 5. D). These results indicated that the low expression of ENPP6 in with pain patients likely acted through above signaling pathways, and pain-mediated through indirect mechanisms such as a variety of cytokines and metabolic enzymes involved in the RAS and AKT pathway.

PPI network and key genes
We used the DEGs PPI networks which constructed by STRING database to predict protein-protein interactions (Fig. 6. A). Using the MCODE in Cytoscape to obtain the Hub genes, it shows the maximum correlation hub gene set, such as GNAI1, CCL1, CCL13 and P2RY14 (Fig. 6. B). Then, we obtained the bar chart of the Top 30 hub gene interactions using MCC analysis in cytoHubba (Fig. 6. C), and the length represents the number of molecular interactions. Next, we show the collection of hub genes obtained through using cytoHubba (Fig. 6. D), with darker colors representing higher related.

ENPP6 is correlated with autophagy phenotype and immunophenotype
A set of 2730 autophagy and immune related genes were obtained from HADb and ImmPort database.
231 DEGs were intersected with autophagy-related genes (232) and immune-related genes (2498) respectively, the result obtained 2 autophagy-related gene and 10 immune-related gene. Then, the Venn graph was generated using VennDiagram package. Subsequently, Pearson's correlation was used to analyze the correlations we analyzed the correlation between ENPP6 and 12 genes related to autophagy and immunophenotype, the phenotype related to differential gene and Pearson's correlation were listed in Table 3. Then, the highly score results were visualized using plot function in R software, displayed as correlation scatter plots (Fig. 7).

cells in ltration, including B cells, Plasma cells, T cells, NK cells and myeloid subsets, between
the high and low ENPP6 expression groups in 16 with pain samples. As shown in the bar chart, each column represents one sample, each color represents one type of immune cell. From this, we can determine the in ltration abundance of each kind of immune cell in each sample (Fig. 8. A). Then we analyzed the expression differences of immune cells between the groups with high and low ENPP6 expression, and found that the activated NK cells were signi cantly highly expressed in the tissues with low ENPP6 expression, which was shown in the violin diagram (Fig. 8. B). Above results suggested that neuropathic pain may be closely related to in ammatory and metabolic pathways, and may change the proportion of immune cells distributed by predicting the targets small-molecule drugs of activated NK cells, thereby improving pain.

The diagnostic value of ENPP6 expression in with pain patients
In order to evaluate the diagnostic performance of ENPP6 for neuropathic pain, the pROC package was used to draw the receiver operating characteristic (ROC) curve to analyze whether the expression of ENPP6 could distinguish between 5 non-pain samples and 16 pain samples, and determined the optimal cut-off value for generating the best likelihood ratio to determine the recognition threshold of ENPP6 for pain. The area under the curve (AUC) value was 0.925 (95% con dence interval [CI] = 0.8-1.0), indicating that the expression of ENPP6 can be of great value in diagnostic in neuropathic pain (Fig. 8. C).

The state of mutation in the ENPP family
Since the results suggested that ENPP6 was associated with the prognosis of NeP in cancer patients, the mutation status of this gene was analyzed furtherly. First, the list of ENPP family members, including ENPP1-7, were obtained from the HGNC database [31]. Then, the mutations of ENPP family members in 32 TCGA pan-cancer data sets was analyzed using the cBioPortal database. The results showed that mutation sites existed in ENPP 1-7 ( Fig. 9. A). Subsequently, the mutations existing in 7 members of the ENPP family were presented, and the results showed that there were somatic, phosphorylation and endonuclease mutation sites of ENPP3, ENPP4 and ENPP5, while most of the mutations in ENPP4, ENPP5, ENPP6 and ENPP7 were phosphorylation sites ( Fig. 9. B).

Discussion
Cancer-related neuropathic pain (NeP) is a common cause of chronic pain, which is one of the most di cult clinical problems [1,2]. Despite an increasing number of available diagnosis and therapies, adequately diagnosis of NeP is often di cult and satisfactory pain control is not always achieved [2,6,7].
Therefore, novel molecular network biomarkers or targets are needed for earlier detection or treatment.
In the present study, we collected malignant tumors of the spine data based on RNA sequencing from published datasets, and demonstrated ENPP6 was signi cantly down-regulated in with pain samples compared with no pain samples. There are few reports on the expression of ENPP6 in cancer, while Yano Y et al. demonstrated that ENPP3 which is the same family as ENPP6 was expressed in the tumor cells of bile duct malignancies and may play a role in tumor in ltration [32]. And ENPP2 was expressed in thyroid cancer cells and may associated with tumor cell motility and tumorigenic capacity [33]. Then, we analyzed the function of ENPP6 in malignant tumors of the spine using GSEA, and the results indicated that the low expression of ENPP6 in with pain tissue may mediate pain through indirect mechanisms such as a variety of cytokines and metabolic enzymes involved in the RAS and AKT pathway. There is evidence indicates that Ras-AKT signaling can promote the progression of glioma [34]. Further analysis using the HPA database showed high expression of ENPP6 in normal brain tissue and low expression in brain tumor tissue (glioma) [35], which consistent with the foregoing results. Above evidence suggested that ENPP6 may be related to cancer and was differentially expressed in cancer.
Subsequently, using GO and KEGG pathway analysis, we found that DEGs was mostly enriched in the immunoin ammation-related and regulation of trans-synaptic signaling GO terms, and In ammation and cancer associated pathways. Accumulating evidence shows that in ammation involved in the development of NeP. It has been demonstrated that ENPP6 is expressed in oligodendrocytes which necessary for transmission of bioelectrical signals and the protection of the normal function of neurons [36]. In addition, ENPP2 have been revealed may play an important role in neural and vascular development, tumour progression and metastasis, as well as in ammation, neuropathic pain and brotic disease [37]. Xiao L et al have reported that ENPP6 may play a role in lipid metabolism during myelin sheath formation and might be required to initiate myelination rapidly in response to differentiation induced signals [36]. ENPP6, ENPP2 and ENPP7 shared recognition of phospholipids with choline.
Evidence supported that ENPP7 is associated with anti-in ammatory and anti-tumorigenic activity by affecting the conversion of sphingomyelin to ceramide [38,39]. Through these data, we found that ENPP6 may be involved in lipid metabolism, in ammation, cancer or nerve signaling, which prompted us to further analyze.
Afterward, the results of CIBERSORT analysis for the proportion of tumor-in ltrating immune cell (TIC) revealed that the expression of ENPP6 was positively correlated with immunophenotype and autophagy phenotype in spinal tumors patients. Then, we analyzed the expression differences of immune cells between the groups with high and low ENPP6 expression, and found that the activated NK cells were signi cantly highly expressed in the tissues with low ENPP6 expression. Natural Killer (NK) cells are lymphocytes with the capacity to target tumor cells via innate and adaptive responses [40,41]. Activated NK cells can rapidly produce cytokines and activate other leukocytes, resulting increased complex uctuations of cytokines observed in the blood and cerebrospinal uid of NeP [42,43]. Immune effector cells, especially NK cells, are associated with NeP [44]. Gao YH et al. reported that electroacupuncture improved NeP by affecting the activity and number of NK cells [45]. Morphine can inhibit the cytotoxic activity of NK cells through opioid receptors and Toll-like receptor-4 (TLR4), which is of great signi cance for the maintenance of immune function during pain [46]. Therefore, it is tempting to speculate that drugs targeted in activated NK cell may improve the proportion of immune cell distribution, and alleviate pain subsequently.
Furthermore, by analyzing the area under the ROC curve (AUC), we surprisingly found that ENPP6 had a good performance in diagnosing NeP. Recent reports indicated that choline plays a role in NeP and neuroin ammatory disorders [47,48], and the role of ENPP6 in supplying choline to the cells have been well demonstrated [11]. Then, we analyzed the mutation status of the ENPP family, and the results showed that there were most of the mutations in ENPP4-7 were phosphorylation sites. Molecular and channel phosphorylation can affect the pathophysiological processes of NeP [49,50]. García G et al reported that with mutations at PKC/PKA phosphorylation sites reversed tactile allodynia in neuropathic rats [51]. These suggested that the mutation of ENPP6 phosphorylation site may be one of the pathogenic links.
There are some limitations of our study. Firstly, the current study was only based on data analysis and additional experiments are needed to demonstrate the biological impact of ENPP6 in NeP. Secondly, the sample size of the data involved was small, and the study failed to cover different regions, which may affect the gene expression in tumors. Thirdly, because our study only focused on those genes that showed signi cant changes in the data set, some biological information might be ignored in our study.
Therefore, further studies about direct mechanisms in NeP are needed.
In conclusion, this is the rst study that uncovers ENPP6 may had a good performance in diagnosing NeP and may had an important role in the regulation of in ammatory and cancer pathways in malignant spinal tumors. The present study may offer new ideas for diagnosis and treatment of malignant spinal tumour patients with NeP.

Declarations
Ethical approvaland consent to participate This study does not contain any studies with human participants or animals performed by any of the authors.

Consent for publication
All the authors have reviewed the nal version of the manuscript and approved it for publication.

Availability of data and material
All data analyzed or generated are included in this article. Green and red tones represent downregulated and upregulated genes. P, with pain samples; NP, no pain samples; DEGs, differentially expressed genes; H, ENPP6 high expression group; L, ENPP6 low expression group.        Supplementary Files