Identification of an Autophagy-Related Prognostic Signature for Clear Cell Renal Cell Carcinoma

Abnormal autophagy is closely related to the development of cancer. Many studies have demonstrated that autophagy plays an important role in biological function in clear cell renal cell carcinoma (ccRCC). This study aimed to construct a prognostic signature for ccRCC based on autophagy-related genes (ARGs) to predict the prognosis of ccRCC. Differentially expressed ARGs were obtained from ccRCC RNA-seq data in The Cancer Genome Atlas (TCGA) database. ARGs were enriched by gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG). The prognostic ARGs used to construct the risk score models for overall survival (OS) and disease-free survival (DFS) were identified by Cox regression analyses. According to the median value of the risk score, patients were divided into a high-risk group and a low-risk group. The OS and DFS were analyzed by the Kaplan-Meier method. The predictive accuracy was determined by a receiver operating characteristic (ROC) curve analysis. Additionally, we performed stratification analyses based on different clinical variables and evaluated the correlation between the risk score and the clinical variables. The differentially expressed ARGs were mainly enriched in the platinum drug resistance pathway. The prognostic signatures based on 11 ARGs for OS and 5 ARGs for DFS were constructed and showed that the survive time was significantly shorter in the high-risk group than in the low-risk group (P < 0.001). The ROC curve for OS exhibited good predictive accuracy, with an area under the curve value of 0.738. In the stratification analyses, the OS time of the high-risk group was shorter than that of the low-risk group stratified by different clinical variables. In conclusion, an autophagy-related signature for OS we constructed can independently predict the prognosis of ccRCC patient, and provide a deep understanding of the potential biological mechanisms of autophagy in ccRCC.


INTRODUCTION
Renal cell carcinoma (RCC) is a malignant tumor originating from the renal tubular epithelium. It is a common malignant tumor of the urinary system, and clear cell renal cell carcinoma (ccRCC) is the most common subtype (1). Currently, surgical resection is the main treatment for ccRCC, but ccRCC has a poor prognosis and is likely to recur (2). Common clinical variables, such as the TNM stage, have good prognostic value (3). However, because of tumor heterogeneity the TNM stage cannot accurately predict the prognosis of patients (4). Therefore, the discovery of new molecular targets in ccRCC is the first requirement for achieving early diagnosis and improving the survival rate of ccRCC patients.
Autophagy is a highly conserved intracellular self-digestion process that maintains cellular homeostasis through lysosomes (5). Autophagy plays a key role in maintaining the balance between the synthesis and degradation of cell components (6). Dysregulation of autophagy is closely related to cancer (7). Initially, autophagy can prevent or delay tumor formation, but once tumors are formed, autophagy can promote tumor progression and protect cancer cells from environmental damage (8).
Many studies have demonstrated the role of autophagy in ccRCC. Autophagy is a therapeutic target for renal cancer (9,10). Studies have shown that promoting autophagy can inhibit the progression of ccRCC (11,12). The lower the level of autophagy is, the higher the stage and grade of ccRCC (13). Many drugs have been developed to promote autophagy in renal carcinoma and have achieved good therapeutic effects (14). However, these studies mainly focused on the influence of autophagy on the progression and treatment of ccRCC, and few researchers have studied the role of autophagy in the prognosis of ccRCC.
In this study, 45 differentially expressed autophagy-related genes (ARGs) were obtained from the expression data of patients in the kidney clear cell carcinoma (KIRC) cohort in The Cancer Genome Atlas (TCGA) database, and the biological functions of these differentially expressed ARGs were analyzed. These analyses provided further insight into the roles of these ARGs in ccRCC. Importantly, we constructed risk score models based on 11 prognostic ARGs for overall survival (OS) and five prognostic ARGs for disease-free survival (DFS) and found that the autophagy-related signature can independently predict the prognosis of ccRCC patients without considering clinical variables, suggesting that those autophagy-related signatures are reliable prognostic marker in ccRCC patients.

Data Sources
We obtained 232 ARGs from the HADb (Human Autophagy Database, http://www.autophagy.lu/) and then downloaded the FPKM-standardized RNA-seq data and the clinical and OS information from the KIRC cohort in the TCGA database (https://portal.gdc.cancer.gov/). Table 1 shows the basic clinical characteristics of the patients with ccRCC in the TCGA database. We obtained the DFS data of 431 ccRCC patients from cBioportal (https://www.cbioportal.org/). A total of 222 ARGs with expression values were obtained.

Enrichment Analysis of Differentially Expressed ARGs
We used a false discovery rate (FDR) < 0.05 and a | log 2 fold Change| (logFC) > 1 as screening criteria to obtain the differentially expressed ARGs. To better understand the role of differentially expressed ARGs, we used the "cluster Profiler" package (15) for enrichment analyses and then used the "GOplot package" (16) for visualization.

Construction of the Autophagy-Related Prognostic Signature
We first obtained the prognostic ARGs in ccRCC by univariate Cox regression analysis and then performed multivariate Cox regression analysis and an optimized risk score model with the step function. The risk score was calculated as follows: where Coef is the coefficient, and x is the expression value of each selected ARG. This formula was used to calculate the risk score for every ccRCC patient. Then we performed Cox regression analyses to demonstrate whether the autophagy-related signature was an independent prognostic factor in ccRCC patients.

Gene Set Enrichment Analysis (GSEA)
Patients were divided into high-and low-risk groups according to the value of the risk scores. GSEA (17) was used to analyze which pathways genes are primarily enriched. GSEA was performed using GSEA3.0 (http://www.broad.mit.edu/gsea/). Differences for which the nominal P < 0.05 and the FDR < 0.25 were considered statistically significant.

Statistical Analysis
Statistical analyses were performed with R software (Version 3.5.1). The Wilcox signed-rank test was used to compare the expression levels of differentially expressed ARGs between cancer tissues and normal tissues and the expression of 11 prognostic ARGs in the high-and low-risk groups. Student's t-test was used to compare the correlation between the risk score and clinicopathological variables. Cox regression analyses were used to screen genes for inclusion in the risk score model. ccRCC patients were divided into a high-risk and a low-risk groups according to the median value of the risk score, and OS and DFS of patients were analyzed by the Kaplan-Meier method and logrank test. Receiver operating characteristic (ROC) curve analysis was performed with the "survivalROC" package. P < 0.05 was considered to be statistically significant.

Differentially Expressed ARGs in Cancer
Tissues and Normal Tissues Figure 1 shows the flow chart of our research process. We first obtained 232 ARGs from the HADb database and then downloaded the RNA-seq and clinical and prognostic data of 530 patients in the KIRC cohort from the TCGA database. We ultimately obtained the expression data of 222 ARGs. With FDR < 0.05 and | log 2 FC| >1as the screening criteria, 45 differentially expressed ARGs were obtained, including 9 downregulated ARGs and 36 upregulated ARGs (Figures 2A,B). The expression of the differentially expressed ARGs between cancer tissues and normal tissues was visualized ( Figure 2C)

Functional Annotation of the Differentially Expressed ARGs
Functional enrichment analysis was performed with the 45 differentially expressed ARGs. In the biological processes, the ARGs were mainly enriched in autophagy, positive regulation of peptidase activity, regulation of apoptotic signaling pathways, regulation of cell growth, autophagy of mitochondrion, autophagosome assembly, etc. In the cellular components, the ARGs were mainly enriched in autophagosomes, autophagosome membranes, cytosolic part, mitochondrial outer membrane, organelle outer membrane, vacuolar membrane, basal plasma membrane, etc. In the molecular functions, the ARGs were mainly enriched in protein phosphatase binding, P53 binding, phosphatase binding, chemokine receptor activity, cytokine activity, integrin binding, peptidase activator activity, receptor ligand activity, etc. ( Table 2 and Figure 3A). In the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways, the ARGs were mainly enriched in platinum drug resistance, epidermal growth factor receptor (EGFR) tyrosine kinase inhibitor resistance, the epidermal growth factor receptor (ErbB) signaling pathway, endocrine resistance, the mitogenactivated protein kinase (MAPK) signaling pathway, calcium signaling pathway, and cytokine-cytokine receptor interaction ( Table 2 and Figures 3B,C). The z scores of these KEGG pathways were >0, indicating that the ARGs were upregulated in these pathways.

Construction of a Prognostic Signature Based on the Prognostic ARGs for OS
After analyzing the expression and functions of the differentially expressed ARGs in ccRCC, we constructed a risk score model for the prediction of the prognosis of patients with ccRCC. After univariate Cox regression analysis, 23 ARGs were associated with the prognosis of ccRCC patients (Figure 4). After multivariate Cox regression analysis, 11 ARGs were identified and used to construct a prognostic signature for OS ( Table 3). The risk score was calculated as follows:

The Correlation Between the Autophagy-Related Signature for OS and Prognosis of ccRCC Patients
To determine the ability of the autophagy-related signature for OS to predict the prognosis of ccRCC patients, Kaplan-Meier analysis was performed to evaluate the OS outcomes in the two groups. The OS rate of patients in the highrisk group was significantly lower than that of patients in the low-risk group (P = 1.221e−15, Figure 5A), and the 5years survival rates of patients in the high-and low-risk groups were 40.1 and 78.8%, respectively. The violin plot shows the expression of the eleven ARGs in the two groups. BID, RGS19, CASP4, IFNG, ATG16L2, and EIF4EBP1 were highly expressed in the high-risk group, and PRKCQ, BAG1, CX3CL1, ERBB2, and BNIP3 were highly expressed in the low risk group ( Figure 5B). The risk score of patients in the high-and low-risk groups were visualized ( Figure 5C). As the risk score increased, an increasing number of patients died ( Figure 5D). These results showed that the risk score accurately reflect the survival of patients and that the autophagyrelated signature for OS accurately predicts the prognosis of patients.
To determine whether the autophagy-related signature for OS is an independent prognostic factor for ccRCC patients, we performed Cox regression analysis. Univariate Cox regression analysis showed that age, stage, grade, T stage, M stage, and risk score were significantly associated with OS in ccRCC patients ( Figure 6A). Multivariate Cox regression analysis showed that age, stage, grade and risk score were independent factors influencing ccRCC prognosis ( Figure 6B). Then, a ROC curve was constructed to determine the predictive accuracy of the autophagy-related signature. The area under the curve (AUC)   of the autophagy-related signature for OS was 0.738, indicating good predictive accuracy ( Figure 6C). Because of the survival differences between the high-risk and low-risk groups, we conducted GSEA to study the functional differences between these groups. The regulation of autophagy and renal cell carcinoma pathways were significantly enriched in the low-risk group (Figures 6D,E), indicating that the regulation of autophagy was mainly involved in low-risk ccRCC patients.

Validation of Prognostic Signature Based on Prognostic ARGs for OS
To validate the applicability of the prognostic signature for OS we constructed based on the entire TCGA data set, we randomly divided the 530 ccRCC patients in the entire data set into a training set (n = 265) and a validation set (n = 265). According to the formula, we calculate the risk score for each patient, and the patients in the training set and the validation set were divided into high-and low-risk groups based on the median value of the risk score. Consistent with the results observed in the entire data set, the OS rate of patients in the high-risk group was lower than that of the low-risk group in the training set (P = 3.023e-12, Figure 7A), and the prognosis of the high-risk group was worse than that of the low-risk group in the validation set (P = 1.341e-05, Figure 7C). The ROC curves of the training set and validation set also show good performance. The AUCs for the 1-year, 3years, and 5-years OS of the training set were 0.775, 0.785, and

Role of the Signature for OS in the Prognosis of ccRCC Patients Stratified by Clinicopathological Variables
To investigate the prognostic value of the signature for OS in ccRCC patients stratified by clinicopathological variables, ccRCC patients were stratified according to age, gender, grade, stage, M stage, and T stage. For all different stratifications, the OS time of the high-risk group was shorter than that of the low-risk group (Figure 8). These results suggest that the autophagy-related signature for OS can predict the prognosis of ccRCC patients without the need to consider clinicopathological variables.

Relationship Between the Prognostic Signature for OS and Clinicopathological Variables
To determine whether the autophagy-related prognostic signature for OS affects the progression of ccRCC, we analyzed the correlations between the autophagy-related prognostic signature for OS and clinicopathological variables. The risk score of G3-4 was higher than that of G1-2 (P = 1.653e−06, Figure 9A), the risk score of M1 was higher than that of M0 (P = 0.002, Figure 9B), the risk score of N1 was higher than that of N0 (P = 0.004, Figure 9C), the risk score of stage III-IV was higher than that of stage I-II (P = 1.102e−05, Figure 9D), and the risk score of T3-4 was higher than that of T1-2 (P = 5.676e−05, Figure 9E). These results suggested that the higher the risk score is, the greater the degree of malignancy of ccRCC. Thus, the prognostic signature for OS could accurately predict the progression of ccRCC.

Construction of a Prognostic Signature Based on Prognostic ARGs for DFS
Considering the significance of DFS in the prognosis of ccRCC, we also established a prognostic signature for DFS. We obtained DFS data for ccRCC from cBioportal, including 431 patients. After univariate Cox regression analysis, we obtained 19 ARGs significantly correlated with DFS in ccRCC patients. After multivariate Cox regression analysis, we obtained 5 ARGs and constructed a prognostic signature: risk score = (0.5163 × BID) + (−0.4748 × BAG1) + (0.1084 × APOL1) + (−0.6522 × NKX2-3) + (0.3866 × EIF4EBP1). The risk score for each patient was calculated according to the formula, and the patients in the entire data set were divided into high-and low-risk groups according to the median value of risk score. K-M analysis showed that the DFS time of the high-risk group was significantly shorter than that of the low-risk group (P = 9.177e−10, Figure 10A), and ROC analysis showed that the AUCs for 1-year, 3-years, and 5-years DFS were 0.745, 0.754, and 0.756, respectively ( Figure 10D). These results showed that the prognostic signature for DFS can also predict the prognosis of ccRCC patients well.
To verify the applicability of the prognostic signature for DFS, 431 patients were randomly divided into a training set (n = 216) and a validation set (n = 215). The risk score for each patient was calculated according to the formula, and the patients were divided into high-and low-risk groups according to the median value of risk score. Consistent with the results obtained from the entire data set, patients in the high-risk group in the training (P = 9.543e-08, Figure 10B) and validation sets (P = 2.277e-05, Figure 10C) had shorter DFS times than the low-risk groups. In the training set, the AUCs for 1-year, 3-years, and 5-years DFS were 0.746, 0.759, and 0.796, respectively ( Figure 10E). In the validation set, the AUCs for 1-year, 3-years, and 5-years DFS were 0.741, 0.748, and 0.711, respectively ( Figure 10F).
To determine whether the prognostic signature for DFS can independently predict the prognosis of ccRCC, we performed Cox regression analysis. Univariate Cox analysis showed that grade, stage, T stage, M stage and risk score were significantly correlated with DFS of ccRCC ( Figure 10G). Multivariate Cox regression analysis showed that grade, stage and risk score were independent factors influencing the DFS of ccRCC ( Figure 10H). These results showed that the autophagy-related signature for DFS could predict the DFS of patients well.

DISCUSSION
ccRCC is a disease that is likely to recur and to have a poor prognosis. Precise diagnostic and therapeutic biomarkers are urgently needed. Many studies have found that autophagy is significantly related to the occurrence and progression of cancer. However, current studies have focused on the influence of ARGs in cancer development and treatment (18)(19)(20), and few have addressed the prognostic value of ARGs in cancer.
Recently, many studies have predicted the prognosis of ccRCC by constructing a prognostic signature based on miRNAs and lncRNAs (21)(22)(23), but few have reported the construction of an autophagy-related prognostic signature to predict the prognosis of ccRCC or a bioinformatic exploration of the possible role of ARGs in ccRCC. Although Wan et al. (24) constructed a prognostic model of autophagy-related genes for ccRCC with partial data sets from TCGA (n = 266). We used the entire TCGA data set (n = 530) for model construction. The differentially expressed ARGs in our two studies are also different, so our two constructed models are completely different. In our study,   in addition to building the prognostic model for ccRCC, we also examined the role of ARGs in ccRCC. First, we analyzed the differentially expressed ARGs in ccRCC and normal kidney tissues, and obtained 45 differentially expressed ARGs. GO and KEGG analyses indicated that the differentially expressed ARGs were mainly enriched mainly in platinum drug resistance. Studies have shown that the induction of autophagy in cancer can increase cisplatin resistance (25,26), consistent with our results, suggesting that these ARGs can promote the progression of ccRCC through platinum drug resistance. However, further experiments are needed to verify the role of autophagy in ccRCC. Treating patients by inducing or inhibiting autophagy remain controversial (27). The overexpression of ARGs in lung cancer tissues can promote the progression of lung cancer (28). ARG Beclin-1 is highly expressed in colorectal cancer, and Beclin-1 high expression is positively correlated with clinicopathological variables and predicts good prognosis (29). However, some studies have found that Beclin-1 is expressed at low levels in bladder cancer (30), salivary gland adenoid cystic carcinoma (31), and pancreatic ductal adenocarcinoma (32) and that patients with low Beclin-1 expression have shorter survival times. ARGs are expressed differently in different cancers, possibly due to tumor heterogeneity.
Autophagy is closely related to the prognosis of cancer patients, so it is important to find a prognostic signature for ccRCC patients. We used univariate Cox regression to analyze ARGs associated with the prognosis of patients with ccRCC. 23 ARGs were found to be significantly associated with the prognosis of ccRCC. Then, multivariate Cox regression analysis was performed, and 11 ARGs (BID, ERBB2, CASP4, PRKCQ, BAG1, IFNG, ATG16L2, EIF4EBP1, CX3CL1, RGS19, and BNIP3) were identified for inclusion in the risk score model for OS. In vivo and in vitro experiments showed that knockdown of CASP4 leads to cell migration and impairs cell-matrix adhesion (33). Silencing BAG1 in breast cancer cells increases resistance to tamoxifen and reduces apoptosis by activating the PI3K/Akt/mTOR signaling pathway (34), and the overexpression of ATG16L2 is related to poor prognosis in epithelial cancer (35). The overexpression of EIF4EBP1 is related to shorter recurrencefree survival in breast cancer patients (36). CX3CL1 is highly expressed in esophageal cancer and can promote its metastasis (37). RGS19 can effectively inhibit Ras-related carcinogenesis in lung cancer (38). BNIP3 has an anticancer effect and is negatively correlated with the expression of the m6A demethylase FTO in breast cancer; BNIP3 can slow down the growth and metastasis of FTO-overexpressing tumors (39). Our GSEA results showed that the regulation of autophagy was mainly enriched in the lowrisk group, indicating that autophagy has a greater regulatory role and influence in the low-risk group than in the high-risk group. However, whether autophagy is negatively regulated in the lowrisk group requires further study. The next study demonstrated that the autophagy-related signature for OS can independently predict the prognosis of ccRCC patients and is a good predictor of ccRCC patient prognosis. The higher the risk score was, the worse the prognosis and the greater the degree of malignancy. Through the internal validation with the training set and the validation set, the autophagy-related signature for OS we constructed was shown to have good predictive performance.
We also found that the autophagy-related signature for OS can predict the prognosis of ccRCC patients without the need to consider clinicopathological variables. In addition, we also established a prognostic signature for DFS and conducted internal validation. The autophagy-related signature for DFS can independently and accurately predict the prognosis of ccRCC patients. However, our research also has some limitations. First, the mechanisms of action of the ARGs in ccRCC need validation in vivo and in vitro experiments. Second, we only used data from the TCGA database for this analysis, and we should validate the results in other databases. Although we have performed internal verification, we still need to perform external validation in other cohorts to test the applicability of the autophagyrelated signature.
In conclusion, we constructed the autophagy-related signature for OS and DFS that can independently predict the prognosis of ccRCC patients and provide new therapeutic targets for ccRCC. We have developed a deep understanding of the biological mechanisms and clinical significance of the identified ARGs in ccRCC, but further experiments are still needed to verify our findings in the future.