Identification of biomarkers and pathways associated with immune infiltration and prognosis in melanoma patients


 Background Skin cutaneous melanoma (SKCM) is one of the most malignancy and aggressive cancers, causing about 72% of deaths in skin carcinoma. Although extensive study has explored the mechanism of recurrence and metastasis, the tumorigenesis of melanoma remains unclear. Exploring the tumorigenesis mechanism may help to identify prognostic biomarkers that could serve to guide cancer therapy . Methods We identified differentially expressed genes (DEGs) between patients with primary melanoma and normal skin in three data sets including GSE46517, GSE15605, GSE114445. Functional annotation including KEGG pathway and gene ontology (GO) enrichment analysis was performed by DAVID. Subsequently, protein ‐ protein interaction network of DEGs was developed and the most significant module was constructed as hub genes for further study. ClueGO was used to investigate significant pathways in hub genes. The Kaplan–Meier method was performed to assess prognostic genes. ROC curves were used to describe diagnosis value of hub genes. Then, genes expression and clinical characteristics were validated by TCGA profiles. In addition, we explored the correlation between prognostic genes and immune cell infiltration. Results A total of 308 DEGs and 12 hub genes were identified. GO enrichment results demonstrated a correlation between DEGs and the immune response, inflammatory response and transcriptional control. Hallmark pathways of hub genes including interleukin-10 signaling, chemokine receptors bind chemokine signaling ， peptide ligand-binding receptors. Survival analysis and ROC curves suggested that CCL4, CCL5, NMU, CXCL9, CXCL10, CXCL13 were independent prognostic factors (except GAL). Furthermore, prognostic genes were highly expressed in tumor tissue and related to pathological stages and Breslow depth. In addition, the expressions of CCL4, CCL5, CXCL9, CXCL10, CXCL13 were positively correlated with six immune cell infiltration (B-cell, CD8+ T cells, CD4+ T cells, macrophages, Neutrophils, Dendritic cells) and 28 types of TILs such as activated CD8 T cells, regulatory T cells, natural killer T cells while NMU and GAL were negatively correlated with it. Conclusion In summary, this study identified significant hub genes and pathways closely associated with immune infiltrations and prognosis in SKCM patients, which might help us evaluate underlying carcinogenesis progress from skin to melanoma.

3 carcinogenesis progress from skin to 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]. It is well known that melanoma is associated with multiple risk factors especially the sun exposure [4]. The general progression models of melanoma are from melanocyte to melanoma in situ, to invasive melanoma [5]. However, extensive research has explored the mechanism of recurrence and metastasis, the tumorigenesis of melanoma remains unclear. Cancer cells act as antigens to induce immune infiltration while immune cells are attracted to tumor environment through chemotaxis to play immune defense roles [6]. In addition, increasing researches showed that immune escape is an important factor in tumorigenesis [7]. Emerging immunotherapy strategies have been applied in melanoma such as anti-PD-1, anti-PDL-1, anti-CTLA4 and MAGE-A3 [8]. However, only a small group of patients benefit from it and 50%-60% showed poor prognosis [9]. The tumor microenvironment (TME) consists of a variety of immune cells and stromal cells, including immune cells, fibroblasts, endothelial cells, extracellular matrix, cytokines, chemokines and receptors [10]. Previous researches have suggested that the increasing of tumor-infiltrating lymphocytes (TILs) was associated with better prognosis, reduced lymph node metastasis and a lengthier disease-free survival [11]. Therefore, exploration of the mechanism of melanoma and identification the immune-related biomarkers with better prognosis will greatly help us evaluate underlying carcinogenesis progress of melanoma, which might improve the precision medicine and immunotherapies.
In the present study, we analyzed gene expression profiles to identify hub genes and pathways associated with immune infiltrations and prognosis. Our results may reveal potential therapeutic targets and provide insights into the molecular mechanisms of melanoma.

Gene expression datasets
Expression profiling of melanoma patients with clinical information was obtained from the Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo) [12]. Among

Identification of DEGs
The DEGs between primary melanoma and normal skin were identified using GEO2R (http://www.ncbi.nlm.nih.gov/geo/geo2r) based on limma package 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 (http://bioinformatics.psb.ugent.be/webtools/Venn/) was applied to detect the overlap DEGs among three datasets.

Function 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 [13]. 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 hub genes selection
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 of DEGs and analyze the functional interactions 5 between proteins [14]. An interaction with a combined score > 0.4 was recognized as statistical significance. The molecular interaction networks were visualized using the Cytoscape (version3.4.0) and the most significant hub genes in the PPI network was narrowed down using MCODE (Cytoscape) with the following criteria: degree cutoff =2, node score cutoff = 0.2, k-core = 2, max depth = 100 [15,16].

Hub genes significant pathways 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 plug-in of Cytoscape [17]. p-value < 0.01 was considered statistically significant.

Prognostic and diagnostic value assessment of hub genes
The Kaplan-Meier method was used to analyze survival differences between groups based on TCGA cohort. The primary endpoint was Overall survival (OS). The follow-up duration was estimated and illustrated by the Kaplan-Meier method with 95% confidence intervals (95%CI) and log-rank test in separating curves. ROC cures were constructed to assess the diagnostic accuracy of the prognostic genes and statistical analysis was performed using SPSS 22.0. A value of P<0.01 was considered statistically significant.

Validation of prognostic genes expression and clinical characteristics
A total of 472 melanoma patients with available RNA sequence data from TCGA database and GEPIA (http://gepia.cancer-pku.cn/index.html) database were used to verify gene expression and clinical characteristics [18]. The overall statistical expression difference of pathological stages or Breslow depth was measured using One-way ANOVA test. Student's t tests were used to compare differential transcriptional expressions levels of prognostic genes between paired tumor and normal samples. Pvalues less than 0.05 were considered significant. The Human Pathology Atlas project (https://www.proteinatlas.org) contains immunohistochemistry (IHC) data using a tissue microarraybased analysis on 44 different normal tissue types, and proteome analysis of 17 major cancer types [19]. In this study, representative proteins expressions of IHC images of prognostic genes were detected in melanoma and normal tissues in Human Protein Atlas.

Immune infiltration analysis
Prognostic gene expression data were used to measure the abundance of six types of infiltrating immune cells (B cells, CD4+T cells, CD8+T cells, neutrophils, macrophages, and dendritic cells) in melanoma patients using the Tumor Immune Estimation Resource (TIMER) algorithm [20]. An integrated repository portal for tumor-immune system interactions (TISIDB, http://cis.hku.hk/TISIDB/index.php) was utilized to examine tumor and immune system interactions in 28 types of TILs across human cancers [21]. Spearman's test was used to measure correlations between prognostic genes and TILs. All hypothetical tests were two-sided and p-values <0.05 were considered significant in all tests.

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 [22]. 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.

Identification of DEGs
After standardization and identification of the microarray results, DEGs (4640 in GSE15065,1096 in GSE46517 and 1895 in GSE114445) were selected. The overlap among three datasets included 308 genes as shown in the Venn diagram ( Figure 2a) (Additional file 1) between primary melanoma and normal skin.

GO and KEGG Enrichment Analysis of the DEGs
Functional and pathway enrichment of DEGs was performed by DAVID, displayed in bubble charts (p<0.05, Figure 3). GO analysis indicated that changes in biologic processes of DEGs were significantly enriched in the inflammatory response, immune response, regulation of transcription. In addition to cellular component, the overlapping DEGs were particularly enriched in extracellular region, extracellular space and extracellular exosome. Changes in molecular function were mostly enriched in chemokine activity, transcriptional activator activity, protein binding and CXCR3 chemokine receptor binding. At last, KEGG pathway analysis demonstrated that DEGs played pivotal roles in pathways in cytokine-cytokine receptor interaction, chemokine signaling pathway and pathways in cancer. These results demonstrated a correlation between DEGs and immune response and inflammatory response.

PPI network construction and hub genes analysis
The PPI network was constructed via Cytoscape, including 247 nodes and 792 edges ( Figure 2b).

Validation of prognostic genes expression and clinical characteristics
In the present study, 472 melanoma patients from the TCGA cohort were used to analyze correlation between prognostic gene expression and pathological stages, Breslow depth, which are useful prognostic indicators for melanoma. As shown in figure 7, CXCL9, 10, 13 and CCL4, 5 were significantly overexpressed in stage Ⅰ and decreased in subsequent stages with P value<0.001.

Immune infiltration analysis
After determining the prognostic value, we performed correlation analysis between prognostic genes and immune infiltration level for melanoma. There was a positive correlation between CXCL9

Discussion 9
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 [23]. The incidence of melanoma continues to increase every year. Once melanoma spreads through the dermis, the prognosis is poor and melanoma is not projected to reduce the death rate in the next few years [24]. A series of therapeutic methods, such as radiotherapies, chemotherapies, targeted therapies and immunotherapies have enhanced advanced patient survival [25]. However, many patients that are being administered these therapies demonstrate a low durable response, drug resistance and poor prognosis [26]. It is imperative to mine the mechanism in melanoma and to detect effective prognostic biomarkers that could serve to guide cancer therapy.
In the present study, 308 DEGs were identified and related to inflammatory response, immune response and regulation of transcription. Moreover, the PPI network were constructed and nodes with a high connectivity degree in the modules were consider as hub genes. Hallmark pathways in hub genes include interleukin-10 signaling, chemokine receptors bind chemokine signaling and peptide ligand-binding receptors. STAT3 is essential for the action of IL-10,IL-6, and IL-6 plays a central role in 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 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 [43]. A recent study found the YAP1-NMU axis is associated with pancreatic cancer progression and poor prognosis [44]. High expression of NMU was found to be related to reginal metastasis of head and neck squamous cell carcinoma [45]. NMU upregulation in the breast cancer tissues was proposed as prognostic biomarker for poor outcome, but only in HER2overexpressing tumor [46]. Moreover, patients with NSCLC (non-small-cell lung carcinoma) and NMUpositive tumor showed significantly shorter survival times than patients with NMU-negative tumor [47]. 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 [46]. 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 [48,49]. 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 [50]. 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 This study has several limitations. First, the microchip data were acquired from a relatively small number of melanomas deposited in a public database, and only 472 patients with corresponding transcriptome data were acquired from the TCGA cohort. Second, further studies, including in vivo and in vitro experiments as well as validated cohorts, are needed to verify the absoluteness of these findings.

Ethics approval and consent to participate
Not applicable.

Consent for publication
Not applicable.     The KEGG pathway analysis of hub genes via ClueGO (P< 0.01). 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 ligand-binding receptors.

Figure 5
Univariate survival analysis of the hub genes was performed using the Kaplan-Meier curve.
The 7 genes of 12 hub genes showed significant difference in overall survival (P<0.01).
28 Figure 9 GEPIA database based on TCGA cohort indicated that CXCL9, 10, 13 and CCL4, 5, NMU were significantly high expressed in tumor tissue compared to normal p<0.01.
29 Figure 10 IHC results revealed CCL4, 5, CXCL9, 10, 13 and NMU were detected in melanoma tissues while not detected in normal tissues using online database.
30 Figure 11 a, There was a positive correlation between CXCL9 expression and the infiltration of B cells

Supplementary Files
This is a list of supplementary files associated with this preprint. Click to download.