Identification of immune-related biomarkers associated with tumorigenesis and prognosis in cutaneous melanoma patients

Background Skin cutaneous melanoma (SKCM) is one of the most malignant 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 cutaneous melanoma remains unclear. Exploring the tumorigenesis mechanism may help identify prognostic biomarkers that could serve to guide cancer therapy. Method Integrative bioinformatics analyses, including GEO database, TCGA database, DAVID, STRING, Metascape, GEPIA, cBioPortal, TRRUST, TIMER, TISIDB and DGIdb, were performed to unveil the hub genes participating in tumor progression and cancer-associated immunology of SKCM. Furthermore, immunohistochemistry (IHC) staining was performed to validate differential expression levels of hub genes between SKCM tissue and normal tissues from the First Affiliated Hospital of Soochow University cohort. Results A total of 308 differentially expressed genes (DEGs) and 12 hub genes were found significantly differentially expressed between SKCM and normal skin tissues. Functional annotation indicated that inflammatory response, immune response was closely associated with SKCM tumorigenesis. KEGG pathways in hub genes include IL-10 signaling and chemokine receptors bind chemokine signaling. Five chemokines members (CXCL9, CXCL10, CXCL13, CCL4, CCL5) were associated with better overall survival and pathological stages. IHC results suggested that significantly elevated CXCL9, CXCL10, CXCL13, CCL4 and CCL5 proteins expressed in the SKCM than in the normal tissues. Moreover, our findings suggested that IRF7, RELA, NFKB1, IRF3 and IRF1 are key transcription factors for CCL4, CCL5, CXCL10. In addition, the expressions of CXCL9, CXCL10, CXCL13, CCL4 and CCL5 were positively correlated with infiltration of six immune cells (B cell, CD8+T cells, CD4+T cells, macrophages, neutrophils, dendritic cells) and 28 types of TILs. Among them, high levels of B cells, CD8+T cells, neutrophils and dendritic cells were significantly related to longer SKCM survival time. Conclusion In summary, this study mainly identified five chemokine members (CXCL9, CXCL10, CXCL13, CCL4, CCL5) associated with SKCM tumorigenesis, progression, prognosis and immune infiltrations, which might help us evaluate several immune-related targets for cutaneous melanoma therapy.


Background
Skin cutaneous melanoma (SKCM) accounts for only 2% of total skin cancers. However, due to its high degree of malignancy and invasiveness, it causes over 72% of deaths in skin carcinoma [1]. The incidence of cutaneous malignant melanoma continues 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 SKCM 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 cutaneous melanoma remains unclear.
In the present study, we analyzed the differentially expressed genes (DEGs) between primary melanoma and normal skin to explore the potential tumorigenesis mechanism of SKCM. Our results mainly identified several chemokine family members (CXCL9, CXCL10, CXCL13, CCL4, CCL5) which were found related to better overall survival (OS) in SKCM patients. Chemokine family members are a group of low-molecular weight cytokines which were involved in many biological processes including angiogenesis, tumor development and metastasis, and the migration of leukocytes [6,7]. In addition, we found that their expression levels were positively associated with infiltration of immune cells (CD4 + T, CD8 + T, B-cell, macrophages, neutrophils, dendritic cells) and tumor infiltrating lymphocytes (TILs). These immune infiltration cells play important roles in tumor microenvironment and can directly or indirectly regulate tumor immunity and modulate tumor immunological for antitumor effects [7,8]. Therefore, our results may identify several immune-related biomarkers that may serve to guide SKCM therapy.

Patients and variables
A total of 46 melanoma and 46 normal tissues were obtained from 92 patients at the Department of Burn and Plastic Surgery, the First Affiliated Hospital of Soochow University (FAHSU, Suzhou, China) from March 2015 to August 2019. None of the patients had received radiotherapy or chemotherapy before operation. Tissue samples, including cutaneous melanoma and normal tissue, were collected during surgery and fixed in 4% paraformaldehyde, available from FAHSU tissue bank. Clinical data was available to obtain from hospital records. This research was supported by the Independent Ethics Committee (IEC) of the FAHSU and all patients were well informed of storing and upcoming use of their resected specimens for further research purposes.

GEO and TCGA datasets
Expression profiling of SKCM patients with clinical information was obtained from the Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo) [9]. Among the inclusion criteria were (a) diagnosis of patients with primary melanoma (PM) and normal skin (NS), (b) detection of gene level in tissue or blood samples. Exclusion criteria included: (a) clinical data without survival time and outcome, and (b) datasets with small sample sizes (n < 15). Finally, three datasets were eligible: accession numbers GSE15605 (46 PM samples and 16 NS samples), GSE46517 (31 PM samples and 8 NS samples) and GSE114445 (16 PM samples and 6 NS samples). For validation, gene expression profiles of 472 melanoma patients were downloaded from the TCGA data portal (https ://tcga-data.nci.nih.gov/tcga/) [10]. Clinical characteristics of patients were also obtained, including age, Clark level, Breslow depth, survival time, and outcome. All of the clinicopathological characteristics of SKCM patients were shown in Table 1.

DAVID database
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 .ncifc rf.gov/,versi on 6.8) online tool [11]. P-value < 0.05 was considered as the cutoff value. GO enrichment and KEGG pathway results were visualized as bubble charts by R software.

STRING and cytoscape
STRING (http://strin g-db.org, version 11.0) database was used to predict the PPI network of DEGs and analyze the interactions between proteins [12]. An interaction with a combined score > 0.4 was recognized as statistical significance. The molecular interaction networks were visualized using the Cytoscape (version 3.7.0) and the most significant model in the PPI network was narrowed down by MCODE, with the following criteria: degree cutoff = 2, node score cutoff = 0.2, k-core = 2, max depth = 100 [13,14].

Functional enrichment
Metascape (https ://metas cape.org) was used to further verify the function enrichment of hub genes [15]. P < 0.05 was set as the cutoff value. Hub genes pathway analysis was performed and visualized by ClueGO (version 2.5.4) and CluePedia (version 1.5.4), the plug-in of Cytoscape [16]. P-value < 0.01 was considered statistically significant.

GEPIA database
GEPIA (http://gepia .cance r-pku.cn/index .html) is an analysis tool containing RNA sequence expression data of 9736 tumors and 8587 normal tissue samples [17]. In this study, GEPIA was used to perform survival analysis based on the data from TCGA database. Student's t test was used to analyze the correlation between the expression and pathological stages. P value < 0.05 was considered statistically significant.

Immunohistochemistry (IHC)
Protein expression levels of hub genes were measured using IHC staining and rabbit polyclonal anti-CCL4 antibody (ab235978), anti-CCL5 antibody (ab9679), anti-CXCL9 antibody (ab9720), anti-CXCL10 antibody (ab9807), anti-CXCL13 antibody (ab272874). Positive or negative staining of a certain protein in one FFPE slide was independently assessed by two experienced pathologists and supervised by a clinician. Based on the staining intensity level (no staining, weak, moderate and strong staining), the score was ranging from 0 to 3, as previous described. The staining extent was graded from 0 to 4 for the coverage percentage of immunoreactive tumor cells (0%, 1-25%, 26-50%, 51-75%, 76-100%). The overall IHC score grading from 0 to 12 was evaluated according to the multiply of the staining intensity and extent score. Negative staining represented grade 0 to 4 and positive staining from 5 to 12 for each sample.

cBioPortal database
The cBioPortal database (www.cbiop ortal .org), a comprehensive web resource, can visualize and analyze multidimensional cancer genomics data [18]. The genetic alterations of prognostic genes were obtained from cBio-Portal based on 475 SKCM samples in TCGA.

TRRUST database
TRRUST (https ://www.grnpe dia.org/trrus t/) is a useful predict tool for human, and mouse transcriptional regulatory networks [19]. The TRRUST database can provide information on how these interactions are regulated, containing 8444 transcription factor (TF)-target regulatory relationships of 800 human TFs.

Immune infiltration analysis
Prognostic gene expression data was 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 SKCM patients using the Tumor Immune Estimation Resource (TIMER) algorithm [20]. Then, an integrated repository portal for tumorimmune system interactions (TISIDB, http://cis.hku.hk/ TISID B/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 statistically significant.

The drug-gene interaction database
DGIdb (version 2.0) is an open-source project that help users mine existing resources and generate assumptions about how genes are therapeutically targeted or prioritized for drug development [22]. The parameters were set as: preset filters: FDA approved; antineoplastic; all the default.

Statistical analysis
The expression heatmap and correlation coefficient were analyzed and visualized by pheatmap and corplot packages in R software. ROC curve of prognostic genes was performed by SPSS 22.0. P-values < 0.05 were considered statistically significant in all tests.

DEGs selection
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 ( Fig. 1a, Additional file 1: Table S1) between primary melanoma and normal skin.

GO and KEGG enrichment analysis of the DEGs
Functional and pathway enrichment of DEGs was performed by DAVID, then displayed in bubble charts by R software (P < 0.05, Fig. 2). GO analysis indicated that changes in biologic processes were significantly enriched in the inflammatory response, immune response, regulation of transcription. As for cellular component, 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. 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 revealed that immune response and inflammatory response play important roles in SKCM tumorigenesis.

Function analysis of hub genes
As expected, functional annotation obtained from Metascape suggested that hub genes were mainly enriched in regulation of T cell chemotaxis, peptide ligand-binding receptors and chemokine receptors bind chemokine signaling (Fig. 3a, b, P < 0.001). Similarly, ClueGO (Fig. 3c, P < 0.01) revealed that the most involved pathways were interleukin-10 signaling, peptide ligand-binding receptors and chemokine receptors bind chemokine signaling.

The prognostic and diagnostic value of hub genes in SKCM patients
Using the Kaplan-Meier method, the prognostic values of the hub genes in SKCM patients were determined. Seven hub genes were significantly associated with OS as shown in Fig. 4.

Validation the aberrant expression of prognostic genes and their clinical characteristics
Based on data from GEPIA database, all of the prognostic genes were found differential expressed in SKCM vs normal skin tissue (P < 0.05, Fig. 6). correlated to SKCM pathological stages (Fig. 7). Among them, CXCL9, CXCL10, CXCL13, CCL4, CCL5 were higher expressed in stage I and decreased in the subsequent stages. Subsequently, all of the samples from TCGA were ordered according to the expression level of CXCL9, and statistical analysis was conducted to compare the head 80 with the tail 80 samples. The heatmap (Fig. 8b) and statistical results suggested that higher expression levels (tail 80 samples) of CXCL9, CXCL10, CXCL13, CCL4, CCL5 were related to more patients of lower Breslow depth (0-1.5 mm, 51.2%), BRAF mutation (60%) and lower Clark levels (I-III 45% and IV-V, 55%), reduced T4 stage (T4, 25%). With the decrease levels of these five chemokines (head 80 samples), there were more patients with higher Breslow depth (0-1.5 mm, 15% and > 1.5 mm, 85%), BRAF WT (wild type, 56.3%), higher Clark levels (I-III 17.5% and IV-V, 82.5%) and T4 stage (T4, 47.5%).
All of the data above suggested that five chemokine family members play critical roles in cutaneous melanoma tumorigenesis and progression.

Genetic alteration, correlation coefficient and key transcription factors analyses of prognostic genes in patients with SKCM
As a result, there were nearly 6% (CXCL9), 5% (CXCL10), 3% (CXCL13), 4% (CCL4), 5% (CCL5), 3% (NMU), 4% (GAL) of SKCM samples had genetic alteration. The most common genetic change among five a b c Fig. 3 The function analysis of hub genes (P < 0.01). a, b The functions of hub genes were mainly enriched in regulation of T cell chemotaxis, peptide ligand-binding receptors and chemokine receptors bind chemokine signaling. c The most significant pathway and related genes chemokines (CXCL9, CXCL10, CXCL13, CCL4, CCL5) was enhanced mRNA expression (Fig. 8a). While genetic change in GAL and NMU were mainly related to amplification. Then, we assessed the correlation coefficients in prognostic genes. As expect, significant positive correlations were observed between five chemokines (Fig. 8c). Among them, CXCL9 and CXCL10 had the highest positive correlation with a Spearman's correlation coefficient of 0.92. CXCL9 exhibited the tightest association with all other chemokines, with a median Spearman's correlation coefficient of 0.85.

Discussion
Skin is the largest organ of human beings. There are about 1500 melanocytes in human epidermis per square millimeter, equivalent to nearly 3 billion melanocytes in an ordinary human skin [23]. The incidence of SKCM 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]. Therefore, more therapeutic targets and prognostic biomarkers must be identified. In our study, 308 DEGs were identified between primary SKCM and normal skin, and the functional annotation indicated that inflammatory response and immune response were closely associated with SKCM tumorigenesis. KEGG pathways in hub genes included IL-10 signaling, chemokine receptors bind chemokine signaling and peptide ligand-binding receptors. STAT3 is essential for the action of IL-10 and IL-6. IL-6 plays a central role in melanoma development and enhances the IL-10 production via STAT3-dependent signaling [27]. IL-6,10/ JAK-STAT3 pathways are found in many carcinomas, and their hyperactivation was generally related to unfavorable clinical prognosis [28]. IL-6 and IL-10 create a favorable environment to support tumorigenesis by increasing the cancer cell proliferation, angiogenesis, metastasis and contribute to enhancing the immune suppression [28]. Smith et al. found that IL-10 directly inhibits CD8 + T cell function by enhancing N-Glycan branching to decrease antigen sensitivity [29]. However, chemokine and its receptors usually attract immune cells such as CD8 + T cell and NKs in the tumor microenvironment for immune activation [30]. Apparently, these two pathways play opposite roles in tumor formation, and the tumorigenesis of SKCM may be related to their abnormality.
Moreover, SKCM patients with high expression of five chemokines (CXCL9, CXCL10, CXCL13, CCL4, CCL5) were associated with better OS. While high levels of NMU and GAL showed poor prognosis. In addition, we found that the high expression of CXCL9, CXCL10, CXCL13, CCL4 and CCL5 were related to better outcomes in Breslow depth, Clark level, T-stages. IHC results from FAHSU cohort verified that CXCL9, CXCL10, CXCL13, CCL4, CCL5 were significantly overexpressed in SKCM tissues. All of the results suggested that chemokines members (CXCL9, CXCL10, CXCL13, CCL4, CCL5) and NMU, GAL were important biomarkers in SKCM tumorigenesis, progression and prognosis.
Chemokines are chemotactic cytokines mediating the migration and localization of immune cells [31]. C-X-C motif ligand 9,10,13 (CXCL9, CXCL10, CXCL13) 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 [30]. Research shows that CD8 + T cell infiltration is CXCR3 dependent [32]. CD8 + T cell plays important role in tumor microenvironment, which inhibits the proliferation and metastasis of tumor cells. A recent research found that the CXCL10-CXCR3 axis may relate to melanoma brain metastasis [33]. Additionally, the CXCR3 ligands enhanced the response rates to immune checkpoint blockade