Abstract

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.

1. 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 autophagy-relevant 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.

2. 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/).

2.2. 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 value < 0.05 as well as [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.

2.3. Development and Verification of the Prognostic Signature in relation to Autophagy

Subsequently, metastatic cutaneous melanoma cases indiscriminately fell to the test set () and the training set () 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 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: . In the calculated formula, stands the coefficient of gene, and 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.

2.4. 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 .

3. Result

3.1. 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 value < 0.05 as well as , 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(b) reveals the expression of Top100 DEGs between the primary and metastatic cancer samples of the TCGA database. Furthermore, we combined the 886 DEGs with 222 ARGs, obtaining 13 DE-ARGs, to carry out the following investigation (Figure 1(c)).

3.2. 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 melanoma cases’ OS (), of which HSPB8 was a risk factor (), and CCR2 was a protective factor () in metastatic cutaneous melanoma. Then, we reported a 2-gene signature based on the Cox regression investigation based on multiple variates (Figure 2(b)). Two ARGs and their corresponding coefficients were utilized for determining the score of risk of the respective metastatic cutaneous melanoma case. The calculated equation in terms of the score of risk is presented as within the training and test sets. Figure 2(c) illustrates the distributions of case risk scores and survivals. The K-M survival investigation revealed that the metastatic cutaneous melanoma cases with high-risk scores had a significantly poorer OS in comparison with that of cases based on small risks (Figure 2(c), ). Besides, according to the ROC depending on time, the AUC of 1-, 3-, and 5-year OS reached 0.733, 0.658, and 0.629, separately (Figure 2(e)). Lastly, we verified the 2-ARGs prognosis signature with the use of OS information according to the test set, complying with the results of the training set (Figures 3(a)3(c)), in which AUC of 1-, 3-, and 5-year OS reached 0.711, 0.627, and 0.683, separately. Taken together, all the results suggested the reliable predicting performance exhibited by the prognosis signature constructed by the two ARGs.

3.3. 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.

3.4. Individual Prognosis Value Achieved by the 2-Gene Signature within Metastatic Cutaneous Melanoma

To further demonstrate whether this risk score acted as an individual factor in terms of the prediction of the course of metastatic cutaneous melanoma according to the clinicopathology characteristics of age, gender, T phase, N phase, M phase, and Pathological Phase. The result of univariate Cox regression investigations indicated that the risk score, age, pathologic phase, T phase, and N phase showed noticeable relations to the prediction of the course of metastatic cutaneous melanoma (Figure 6(b), ). Based on these significant clinicopathological features, we further performed Cox regression investigation based on multiple variates. As revealed from the result, the risk score, T phase, N phase, and age are significantly correlated with the OS (Figure 6(b), ).

3.5. 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(a), many immune-related GO terms, including GO_TOLL_LIKE_RECEPTOR_7_SIGNALING_PATHWAY, GO_REGULATION_OF_NATURAL_KILLER_CELL_MEDIATED_IMMUNITY, GO_NATURAL_KILLER_CELL_CHEMOTAXIS, GO_T_CELL_ACTIVATION_VIA_T_CELL_RECEPTOR_CONTACT_WITH_ANTIGEN_BOUND_TO_MHC_MOLECULE_ON_ANTIGEN_PRESENTING_CELL, and GO_POSITIVE_REGULATION_OF_TYPE_2_IMMUNE_RESPONSE, showed enrichment within the groups with the score based on small risks. Besides, the KEGG pathway analyses also indicated the KEGG_INTESTINAL_IMMUNE_NETWORK_FOR_IGA_PRODUCTION, KEGG_PRIMARY_IMMUNODEFICIENCY, KEGG_ANTIGEN_PROCESSING_AND_PRESENTATION, KEGG_AUTOIMMUNE_THYROID_DISEASE, KEGG_GRAFT_VERSUS_HOST_DISEASE, KEGG_ALLOGRAFT_REJECTION, and KEGG_TYPE_I_DIABETES_MELLITUS were enriched in high-risk groups (Figure 8(b)), further suggesting that these prognostic genes might participate in the progression of cutaneous melanoma metastasis.

3.7. 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), ). Moreover, we analyzed the ESTIMATE of the two groups and obtained the same trends (Figure 9(c), ). 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 , Figure 9(d)). Accordingly, the immune infiltration in metastatic cutaneous melanoma may act as targets for immunotherapy and may have potential clinical implications.

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 special 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 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.

4. Discussion

Autophagy and autophagy-related genes play an important role in metastatic cutaneous melanoma. Ryabaya et al. reported that autophagy inhibitor integration chloroquine 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 (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 autophagy-related 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_NATURAL_KILLER_CELL_CHEMOTAXIS, and GO_POSITIVE_REGULATION_OF_TYPE_2_IMMUNE_RESPONSE. These pathways are associated with the tumor immunity of melanoma [3032]. The KEGG pathway analyses also indicated the KEGG_INTESTINAL_IMMUNE_NETWORK_FOR_IGA_PRODUCTION, KEGG_PRIMARY_IMMUNODEFICIENCY, KEGG_ANTIGEN_PROCESSING_AND_PRESENTATION, KEGG_AUTOIMMUNE_THYROID_DISEASE, KEGG_GRAFT_VERSUS_HOST_DISEASE, KEGG_ALLOGRAFT_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 antigen-specific 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 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.