Development of an autophagy-related gene expression signature for prognosis prediction in prostate cancer patients

Background Prostate cancer (PCa) is one of the most prevalent cancers that occur in men worldwide. Autophagy-related genes (ARGs) may play an essential role in multiple biological processes of prostate cancer. However, ARGs expression signature has rarely been used to investigate the association between autophagy and prognosis in PCa. This study aimed to identify and assess prognostic ARGs signature to predict overall survival (OS) and disease-free survival (DFS) in PCa patients. Methods First, a total of 234 autophagy-related genes were obtained from The Human Autophagy Database. Then, differentially expressed ARGs were identified in prostate cancer patients based on The Cancer Genome Atlas (TCGA) database. The univariate and multivariate Cox regression analysis was performed to screen hub prognostic ARGs for overall survival and disease-free survival, and the prognostic model was constructed. Finally, the correlation between the prognostic model and clinicopathological parameters was further analyzed, including age, T status, N status, and Gleason score. Results The OS-related prognostic model was constructed based on the five ARGs (FAM215A, FDD, MYC, RHEB, and ATG16L1) and significantly stratified prostate cancer patients into high- and low-risk groups in terms of OS (HR = 6.391, 95% CI = 1.581– 25.840, P < 0.001). The area under the receiver operating characteristic curve (AUC) of the prediction model was 0.84. The OS-related prediction model values were higher in T3-4 than in T1-2 (P = 0.008), and higher in Gleason score  > 7 than  ≤ 7 (P = 0.015). In addition, the DFS-related prognostic model was constructed based on the 22 ARGs (ULK2, NLRC4, MAPK1, ATG4D, MAPK3, ATG2A, ATG9B, FOXO1, PTEN, HDAC6, PRKN, HSPB8, P4HB, MAP2K7, MTOR, RHEB, TSC1, BIRC5, RGS19, RAB24, PTK6, and NRG2), with AUC of 0.85 (HR = 7.407, 95% CI = 4.850–11.320, P < 0.001), which were firmly related to T status (P < 0.001), N status (P = 0.001), and Gleason score (P < 0.001). Conclusions Our ARGs based prediction models are a reliable prognostic and predictive tool for overall survival and disease-free survival in prostate cancer patients.

Several studies reported that autophagy could play a role in tumor progression or tumor suppression in different stages of cancers [3,4]. However, the role of autophagy in tumorigenesis is still rudimentary.
Prostate cancer (PCa) is a common malignancy of the urinary system and the second cause of cancerrelated death of males in western developed countries [5]. In China, the annual incidence of PCa was more than 60,000 cases and 26,600 patients who succumbed to PCa in 2015 [6]. The majority of early-stage PCa patients have an excellent prognosis with a low mortality rate [7]. However, there are still a large number of PCa patients who develop the resistance to androgen deprivation therapy (ADT) and become castrationresistant PCa (CRPC), which results in a short survival time [8].
The relationship between autophagy and multiple biological processes of prostate cancer has been previously reported [9]. For instance, in the early stage of PCa, autophagy may increase tumor cell death. However, elevated autophagy promotes prostate cancer invasion and progression and reduces the damage of chemotherapy drugs in the late stage. Cao et al. showed that the induction of autophagy might increase susceptibility to radiation in prostate cancer cell lines [10,11]. However, large-scale gene expression signature has rarely been used to investigate the association between autophagy and prognosis in prostate cancer. To better understand the impact of tumor genetic composition on clinical outcomes, the Cancer Genome Atlas (TCGA) database has been established for discovering gene signatures.
There are many research and prognostic models based on gene expression profiles in prostate cancer, such as lncRNAs and miRNAs [12]. Nevertheless, prognostic models for autophagy-related genes have not been reported. In this study, we used gene expression microarray data obtained from TCGA to develop autophagy-related gene expression signature and develop a prognostic model as an independent index for overall survival (OS) and disease-free survival (DFS).
RNA-seq data for prostate cancer patients were downloaded from the TCGA data portal (https ://tcga-data.nci. nih.gov/tcga/), which contains 485 prostate cancer and 51 adjacent non-tumor tissues. We searched the cBio Cancer Genomics Portal (http://cbiop ortal .org) to identify the clinical data, including OS and DFS.

Differentially expressed ARGs and enrichment analysis
Data analysis of differential expression of ARGs between PCa and their non-tumor counterparts was performed using package limma in R, with thresholds of |log 2 fold change (FC)|> 2 and adjusted P-value < 0.05. Then, we performed gene ontology (GO) enrichment analyses to find the major biological attributes of differential expression ARGs. The visual GO enrichment maps of annotation analysis results were performed by R with the "ggplot2" and "GOplot" packages.

Construction of prognostic signature based on ARGs
Univariate Cox and multivariate Cox regression analyses were performed to find out the OS-related and DFS-related ARGs in PCa. Then, the OS-related and DFS-related prediction formulas were applied to build prognostic models using package "glmnet" based on the multivariate Cox regression. The survival analysis was assessed by Kaplan-Meier (K-M) methods to compare the high-risk and low-risk groups according to predictive signatures. Finally, the predictive value of prognostic prediction models was evaluated by areas under the curve (AUC) of the receiver-operator characteristic (ROC) curve using package "survivalROC" in R.

Statistical analysis
All of the statistical tests were done with R 3.3.1 (https ://www.r-proje ct.org/) and GraphPad Prism 7 (San Diego, CA, USA). All analyses performed were twosided, and statistical significance was defined as a P-values < 0.05.

Differentially expressed ARGs between prostate cancer and adjacent non-tumor tissues
A total of 485 primary PCa patients with RNA-seq data and clinical follow-up information were involved in the present study. Among 234 autophagy-related genes, there were 13 differentially expressed ARGs, including 5 up-regulated (ATG9B, BIRC5, CAMKK2, CDKN2A, and NKX2-3) and 8 down-regulated ARGs (DNAJB1, FAM215A, HSPB8, ITGB4, ITPR1, NRG1, NRG2, and TP63), with thresholds of |log 2 fold change (FC)| > 2 (Fig. 1a). Then, the volcano plot and box plot were visualized to show the expression pattern of the differentially expressed ARGs between PCa and non-tumor tissues (Fig. 1b, c). GO enrichment analysis of differentially expressed autophagy-related genes GO enrichment analysis was performed according to the differentially expressed ARGs. According to the results of DAVID, we found that the top enriched GO terms for biological processes were autophagy, process utilizing autophagic mechanism, and odontogenesis of dentin-containing tooth. The heatmap of the relationship between ARGs and GO enrichment analysis was also displayed (Fig. 2a, b). We divided the 485 prostate cancer cases into highand low-risk groups according to the median values of the OS-related prediction model. Kaplan-Meier survival curves showed that low-risk group had a lower mortality rate than high-risk group (HR = 6.391, 95% CI = 1.581-25.840, P < 0.001) (Fig. 3a). The ROC curves of OS-related predictive signatures were demonstrated in Fig. 3b, with
Among 485 primary PCa patients in the present study, only two of them had distant metastasis. Therefore, the relationship between the M status and prediction model has not been analyzed.

Discussion
Much evidence has indicated that autophagy participates in multiple signaling pathways to play a role in the proliferation and invasion of prostate cancer [9,13]. Additionally, Blessing et al. demonstrated that four core autophagy genes (ATG4B, ATG4D, ULK1, and ULK2) regulate androgen receptor (AR) activity, thereby affecting the biological behavior of prostate cancer [14].
In the present study, we mined the expression profiles of ARGs from the TCGA database and aimed to analyze the association between ARGs and the prognosis of prostate cancer patients. Firstly, we screened 13 differentially expressed ARGs between prostate cancer and non-tumor tissues, that many of them play a role in the biological processes by GO term analysis. Then, a total of 14 OS-related ARGs were found in the univariate Cox regression analysis. Further multivariate Cox regression analysis was performed to determine five OS-related ARGs (FAM215A, FDD, MYC, RHEB, and ATG16L1) and construct the OS-related prediction model, which could be an independent prognostic indicator for PCa patients. Handle et al. found that MYC activity is closely related to AR, which regulates the growth of anti-androgen resistant cell lines [15]. Kobayashi et al. demonstrated that RHEB mRNA and protein expression was higher in more aggressive prostate cancer cell lines (PC3 and DU145) compared with the less aggressive LNCaP. Moreover, inhibition of RHEB can lead to the suppressed proliferation of prostate cancer cell lines [16]. Previous research analyzed the relationship between genetic variants of the autophagy pathway and clinical outcomes in 458 prostate cancer patients, which indicated that high expression of ATG16L1 was correlated with lower tumor aggressiveness and favorable prognosis [17]. Fu et al. reported that high expression of FAM215A was associated with low tumor grades, early disease stages, and favorable overall survival in epithelial ovarian cancer [18]. FDD is a component of FMNL3, and high expression of FMNL3 associated with cancer cell migration, invasion, and unfavorable prognosis in tongue squamous cell carcinoma [19]. However, the function of FAM215A and FDD gene has not been reported in prostate cancer, indicating that functional studies on these genes may help us to understand the prognosis-related biological behavior of bladder cancers more accurately.
In the present study, a total of 22 ARGs were significantly associated with DFS of PCa in multivariate Cox regression analysis, including ULK2, NLRC4, MAPK1, ATG4D, MAPK3, ATG2A, ATG9B, FOXO1, PTEN, . Previous research has shown that ULK2 and ATG4D were hub autophagy genes, which are necessary for maximal androgen-mediated autophagy and cell proliferation and also associated with poor prognosis in PCa [14]. Li et al. demonstrated that MAPK1 plays an important role in regulating cancer cell invasion and metastasis in vitro and in vivo [20]. Forkhead box transcription factor-1 (FOXO1) is a tumor suppressor that is downregulated in human prostate cancer, which acts as a repression target of EZH2 and an essential mediator of EZH2 inhibition-induced cell death [21]. PTEN is one of the most commonly altered tumor suppressor genes in prostate cancer, which negatively regulates the PI3K/AKT/mTOR signaling pathway. PTEN deletion is associated with poorer cancer-specific outcomes, increasing stage, and higher Gleason score [22]. Chuang et al. suggested that HDAC6 has anti-cancer activity in prostate cancer, which participates in regulating the cRaf-PP1-ERK signaling pathway and contributing to M phase cell-cycle transition [23]. Many studies have shown that multiple oncogenes promote PCa cell proliferation, migration, invasion, and inhibiting apoptosis through activating the PI3K-AKT-mTOR signaling pathway [24]. Chen et al. found that the high levels of NPRL2 gene expression in prostate cancer cells promote resistance to EVS (an inhibitor of the mTOR) by enhancing autophagy [25]. In addition, TSC1 was significantly associated with DFS in PCa, which is an essential component of the PI3K/ AKT/mTOR signaling pathway [26]. Among the 22 DFS-related ARGs, except for those mentioned above, other ARGs are either poorly investigated or have not been reported, which means our findings suggested further research for them is imperative. The OS and DFS-related prediction model values were both associated with the T stage and Gleason score in PCa patients, higher in T3/4 than in T1/2, and higher in Gleason score > 7 than ≤ 7. Patients with T3 or T4 stage are also known as locally advanced prostate cancer. Krimphove et al. reported that PCa with T3 or T4 had a worse overall survival [27]. Meanwhile, the Gleason score is the sum of the two most common grade patterns in PCa, which act as the single most potent predictor of PCa outcomes [28]. DFS-related prediction model values were higher in N1 than in N0. N1 is defined as regional lymph node metastasis in AJCC/UICC N category. Jin