Identification of new biomarkers and pathways associated with treatment and prognosis in melanoma patients

Skin cutaneous melanoma (SKCM) is known as the most malignancy and treatment-resistant in human tumor, causing about 72% of deaths in skin carcinoma. However, the potential mechanism and new effective targets remain to be further elucidated. Available datasets such as Gene Expression Omnibus (GEO) can be utilized to search for novel therapeutic targets and prognostic biomarkers. Three data sets were downloaded from GEO database . The differentially expressed genes (DEGs) were identified via Venn software. Protein‐protein interaction network of DEGs was developed and the module hub genes analysis was constructed by Cytoscape. Subsequently, multiple online tools and Kaplan-Meier survival curves were analyzed to detect underlying signaling pathways, gene expression, drug-gene interaction and prognostic value of hub genes. In addition, we explored the correlation between hub genes and immune cell infiltration. At last, the related miRNA, lncRNA networks were constructed by R software. enrichment CXCL9, CXCL10 in the occurrence and development of melanoma. Our results revealed that CXCL9, CXCL10, CXCL13 are positively associated with immune cell infiltration.


Conclusion
In conclusion, CCL4, CCL5, NMU, GAL, CXCL9, CXCL10, CXCL13 were of high prognostic value and may be potential targets for the diagnosis and therapy of patients with melanoma.

Background
Skin cutaneous melanoma (SKCM) accounts for only 2% of total skin cancer. However, due to its high degree of malignancy and invasiveness, it causes about 72% of deaths in skin carcinoma [1]. The incidence of cutaneous malignant melanoma continue to increase annually [2]. Melanoma has become a serious public health problem, bringing great economic burden for society [3]. Early diagnosis and treatment can improve the survival rate of patients with SKCM [4]. Melanoma is associated with multiple risk factors especially the sun exposure [5]. The general progression models of melanoma are from nevus to dysplastic nevus, melanoma in situ to invasive melanoma [6]. However, there is a lack of comprehensive research on melanoma in different periods, which limits people's understanding of the key genes in these periods, such as the initial stage and metastasis stage. Therefore, it is extremely important to explore the potential molecular mechanism and immune microenvironment of melanoma for better treatment and prognosis. During the last decades, microarray technology and bioinformatic analysis have been widely used 4 to screen genetic alterations at the genome level, which have helped us identify the differentially expressed genes (DEGs) and functional pathways involved in the tumorigenesis and progression of melanoma. In the present study, we chose and analyzed 3 mRNA microarray datasets from Gene Expression Omnibus (GEO) between primary melanoma (PM) and normal skin (NS). Subsequently, Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis and proteinprotein interaction (PPI) network analyses were performed to help us understand the molecular mechanisms underlying carcinogenesis and progression. In conclusion, a total of 308 DEGs and 12 hub genes were identified. We assessed the diagnostic, treatment and prognostic value of hub genes, which may be new biomarkers for SKCM.

Methods
Microarray data GEO (http://www.ncbi.nlm.nih.gov/geo) is a public functional genomics data repository of high throughout gene expression data, chips and microarrays [7]. Three gene expression profiles in primary melanoma and normal skin were downloaded from GEO, including GSE15605, GSE46517 and GSE114445 (Affymetrix GPL570 platform, Affymetrix GPL96 platform, Affymetrix GPL570 platform, respectively Identification of DEGs GEO2R (http://www.ncbi.nlm.nih.gov/geo/geo2r) is regarded as an interactive online tool designed to compare two or more datasets in a GEO series for the purpose of DEGs 5 identification across experimental conditions. The DEGs between primary melanoma and normal skin were identified using GEO2R with the threshold of |logFC|>1 and P value < 0.05 which were considered of statistically significance. For the next step, the online Venn software was applied to detect the intersection DEGs among three datasets.

Enrichment analysis for DEGs
To elucidate the biological functions of the overlapping DEGs, we performed functional annotation and pathway enrichment analysis via DAVID (The Database for Annotation, Visualization and Integrated Discovery, http://david.ncifcrf.gov/,version 6.8) online tool [8]. P value < 0.05 was considered statistically significant. GO enrichment and KEGG pathway of DEGs was analyzed and displayed using bubble charts.
PPI network construction and module analysis STRING (Search Tool for the Retrieval of Interacting Genes, http://string-db.org, version 11.0) online database was used to predict the PPI network which may further explain the mechanisms of the occurrence and progression of diseases [9]. By using STRING database, PPI network of DEGs was analyzed and an interaction with a combined score > 0.4 was recognized as statistical significance. The plug-in MCODE (Molecular Complex Detection) app of Cytoscape (an public bioinformatics software, version 3.7.2) is constructed for clustering a network based on topology to determine intensively connected regions [10,11]. The PPI network was plotted with the application of Cytoscape and the most significant module in the PPI network was narrowed down using MCODE with the following criteria: degree cutoff = 2, node score cutoff = 0.2, k-core = 2, max depth = 100.

Hub genes selection and analysis
The pathway analysis of hub genes was performed and visualized by ClueGO (version 2.5.4) and CluePedia (version 1.5.4), which was a functional extension of ClueGO and a 6 plug-in of Cytoscape [12]. GEPIA (http://gepia.cancer-pku.cn/index.html) is a web-based tool to provide key interactive and customizable functions based on TCGA (The Cancer Genome Atlas) and GTEx (Genotype-Tissue Expression) data [13]. Furthermore, GEPIA was applied for validating the gene expression. The overall survival of hub genes was analyzed using Kaplan Meier-plotter online tool in GEPIA which is commonly applied for assessing the effect of genes on survival based on TCGA database.

The correlation between hub genes and immune cell infiltration
The TIMER (https://cistrome.shinyapps.io/timer/) database is a practical platform for the systematical analysis of immune infiltrates in diverse cancer types [14]. We used the TIMER Gene module to explore the correlation between the abundances of six immune infiltrates and the level of hub genes.

Drug-gene interaction predictive analysis
The Drug-Gene Interaction database (DGIdb2.0) mines existing resources and generates assumptions about how genes are therapeutically targeted or prioritized for drug development [15]. In this study, DGIdb2.0 was used to predict drugs based on the module hub genes. The parameters were set as: preset filters: FDA approved; antineoplastic; all the default. All the drug-gene relationship pairs related to the hub genes were predicted, and Cytoscape was used to construct the network map.
Transcription factor analysis of seven significant hub genes Transcription factor regulation network(including related lncRNA, target miRNA and protein-protein interactions)were constructed among CCL4, CCL5, NMU, GAL, CXCL9,

Hub gene analysis
The KEGG pathway analysis of hub genes by using ClueGO (Fig. 3). The most important pathways include interleukin-10 signaling, peptide ligand-binding receptors, chemokine receptors bind chemokine, toll-like receptor signaling pathway and cytosolic DNA-sensing pathway. Kaplan-Meier survival curves showed significant overall survival outcomes of seven significant hub genes (P values < 0.05) in Fig. 4. Subsequently, we chose these seven genes (CCL4, CCL5, NMU, GAL, CXCL9, CXCL10, CXCL13) for further analysis.
GEPIA database verified that the expression of CCL4, CCL5, CXCL9, CXCL10 and CXCL13 presented significant difference (P values < 0.05) between melanoma samples and normal samples (Fig. 5). Datasets from TCGA was applied to validate the correlations between the seven hub genes and clinical stages. Five of those seven hub genes were found to be related to clinical stages (P values < 0.05), including CCL4, CCL5, CXCL9, CXCL10 and CXCL13 (Fig. 6).

The correlation between hub genes and immune cell infiltration
We analyzed the correlation between hub genes and immune cell infiltration via TIMER.
Our results revealed a positive correlation between gene expression (CCL4, CCL5, CXCL9, Transcription factor analysis of seven significant hub genes Transcription factor regulation network(including related lncRNA, target miRNA and protein-protein interactions) were constructed among CCL4, CCL5, NMU, GAL, CXCL9, CXCL10 and CXCL13 in Fig. 9. Significant nodes were marked in different colors in line with significant hub genes.

Discussion
Skin is the largest organ of human beings. There are about 1,500 melanocytes in human epidermis per square millimeter, equivalent to nearly 3 billion melanocytes in an ordinary human skin [16]. The incidence of melanoma continues to increase every year. Because of its highly aggressive and treatment-resistant, melanoma is not projected to reduce the death rate in the next few years [17]. Melanoma is associated with a huge burden of genetic alterations [18]. Research on genetic characterizations of melanoma has made a great progress [19,20]. BRAF, NRAS and NF1 are related to initiation and proliferation of tumor cells [21]. CDKN2A and TERT are associated with progression and cell cycle control [19]. While PTEN and KIT play important roles in advanced progression and metabolism [22]. Targeted therapeutic drugs for the above mutations have been rapidly developed and used for the treatment of melanoma with an improvement in overall survival [23]. The classical progression models of melanoma are from nevus to primary melanoma, to invasive melanoma [6]. However, the lack of research on different periods limits people's understanding of the key genes of melanoma in these periods, especially the initial stage and metastasis stage. Most of the existing targeted drugs are aimed at patients with advanced melanoma, which has the extensive genetic heterogeneity and genomic instability [24]. These two reasons may lead to the emergence of drug resistance,as seen in targeted therapeutic approaches for BRAF inhibitors [25]. Cancer molecular prevention suggests that we can interrupt the prime driver at an early stage of tumor [26]. The shortage of intermediate biomarkers has been a barrier for the development of new molecular target drugs. In addition, melanoma is a highly immunogenic cancer. It is of great value to clarify the potential mechanism of immune activation and immunosuppression in tumor microenvironment. Microarray technology enables exploration of immune mechanism and new targets for melanoma.
In the present study, a total of 308 DEGs and 12 hub genes were identified from three microarray datasets. GO analysis indicated that changes in DEGs were significantly Subsequently, we chose these seven genes for further analysis. GEPIA verified that the expression these genes except NMU and GAL has significant difference (P values < 0.05) between tumor tissues and normal. In addition, five of these seven hub genes were found to be related to clinical stages (P values < 0.05), including CCL5, CXCL13, CCL4, CXCL9, CXCL10. In recent years,a variety of new biomarkers, such as microRNAs (miRNAs) and lncRNA, have become powerful new tools for diagnosis and monitoring in cancer patients 11 [27]. Therefore, we constructed the related miRNA and lncRNA networks.
C-X-C motif ligand 9,10 (CXCL9, CXCL10) belong to CXC family chemokine. They are predominantly expressed by immune cells such as macrophages and T cells. CXCL9, CXCL10 exert their function through binding to C-X-C motif receptor CXCR3 that was highly expressed on the activated T cells [28]. Research shows that CD8 + T cell infiltration is CXCR3 dependent [29]. CD8 + T cell plays important role in tumor microenvironment, which inhibits the proliferation and metastasis of tumor cells. In our study, high expression of CXCL9 and CXCL10 associated with better overall survival in melanoma patients. What's more, CXCL9 and CXCL10 were found to be highly expressed in stage I melanoma compared to other stages. A recent research found that the CXCL10-CXCR3 axis may relate to melanoma brain metastasis [30]. Additionally, the CXCR3 ligands enhanced the response rates to immune checkpoint blockade (anti-PD-1, anti-CTLA-4) by recruiting T cells [31]. Kim J et al. found that CXCR3-deficient natural killer cells fail to migrate to B16F10 melanoma cells [32]. These results reveal the importance of manipulating this axis, which provide us opportunities to manipulate chemokine expression for tumor therapy and prevention.
C-X-C motif ligand 13 (CXCL13) also belongs to CXC family. There are very few studies on the relationship between CXCL13 and melanoma. One research suggested that CXCL13 is closely interact with NRAS, BARF and MITF in MAPK signaling pathway, which have implications for melanoma genesis and metastatic progression [33]. In our study, CXCL13 showed similar outcome in overall survival and expression with CXCL9, CXCL10. Pathway analysis suggested they were involved in peptide ligand-binding receptors and chemokine receptors bind chemokine. Therefore, CXCL13 may play the same role as CXCL9, CXCL10 in the occurrence and development of melanoma. Our results revealed that CXCL9, CXCL10, CXCL13 are positively associated with immune cell infiltration. Interestingly, 12 CXCL10 and CXCL2 were found in the predicted drug-gene interaction, providing new thoughts for melanoma treatment. Whether to detect the level of chemokine expression in peripheral blood can help us diagnose melanoma and its stage is unknown. Further external experiments are needed to verify this conjecture.
C-C motif ligand CCL4, CCL5 belong to CC family chemokine. CCL4, CCL5, CXCL9, CXCL10 were confirmed to be preferentially expressed in tumors that contained T cells [34]. This suggests that CCL4, 5 and CXCL9,10 may act in the same pathway in T-cell recruitment.
The experimental results carried out by Michelle H et al. showed that expression of CXCR3 ligands and CCL5 in chemotherapy tumors did lead to a synergistic increase in T-cell infiltration [35]. In our study, these four genes were mainly enriched in cytokine-cytokine receptor interaction, chemokine signaling pathway and toll-like receptor signaling pathway. The expression of CCL4, CCL5 showed positive correlation with multiple immune cells infiltration (B-cell, CD8 + T cells, CD4 + T cells, macrophages, Neutrophils, Dendritic cells). It is also been reported that CCL5 released by tumor cells acts its function via paracrine signaling to attract NKs to the tumor bed [36]. However, previous studies showed that tumor growth is delayed in CCR5-/-mice, associated with reduced tumor regulatory T cells infiltration [37]. Edwin L et al. found that a selective increase in CXCL10, CCL4, and CCL17 may contribute to aggressive development of brain metastasis [38]. In the present study, high expression of CCL4 and CCL5 associated with better overall survival and were found to be closely related to stage I melanoma patients compared to other stages. Taken together, we assumed that CCL4 and CCL5 may be potent biomarkers and targets for melanoma.
Neuromedin U (NMU), a neuropeptide belonging to the neuromedin family, is involved in stress, obesity, immune, and energy metabolism [39]. A recent study found the YAP1-NMU axis is associated with pancreatic cancer progression and poor prognosis [40]. High 13 expression of NMU was found to be related to reginal metastasis of head and neck squamous cell carcinoma [41]. NMU upregulation in the breast cancer tissues was proposed as prognostic biomarker for poor outcome, but only in HER2-overexpressing tumor [42]. Moreover, patients with NSCLC (non-small-cell lung carcinoma) and NMUpositive tumor showed significantly shorter survival times than patients with NMUnegative tumor [43]. All these findings suggest that NMU is associated with cancer development and NMU may serve as a new biomarker of poor prognosis. Furthermore, NMU has been proved to be an important resistance-enhancing factor. Overexpression of NMU in sensitive cells makes it resistant to the tested drugs, while NMU silencing sensitized resistant cells [42]. Therefore, silencing NMU can be an effective way for cancer treatment. Researchers found the results in proliferation and survival of NMU appeared to be cancer-type dependent [44,45]. However, research on melanoma and NMU is still lacking. In our study, high expression of NMU result in poor overall survival. NMU receptors are expressed on macrophages and endothelial cells, indicating that NMU is highly related to the tumor microenvironment [46]. Therefore, NMU is a potential tumor growth and progression biomarker in melanoma and silencing NMU may be an effective way for melanoma treatment. More research is needed to reveal its role in melanoma microenvironment.
Galanin and GAMP prepropeptide (GAL) is a neuropeptide, which is found to be expressed in several cancers such as gastrointestinal cancer, head and neck squamous cell carcinoma, brain tumor and salivary duct carcinoma [47]. Galanin acts its function through binding to three receptors (GALR1, GALR2, GALR3). The activation of GALR1 is usually anti-proliferation, while the activation of GALR2 may have the effect of promoting proliferation or anti-proliferation [48]. GALR1 and GALR2 inhibit proliferation of head and neck squamous cell carcinoma (HNSCC) cells. Hypermethylation of GALR1, GALR2 is associated with poor survival and high recurrence in HNSCC [49]. Therefore, galanin and galanin receptors (GALRs) are gradually used as new biomarkers for prognosis in HNSCC.
Another research found that the down-regulation of Galanin and GALR1 is related to development of gastric cancer [50]. Besides, Galanin and GAL receptors have been found to be expressed in human skin [51]. Few studies have proved the role of GAL in melanoma. In our results, overall survival of samples with high expression of GAL showed in poor outcome. Four drugs interacted with GAL in the predicted drug-gene interaction.
The relationship between GAL, GALRs and melanoma remains unknown, and more research is needed to verify it.
In our study, we used ClueGO to investigate pathway in hub genes. The most important pathways include interleukin-10 signaling, peptide ligand-binding receptors, chemokine receptors bind chemokine, toll-like receptor signaling pathway, cytosolic DNA-sensing pathway. Interleukin-10 signaling mediates immunosuppressive via its immune-regulatory function, while chemokine and its receptors usually attract immune cells such as CD8 + T cell in tumor microenvironment for immune activation [52]. Melanoma is a highly immunogenic cancer. It is valuable to clarify the potential mechanism of immune activation and immunosuppression in tumor microenvironment. Our results revealed a strong positive correlation between gene expression (CCL4, CCL5, CXCL9, CXCL10 and  [53]. Altogether, these hub genes and pathways may play critical role in formation, prognosis and therapy of melanoma, while further investigation is necessary. 15

Conclusion
In conclusion, bioinformatic analysis was used to identify DEGs and hub genes between melanoma and normal skin that may help us understand the potential pathogenesis and prognosis factors in melanoma. CCL4, CCL5, NMU, GAL, CXCL9, CXCL10 and CXCL13 were significantly associated with prognosis and tumor microenvironment. Function and pathway analysis of hub genes revealed the part mechanism in formation, prognosis and therapy of melanoma. However, whether these hub genes can be regard as new targets or prognosis biomarkers for melanoma still need further external experiments to verify.

Consent for publication
Not applicable.

Availability of data and materials
The datasets supporting the conclusions of this article are available in the GEO (https://www.ncbi.nlm.nih.gov/geo/) and TCGA repository (https://www.cancer.gov/aboutnci/organization/ccg/research/structural-genomics/tcga/studied-cancers). The KEGG pathway analysis of hub genes via ClueGO (P< 0.05). a, the most significant pathway and related genes. b, the number of genes involved in each pathway. c, the classification pie charts were listed as follows: 80.00% terms belong to interleukin-10 signaling and 20.00% terms belong to peptide ligandbinding receptors.