Identication of Biomarkers with Therapeutic And prognostic Value in Skin Cutaneous Melanoma

The morbidity and mortality of skin cutaneous melanoma (SKCM), the most deadly type of skin cancer, are on the rise worldwide. Through in-depth study of the tumor microenvironment (TME) of SKCM, this study further identied biomarkers with therapeutic and prognostic value. The gene expression proles and clinical data of SKCM patients were downloaded from The Cancer Genome Atlas database. Then we calculated the immune score and stromal score of patients with skin melanoma by using the estimate algorithm, and divided all patients into the high/low immune/stromal score groups and discussed the correlation between them and clinical characteristics.Then, limma R package was used to screen out the differential genes in the high/low immune/ stromal score groups, and the heat map of the differential genes was drawn. At the same time, the functional enrichment analysis of the differential genes was carried out. The protein ‒ protein interaction (PPI) network of different genes was constructed by using STRING and Cytoscape, and the key genes related to the prognosis of cutaneous melanoma were further selected. Finally, kinase target, co expression genes and immune inltrating cells of key genes were discussed. leukocytes, T-cell activation, T-cell aggregation, lymphocyte aggregation, leukocyte aggregation, adaptive immune response, regulation of leukocyte activation, regulation of lymphocyte activation, cytokine activity, cytokine receptor activity, chemokine activity, cytokine receptor binding, chemokine receptor binding, CCR chemokine receptor binding, antigen binding, carbohydrate binding, G-protein coupled chemoattractant receptor activity, and chemokine receptor activity. Our results conrm previous reports describing the critical role of immune cells and ECM molecules in the construction of the TME and the relationship between this microenvironment and SKCM progression[17, 18].We also studied the role of DEG in the prognosis of SKCM and identied 15 genes related to the TME. The expression of TLR7, AQP9, CCR8, PATL2, RGL4, OR9G1, CXCR5, IPCEF1, FCRL3, RCSD1, CP, NEXN, MAFB, SLC40A1, and TNFRSF13B is related to the survival of patients with SKCM.

Although skin cutaneous melanoma (SKCM) accounts for only 2% of all skin cancers [1], it is the most deadly type of skin cancer because of its high degree of malignancy and invasiveness. Indeed, the higher the SKCM stage, the higher the 5-year mortality. The I-IV survival rates were 97% (IA), 84% (IB), 68% (II), 55% (III), and 17% (IV) [2]. While continuous progress and medical technological innovations have prolonged the survival of patients and decreased mortality rates, the incidence of malignant SKCM continues to increase annually at an alarming rate [3]. Previous reports indicated that early diagnosis and accurate prognosis of patients with SKCM are very important for formulating optimal treatment plans and improving quality of life. However, no clinical biomarker that can accurately predict the prognosis of patients with SKCM is yet available. Thus, exploring the biomarker of SKCM is necessary to improve the treatment and prognosis of patients.
The tumor microenvironment (TME), which refers to the tumor cell environment composed of extracellular stromal molecules, in ammatory mediators, and various cells (e.g., immune cells, interstitial cells, and endothelial cells) [4,5], plays an important role in SKCM tumor metastasis, cell proliferation, invasion, and other biological behaviors [6]. Immune and stromal cells, as important non-tumor components of the TME, participate in tumor information transmission, immune monitoring, and niche formation. Previous studies demonstrated the diagnostic and prognostic value of immune and stromal cells in tumors [7,8]. This algorithm is based on gene expression features and can re ect the proportion of stromal and immune cells in tumor tissues to evaluate tumor purity [9]. An increasing number of researchers have used this algorithm to study different tumor types, such as glioblastoma [5], breast cancer [10], prostate cancer [11], and colon cancer [12]. However, applying ESTIMATE algorithm is necessary to explore the function of stromal and immune cells in SKCM.
In this study, the RNA expression pro les and clinical data (age, sex, tumor stage, survival status and prognosis) of patients with SKCM were extracted from The Cancer Genome Atlas (TCGA) database, which was established by the U.S. government to identify and catalog oncogene changes [13]. We then used the ESTIMATE algorithm [9] and tumor immune evaluation resource (TIMER) [14] to analyze the immune cells and genes in the SKCM TME. We not only studied the correlation between immune/stromal scores and clinical symptoms but also examined the relationship between TME-related genes and the prognosis of patients with SKCM. Finally, we identi ed the hub genes of SKCM and identi ed their kinase targets, co-expressed genes, and correlations with immune cell in ltration. Our results show that the composition of the TME affects the clinical prognosis of patients with SKCM. This nding may provide a basis for the future development of new prognostic biomarkers and treatment methods, especially immunotherapy.

Gene expression datasets
The gene expression pro les (level 3 data) of skin cutaneous melanoma patients were extracted from the TCGA database (https://tcga-data.nci.nih.gov/tcga/).The RNA expression patterns of SKCM patients were analyzed using the Illumina HiSeq. Then the clinical data of SKCM patients were extracted, including gender, age, grade, stage and survival status. After downloading and transposing the data, immune and stromal scores were calculated by performing the ESTIMATE algorithm .

Survival rates of immune/gene scores and their association with clinical characteristics
According to the ESTIMATE algorithm results, all samples were divided into high/low immune score groups and high/low stromal score groups to select crossover genes. Kaplan-Meier method was used to assess the prognostic value of the scores in skin cutaneous melanoma, and the overall survival curve was generated by R software. In addition, we studied the correlation between immune/stromal scores and clinical characteristics.

Identi cation of differentially expressed genes
Evidence of differential gene expression was analyzed using the "limma" package. Differentially expressed genes (DEGs) were screened by|logFC|>2 and adj.P<0.05.Clustering analysis and the construction of immune and stromal heatmaps were done using the pheatmap R package .

Enrichment analysis of DEGs
Functional annotation of DEGs was carried out by means of gene ontology (GO) analysis and Kyoto encyclopedia of genes and genomes (KEGG).Go analysis is primarily used to study cellular composition (CC), molecular function (MF), and biological processes (BP).GO and KEGG pathway analyses results were processed by the "clusterPro ler", "encichplot", and "ggplot2" packages. The P-value and q-value cut-offs were set at 0.05 each.

Construction of PPI network
Protein-protein interaction (PPI) networks can be generated using the STRING database [15]. The required minimum interaction score was set at 0.4. In addition, we used CytoHubba plug-in in Cytoscape to detect each module of PPI network and explore the key genes.

Overall survival curve of all DEGs
The prognostic value of DEGs in SKCM was evaluated using the Kaplan-Meier method, and overall survival curves were generated by the R software. The log-rank test was used to assess the prognostic value of the identi ed DEGs.

LinkedOmics
The Linkedomics (http://www.links.org/) database is a platform for accessing, analyzing, and comparing internal recombination data for different tumor types [16].In this study, the LinkInterpreter module of LinkedOmics was used to explore the kinase target of SKCM hub genes. At the same time, the coexpression of key genes was also studied.

TIMER
The TIMER (https://cistrome.shinyapps.io/timer/) database is a practical platform for the systematical analysis of immune in ltrates in diverse cancer types. The current study used the TIMER Gene module to explore the correlation between the abundances of six immune in-ltrates and the level of hub genes.

Correlation between immune/stromal scores and clinicopathologic variables of skin cutaneous melanoma
In this study, complete gene expression pro les and clinicopathological information from the TCGA database were used to analyze the data of all cases. According to the data of the ESTIMATE algorithm, stromal scores ranged from −1806.74 to 1892.23, immune scores ranged from −1481.02 to 3769.31, and ESTIMATE scores ranged from −2632.50 to 5077.47. The mean immune (P=0.009), stromal (P=0.012), and ESTIMATE (P=0.006) scores of patients with SKCM in the older-age group were consistently lower than those of patients in the lower-age group (Fig. 1A). Moreover, the mean immune (P=1.78e−04), stromal (P=3.69e−03), and ESTIMATE (P=3.95e−04) scores differed among stages I, II, III, and IV (Fig. 1B); here, the overall trend showing an initial decrease followed by a gradual increase. The mean immune (P=5.14e−04), stromal (P=0.017), and ESTIMATE (P=1.83e−03) scores also varied according to tumor stage (Fig. 1C). These results suggest that stromal, immune, and assessment scores are correlated with certain clinicopathologic variables.
Kaplan-Meier survival curves were drawn to assess the potential associations between overall survival and immune, stromal, and assessment scores. First, the 471 patients with SKCM were divided into high and low groups according to their scores. Immune score results showed that, compared with the highscore group, the overall survival time of patients with SKCM in the low-score group was longer ( Fig. 2A,  P=0.001). ESTIMATE score results showed that the overall survival time of patients with SKCM in the lowscore group was longer than that in the high-score group (Fig. 2, P=0.001). Finally, stromal scores (Fig.  2B, P=0.069) indicated that the overall survival time of patients in the low-score group was longer than that in the high-score group, but the difference was not statistically signi cant.
Next, we analyzed Affymetrix chip data from our 471 patients to reveal differences in gene expression pro les based on immune and stromal scores. A heat map drawn on the basis of immune scores is shown in Fig. 3A, which showed 1049 and 94 up-regulated genes in the high-and low-score groups, respectively (|fold change|>2, P<0.05). Fig. 3B shows the heatmaps drawn on the basis of stromal scores; here, 1369 and 30 genes were up-regulated in the high-and low-score groups, respectively (|fold change|>2, P<0.05). Differentially expressed genes (DEGs) in the immune and stromal score groups were identi ed by Venn map. Our results showed that 909 genes were up-regulated in the high-score group (Fig. 4A) while 5 genes were down-regulated in the low-score group (Fig. 4B).  (Fig. 5A). KEGG pathway analysis also indicated that the cytokine-cytokine receptor interaction and chemokine signaling pathways are involved in the tumor immune response and SKCM development.

Role of individual DEGs in overall survival in skin cutaneous melanoma
In this study, Kaplan-Meier survival curves were drawn to determine the potential function of 909 upregulated and 5 down-regulated genes in the overall survival of SKCM patients. The log-rank test showed that, among 1014 DEGs, AQP9, CCR8, PATL2, RGL4, OR9G1, CXCR5, IPCEF1, FCRL3, RCSD1, CP, NEXN, MAFB, SLC40A1, and TNFRSF13B have signi cant predictive effects on the overall survival of SKCM (  Fig. 7, all P<0.05).

Immune landscape of the tumor microenvironment in melanoma
We analyzed the correlation between hub genes and immune cell in ltration. Our results revealed a positive correlation between CCR8 expression and CD8+ T cells

Discussion
SKCM is one of the most deadly skin diseases worldwide. Although the overall survival of patients with SKCM has improved, mortality remains high mainly because the prognosis and therapeutic biomarkers of the disease have not been identi ed. In this study, SKCM data were extracted from the TCGA database to investigate the TME genes affecting the overall survival and tumorigenesis of SKCM. We then identi ed potential biomarkers for the prognosis and treatment of this lethal form of melanoma.
We used the ESTIMATE algorithm to obtain immune and stromal scores to understand the TME of SKCM. The relationship between the clinicopathological variables of SKCM and immune and stromal scores was then evaluated. Our results showed that immune and stromal scores were correlated with speci c clinicopathologic variables, namely, age, tumor stage, and T stage. In addition, the immune and stromal scores could predict the prognosis of patients with SKCM accurately. We divided the samples into highand low-immune/stromal score groups to understand the correlation between TME and survival rate. Our results suggest that patients with SKCM and low-immune, stromal, and ESTIMATE scores have longer overall survival times than those with high scores.
We identi ed 914 DEGs after comparing the high-and low-immune/stromal score groups. GO and KEGG pathway analysis revealed that many DEGs are involved in the tumor immune response, tumorigenesis, and TME maintenance. For example, GO analysis showed that DEGs are enriched in leukocyte cell-cell leukocytes, T-cell activation, T-cell aggregation, lymphocyte aggregation, leukocyte aggregation, adaptive immune response, regulation of leukocyte activation, regulation of lymphocyte activation, cytokine activity, cytokine receptor activity, chemokine activity, cytokine receptor binding, chemokine receptor binding, CCR chemokine receptor binding, antigen binding, carbohydrate binding, G-protein coupled chemoattractant receptor activity, and chemokine receptor activity. Our results con rm previous reports describing the critical role of immune cells and ECM molecules in the construction of the TME and the relationship between this microenvironment and SKCM progression [17,18].We also studied the role of DEG in the prognosis of SKCM and identi ed 15 genes related to the TME. The expression of TLR7, AQP9, CCR8, PATL2, RGL4, OR9G1, CXCR5, IPCEF1, FCRL3, RCSD1, CP, NEXN, MAFB, SLC40A1, and TNFRSF13B is related to the survival of patients with SKCM.
Next, we selected CCR8 and CXCR5 as key genes by combining signi cant DEGs and prognosisassociated genes. Loss of CXCR4, CXCR5, and CCR9 in circulating T-cell subsets has been reported to be associated with pulmonary involvement [21]. As a skin homing receptor, CCR8 is expressed in most T cells in normal skin and produces proin ammatory cytokines [22]. We then used the LinkedOmics database to explore kinase targets for the CCR8 and CXCR5 hubs. PRKCA, ATM, PLK1, CHEK1, and ATR are the ve most signi cant kinase targets of CXCR5. Kinases LYN, FYN, AURKB, PAK1, and LCK are the ve most prominent target networks for CCR8. We continued to study the coexpression genes of the CCR8 and CXCR5. Studies have shown that MAPK1 and other kinase targets can regulate the cell cycle, epithelialmesenchymal transformation, and tumor cell invasion and metastasis [23]. The results also suggest that CCR8 and CXCR5 may regulate the cell cycle, epithelial-mesenchymal transformation, and tumor cell invasion and metastasis of SKCM. Increasing evidence of a signi cant association between immune cell in ltration and patient prognosis in many tumor types has been reported [24].
Finally, we investigated the immune cell in ltration of CCR8 and CXCR5. The results showed that CCR8 and CXCR5 are positively correlated with the in ltration of neutrophils, macrophages, and CD8+ T, CD4+ T, B, and dendritic cells. Studies have shown that enhanced CD8+T-cell tra cking is very important in the anti-PD-1/anti-CTLA-4 e cacy of melanoma brain metastases, which may be an immunotherapeutic enhancement strategy [25]. The interaction between melanoma and immune cells has also been studied, and the key role of these cells in tumor growth and the immune response has been emphasized [26]. Thus, these genes may provide additional information on the correlation between immune cell in ltration and clinical outcomes in SKCM patients.
The TME is involved in the development of cancer and interacts with SKCM, thus affecting the progression, metastasis, drug resistance, and prognosis of the disease. This study evaluated the relationship between the immune/stromal scores and clinical characteristics of patients with SKCM. The function and prognostic value of TME-related genes were also investigated. Our ndings add to the existing body of knowledge on the complex interactions between SKCM and its tumor environment. However, the results of our study are very limited because it was carried out in a small population. Therefore, a more comprehensive analysis of a larger corpus of patient data is necessary to reveal the complex relationship between SKCM and its TME.

Conclusions
In summary, we evaluated the immune score and stromal score of SKCM by using the ESTIMATE algorithm, and studied the genes related to the tumor microenvironment by functional enrichment analysis.Our study provides new biomarkers for the treatment and prognosis of SKCM. 6. Declarations

Ethics approval and consent to participate
The data in this study are from public databases on the Internet and do not involve ethical issues.

Consent for publication
All authors read and approved the nal manuscript.

Availability of data and materials
The data of this study includes two aspects.Protein data can be downloaded from the TCPA database (https://www.tcpaportal.org), while clinicopathological data can be downloaded from the TCGA database (https://www.tcpaportal.org/tcpa/).

Competing interests
The authors declare that there are no con ict of interests.

Authors' contributions
Guangwen Wang and Yonghuo Ling have contributed equally to this work. 6.7 Acknowledgements  Figure 1 Association between immune/stromal scores and clinicopathologic variables in SKCM. (A) A signi cant association was revealed between age and the level of immune scores (P = 0.009) stromal scores (P = 0.012), and ESTIMATE score (P = 0.006). (B) A signi cant association was revealed between T stage and the level of immune scores (P = 5.146e-04) stromal scores (P = 0.017), and ESTIMATE score (P = 0.002).

Figure 2
Association between immune/stromal scores and overall survival in SKCM. (A) Overall survival in SKCM patients with high immune scores is poorer than of patients with low immune scores (P = 0.001). (B) Overall survival in SKCM patients with high stromal scores is poorer than that of patients with low stromal immune scores (P = 0.069). (C) Overall survival in SKCM patients with high ESTIMATE score is poorer than that of patients with low ESTIMATE scores (P = 0.001). The logrank test was used to calculate the P-value.   Commonly DEGs PPI networks constructed by STRING tool and Cytoscape. (A) There were a total of 162 DEGs in the PPI network complex. The nodes meantgenes; the edges meant the interaction of genes; (B) The signi cant genes in PPI networks via CytoHubba plug-in (degree cutoff = 2, node score cutoff = 0.2, kcore = 2, and max. Depth = 100).