A Novel Four Genes of Prognostic Signature for Uveal Melanoma

Autophagy and immunity play critical roles in various cancers, but the prognostic impact of autophagy and immunity for uveal melanoma (UM) remains lacking. Therefore, the RNA sequencing of data in the TCGA-UVM dataset was downloaded from UCSC Xena database. The prognostic autophagy- and immunity-related genes (AIRGs) were selected via univariate Cox regression. Next, we applied LASSO method to construct four genes of signature in the TCGA-UVM and verified in another two GEO datasets (GSE84976 and GSE22138). This signature intimately associated with overall survival (OS) time and metastasis-free survival (MFS) time of UM, which could be considered as a prognostic indicator. Besides, by applying risk assessment, the patients of UM can be divided into two subgroups (high/low risk) with different survival time, distinct clinical outcomes, and immune microenvironments. Gene set enrichment analysis (GSEA) manifested that cancer hallmark epithelial-mesenchymal transition and KRAS pathways were positively activated in the high-risk group. Moreover, the high-risk group could be more sensitive to chemotherapies than the low-risk group. Thus, our finding suggested that the four genes of signature closely linked with UM risk and survival can afford more accurate survival prediction and potential therapeutic targets for clinical application.


Introduction
Uveal melanoma (UM) is the most common primary intraocular malignancy which accounts for 70% of all ocular cancers. It has become a growing global public health concern with about 50% of UM patients dying of metastatic disease [1][2][3]. Despite there are certain advances in treatment of UM, the prognosis of UM patients is still poor [4,5]. Thus, it is urgently required to discovery the novel prognostic biomarkers of UM to unveil the underlying potential molecular mechanism as well as therapeutic targets.
Autophagy is an important process mediating intracellular degradation which regulates cellular and biological homeostasis. Deregulation of autophagy has been proven in various pathological and disease processes, including differentiation, development, and tumorigenesis [6]. As for malignant tumors, the function of autophagy in tumorigen-esis seems to be a two-edged sword, which may change at different periods [7]. Recently, many researches revealed that autophagy can reduce cell damage and maintain chromosomal stability by eliminating damaged protein and organelles in early phases of cancer progression [8]. In contrast, autophagy seems to promote tumor cell proliferation and escape from immune surveillance at the late stage of cancer [9]. Moreover, autophagy is intimately associated with immune response, inflammation, and therapeutic resistance [10,11]. For instance, recent researches have showed that therapeutic targets on autophagy will enhance the immune responses and antitumor effects and overcome drug resistances [12,13]. Thus, it is undoubted that autophagy and immunity play crucial functions in diverse physiological and pathophysiological processes, especially in cancers. In addition, growing researches have demonstrated that autophagy-and immunity-related gene (AIRG) biomarkers are closely correlated with the prognosis of various types of cancer, including glioma, breast cancer, pancreatic ductal adenocarcinoma, and hepatocellular carcinoma [14][15][16]. In spite of the significant role of AIRGs in various cancers have been proven, the prognostic value for uveal melanoma is poorly understood.
Luckily, large-scale transcriptome datasets are deposited in accessible repositories, such as UCSC Xena website and GEO database, which afforded valuable resources to explore potential signatures in various tumors for different organisms and biological conditions [17].
Therefore, in our research, AIRGs were firstly downloaded from the ImmPort Database website and the Human Autophagy Database. Next, the univariate Cox analysis and LASSO model were performed to construct a robust four genes of signature, which can distinguish different clinical outcomes, immune microenvironment, and chemotherapy response of UM. The present study indicated that the expression of AIRGs takes pivotal roles in the prognosis of UM and could be considered as a prognostic marker for UM therapy.

Gene Expression Profile and Clinical Information.
The RNA sequencing data of 80 UM samples as well as clinical characteristics in the TCGA-UVM dataset were downloaded from UCSC Xena website, which was used for training set. Another two gene expression datasets (GSE84976 and GSE22138) contained 91 UM patients, which was acquired from GEO database and used as outside validation sets. The workflow of overall bioinformatic analysis in this study is shown in Figure 1.

2.2.
Autophagy-and Immunity-Related Genes (AIRGs). The autophagy-related genes were collected from the Human Autophagy Database (http://www.autophagy.lu/), which has been demonstrated to take part in the process of autophagy according to the previous studies. Moreover, the list of immunity-related genes was derived from the ImmPort Database (https://immport.niaid.nih.gov), which was the largest accessible human immunology database.
2.3. Identification of Prognostic AIRG Model. The overall survival (OS) time was applied to assess the prognostic value of AIRGs in the TCGA-UVM. The univariate cox regression and Benjamini and Hochberg (BH) correcting methods were used to select prognostic AIRGs (adjust p values < 0.05 and HR > = 1:10 | HR< = 0:90Þ: Next, we used LASSO algorithm to develop risk model via the prognostic AIRGs. The qualified prognostic AIRGs were selected out to construct the risk system and generate risk scores for UM patients based on the corresponding coefficients. Afterwards, the patients of UM were classified into two subgroups (high or low risk) by the median cutoff value of risk score. The Kaplan-Meier survival curve (log-rank test) was applied to estimate the prognosis of the high-risk and low-risk groups. The sensitivity and specificity of risk model were assessed by receiver operating characteristic (ROC) analysis with area under the curve (AUC).

Associations between Risk Score with Clinical Features
and Immune Microenvironment. To evaluate the association between the risk score and clinical features, the risk scores were firstly separated by subgroups of clinical features including age, stage, chromosome 3 status, metastasis, histological type, and vital status. Next, the Kruskal-Wallis or Wilcoxon test was applied to determine significant associations. Moreover, the uni-and multivariable Cox regression were applied to explore the possible prognostic factors of the risk score and clinical features. Besides, several deconvolution methodologies such as CIBERSORT [18], xCell [19], and ESTIMATE [20] were performed to decode tumor microenvironment (TME) contexture. Then, the correlation analysis was conducted to investigate the relationship between risk score and tumor microenvironment.
2.5. Gene Set Enrichment Analysis (GSEA). The GSEA was performed to discover the pathway enrichment and significant molecular mechanisms of the low-vs. high-risk groups by using "clusterProfiler" package [21]. All genes were firstly assessed by "Limma" differential analysis for the low-vs. high-risk groups and then preranked using the log2 fold change of the expression values. Then, the GSEA analysis of the low-vs. high-risk groups in UM was performed; the significant pathways were screen by FDR < 0:05. The cancer hallmarks (h.all.v7.0.symbols) set was used for the GSEA.
2.6. Chemotherapeutic Response Prediction. Drug information data was extracted from the Genomics of Drug Sensitivity in Cancer (GDSC) database (https://www.cancerrxgene .org). To explore the likelihood of chemotherapeutic drugs, the "pRRophetic" algorithm was used to predict the chemotherapeutic response for each sample, and the significant chemo drugs were selected by p < 0:05 [22]. 2.7. Statistical Analysis. Every statistical analysis was executed with R software (v.3.5.2) and corresponding statistical packages.

Prognosis Associated AIRGs.
A number of 222 autophagy-related genes and 1793 immunity-related genes were collected for univariable Cox regression in the TCGA-UVM dataset. According to the selection standard, 44 autophagy-related genes were significantly associated with survival of UM patients, which contained 27 risk genes and 17 protect genes (Figure 2(a)). In addition, 365 immunity-related genes were survival-associated, which consisted of 301 risk genes and 64 genes (Figure 2(b)). Apart from the overlapped genes, a total of 409 AIRGs were identified for the subsequent analysis.

Discovery and Validation of the Four Genes of Prognostic
Signature. Firstly, 1,000 iterations of LASSO modeling were performed to evaluate the qualified variables from the 409 AIRGs for constructing prognostic biomarker. Across 2 Journal of Oncology 1,000 iterations, 12 survival-associated AIRGs were selected out to build biomarker for appearing more than 990 times in LASSO modeling (Figure 2(c)). Moreover, based on the 5-year AUCs of different gene combinations, we found that 4 genes of signature (PRKCD, MPL, EREG, and JAG2) were arrived the max value of 0.829 (Figures 2(d) and 2(e)). The significant correlations were also found among these AIRGs ( Figure 2(f)). Finally, a 4-gene signature was identified, and the formula was computed as follows: risk score = ∑ N i=1 ð coef i × expr i Þ. The coef i means the Cox coefficient of gene and N is the number of gene. The expr i represents the relative expression of the gene in risk model. The risk scores of UM patients were calculated by risk formula and next scaled range from 0 to 1. Then, UM patients were divided into high-risk (n = 40) and low-risk groups (n = 40) by the median cutoff value. The curves of the Kaplan-Meier (KM) survival analysis manifested that high-risk patients have worse prognosis than those in the low-risk with log-rank test p value < 0.001 (Figure 3(a)). The ROC curves manifested that the 3-and 5-year AUCs were 0.916 and 0.829, respectively ( Figure 3(b)). The distribution risk scores, overall survival, vital status, and corresponding expression of AIRGs in the TCGA-UVM were illustrated in Figure 3(c). To verify the robustness and applicability of signature, the validated analyses were applied in another two GEO datasets (GSE22138 and GSE84976). These datasets were classified into subrisk groups accordingly. The curves of KM revealed the similar results which indicated that the low-risk patients had significantly longer survival time than patients in the high-risk group (GSE22138 log-rank p value = 0.0019 and GSE84976 log-rank p value < 0.001, respectively) (Figures 4(a) and 5(a)). The 3-and 5-year AUCs in GSE22138 were 0.746 and 0.713 (Figure 4(b)). Higher AUCs (3-year: 0.836 and 5-year: 0.872) were also observed in GSE84976 ( Figure 5(b)). The distribution risk scores, metastasis-free survival, vital status, and corresponding expression of AIRGs in GSE22138 were illustrated in Figure 4(c). The distribution of risk scores, overall survival, vital status, and corresponding expression of corresponding AIRGs in GSE84976 were illustrated in Figure 5(c).

Performance Evaluation of Four Genes Signature.
To evaluate the performance of four genes signature and other clinical features for prognostic prediction, uni-and multivariable cox regression were conducted by the overall or metastasis-free survival (MFS or OS) time in multiple datasets ( Table 1). The results of univariable Cox revealed that stage, age, metastasis, histological type, and risk sore were significantly 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 the TCGA-UVM 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). Compared to the time-dependent AUC estimation of traditional clinical variables, the risk score of four genes signature was close to chromosome 3 status and even superior than other variables regardless in the TCGA-UVM ( Figure 6  3 Journal of Oncology The number of risk gene is 27 The number of protect gene is 17 The number of risk gene is 301 The number of protect gene is 64 Univariate cox coefficient  Gene  EREG  JAG2  MPL  PRKCD  TNFRSF19  TUSC1  GABARAPL1  IL11  CTF1  CCL24  GAST  IL36B  ULBP1  PTGS2  CCL20  PGR  S100A2   CXCL2  IL2  PGLYRP2  SEMA3E  ESR2  LCN1  ADM2  IFNL1  CSF3  EDNRA  HDAC6  NPR1  ISG20  NR6A1  ROBO3  TCF7L2 (c)

Correlation between Signature and Clinical Features.
The correlation between risk score of signature and clinical features was evaluated, and the results of stratified analyses showed that metastasis (Figure 7(a)), vital status (Figure 7(b)), chromosome 3 status (Figure 7(c)), tumor stage (Figure 7(d)), and histological type (Figure 7(e)) were significantly correlated with the risk score. Other clinical feature, such as age (Figure 7(f)), has no relationships with the risk score. Thus, the signature was intimately correlated with these traditional clinical features.

Associations between Risk Score and Immune
Microenvironment. Through the CIBERSORT algorithm, 22 kinds of immunity cells were selected for the stratifying analysis. The box plots showed that the fractions of T cell CD8, T cells CD4, T cells regulatory, NK cells activated, monocytes, macrophages M1, and mast cells in the highrisk subgroup significantly different from those in the lowrisk group (Figure 8(a)). Next, the microenvironment score, immune score, and stromal score were calculated by applying xCell and ESTIMATE algorithms. The heat map of microenvironment score, immune score, and stromal scores between high-and low-risk groups was illustrated in Figure 8(b). The box plots suggested that the high-risk group has a higher microenvironment, immune score, and stromal score than the low-risk group (Figure 8(c)). Moreover, some important immunity-related pathways (such as MHC class-II, MHC class-I, immune checkpoint, CD8 T effector, and ICB resistance) and tumor-associated pathways (like cell cycle, EMT, and ferroptosis) were assessed by single sample GSEA (ssGSEA) algorithm [23]. The heat map (Figure 8(d)) and subgroup (Figure 8(e)) analyses revealed that the scores of these pathways in the high-risk subtype were generally higher than those in the low-risk subtype.
3.6. GSEA. According to selected criterion, we confirmed that six positive correlated pathways were enriched in the high-risk group, including epithelial-mesenchymal transition (EMT), 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 inflammatory response actively associated with the low-risk group (Figure 8(f)).

High-Risk Subgroup More Sensitive to Chemo Drugs.
Due to the resistance for regular chemotherapeutics,    Journal of Oncology scientists are struggling to discover new potential compounds for UM. Therefore, the GDSC database was used to predict sensitivity of chemo drugs. Based on the selecting criteria, 45 kinds of drugs were screened out in the TCGA-UVM set. 34 and 14 kinds of drugs were identified in the GSE22138 and GSE84976, respectively ( Table 2). The Venn plot analysis suggested that AMG.706 and JNK.Inhibitor.-VIII were the common drugs among the three datasets (Figure 9(a)). We observed significant differences for the estimation of IC50 between the high-and low-risk group, and these drugs in the high-risk group could be more sensitive to chemotherapies (Figures 9(b)-9(d)).

Discussion
Like many kinds of tumor, the prognosis of uveal melanoma is largely dependent on the early diagnosis and treatment [3]. With the growing clinical application of new therapeutic targets, the present studies are mainly focused on the detection of novel prognostic biomarkers which are often used to assess disease risk and early diagnosis [24][25][26][27]. Clinical application of prognostic biomarkers can early guide highrisk patients who undergone individual treatment and management, even prevent life-threatening metastases [28][29][30][31][32]. For example, massive previous researches were struggled to identified potential genes, miRNAs, or DNA methylations combination biomarkers via bioinformatics to predict survival of UM [33][34][35]. Hence, as far as we known, our study is the first time to investigate the prognostic role of AIRG biomarker in UM by applying multiple public datasets. Then, we constructed a robust 4 genes of signature in UM prognosis and validated in multiple independent datasets. Our prognostic biomarker can subsequently classify patients into subgroups with different survival events. Compared to traditional clinical features, our AIRG biomarker achieved higher accuracy than most of clinical indicators (e.g., stage, age, metastasis, and histology type). Notably, the results of multivariate cox regression manifested that the risk score of 4-gene signature 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 classified the UM patients into high-and low-risk groups. The Kaplan-Meier curves revealed that patients in the high-risk group have a poor survival. Among the four AIRGs, almost all of the genes have been reported to be associated with prognosis 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 [36][37][38]. Additionally, PRKCD has been regarded as a novel prognostic biomarker for ovarian cancer patient response to overall disease-specific survival [39]. MPL (thrombopoietin receptor) is a member of homodimeric class I receptor family, which mainly exerted its influence on regulation of megakaryopoiesis, formation of immune synapses, and enhancing antitumor function [40,41]. Nishimura et al. recently reported that C-MPL can work as a novel immunotherapeutic target to promote antitumor activity of T cells by mediation of cytokine pathways [42]. Therefore, MPL could be considered as a potential therapeutic target for UM patients in the future. EREG (epiregulin) belongs to the epidermal growth factor family and benefits to wound healing, tissue repair, and inflammation in normal physiology, while dysfunction of epiregulin will increase the progression of different malignancies. Previous studies revealed that the     15 Journal of Oncology overexpressed EREG closely correlates with many cancers, including prostate cancer, colon cancer, breast cancer, and bladder cancer [43][44][45][46]. What is 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 [47].
To better understand the underlying molecular functions in the subtypes of UM, the GSEA analyses were conducted to find the potential 4 genes of signature-related pathways in cancer hallmarks database. Cancer hallmark pathways such as epithelial-mesenchymal transition (EMT) and KRAS signaling were positively enriched in the high-risk expression subgroup. It is generally recognized that EMT pathway is a crucial factor to control 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 finally lead to poor prognosis [48]. Besides, Kirsten rat sarcoma viral oncogene homolog (KRAS) regard as the most lethal cancer-related proteins takes a crucial role in human aggressive cancers [23]. Previous study demonstrated that the activation of KRAS-ERK signaling can promote development and migration of uveal melanoma cells [49]. Meanwhile, we observed that interferon gamma, interferon alpha, and inflammatory response were actively enriched in the low-risk subgroup. It is widely acknowledged that interferon not only has antiviral activities but also powerful antitumor activities via various inflammatory response pathways [50]. Substantial evidences suggested that interferon alpha and gamma can inhibit the proliferation of uveal melanoma cells in vitro and indicated a potential adjuvant treatment for patients with uveal melanoma [51,52]. Therefore, we can reasonably believe that our results were consistent with the previous conclusions indicated that the high-risk subgroup has a poor survival in UM due to the positive enrichment of EMT and KRAS as well as the negative association with interferon gamma and interferon alpha. Heat map of microenvironment score, immune score, and stromal scores between the high-and low-risk groups. Red color indicates higher score, and blue color refers to lower score. (c) Box plots of microenvironment scores, immune score, and stromal scores between the high-and low-risk groups. (d) Heat map of scores for immunity-and tumor-related pathway between the high-and low-risk groups. Red color indicates higher score, and blue color refers to lower score. (e) Box plots of scores for immunity-and tumor-related pathway between the high-and low-risk groups. * means p < 0:05; * * means p < 0:01; * * * means p < 0:001; and * * * * means p < 0:0001. (f) GSEA analysis of high-vs. low-risk score groups. 16 Journal of Oncology By using various TME-decoding algorithms, we firstly observed that the high-risk group has a significant higher microenvironment score, immune score, and stromal score than the low-risk group, which indicated that a distinct TME contexture exists between the two subtypes. Next, via CIBERSORT assessment, we found that the infiltration of CD8+ T cells and M1 macrophages were highly expressed in the high-risk group. What is more, immunity-related pathways including MHC class-II, MHC class-I, immune checkpoint, CD8 T effector, ICB resistance, and tumorassociated pathways containing cell cycle, EMT, and ferroptosis were all positively enriched in the high-risk group. These observations appear counterintuitive because high infiltration of immune cells and activation of immunerelated pathways are usually correlated with good prognosis in many solid tumors [53]. However, the increasing researches revealed that an immunosuppressive environment definitely exists in UM. For eye as an immuneprivileged organ, the high infiltration of immune cells and activation of immune immune-related pathways will promote ocular tumor immune evasion. Notably, recent studies also demonstrated that the increased level of CD8+ T cells is a predictive factor for poor prognosis in UM, and the accumulation of macrophages is poorly associated with the prognosis of melanoma patients [54][55][56]. In addition, Valdor and Macian reviewed the articles and proposed that the activa-tion of autophagy could lead to strong MHC class antigen presentation and increased CD8+ T cell activity [57]. Therefore, it is logical to speculate that our signature could become an effective therapy target in UM due to the correlation with autophagy-regulated immune response.
Furthermore, in order to explore the different chemotherapeutic response between the high-and low-risk subgroups, GDSC drug database was performed to explore the candidate drugs, and the results showed that the high-risk subgroup is more sensitive to chemo drugs. Moreover, the three datasets showed that the AMG.706 and JNK.Inhibitor.VIII were the commonly enriched drugs. The estimated IC50 of the AMG.706 and JNK.Inhibitor.VIII in the highrisk group is significantly lower than the low-risk groups. The AMG.706 (motesanib) is a multikinase inhibitor which significantly inhibits VEGF-induced vascular permeability and angiogenesis and induces regression in tumor xenografts [58]. The Rosen et al. study showed that the AMG.706 has favorable pharmacokinetics and good tolerability in advanced refractory solid tumors [59]. A phase III clinical trial suggested that the AMG.706 can significantly improve progression-free survival of patients with stage IV lung cancer [60]. The JNK.Inhibitor.VIII, that is, a JNK inhibitor, can inhibit the cell apoptosis and block the phosphorylation of c-Jun and JNKAR1 which promote tumorigenesis and metastasis of cancers [ 17 Journal of Oncology proven that the JNK inhibitors could be regarded as potential targeting therapy in various type of tumors [62]. Hence, it is reasonable to believe that the AMG.706 and JNK.Inhibitor.VIII could be considered as potential chemo drugs for further clinical treatment of UM.
To sum up, our research constructed a robust 4 genes of signature prognostic biomarker in UM which could be regard as an independent prognostic model in clinical application. The patients with the high risk score could benefit more chemotherapy. However, the discoveries are obtained from the bioinformatics analysis and needed to validate in further basic and clinical studies.

Conflicts of Interest
The authors declare that they had no conflicts of interests.

Authors' Contributions
HBD was responsible for writing-original draft preparation; YL was responsible for revision and polish; YH was responsible for writing-review and editing; WL and QW were responsible for data curation; WHW and XHL was responsible for project administration and funding acquisition. All the authors commented and approved the text. Yan Liu and Huibin Du have contributed equally to this work and share first authorship.