Prognostic Impact of Autophagy and Immunity Related Biomarker for Uveal Melanoma

Background: Autophagy and immunity related genes serve crucial roles in carcinogenesis, but little is known about the prognostic impact for uveal melanoma (UM). Methods: Autophagy related and immunity related genes (AIRGs) expression data of 80 UM patients were obtained from the cancer genome atlas project (TCGA) database. Next, univariate cox regression analysis and the least absolute shrinkage and selection operator (LASSO) algorithms were applied to build a robust AIRGs signature in TCGA and validated in another two independent datasets. Besides, UM patients classied into two subgroups based on the risk model. Differences of clinical outcome, tumor microenvironment and the likelihood of chemotherapeutic response were further explored. Results: In total, a 4-AIRGs signature was constructed and validated in various datasets, which can robustly predict patients’ metastasis-free survival (MFS) and overall survival (OS) and is an independent prognostic factor in UM. The UM patients can be classied into high and low risk subgroups by applied risk score system. The high risk group have poor clinical outcomes, higher CD8+ T cell and macrophage immune-inltrating and more sensitive to chemotherapies. In addition, Gene Set Enrichment Analysis (GSEA) analysis revealed that hallmark epithelial-mesenchymal transition and KRAS pathways are commonly enriched in high-risk expression phenotype. Conclusion: Thus, our ndings provide a new clinical strategy for the accurate diagnosis and identify a novel prognostic autophagy and immunity associated biomarker for the treatment of uveal melanoma.

risk with a log-rank test p<0.001 ( Figure 3A). The ROC curve was drawn and 3,5 years of AUCs were 0.916 and 0.829 respectively ( Figure 3B). The distribution of risk score, OS, vital status, and expression of corresponding AIRGs in TCGA-UVM were shown in Figure 3C-E. In order to verify the robustness of the result, validated analyses were applied in another two independent datasets including GSE22138 and GSE84976. Patients in GSE22138 and GSE84976 datasets were accordingly classi ed into high-risk and low-risk groups. KM analyses manifested a similar result that patients in low risk group had signi cantly longer survival time than those in high-risk group (GSE22138 log-rank p=0.0019 and GSE84976 log-rank p<0.001 respectively) ( Figure 4A and Figure 5A). The 3,5 years of AUCs in GSE22138 were 0.746, 0.713 ( Figure 4B). Higher AUCs (3-year:0.836; 5-year:0.872) also observed in GSE84976 ( Figure   5B). The distribution of risk score, MFS, vital status, and expression of 4 AIRGs in GSE22138 was shown in Figure 4C-E. Meanwhile, the same risk score heatmap of 4 AIRGs in GSE84976 was shown in Figure 5C-E.

Prognostic value of biomarker
To explore prognostic factors for OS or MFS in multiple datasets, the risk score of AIRGs biomarker and clinical variables were conducted by the univariate and multivariate cox regression analyses ( Table 1). The results of univariate cox regression revealed that stage, age, metastasis, histological type and risk sore were signi cantly associated with MFS or OS. According to the multivariate cox regression analyses, only the risk score was stably remained independent prediction for OS or MFS in TCGA-VUM dataset (HR= 34.951, 95% CI = 4.891-249.759, P < 0.001), GSE22138 dataset (HR = 45.623, 95%CI = 10.025-564.298, P < 0.001) and GSE84976 dataset (HR = 3.532, 95%CI = 1.783-9.276, P =0.042). The 5 years AUCs of age, stage, metastasis, histological type and risk sore in TCGA-VUM set were 0.741, 0.778, 0.658, 0.424 and 0.829 respectively ( Figure 6A). In GSE22138, the 5 years AUCs of Chromosome 3 status and risk score were 0.762 and 0.713 ( Figure 6B). In addition, the 5 years AUCs of Chromosome 3 status, metastasis and risk sore in GSE84976 set were 0.912, 0.858, and 0.872 respectively ( Figure 6C).

Associations between risk model with clinical characteristics and immune in ltration
The relationships between the risk score signature and various clinical characteristics were investigated in TCGA-VUM dataset. The box plot manifested that UM patients in stage IV, metastasis, epithelioid histological type and dead status have higher risk score levels than other subtypes ( Figure 7A). Other clinical characteristics (age and gender) have no relationships with risk score. Thus, the risk model was intimately associated with tumor stage, metastasis, histological type and vital status. Through CIBERSORT analysis, the proportion of 22 different in ltrating immune cells was estimated in TCGA-VUM dataset. The stratifying analysis showed that the fractions of T cell CD8, T cells CD4, T cells gamma delta, B cells naïve, NK cells activated, Monocytes, Macrophages M1 and Neutrophils in high risk subgroup signi cantly different from those in low risk group ( Figure 7B).

Gene set enrichment analysis (GSEA)
In order to explore the different hallmark pathways enriched in high and low risk group, GSEA analyses were performed. According to selected criterion, we con rmed that six positive correlated pathways were enriched in high risk group, which including epithelial mesenchymal transition, estrogen response late, myogenesis, estrogen response early, apical junction and KRAS signaling. moreover, six cancer hallmark pathways contained interferon gamma response, protein secretion, interferon alpha response, MTORC1 signaling, JAK STAT3 signaling and in ammatory response were actively associated with low risk group ( Figure 7C).
High-risk subgroup more sensitive to chemotherapies At present, chemotherapy is a common treatment for UM, therefore GDSC database was also applied to predict sensitivity to chemo drugs. Assessment by selecting criteria, 45 kinds of drugs were screened out in TCGA-UVM set. 34 and 14 kinds of drugs were sensitive to high risk group in GSE22138 and GSE84976 respectively ( Table 2). The Venn plot of three datasets suggested that AMG.706 and JNK.Inhibitor.VIII were the common drugs ( Figure 8A). We could observe a signi cant difference in the estimated IC50 between high and low risk group for these chemo drugs where high-risk group could be more sensitive to chemotherapies ( Figure 8B-D).

Discussion
Like many kinds of tumor, the prognosis of uveal melanoma is largely depend on early diagnosis and treatment [3]. The growing clinical application of new therapeutic targets and the present studies are mainly focused on detection of novel prognostic biomarkers which are often used to assess disease risk and early diagnosis [19,20]. Clinical application of prognostic biomarkers can early guide high-risk of patients undergone individual treatment and management, even prevent life-threatening metastases [21][22][23][24]. Hence, as far as we known, our study is the rst time to investigate the prognostic role of AIRGs biomarker in UM. We also constructed a robust 4-AIRGs biomarker in UM prognosis and validated in multiple independent datasets. Our prognostic biomarker can subsequently classify patients into subgroups with different survival events. Compared with traditional clinical characteristics, our AIRGs biomarker achieved higher accuracy than most of clinical indicators (eg. stage, age, metastasis and histology type). Notably, the results of multivariate cox regression manifested that the risk score of 4-AIRGs biomarker can afford a robustly accurate prediction of OS or MFS in UM and could be considered as an independent prognostic model.
In this research, we developed a prognostic biomarker with four selected AIRGs (PRKCD, MPL, EREG and JAG2), which classi ed the OS or MFS of UM into high and low risk groups. Kaplan-Meier curves revealed that patients in high risk group have a poor survival. Among the four AIRGs, almost all of genes have been identi ed to be prognostic biomarker in other malignancies. For instance, PRKCD (Protein kinase C delta) is a member of protein kinase C (PKC) family which is a critical regulator of the chemosensitivity in cancers such as non-small cell lung cancer, ovarian cancer and prostate cancer [25][26][27]. Additionally, formation of immune synapses and enhancing anti-tumor function [29,30]. Christopher et al recently reported that C-MPL can work as a novel immunotherapeutic target to promote anti-tumor activity of T cells by mediation of cytokine pathways [31]. Therefore, MPL will be regarded as a potential therapeutic target for UM patients in future. As for EREG (Epiregulin), which belongs to the epidermal growth factor family and bene ts to wound healing, tissue repair, in ammation in normal physiology. While dysfunction of Epiregulin will increase the progression of different malignancies. Previous studies revealed that the overexpressed EREG closely correlates with many cancers, including prostate cancer, colon cancer, breast cancer and bladder cancer [32][33][34][35]. What's more, Asnaghi et al reported that JAG2 plays an important role in promoting tumor cells growth and metastasis in uveal melanoma, which could serve as a new therapeutic target to further investigation [36].
To better understand the underlying molecular functions, GSEA analyses were conducted to nd the potential 4-AIRGs related pathways in cancer hallmarks database. The results showed that hallmark epithelial-mesenchymal transition (EMT) pathway positively enriched in high risk expression subgroup. It's generally recognized that EMT pathway is a crucial factor controlling the clinical outcome of patients with various tumors. Asnaghi et al proven that EMT related factors will promote invasive properties of uveal melanoma cell and nally lead to poor prognosis [37]. Interestingly, our result consistent with this conclusion indicated that high risk subgroups positive enriched in epithelial-mesenchymal transition with a poor survival in UM. Additionally, we also identi ed that KRAS signaling was also positively associated high risk group. Kirsten rat sarcoma viral oncogene homolog (KRAS) regards as the most lethal cancer-related proteins which takes a crucial role in human aggressive cancers. Although the current study about relationship between KRAS and UM is few, KRAS has been recognized an induction factor of lung tumorigenesis for more than 20 years[38]. By CIBERSORT assessment, we also observed that high risk group have a higher in ltration of CD8 + T cells and M1 macrophage. Notably, recent studies have demonstrated that the increased level of CD8 + T cells is a predictive factor for poor prognosis in UM and the accumulation of macrophages are poorly associated with the prognosis of melanoma patients [39][40][41].
Therefore, it is logical to speculate that various malignant hallmarks and tumor in ltrating immune cells lead to the poor outcome of high risk patients.
Furthermore, in order to explore the different chemotherapeutic response between high and low risk subgroups, GDSC drug database was performed and the results showed that high risk subgroup is more sensitive to chemotherapies. Moreover, the three datasets showed that AMG.706 and JNK.Inhibitor.VIII were the commonly enriched drugs. The estimated IC50 of AMG.706 and JNK.Inhibitor.VIII with high risk group is signi cantly lower than the low risk groups. AMG.706 (Motesanib) is a multikinase inhibitor which signi cantly inhibits VEGF-induced vascular permeability, angiogenesis and induces regression in tumor xenografts [42]. As for JNK.Inhibitor.VIII, that is a JNK inhibitor which can inhibit the cell apoptosis and block the phosphorylation of c-Jun and JNKAR1 which promote tumorigenesis and metastasis of cancers [43]. Hence, it's reasonable to believe that AMG.706 and JNK.Inhibitor.VIII could be considered as potential chemo drugs for further clinical treatment.

Conclusion
Sum up, our research constructed a robust 4-AIRGs prognostic biomarker in UM, which could be regard as an independent prognostic model in clinical application. The patients in the high risk score group could bene t more chemotherapy. Furthermore, basic and clinical studies are also needed to validate the accuracy for predicting prognosis and managing patients.

Gene expression pro le and clinical information
Gene expression pro les of 80 UM and clinical information were obtained from TCGA-UVM database for training set. Another two datasets GSE84976 and GSE22138 consisting of 91 UM samples from Gene Expression Omnibus (GEO) were used as outside validation sets. The complete analysis work ow in this study is illustrated in Figure 1.

Autophagy & Immunity Related Genes (AIRGs)
The AIRGs were obtained from the Human Autophagy Database (http://www.autophagy.lu/) which contains 222 genes directly or indirectly associated with the autophagy process according to previous studies [18]. Moreover, the latest version of immunity related genes was acquired from the ImmPort Database (https://immport.niaid.nih.gov).

Identi cation of prognostic AIRGs model
The prognostic values of AIRGs were assessed by using overall survival (OS) time of patients in TCGA-UVM. Univariate cox-regression analysis was used to select AIRGs that were signi cantly correlated with uveal melanoma's OS (p values <0.05&HR <=0.90|HR>= 1.10). Next, we used LASSO algorithm to develop risk model via the prognostic AIRGs. Based on LASSO approach, the quali ed prognostic AIRGs were selected out to construct the risk system and generate risk score for each patient of UM according to the corresponding coe cients. Afterwards, the patients were accordingly classi ed into high-and lowrisk group by median cutoff. Kaplan-Meier survival curves (Log-rank) were applied to estimate the prognosis of high-and low-risk group. Receiver operating characteristic (ROC) analysis and area under the curve (AUC) were calculated to assess the sensitivity and speci city of risk model.

Relationships between risk score with clinical characteristics and immune microenvironment
To assess the relationship between risk score and clinical variables, the risk score distribution for subgroup of age, sex, stage, histological type, chromosome 3 status, metastasis and vital status was performed. Moreover, univariate and multivariate cox regression analysis for OS or metastasis-free survival (MFS) were carried out to evaluate the prognostic value of risk score and clinical parameters. The hazard ratios (HR) and 95% con dence intervals (95% CI) of the prognostic factors were calculated. Besides, to analyze the association between risk score and immune cell in ltration, "CIBERSORT" algorithm was performed to quantify the proportion of 22 different in ltrating immune cells in tumor immune microenvironment.

Gene set enrichment analysis
To illustrate the potential mechanism between the low-and high-risk groups, we performed Gene Set Enrichment Analysis (GSEA) via "clusterPro ler" R package. Firstly, the differentially expressed genes between low-and high-risk groups were ranked by the value of log2 fold change. Then the cancer hallmark database (h.all.v7.0.symbols) was applied in GSEA to explore the signaling pathways correlated with different subgroups of UM. The FDR q-value < 0.05 was used to screen signi cant pathways in each phenotype.

Chemotherapeutic response prediction
To explore the likelihood of chemotherapeutic drugs, the Genomics of Drug Sensitivity in Cancer (GDSC) database (https://www.cancerrxgene.org) was used to predict the chemotherapeutic response for each sample and the signi cant chemo drugs were selected (P<0.05). This analysis process was conducted by R package "pRRophetic".

Statistical analysis
Every statistical analysis was executed with R software(v.3.5.2) and corresponding statistical packages. The datasets used and analysed during the current study available from the corresponding author on reasonable request

Competing interests
All authors declare that they have no competing interests.

Funding
There is no sponsorship or funding arrangements relating to our research Author contributions HBD was responsible for writing-original draft preparation; YH was responsible for writing-review and editing; WL were responsible for data curation; QW was responsible for project administration and funding acquisition. All the authors commented and approved the text.   Figure 1 The complete work ow of the analysis in this study.