Screening of Autophagy-Related Prognostic Genes in Metastatic Skin Melanoma

Cutaneous melanoma refers to a common skin tumor that is dangerous to health with a great risk of metastasis. Previous researches reported that autophagy is associated with the progression of cutaneous melanoma. Nevertheless, the role played by genes with a relation to autophagy (ARG) in the prediction of the course of metastatic cutaneous melanoma is still largely unknown. We observed that thirteen ARGs showed relations to overall survival (OS) in the Cox regression investigation based on a single variate. We developed 2-gene signature, which stratified metastatic cutaneous melanoma cases to groups at great and small risks. Cases suffering from metastatic cutaneous melanoma in the group at great risks had power OS compared with cases at small risks. The risk score, T phase, N phase, and age were proved to be individual factors in terms of the prediction of OS. Besides, the risk scores identified by the two ARGs were significantly correlated with metastatic cutaneous melanoma. Receiver operating characteristic (ROC) curve analysis demonstrated accurate predicting performance exhibited by the 2-gene signature. We also found that the immunization and stromal scores achieved by the group based on large risks were higher compared with those achieved by the group based on small risks. The metastatic cutaneous melanoma cases achieving the score based on small risks acquired greater expression of immune checkpoint molecules as compared with the high-risk group. In conclusion, the 2-ARG gene signature indicated a novel prognostic indicator for prognosis prediction of metastatic cutaneous melanoma, which served as an important tool for guiding the clinical treatment of cutaneous melanoma.


Introduction
Cutaneous melanoma refers to one type of skin malignant tumor exhibiting high malignancy and ineffective prediction of disease courses [1,2]. Occult onset and easy invasion and metastasis are important clinical features of cutaneous melanoma [3]. According to the statistics, cutaneous melanoma's incidence in China rose to about 3-5% [4]. Although cutaneous melanoma's incidence within China and other Asian countries is relatively low compared with those in Europe and America, cutaneous melanoma's incidence within China is increasing rapidly [5]. Once cutaneous melanoma cases have distant metastasis, they are diagnosed as advanced or metastatic cm, so the survival time of cutaneous melanoma cases is often short [6]. In the current treatment strategies of metastatic cutaneous melanoma, targeted therapy and immunotherapy play an important role. Metastatic tumor surgery and radiotherapy can also be used selectively [7].
Due to the high metastasis rate of cutaneous melanoma, it is necessary to find a new prognosis model to provide theoretical guidance for the treatment of metastatic cutaneous melanoma.
Autophagy is a process in which cytoplasmic components, or organelles are encapsulated and transported to lysosomes for degradation by forming double membrane autophagosomes [8,9]. Autophagy can be induced by DNA damage, chemical drugs, ion irradiation, reactive oxygen species, and abnormal growth of tumor cells [10]. According to existing works, autophagy refers to a barrier against malignant transformation of carcinoma cells [11,12]. Some major oncogenes, such as mTOR and Akt, are considered to be negative regulators of autophagy [13,14]. According to considerable works, mutant tumor suppressors such as PTEN and TSC1/2 can activate autophagy [15]. It is controversial whether autophagy has a tumor promoting or antitumor effect on the occurrence and development of cancer. At present, according to some studies, autophagy impacts tumor inhibition in the early stage of cancer, but it plays a role in promoting cancer in the formed tumor and contributes to the generation of drug resistance of cancer cells [16,17]. Therefore, autophagy is considered to help promote the survival of cancer cells in the advanced phase. At present, there are few comprehensive studies on exploring autophagyrelevant genes within the disease course prediction and immunotherapy of metastatic cutaneous melanoma.
Here, we obtained RNA-seq and clinic information regarding cutaneous melanoma cases from The Cancer Genome Atlas (TCGA) database. By several bioinformatics investigations, a multigene signature of ARG was constructed. Relationships of risk model and clinicopathological features of metastatic cutaneous melanoma were confirmed. Then, we carried out Cox regression investigations based on single and multiple variates for the identification of individual factors for the OS of metastatic cutaneous melanoma. A nomogram containing independent prognostic factors was built using "rms" package. We carried out GSVA for exploring the biological processes and pathways involved in the groups based on great and small risks. Furthermore, the analysis was conducted on the landscape of immune infiltration and the expression of immune checkpoint molecules in metastatic cutaneous melanoma. This work might provide a new idea for prognosis and immunotherapy of later phase metastasis melanoma.

Material and Method
2.1. Data Acquisition and Preprocessing. We retrieved the RNA-seq and clinic data regarding the cases with cutaneous melanoma according to The Cancer Genome Atlas (TCGA) cohort (https://portal.gdc.cancer.gov/), which contained 103 primary and 367 metastatic cancer cases. Totally, 232 ARGs were acquired in the Human Autophagy Database (http://autophagy.lu/).

Differentially Expressed Analysis.
For screening the gene that achieves different expressions (DEGs) in the metastatic and primary carcinoma samples, the "limma" package was used with regulated P value < 0.05 as well as |log 2 ðfold changeÞ | >0:5 [18]. The expression of Top100 genes that achieved different expressions between the metastatic and primary carcinoma samples was shown in a heat map. Then, the DEGs were overlapped with the ARGs to obtain the ARGs that achieved different expressions (DE-ARGs), which were chosen to conduct the subsequent investigation.

Development and
Verification of the Prognostic Signature in relation to Autophagy. Subsequently, metastatic cutaneous melanoma cases indiscriminately fell to the test set (n = 93) and the training set (n = 217) at 3 : 7. To explore whether each DE-ARG is related to overall survival (OS), we performed Cox regression investigation based on a single variate in the training set. The DE-ARGs with the P value < 0.05 were identified, followed by the subsequent analysis based on Cox regression investigation based on multiple variates to obtain the best risk model. In the Cox regression investigation based on multiple variates, this study applied the stepwise regression function and set the "direction" as "both." Based on the risk model, the score of risk of the respective metastatic cutaneous melanoma case was obtained by: risk score = ðβ 1 G 1 + β2 G 2 + β3 G 3+⋯+β n G nÞ. In the calculated formula, β stands the coefficient of gene, and G stands the expression level of each gene. The metastatic cutaneous melanoma cases were then stratified into the group based on small risks and group based on large risk in accordance with the mean value of risk score. Furthermore, the OS of these groups was compared using the Kaplan-Meier (K-M) approach on the basis of the log-rank test. In addition, using "survivalROC" R package, we obtained the curves of receiver operating characteristic (ROC) of 1, 3, and 5 years [19]. To be specific, we obtained the area under the curve (AUC) for assessing the risk model's effectiveness. The model of risk was further verified based on the test set.

Functional Enrichment
Analysis. According to the gene sets files, we carried out GSVA (PMID: 23323831) for the exploration of the potential biology process and pathways most relevant to groups based on large risk and small risk of metastatic cutaneous melanoma cases. Using the "gsva" package of R software, we carried out a single sample gene set enrichment investigation (ssGSEA) for calculating infiltrating immune cells' proportion in cases with metastatic cutaneous melanoma [20].
2.5. Evaluation of Immune Microenvironment. We implemented "ESTIMATE" R package for obtaining the immunization and stromal scores of metastatic cutaneous melanoma cases within TCGA database [21]. Furthermore, the expressions of immunization checkpoint molecules were examined in the metastatic cutaneous melanoma samples.
2.6. Statistical Analysis. We carried out the statistical investigations with R software (Version 3.5.3). We investigated various groups' OS on the basis of K-M investigation and compared OS by the log-rank test. Cox regression investigations based on single and multiple variates were applied to investigate the individual prognostic factors for OS. The nomogram containing clinicopathological features was constructed by "rms" package. The differences between two groups were compared using Wilcox.test. Differences were considered statistically significant when P < 0:05.

Identification of Autophagy-Related Genes with Different Expressions (DE-ARGs) in Metastatic Cutaneous Melanoma.
To seek the ARGs in relations to the disease course prediction of cutaneous melanoma, we first analyzed the DEGs between the primary and metastatic cancer samples of TCGA database using "limma" package. Under the threshold of regulated P value < 0.05 as well as |log 2 ðfold changeÞ | >0:5, we identified 886 DEGs in total, covering 554 significantly upregulated and 332 significantly decreased genes within metastatic cancer samples in comparison with the primary cancer samples (Figure 1(a)). Figure 1

Establishment and
Validation of the Prognostic Signature in relation to Autophagy. For more specifically assessing whether the DE-ARGs are related to the survival of metastatic cutaneous melanoma cases, the Cox regression investigation based on single variate was performed within the training set (Figure 2(a)). The result indicated that two genes had significant relations to the metastatic cutaneous

Disease Markers
Taken together, all the results suggested the reliable predicting performance exhibited by the prognosis signature constructed by the two ARGs.

The Associations between the Risk Score and the Clinicopathological Features in Cases with Metastatic
Cutaneous Melanoma. The expression of the two screened ARGs of the small-and great-risk samples within the TCGA dataset is illustrated by heatmaps (Figure 4(a)). We observed differences with statistical significance in these groups within the training and test sets. To further investigate the associations of the risk scores and clinicopathology characteristics, this study quantitatively analyzed the risk score in metastatic cutaneous melanoma (Figures 4(b)-4(g)). As a result, the risk scores and the survival probability were both significantly different in these groups with the TCGA dataset compartmentalized by Phase and T phase. However, the risk scores were no significant differences in these groups divided by age, gender, M, and N phase. Moreover, the survival probability was significantly different in these groups classified by age and N phase.
By the same taken, we performed the stratified survival investigation on the clinicopathological features. According to Figure 5(a), greater risk scores showed relations to lower survival according to male cases, whereas female cases showed an insignificant difference. Furthermore, high-risk score noticeably caused a poorer OS in metastatic cutaneous melanoma cases with Phase I-II, Phase III-IV ( Figure 5(b)), M0, N0, and N1, whereas it was not found to be risk factors in terms of metastatic cutaneous melanoma cases with phase M1 (Figures 5(c) and 5(d)). The results demonstrated that the risk scores identified by two ARGs were significantly correlated with metastatic cutaneous melanoma.

Development and Evaluation of the Nomogram for OS in
Metastatic Cutaneous Melanoma. Then, we established the nomogram using the five clinicopathological features including risk score, age, pathologic phase, T phase, and N phase (Figure 7(a)). The accurate prediction efficiency of 1-year survival and 3-year survival in the TCGA database was investigated by the calibration curve (Figure 7(b)). Moreover, according to the analysis in terms of decision curve (DCA), the risk model with the addition of clinicopathological features showed better net benefit than the risk only model (Figure 7(c)), which suggested the ability of the nomogram in the accurate prediction of the prognosis of metastatic cutaneous melanoma cases.
3.6. Functional Analyses in the TCGA Database. Furthermore, GSVA was performed for elucidating the biology process and channels related to the risk score. According to Figure 8 (Figure 8(b)), further suggesting that these prognostic genes might participate in the progression of cutaneous melanoma metastasis.

The Landscape of Immune Infiltration within Metastatic
Cutaneous Melanoma. With the use of the ESTIMATE, the content of stromal and immunization cells within metastatic cutaneous melanoma tumor tissues was calculated. We found that the immunization and stromal scores achieved by the group based on great risks exceeded those achieved by the group based on low risks (Figures 9(a) and 9(b), P < 0:05). Moreover, we analyzed the ESTIMATE of the two groups and obtained the same trends (Figure 9(c), P < 0:05 ). To investigate relations of the score of risk and immunization state, this study determined the enrichment score of immunization gene sets. Interestingly, the score of Th1 cells, TFH, Tem, Tcm, T helper cells, T cells, pDC NK CD56dim cells, neutrophils, mast cells, macrophages, iDC, iDC, eosinophils, DC, cytotoxic cells, CD8 T cells, B cells, aDC, Th17 cells, Th2 cells, and TReg showed noticeable distinctions in the groups based on small and great risks in the TCGA group (all P < 0:05, Figure 9(d)). Accordingly, the immune   3.8. The Expressions of Immune Checkpoint Molecules. Programmed cell death receptor ligand 1 (PD-L1) and blocking programmed cell death 1 receptor (PD-1) have been a spe-cial interest in developing antibodies for a subset of cancer cases (PMID: 31488176). Therefore, immune checkpoint proteins have diverse clinical implications in the immunotherapy of cancers. We then investigated any potential relation of the score of risk and the expressions achieved by immunization checkpoint molecules. According to Figure 10, the metastatic cutaneous melanoma cases  Pathologic_n   10  20  30  40  50  60  70  80  90  100   10 20 30 40 50 60 70 80 90   T0  T3   T1  T2  T4   N0  12 Disease Markers achieving low-risk score had greater expression of immune checkpoint molecules than the high-risk group. Accordingly, the low-risk cases suffering from metastatic cutaneous melanoma might have a more promising treatment to respond for immunotherapies. 13 Disease Markers or LY294002 and TMZ could enhance the cytotoxicity of alkylating agents on human melanoma cell lines [22]. Zhang et al. investigated that CX-F9, a novel Ribosomal S6 Kinase 2 (RSK2) inhibitor, could significantly suppress the proliferation, invasion, and autophagy of melanoma in vitro and in vivo [23]. In the present study, thirteen ARGs showed correlations to OS in the Cox regression investigation based on a single variate. A 2-gene signature was developed, which stratified metastatic cutaneous melanoma cases into the groups based on great and small risks. Cases with metastatic cutaneous melanoma in the high-risk group had worse OS than that of the group based on small risks. The risk score, T phase, N phase, and age were proved to be individual factors for predicting OS. Besides, the risk scores identified by the two ARGs were significantly correlated with metastatic cutaneous melanoma. Receiver operating characteristic  14 Disease Markers

Discussion
(ROC) curve analysis demonstrated accurate predictive performance of the 2-gene signature. Functional enrichment analysis indicated that immune-related biological processes and channels were significantly enriched. The infiltrating immune cell content was different between the two risk groups. We also found that the immune scores and stromal scores of the high-risk group were higher compared with that of group based on low risks. The metastatic cutaneous melanoma cases achieving low-risk scores had greater expression of immune checkpoint molecules as compared with the high-risk group. Although autophagy-related genes play a crucial role in some diseases, there is no comprehensive study on the prognosis and immunotherapy of autophagy-related genes in cases suffering from metastatic cutaneous melanoma. This study reports for the first time the role of autophagyrelated genes in the prognosis and immunotherapy of cases with metastatic cutaneous melanoma. The risk models of HspB8 and CCR2 were established by univariate Cox and multivariate Cox analysis. HspB8 acts as an oncogene in several cancers. Shen et al. reported that HSPB8 promoted cancer cell growth by activating the ERK-CREB pathway and predicted a poor prognosis in gastric cancer cases [24]. The expression of HSPB8 is investigated to correlate with breast cancer progression [25]. In addition, HSPB8 is responsible for the rug resistance of breast cancer cells. The mTOR inhibitor (AZD8055) could inhibit the tamoxifen resistance in breast cancer cells by suppressing the expression of HSPB8 [26]. Several studies have shown that CCR2 acts as a novel biomarker in metastatic cutaneous melanoma [27,28]. Furthermore, Zhang et al. demonstrated that Toll-like receptors 7 and 8 expression correlated with the expression of immune biomarkers (CCR2, CCR5, CCL3, and CCL5) and positively predicted the clinical outcome of cases with melanoma [29]. This work is the first time to stratify cases with metastatic cutaneous melanoma based on autophagy-related genes, which provides new insights for predicting the efficacy of immunotherapy and possible differentiation targets.
We further performed functional analyses of ARGs in the TCGA database. The results of GSVA elucidated that autophagy-related genes may be closely related to tumor immunity. For example, the ARGs were enriched in GO_TOLL_ LIKE_RECEPTOR_7_SIGNALING_PATHWAY, GO_NAT-URAL_KILLER_CELL_CHEMOTAXIS, and GO_POSI-TIVE_REGULATION_OF_TYPE_2_IMMUNE_RESPONSE. These pathways are associated with the tumor immunity of melanoma [30][31][32]. The KEGG pathway analyses also indicated the KEGG_INTESTINAL_IMMUNE_NETWORK_FOR_ IGA_PRODUCTION, KEGG_PRIMARY_IMMUNODEFI-CIENCY, KEGG_ANTIGEN_PROCESSING_AND_PRESEN-TATION, KEGG_AUTOIMMUNE_THYROID_DISEASE, KEGG_GRAFT_VERSUS_HOST_DISEASE, KEGG_ALLO-GRAFT_REJECTION, and KEGG_TYPE_I_DIABETES_ MELLITUS were enriched in high-risk groups. These pathways participated in the process of immune escape in cutaneous melanoma metastasis [33,34]. The results revealed that these prognostic genes might participate in the progression of cutaneous melanoma metastasis.
Since the discovery of immune checkpoint proteins, the development of antibodies against programmed cell death receptor-1 (PD-1) and programmed cell death receptor ligand-1 (PD-L1) has aroused special interest in the treatment of some cancer cases [35,36]. PD-1 signal carries out the negative regulation of T cell-mediated immune response, which is one of the mechanisms of tumor escaping antigenspecific T cell immune response [37]. It facilitates tumor development and progression by improving the survival rate of tumor cells [38]. In this context, PD-1 signaling is a valuable new and effective target for cancer immunotherapy. Javed et al. reported that significant differences existed in PD-L1 expression between metastatic uveal melanoma and  TNFSF18  CD28  CD47  LAG3  ICOS  TNFRSF14  TNFSF4  CD40  TIGIT  CD70  CD86  CD160  CD274  CD80  PDCD1  TNFSF14  CEACAM1  CD48  PVR  CD226  TNFRSF18 LGALS9 Genes expression log2 (FPKM+1)
15 Disease Markers metastatic cutaneous melanoma. The higher PD-L1 expression was observed in metastatic cutaneous melanoma [39]. This work reported that the expressions achieved by the mentioned key immune checkpoints increased in the group based on small risks. The key immune checkpoints including BTLA, CD86, CD244, and PDCD1 are recognized as predictors of sentinel lymph node metastasis in cutaneous melanoma [40,41]. Our results indicated that the low-risk cases with metastatic cutaneous melanoma might have a more promising treatment to respond for immunotherapies.
In conclusion, the 2-ARG gene signature indicates a novel prognostic indicator for prognosis prediction of metastatic cutaneous melanoma, which served as an important tool for guiding the clinical treatment of cutaneous melanoma.

Data Availability
The datasets used and analyzed during the current study are available from the corresponding author on reasonable request.

Ethical Approval
Ethical approval is not necessary.

Conflicts of Interest
The authors declare that there is no conflict of interest regarding the publication of this paper.

Authors' Contributions
Cao-Jie Chen, Shigeki Sakai, and Kazuo Kishi conceived and designed the study. Cao-Jie Chen, Hiroki Kajita, and Noriko Aramaki-Hattori did bioinformatics analysis and wrote the manuscript. Shigeki Sakai and Noriko Aramaki-Hattori had supervised this project and contributed to writing and revision of the manuscript. All authors contributed to the article and approved the submitted version.