Identification of an independent autophagy-gene prognostic index for papillary renal cell carcinoma.

BACKGROUND: Autophagy was a significant catabolic process which played a critical role in the maintenance of cellular homeostasis and viability in a stressed state. The dysregulation of autophagy was correlated with various diseases. The aim of our study was to develop a prognostic signature for papillary renal cell carcinoma (RCC).
METHODS: First, 40 differently expressed genes related with autophagy (ARGs) were examined via high-throughput sequencing and large-scale databases. Then, functional enrichment analysis was performed to explore the biological attributes of these ARGs. The Cox proportional hazard regression hinted that four ARGs (P4HB, BIRC5, NGR1 and PRKN) were significantly correlated with overall survival (OS). Thus, we got genes with prognostic value. Finally, a prognostic index (PI) was constructed.
RESULTS: After identifying the 4 ARGs, we profiled our risk signature. Based on the PI we developed, papillary RCC patients were stratified into high-risk and low-risk groups. High-risk patients had significant shorter OS than low-risk patients (P<0.001) and the mortality of high scoring patients was higher than low scoring patients. Additionally, we explored the relationship between the 4 ARGs and clinical parameters and found that the expression of P4HB, BIRC5 and NGR1 was correlated with clinicopathological features.
CONCLUSIONS: Our study suggested that the four-gene signature was an independent prognostic factor which could act as a novel indicator for the prognosis of papillary RCC. 2020 Translational Andrology and Urology. All rights reserved.

Introduction

Renal cell carcinoma (RCC), as the most common malignancy in renal neoplasia, caused nearly 12,000 deaths annually worldwide (1). Papillary RCC, accounting for 15% in RCC, was one of the common subtypes of RCC. With lower histological grade and distant organ metastasis, papillary RCC was generally considered to have better disease-free survival (DFS) and overall survival (OS) than clear RCC. Recommended therapy for papillary included partial nephrectomy, nephrectomy, radiofrequency ablation, molecular target therapy and administration of immune checkpoint inhibitor (2). Despite of great progress with regard to therapy and pathogenesis of papillary achieved recent years, prognosis of papillary RCC was still worse. So, the exploration of underlying molecular mechanisms of tumorigenesis and metastasis in papillary RCC was still urgent nowadays.

Autophagy, a highly conserved and evolutionarily ancient biological process in yeast and eukaryotic cells, played a significant role in maintaining the homeostasis of cell in stressed states. The fundamental process of autophagy contained the activation of PI3KC3 complex and initiation of isolated membrane, elongation and closure of autophagosome membrane and formation of autolysosome (3,4). The dysregulation of autophagy was now found accounted for many pathological neurodegenerative diseases (5), cardiomyopathy (6). Some researched have pointed out the momentous role of autophagy in the tumorigenesis and metastasis in papillary RCC (7-9). Besides, the administration of chloroquine (CQ), an inhibitor of formation of autophagosomes, have been found to increase the survival of patients and have synergetic effects when combining with other molecular target drugs (10). However, the concrete role of autophagy in RCC still remained controversial and mysterious.

In this study, we tried to illustrate a comprehensive relation between autophagy and papillary RCC. At first, we identified 40 autophagy related genes (ARGs) which differently expressed genes between adjacent tumor tissues and tumor tissues. Then, by applying COX regression model, 4 core ARGs (P4HB, BIRC5, NGR1 and PRKN) were selected for the further analysis. Finally, we developed prognostic index (PI) based on these 4 core ARGs as an independent index for predicting the OS in papillary RCC and clinical parameters were extracted to reveal the relationship between these 4 core ARGs and clinical status. We present the following article in accordance with the MDAR reporting checklist (available at http://dx.doi.org/10.21037/tau-20-906).

Materials

Data extraction

We identified 234 ARGs from The Human Autophagy Database (HADb, http://www.autophagy.lu/index.html), an autophagy-related database which recorded ARGs that have been described in literature. Gene profile of level 3 RNA-sequence in The Cancer Genome Atlas (TCGA) cohort was downloaded by using R Studio and only data regarding papillary RCC was extracted for further analysis by checking pathological examination. Then, some clinical parameters like gender, Stage, pathological grade, TMN stage, vital status and survival time were also acquired.

Differently expressed gene enrichment analysis

By using EdgeR package in R Studio, 40 differently expressed ARGs were estimated by comparing tumor tissues and adjacent tumor tissues, with the cut-off criterion of fold change >2 and P value <0.01. Then, in order to explore the major biological process of these 40 ARGs, these genes were utilized to conduct the gene functional enrichment analyses including Kyoto Encyclopedia of Genes and Genomes (KEGG) and gene ontology (GO). We utilized Database for Annotation, Visualization, and Integrated Discovery (DAVID, https://david.ncifcrf.gov/) to identify enriched KEGG and GO themes.

Establishment of a PI based on ARGs

In order to find out ARGs which gene expression were associated with OS in patients with papillary RCC, several genes were identified by using univariate Cox regression analyses. Then, we conducted multivariate Cox regression analysis to obtain genes that might be an independent indicator. Totally, 4 ARGs (P4HB, BIRC5, NGR1 and PRKN) were selected for the construction of PI. The construction of PI based on a linear combination of the relative expression level of genes multiplied regression coefficients. The risk coefficients of ARGs were calculated by utilizing a multivariable Cox proportional hazards model. According to the median PI value as the risk cutoff value, all the patients with papillary RCC were divided into low- and high- risk groups. OS of these two groups were analyzed by Kaplan-Meier method and the difference between low and high groups were tested by log-rank test. Finally, we plotted receiver operating characteristic curves (ROCs) to verify the prediction value of the model.

Statistical analysis

The relationship between clinical parameters and gene expression profile was analyzed by Fisher’s exact test or chi-square test. Continuous variables were compared by Student’s t-test. The OS analyses were performed by the Kaplan-Meier method and log-rank test. We carried out all statistical analyses using R Studio 3.6.1 (https://www.r-project.org/), SPSS 24.0 (Chicago, IL, USA) and GraphPad Prism 8 (San Diego, CA, USA). Probability P value <0.05 was considered to be statistically significant.

Ethical statement

All data from TCGA are publicly available and their use do not require the approval of a local ethics committee. All procedures performed in this study were in accordance with the Declaration of Helsinki (as revised in 2013).

Results

Validation of ARGs in papillary RCC

Totally, after checking the pathological examination of all the RCC tissues, 289 papillary RCC tissues and 32 adjacent tumor tissues were downloaded from TCGA (). Among them, 45 patients with gene expression profile and clinical information were included for the further analysis. We compared the differently expressed genes between tumor tissues and adjacent tumor tissues by applying EdgeR package in R Studio. With P value <0.01 and fold change >2 as the cut-off criteria, 31 up-regulated ARGs and 9 down-regulated ARGs were obtained in the end. Volcano map, heat map and box plots were carried out to visualize the differently expressed ARGs ().

Table 1

Clinicopathological characteristics of patient samples

Characteristics	Number of cases (%)
Age (years)
≥60	167 (57.4)
<60	121 (41.6)
Unknow	3 (1.0)
Gender
Male	214 (73.5)
Female	77 (26.5)
Pathologic grade
Unknow	291 (100.0)
Clinical stage
I	173 (59.5)
II	21 (7.2)
III	52 (17.9)
IV	15 (5.2)
Unknow	30 (10.3)
T classification
T1	194 (66.7)
T2	33 (11.3)
T3	60 (20.6)
T4	2 (0.7)
Tx/unknow	2 (0.7)
N classification
N0	50 (17.2)
N1	24 (8.2)
N2	4 (1.4)
Nx/unknow	213 (73.2)
Metastasis
No	95 (32.6)
Yes	9 (3.1)
Unknow	187 (54.3)
Vital states (at follow-up)
Alive	251 (86.3)
Dead	40 (13.7)

Figure 1

Differentially expressed autophagy-related genes (ARGs) between papillary renal cell carcinoma (papillary RCC) and normal tissues. (A) The volcano plot for the 40 ARGs from the TCGA data. Red indicated high expression and green low expression. Black showed those genes showing no difference between papillary RCC and normal tissues. (B) Hierarchical clustering of differentially expressed ARGs expression levels. (C) The expression profiles of ARGs in papillary RCC and normal tissues. Red represented samples of tumor and blue samples of normal tissues.

Functional annotation of these 40 differentially expressed ARGs

In order to fully understand the biological attributes of these 40 differentially expressed ARGs, we conducted KEGG and GO analysis. Based on the results of DAVID, the top ten enriched GO terms were: selective autophagy, response to toxic response, response to starvation, response to oxygen level, regulation of release of cytochrome c from mitochondria, regulation of peptidase activity, regulation of autophagy, process utilizing autophagic mechanism, regulation of endopeptidase activity and positive regulation of tumor necrosis factor-mediated signaling pathway (). The top biological pathway enriched were: human cytomegalovirus infection, hepatitis B, Kaposi sarcoma-associated herpesvirus infection ().

Figure 2

The GO and KEGG analysis of these ARGs. The length of bars represented the number of genes and the color corresponds to the P value. (A,B,C) Enriched GO terms. (D,E,F) Enriched KEGG pathways. GO, gene ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; ARGs, autophagy-related genes.

Identification of prognostic ARGs and construction of the PI

For the purpose of revealing the relationship between these 40 ARGs and OS, 12 prognostic ARGs were identified by utilizing univariate Cox regression analyses ( and ). Then, multivariate Cox regression analysis was applied to increase the robustness and select the independent indicators for the OS and finally we got 4 ARGs (P4HB, BIRC5, NGR1 and PRKN) for the further analysis [P4HB: hazard ratio (HR) =2.686, 95% confidence interval (CI): 1.394–5.172, P=0.003; BIRC5: HR =2.290, 95% CI: 1.604–3.272, P<0.001; PRKN: HR =0.502, 95% CI: 0.255–0.986, P=0.046; NGR1: HR =1.652, 95% CI: 1.095–2.494, P=0.017] (). After the construction of PI, patients were separated into high- and low- risk (). Heap map was utilized to visualize the difference of gene expression profile in low- and high- risk patients in papillary RCC (). The results from Kaplan-Meier analysis indicated that high-risk patients had lower OS than low-risk patients in papillary RCC (P<0.001) (). Further, we applied this prognostic model to predict the vital status of patients and we found that the mortality of high scoring patients was higher than low scoring patients ().

Figure 3

Identification of prognostic ARGs and construction of the prognostic index. (A) Top 12 genes significantly associated with the survival time of patients in the training dataset. (B) The four-gene signature risk score distribution. (C) The heat-map of the four genes expression profiles in high-risk and low-risk groups. Red represented a higher expression and blue a lower expression. Red bar: low-risk group. Blue: high-risk group. (D) Kaplan-Meier plot showed that patients in the high-risk group had significantly shorter overall survival time than those in the low-risk group. (E) The mortality of high scoring patients was higher than low scoring patients. ARGs, autophagy-related genes.

Table 2

Results of multivariate Cox regression analysis for autophagy genes in papillary renal cell carcinoma

ID	HR	HR 95% lower limit	HR 95% upper limit	P value	Coef.
NRG1	1.652431834	1.094811831	2.494064175	0.016791909	0.502248042
PRKN	0.501990075	0.255444139	0.986493703	0.04556346	−0.689174931
BIRC5	2.290790468	1.60371499	3.27222792	5.21E-06	0.82889694
P4HB	2.685524806	1.39446036	5.171924345	0.003133078	0.987876167

Based on the Kaplan-Meier analysis, we found that patients with higher P4HB expression level showed an inferior OS and DFS (HR =3.1, P=0.00055, ; HR =1.9, P=0.029, ). Also, overexpression BIRC5 and NRG1 indicated worse OS and DFS (HR =2.4, P=0.0047, ; HR =4.4, P=2.1e-06, ; HR =2.9, P=0.00072, ; HR =1.8, P=0.041, ). On the contrary, patients with higher PARK2 expression had a better DFS and OS (HR =0.39, P=0.0029, ; HR =0.5, P=0.017, ). This might reveal that PARK2 was positive to better clinical outcomes whereas P4HB, BIRC5, NRG1 were negative.

Figure 4

The relationship between the four ARGs and papillary RCC patients disease-free survival (DFS) and overall survival (OS). Kaplan-Meier plots summarized the correlation between P4HB (A,B), BIRC5 (C,D), NRG1 (E,F) and PARK2 (G,H) expression and DFS and OS. RCC, renal cell carcinoma; ARGs, autophagy-related genes.

Relationship between 4 key genes and clinical parameters

We tried to illustrate the association between some clinical parameters and the gene expression level of these 4 key ARGs (P4HB, BIRC5, NGR1 and PRKN) (). Clinical parameters including T stage, N stage, M stage, age, gender and Stage were extracted in the current study. We found that the expression level of P4HB was correlated with age of patients. Patients <60 years old had higher P4HB expression level than patients >60 years old (P=0.033) (). As for BIRC5, its gene expression level was associated with T stage, M stage, N stage and Stage (P=0.003, P=0.004, P=0.008, P=1.862e-04, respectively) (). Patients with higher stage (T3–4, M1, N1–2, Stage III–IV) had higher BIRC5 expression than patients with lower stage (T1–2, N0, M0, Stage I–II). NRG1 was found to have relationship with T stage and Stage (P=0.002, ; P=0.002, ). Again, patients with higher NRG1 expression level showed a higher T stage and Stage (T3–4, Stage III–IV). Besides, no clinical parameters were illustrated to correlate with PRKN. Finally, we found that high scoring patients had higher T stage and Stage (T3–4, Stage III–IV) than low scoring patients (T1–2, Stage I–II) while other clinical parameters showed no significance (P=0.018, ; P=0.003, ).

Table 3

Correlation between four autophagy-related genes & risk score and clinical factors in papillary renal cell carcinoma (correlation value & P value)

Id	Age	Gender	Stage	T	M	N
NRG1	−0.102 (0.919)	1.359 (0.185)	−3.309 (0.002)	−3.381 (0.002)	−0.902 (0.387)	−1.31 (0.203)
PRKN	−0.727 (0.472)	0.84 (0.415)	−0.022 (0.982)	−0.343 (0.733)	0.762 (0.456)	0.36 (0.724)
BIRC5	1.336 (0.191)	−0.532 (0.599)	−4.102 (1.862e−04)	−3.156 (0.003)	−3.822 (0.004)	−2.959 (0.008)
P4HB	2.225 (0.033)	1.032 (0.317)	−1.778 (0.083)	−0.943 (0.351)	−1.172 (0.260)	−1.219 (0.236)
Risk Score	2.075 (0.052)	0.797 (0.439)	−3.262 (0.003)	−2.524 (0.018)	−1.774 (0.107)	−2.118 (0.054)

Figure 5

The relationship between 4 key genes and clinical parameters. (A) Patients <60 years old had higher P4HB expression level than patients >60 years old (P=0.033). (B,C,D,E) BIRC5 expression level was associated with T stage, M stage, N stage and Stage. (F,G) Patients with T3–4 and Stage III–IV had higher NRG1 expression than patients with T1–2 and Stage I–II. (H,I) High scoring patients had higher T stage and Stage than low scoring patients.

The PI was an independent prognostic factor for papillary RCC

In order to investigate whether the risk signature was independent from clinical parameters such as patients’ age, gender, Stage, lymphatic invasion, tumor size and metastasis, we conducted a univariate Cox and multivariate Cox regression analysis and explored the area under curve (AUC) (). According to the results, we found that metastasis stage and risk-score were significantly correlated with OS. Finally, we came to the conclusion that our risk signature was an independent prognostic factors for papillary RCC. To further provide a more quantitative and comprehensive prediction for the 3- and 5-year survival probability, we integrated out risk signature and some clinicopathological factors to construct a nomogram ().

Figure 6

Identification of the independence of risk signature and some other clinical parameters predictive of OS for papillary RCC. (A) Univariate Cox regression analysis. (B) Multivariate Cox regression analysis. (C) ROC curve analysis. The higher the area under the curve (AUC) value was, the more predictive accuracy the prognostic factors had. (D) Prognostic nomogram of survival probability for papillary RCC patients. OS, overall survival; RCC, renal cell carcinoma; ROC, receiver operating characteristic.

Discussion

While RCC was one of the most common malignancy in renal neoplasia, papillary RCC was a significant subtype of it. Thus, it was urgent to find molecular biomarkers that could aid in a better prognosis and OS for papillary RCC patients. It seemed that exploration of autophagy mechanism had shed light on this new perspective. To date, numerous studies have focused on the role that autophagy played in tumorigenesis and prognosis. Although autophagy has been a hotspot, it still remained a controversy whether autophagy suppressed or promoted tumorigenesis and metastasis. For example, several pathways such as p53 status, RAS family status, activation of JAK-STAT and PI3K signaling might all impose influence on the determination of autophagy dependence within cancer cells (11). Researches have shown that autophagy suppressed the initiation of tumor while promoting the progression of tumor (12). This implied that ARGs could act as a potential anticancer target in the therapeutic strategies and the prognosis of patients. However, while many researches of autophagy-related genes relied on the data retrieved from cell lines or animal models, our current study took use of high-throughput data of ARGs to investigate the prognosis and clinical outcomes of patients with papillary RCC. What’s more, we also explored the relationship between 4 key genes and clinical parameters, which could contribute to a better understanding of our prognostic signature.

Given the great development of high-throughput sequencing and the emergence of large-scare databases, we were able to select gene signatures. First, to investigate the differentially expressed ARGs, we compared the different expression level of autophagy between tumor tissues and adjacent tumor tissues and obtained 40 autophagy-related genes, including some popular oncotherapeutic target genes such as BIRC5 and NRG1. In order to figure out the functional attributes of these 40 ARGs, we then performed GO and KEGG analysis to figure out the biological attributes of these genes. Results demonstrated that these enriched genes were often associated with the following pathways: p53 signaling pathway, apoptosis, cellular senescence and so on. Relationship between p53 and autophagy has been observed. It seemed that autophagy could promote cancer by repressing p53, which was an important mechanism of tumorigenesis. What’s more, p53 also activated the transcription of ARGs and thus promoted autophagy (13). The functional association between autophagy and apoptosis was rather complex. Frequently, autophagy could serve as a trigger of apoptosis and aided in the process of the cellular suicide (14). For cellular senescence, however, it remained in controversy whether autophagy influenced senescence positively or negatively. Kang et al. claimed that selective autophagy functioned as an antisenescence mechanism, whereas general autophagy functioned as a pro-senescence mechanism (15).

Based on the results of univariate Cox regression analysis, 12 ARGs associated with OS were picked out. To further increase the robustness, multivariate Cox regression analysis was performed and finally 4 independent prognostic ARGs (P4HB, BIRC5, NGR1 and PRKN) were screen out to develop the PI. After constructing the 4-gene signature, patients were separated into high-risk and low-risk groups. Our results showed that high risk-score patients had an inferior survival probability and clinical outcomes.

BIRC5 (also named as survivin) was a well-known cancer therapeutic target which inhibited apoptosis and regulated mitosis both in embryonic cells during embryogenesis and in cancer cells during tumorigenesis (16). Studies have found that the mammalian target of rapamycin (mTOR) pathway was one of important pathways of autophagy, the inhibition of which was correlated with the promotion of autophagy (17). The post-transcription regulation of BIRC5 might be mediated by mTOR to increase the mRNA stability and translation of it (18). This might reveal the relationship between autophagy and BIRC5. Besides, the overexpression of BIRC5 was observed in almost all human malignancies, and the increased expression of BIRC5 was correlated with poor clinical outcomes, tumour recurrence and drug resistance in cancer patients (19). Our results showed that enriched gene expression of BIRC5 was associated with worse clinical index, which meant the coefficient of BIRC5 in the PI was positive and was in accordance with the accepted opinion of BIRC5 function in cancer. P4HB was a multifunctional protein which encoded the beta subunit of prolyl 4-hydroxylase and catalyzed the formation and rearrangement of disulfide bonds. Zhou et al. found that knockdown of P4HB could inhibit proliferation and promote apoptosis of human HT29 colon cancer cells via accumulation of reactive oxygen species and inhibition of STAT3 signaling (20). Several studies have reported that P4HB was significantly increased in several solid tumors including bladder cancer, brain and central nervous system (CNS) cancer, lung cancer, prostate cancer (21) and was correlated with poor prognosis. Xie et al. and Zhu et al. had shown that P4HB overexpressed in human clear cell RCC and was correlated with a poor prognosis, confirming that P4HB might be a novel biomarker (21,22). According to our results, P4HB showed a higher expression level in samples under 60 years old than those over 60 (P=0.03). Researches demonstrated that aging was accompanied by reduced autophagy level in various organisms, which was in accordance with our results. This phenomenon might result from an aging-dependent reduced expression of genes important for autophagosome-lysosome fusion, such as lysosome-associated membrane protein 2 (LAMP2A) (23). NRG1, a member of the NRG family, was a ligand for human epidermal growth factor (HER) 3 and HER4 which could activate cell signaling pathways to promote tumorigenesis and metastasis (24). Studies have indicated that NRG1 was overexpressed in various cancer such as breast cancer and gastric cancer (24,25). Han et al. found that increased NRG1 expression was related to advanced pathological stage, lymphnode metastasis and poor prognosis (26). Our research showed that expression of NRG was higher in T3–4 stage than in T1–2 stage, which was in accordance with Han et al. This indicated that NRG was associated with unfavourable clinicopathologic features in papillary RCC and was a negative prognostic factor. Research that targeted NRG1 pathway in lung cancer has succeeded, which meant NRG1 was an independent and effective biomarker for treatment and provided a novel strategy (27).

The advantage of our study was that instead of using cell lines and animal model, we turned to high-throughput data and large-scale databases to conduct a 4-gene prognostic signature statistically, which catered to the urgent need of an effective index for papillary RCC. Besides, our study could provide a better understanding of the role of autophagy in papillary RCC. However, our study did have some limitations. First of all, clinical parameters such as age, pathological stage were not integrated into our PI formula. We mainly focused on molecular markers while clinical perspective should also be explored. Second, our risk signature was composed of only ARGs and did not represent other potential gene transcription expression profiles correlated to OS in papillary RCC.

To date, the use of prognostic signature was proposed due to rapid development of RNAseq. Our study developed an independent 4-gene signature to reflect the clinical outcomes of papillary RCC. Further work will be done to deeply illustrate the role of autophagy in tumor and to improve our prognostic signature. In all, this novel therapy could provide a promising future for both diagnosis and prognosis for not only papillary RCC but also various cancer.

The article’s supplementary files as