Abstract

   Purpose: Autophagy is a major catabolic system by which eukaryotic
   cells undergo self-degradation of damaged, defective, or unwanted
   intracellular components. An abnormal autophagic level is implicated in
   the pathogenesis of multiple diseases, including cancers. The aim of
   this study is to explore the prognostic value of autophagy in bladder
   cancer (BC), which is a major cause of cancer-related death globally.

   Patients and methods: First, 27 differentially expressed
   autophagy-related genes (ARGs) were identified in BC patients based on
   The Cancer Genome Atlas (TCGA) database. Functional enrichment analyses
   hinted that autophagy may act in a tumor-suppressive role in the
   initiation of BC. Then, the Cox proportional hazard regression model
   were employed to identify three key prognostic ARGs (JUN, MYC, and
   ITGA3), which were related with overall survival (OS) significantly in
   BC. The three genes represented important clinical significance and
   prognostic value in BC. Then a prognostic index (PI) was constructed.

   Results: The PI was constructed based on the three genes, and
   significantly stratified BC patients into high- and low-risk groups in
   terms of OS (HR=1.610, 95% CI=1.200–2.160, P=0.002). PI remained as an
   independent prognostic factor in multivariate analyses (HR=2.355, 95%
   CI=1.483–3.739, P<0.001). When integrated with clinical characteristics
   of age and stage, an autophagy-clinical prognostic index (ACPI) was
   finally validated, which had improved performance in predicting OS of
   BC patients (HR=2.669, 95% CI=1.986–3.587, P<0.001). The ACPI was
   confirmed in datasets of [40]GSE13507 (HR=7.389, 95% CI=3.645–14.980,
   P<0.001) and [41]GSE31684 (HR=1.665, 95% CI=0.872–3.179, P=0.122).

   Conclusion: This study provides a potential prognostic signature for
   predicting prognosis of BC patients and molecular insights of autophagy
   in BC.

   Keywords: autophagy-related genes, prognostic index, bladder cancer,
   The Cancer Genome Atlas

Introduction

   Autophagy, also known as type II programmed cell death, is a major
   catabolic system by which eukaryotic cells undergo self-degradation of
   damaged, defective, or unwanted intracellular components.[42]^1^,[43]^2
   This is, in part, to quality control of intracellular organelles by
   continually renewing fresh, better-quality ones. Therefore, stability
   of cellular renovation, homeostasis, and maintaining physiological
   level are inseparable from autophagy. An abnormal autophagic level
   implicated in the pathogenesis of multiple diseases, including
   inflammation, neurodegenerative diseases, and tumors.[44]^3^–[45]^7
   However, the knowledge of autophagy-related mechanism in cancer is
   still rudimentary and inconclusive. Due to the complex function of
   autophagy in cancer, the further research on the relation of autophagy
   and tumors, underlying biological process, and then to apply this
   knowledge in well-designed therapeutic strategy could be valuable in
   the new route of cancer therapy. Even whether autophagy is a friend or
   a foe for cancers cannot draw reliable conclusions for
   now.[46]^8^–[47]^10

   Bladder cancer (BC) is a major cause of cancer-related death globally,
   causing 165,100 deaths per year.[48]^11 In the United States, thee were
   an estimated 79,030 newly-diagnosed cases in 2017 and 16,870 patients
   who succumbed to BC.[49]^12 Recently, several studies reported that
   autophagy could be an indispensable mechanism of the onset and
   progression of BC, which provided a new direction for the clinical
   management of BC.[50]^13^–[51]^17 Su et al[52]^18 observed increased
   autophagic proteins in high grade urothelial bladder carcinoma, which
   were regulated via AMPK activation and mTOR inhibition for tumor cells
   survival, and inhibition of autophagy led to cancer cell death. Some
   studies have also suggested that targeting autophagy could improve
   sensitivity to anti-bladder cancer chemotherapy agents.[53]^14^,[54]^15
   Thus, exploring the appropriate molecular biomarkers focused on
   autophagy has attractive value in estimating the deterioration of BC
   reliably, and may be an important means of fighting BC.

   Here we examined the correlation between expression profiles of
   autophagy-related genes (ARGs) and clinical outcome in 412 BC patients
   from The Cancer Genome Atlas (TCGA) and developed prognostic index (PI)
   as an independent index for overall survival (OS) prognosis based on
   ARGs. To leverage the complementary value of molecular and clinical
   characteristics, we integrated the PI with clinical factors to build a
   composite autophagy-clinical prognostic index (ACPI), which allowed us
   to improve the prognostic efficiency of BC patients. Further validation
   based on other databases evidently support our risk score model. These
   findings could also provide an effective multi-dimensional biomarker
   strategy that would be effective in monitoring autophagy and predicting
   the prognosis in BC patients.

Materials and methods

Data acquisition

   The Human Autophagy Database (HADb,
   [55]http://www.autophagy.lu/index.html) is an autophagy-dedicated
   database aiming to reserve human genes involved in autophagy. A variety
   of ARGs were obtained from the database. RNA-sequencing (RNA-seq) data
   of ARGs and the clinical information of the bladder urothelial cancer
   (BLCA) cohort were downloaded and extracted from the TCGA data portal.

Differentially expressed ARGs enrichment analysis

   EdgeR package in R statistical software was applied to estimate
   differentially expressed ARGs between BC and non-tumor samples. Genes
   exhibiting at least 2-fold changes corresponding to an adjusted P-value
   less than 0.05 were selected as the significantly differentially
   expressed ARGs. Then, we performed a series of gene functional
   enrichment analyses to find the major biological attributes of these
   genes, including gene ontology (GO) and Kyoto Encyclopedia of Genes and
   Genomes (KEGG). The Database for Annotation, Visualization, and
   Integrated Discovery (DAVID, [56]https://david.ncifcrf.gov/), a widely
   used functional annotation tool, was used to identify enriched GO and
   KEGG themes. To provide high-dimensional information, the GOplot
   package of R was performed to concentrate on the visualization of
   enrichment terms.

Construction of an individualized prognostic index based on ARGs

   ARGs expression profiles downloaded from TCGA were normalized by
   [log2(count+1)] transformed. Univariate Cox regression analyses were
   performed to select the ARGs whose expression profiles were
   significantly associated with BC patients’ OS. Subsequently, these
   survival-related genes were subjected to a multivariate Cox regression
   analysis to remove the genes that might not be an independent indicator
   in prognosis monitoring. Finally, several prognostic ARGs were obtained
   and the PI composed of these genes was developed. The formula of PI
   based on a linear combination of the relative expression level of genes
   multiplied regression coefficients, which represented the relative
   weight of genes in the multiple Cox analysis. BC patients were
   separated into high- and low-risk groups by the median PI value as the
   risk cutoff value. The survival curves were plotted by Kaplan–Meier
   (K–M) method, and differences in the survival rates between high- and
   low-risk groups were assessed using the log-rank test.

   To investigate if the autophagy-related PI could be an independent
   predictor of OS in the TCGA cohort of BC patients, the multivariate Cox
   regression analysis was conducted. The PI, age, gender, tumor subtype,
   pathological stages, and histological grades were used as covariates.
   Age, stage, and PI were coded as continuous variables. Specifically,
   stage was coded as I=1, II=2, III=3, and IV=4. The risk factors of
   gender, subtype, and histologic grade are male, non-Papillary, and high
   grade.

Statistical analysis

   All statistical analyses were conducted using SPSS 24.0 (Chicago, IL,
   USA) and R 3.3.1 ([57]https://www.r-project.org/). R, GraphPad Prism 5
   (San Diego, CA, USA), and OriginPro 2017 (Northampton, MA, USA) were
   performed to draw plots. Univariate Cox regression analyses were used
   to evaluate the association between expression profiles and OS. The
   Multivariate Cox proportional hazards regression model was used to
   construct the PI and ACPI model based on the factor correlated with
   survival. Receiver operating characteristic (ROC) curve and the
   corresponding area under the ROC curve (AUC) for each dataset to
   measure the prognostic value of ACPI were performed by the package of
   “survivalROC” in R. All statistical significance was defined as a
   P-value less than 0.05. Meta-analysis of the selected genes for ACPI
   was performed based on the expression data of BC and non-BC tissues
   from Oncomine ([58]https://www.oncomine.org/) with softwares of STATA
   (version 12.0) and Meta-DiSc (version 1.4).

Results

Differentially expressed ARGs

   Altogether RNA-seq and clinical data of 414 BC tissue samples and 19
   non-tumor samples were downloaded from TCGA. Among these patients, a
   total of 408 primary BC patients with gene expression data and clinical
   follow-up information was involved in the current study. Expression
   values of 234 ARGs were extracted. Considered as the criteria of a FDR
   <0.05 and |log2(Fold Change)|>1, we finally obtained nine up-regulated
   and 18 down-regulated ARGs ([59]Figure 1). Furthermore, scatter plots
   were visualized to display the expression pattern of the 27
   differentially expressed ARGs between BC and non-tumor tissues
   ([60]Figure 2). Scatter plots displayed expression patterns of 18
   down-regulated genes (FOS, JUN, HSPB8, CDKN1A, ITPR1, TP53INP2,
   PPP1R15A, DLC1, BAG3, MYC, GABARAPL1, BLC2, CCL2, PRKN, NAMPT, CXCR4,
   NRG2, and CX3CL1) and nine up-regulated genes (BIRC5, CDKN2A, BID,
   EVA1A, TP73, RGS19, EIF4EBP1, ITGB4, and ITGA3).

Figure 1.

   [61]Figure 1
   [62]Open in a new tab

   Differentially expressed autophagy-related genes (ARGs) between bladder
   cancer (BC) and normal bladder tissues. (A) The volcano plot for the
   234 ARGs from the TCGA data portal. Red indicates high expression and
   green low expression. Blue shows those genes showed no difference
   between BC and normal bladder tissues. (B) Hierarchical clustering of
   differentially expressed ARGs expression levels.

   Abbreviations: FDR, false discovery rate; TCGA, The Cancer Genome
   Atlas.

Figure 2.

   [63]Figure 2
   [64]Open in a new tab

   The expression patterns of 27 autophagy-related genes (ARGs) in bladder
   cancer types and paired non-tumor samples. Each red dot represents a
   distinct tumor sample and blue a non-tumor sample. The red bar above
   the gene name shows a significantly high expression and the blue bar a
   low expression.

Functional annotation of the differentially expressed ARGs

   Functional enrichment analysis of the 27 differentially expressed ARGs
   offered that the biological understanding of these genes. The GO terms
   function and KEGG pathway enrichment of these genes were summarized in
   [65]Table 1. According to the results of DAVID, we found that the top
   enriched GO terms for biological processes were: response to drug,
   response to gamma radiation, and apoptotic process; and for cellular
   components were: cytosol, protein complex, and mitochondrion. On the
   basis of molecular function, genes were mostly enriched in terms of
   transcription factor binding, ubiquitin protein ligase binding, and
   protein heterodimerization activity. The overview schematic of the
   analysis results is displayed in [66]Figure 3. Besides, in the KEGG
   pathway enrichment analysis for the differentially expressed ARGs,
   these genes were shown to be notably associated with Pathways in
   cancer, colorectal cancer, Hepatitis B, ErbB signaling pathway, p53
   signaling pathway, and so on. As shown in [67]Figure 4A, the Z-score of
   enriched pathways less than zero indicated that most of the cancer
   pathways were more likely to be decreased. The heatmap of the
   relationship between ARGs and pathways was also displayed ([68]Figure
   4B).

Table 1.

   GO and KEGG analysis of differentially expressed autophagy-related
   genes
   Category ID Term P-value Genes
   Biological Process GO:0042493 Response to drug 4.44E−06 FOS, CDKN1A,
   JUN, BCL2, ITGA3, MYC, TP73
   Biological Process GO:0010332 Response to gamma radiation 1.27E−05
   CCL2, BCL2, MYC, TP73
   Biological Process GO:0006915 Apoptotic process 1.36E−05 DLC1, EVA1A,
   CDKN2A, CXCR4, BCL2, BIRC5, PPP1R15A, TP73
   Biological Process GO:0007346 Regulation of mitotic cell cycle 2.78E−05
   CDKN1A, BIRC5, MYC, TP73
   Biological Process GO:0043524 Negative regulation of neuron apoptotic
   process 4.06E−05 CCL2, JUN, BCL2, BIRC5, TP73
   Biological Process GO:0007050 Cell cycle arrest 5.25E−05 CDKN1A,
   CDKN2A, PPP1R15A, MYC, TP73
   Biological Process GO:0006974 Cellular response to DNA damage stimulus
   2.36E−04 CDKN1A, BCL2, PPP1R15A, MYC, TP73
   Biological Process GO:0001836 Release of cytochrome c from mitochondria
   5.28E−04 BID, JUN, BCL2
   Biological Process GO:0045893 Positive regulation of transcription,
   DNA-templated 8.50E−04 FOS, CDKN2A, JUN, MYC, TP53INP2, TP73
   Biological Process GO:1900740 Positive regulation of protein insertion
   into mitochondrial membrane involved in apoptotic signaling pathway
   9.02E−04 BID, BCL2, TP73
   Cellular Component GO:0005829 Cytosol 3.94E−06 BID, DLC1, NAMPT,
   GABARAPL1, BIRC5, TP73, FOS, CDKN1A, EIF4EBP1, CDKN2A, JUN, BCL2, BAG3,
   MYC, PPP1R15A, TP53INP2
   Cellular Component GO:0043234 Protein complex 0.002235 EIF4EBP1,
   CDKN1A, CDKN2A, MYC, ITPR1
   Cellular Component GO:0005739 Mitochondrion 0.007959 BID, GABARAPL1,
   CDKN2A, BCL2, PPP1R15A, MYC, TP73
   Cellular Component GO:0005654 Nucleoplasm 0.008942 NAMPT, FOS,
   EIF4EBP1, CDKN1A, CDKN2A, JUN, HSPB8, BIRC5, MYC, TP73
   Cellular Component GO:0005741 Mitochondrial outer membrane 0.017612
   BID, BCL2, PPP1R15A
   Cellular Component GO:0030054 Cell junction 0.024203 NAMPT, CXCR4,
   ITGB4, TP73
   Cellular Component GO:0005783 Endoplasmic reticulum 0.025005 FOS,
   GABARAPL1, BCL2, PPP1R15A, ITPR1
   Cellular Component GO:0005667 Transcription factor complex 0.028532
   FOS, JUN, TP73
   Cellular Component GO:0005794 Golgi apparatus 0.028578 GABARAPL1,
   HSPB8, RGS19, PPP1R15A
   Cellular Component GO:0008305 Integrin complex 0.036412 ITGB4, ITGA3
   Cellular Component GO:0005515 Protein binding 3.97E−05 BID, DLC1,
   NAMPT, GABARAPL1, ITGB4, RGS19, BIRC5, ITGA3, CX3CL1, TP73, ITPR1, FOS,
   CDKN1A, EIF4EBP1, CDKN2A, CXCR4, JUN, HSPB8, BCL2, BAG3, MYC, PPP1R15A,
   TP53INP2
   Molecular Function GO:0008134 Transcription factor binding 4.25E−05
   FOS, CDKN2A, JUN, BCL2, MYC, TP73
   Molecular Function GO:0031625 Ubiquitin protein ligase binding 6.65E−04
   BID, GABARAPL1, CDKN1A, CXCR4, BCL2
   Molecular Function GO:0046982 Protein heterodimerization activity
   0.003898 FOS, JUN, BCL2, BIRC5, ITGA3
   Molecular Function GO:0001077 Transcriptional activator activity, RNA
   polymerase II core promoter proximal region sequence-specific binding
   0.004395 FOS, JUN, MYC, TP73
   Molecular Function GO:0043565 Sequence-specific DNA binding 0.005716
   FOS, JUN, BCL2, MYC, TP73
   Molecular Function GO:0000978 RNA polymerase II core promoter proximal
   region sequence-specific DNA binding 0.013453 FOS, JUN, MYC, TP73
   Molecular Function GO:0004861 Cyclin-dependent protein serine/threonine
   kinase inhibitor activity 0.015533 CDKN1A, CDKN2A
   Molecular Function GO:0019901 Protein kinase binding 0.015684 CDKN2A,
   RGS19, PPP1R15A, TP73
   Molecular Function GO:0042802 Identical protein binding 0.020127 JUN,
   HSPB8, BCL2, BIRC5, TP73
   KEGG PATHWAY hsa05200 Pathways in cancer 3.16E−07 BID, FOS, CDKN1A,
   CDKN2A, CXCR4, JUN, BCL2, BIRC5, ITGA3, MYC
   KEGG PATHWAY hsa05210 Colorectal cancer 2.06E−05 FOS, JUN, BCL2, BIRC5,
   MYC
   KEGG PATHWAY hsa05161 Hepatitis B 3.48E−05 FOS, CDKN1A, JUN, BCL2,
   BIRC5, MYC
   KEGG PATHWAY hsa04012 ErbB signaling pathway 7.87E−05 EIF4EBP1, CDKN1A,
   JUN, NRG2, MYC
   KEGG PATHWAY hsa04115 p53 signaling pathway 7.56E−04 BID, CDKN1A,
   CDKN2A, TP73
   KEGG PATHWAY hsa04151 PI3K-Akt signaling pathway 0.001961 EIF4EBP1,
   CDKN1A, BCL2, ITGB4, ITGA3, MYC
   KEGG PATHWAY hsa04668 TNF signaling pathway 0.002843 FOS, CCL2, JUN,
   CX3CL1
   KEGG PATHWAY hsa05166 HTLV-I infection 0.004594 FOS, CDKN1A, CDKN2A,
   JUN, MYC
   KEGG PATHWAY hsa05219 Bladder cancer 0.00551 CDKN1A, CDKN2A, MYC
   KEGG PATHWAY hsa04921 Oxytocin signaling pathway 0.008685 FOS, CDKN1A,
   JUN, ITPR1
   [69]Open in a new tab

   Abbreviations: GO, gene ontology; KEGG, Kyoto Encyclopedia of Genes and
   Genomes.

Figure 3.

   [70]Figure 3
   [71]Open in a new tab

   The bubble plot of enriched gene ontology (GO) terms. The z-score is
   assigned to the x-axis, and the negative logarithm of the P-value to
   the y-axis, as in the barplot (the higher the more significant). The
   size of the displayed circles is proportional to the number of genes
   assigned to the term. Greed circles correspond to the biological
   process, red indicates the cellular component, and blue shows the
   molecular function category.

Figure 4.

   [72]Figure 4
   [73]Open in a new tab

   Kyoto Encyclopedia of Genes and Genomes analysis of differentially
   expressed autophagy-related genes (ARGs). (A) The outer circle shows a
   scatter plot for each term of the logFC of the assigned genes. Red
   circles display up-regulation, and blue ones down-regulation. (B) The
   heatmap of the relationship between ARGs and pathways. The color of
   each block depends on the logFC values.

Identification of prognostic ARGs

   The relationships between the expression profiles of 27 differentially
   expressed ARGs and OS were assessed based on the data obtained from
   TCGA, resulting in four prognosis-related ARGs. In order to improve the
   robustness, four prognosis-related ARGs (JUN, MYC, ITGA3, and NAMPT)
   were selected for further multivariate Cox regression model by SPSS
   24.0 ([74]Table 2). However, the gene of NAMPT showed no significant
   prognostic value with P>0.05. Finally, three genes including JUN, MYC,
   and ITGA3 were identified to develop the PI model ([75]Table 2). The
   results from K-M analysis indicated that the up-regulation of JUN was
   strongly correlated with the inferior OS of BC patients (HR=1.925, 95%
   CI=1.325–2.798, P<0.001; [76]Figure 5A). Also, MYC overexpression lead
   to worse OS (HR=1.931, 95% CI=1.426–2.614, P<0.001; [77]Figure 5C). On
   the contrary, up-regulated ITGA3 indicated BC patients has a longer
   survival time (HR=0.659, 95% CI=0.487–0.893, P=0.007; [78]Figure 5E).
   According to the median value of the three genes to group, the same
   trend was obtained ([79]Figures 5B, [80]D, and [81]F).

Table 2.

   Expression and Cox regression analysis data of the prognosis-related
   ARGs in bladder cancer by TCGA
   Gene Expression Univariate Cox Multivariate Cox
   mean±SD P-value HR z P-value HR z P-value
   JUN Tumor 13.002773±1.164031 <0.0001 1.198762 2.651598 0.008011 1.1786
   2.38 0.017
   Non-tumor 15.027896±1.185059
   MYC Tumor 11.541411±1.554496 <0.0001 1.14059 2.606545 0.009146 1.1682
   2.97 0.003
   Non-tumor 13.287664±1.265550
   ITGA3 Tumor 13.822128±1.515658 <0.0001 0.907077 −2.0139 0.04402 −0.1505
   −2.97 0.003
   Non-tumor 12.975430±1.091134
   NAMPT Tumor 12.330813±1.009040 <0.0001 1.172207 2.001379 0.045352 0.340
   Non-tumor 13.308780±1.459465
   FOS Tumor 13.343147±1.608937 <0.0001 1.063323 1.267733 0.204893
   Non-tumor 16.672991±1.404332
   HSPB8 Tumor 10.106659±2.154647 <0.0001 1.068919 1.911278 0.055969
   Non-tumor 13.749171±1.702493
   CDKN1A Tumor 13.124746±1.128535 <0.0001 1.01594 0.22579 0.821364
   Non-tumor 14.389922±1.347070
   ITPR1 Tumor 9.5318665±1.169013 <0.0001 0.997873 −0.03415 0.972754
   Non-tumor 11.929274±1.585768
   TP53INP2 Tumor 10.462387±0.976958 <0.0001 1.119264 1.518161 0.128974
   Non-tumor 12.653637±1.109542
   BIRC5 Tumor 11.387251±1.094559 <0.0001 1.06901 0.897696 0.369347
   Non-tumor 7.9468232±2.565110
   PPP1R15A Tumor 12.441886±0.891543 <0.0001 0.940067 −0.71379 0.475356
   Non-tumor 13.738993±1.001618
   DLC1 Tumor 9.5131164±1.257221 <0.0001 1.021996 0.347978 0.727857
   Non-tumor 11.242183±0.995269
   BAG3 Tumor 11.719572±0.731981 <0.0001 1.192655 1.702932 0.088581
   Non-tumor 12.704785±1.108736
   GABARAPL1 Tumor 11.058632±0.792235 <0.0001 1.051224 0.525272 0.599394
   Non-tumor 12.255963±0.614109
   BCL2 Tumor 8.3027604±1.178505 <0.0001 1.051513 0.768008 0.442482
   Non-tumor 10.256868±0.425816
   CDKN2A Tumor 8.9951212±3.453185 <0.0001 0.990449 −0.44257 0.658076
   Non-tumor 6.6235955±1.685343
   CCL2 Tumor 9.0438964±1.821760 <0.0001 1.053817 1.273601 0.202805
   Non-tumor 11.704824±2.040782
   PRKN Tumor 5.2978598±1.521700 <0.0001 1.009973 0.201647 0.840193
   Non-tumor 8.1483213±1.213423
   BID Tumor 11.340986±0.737984 <0.0001 0.839371 −1.70367 0.088442
   Non-tumor 10.105586±0.651137
   EVA1A Tumor 6.5437907±2.147278 <0.0001 1.057202 1.528728 0.126332
   Non-tumor 3.8921911±1.267249
   CXCR4 Tumor 10.259886±1.380760 <0.0001 1.013263 0.261675 0.793572
   Non-tumor 11.447608±1.958345
   NRG2 Tumor 4.7761532±2.145474 <0.0001 0.985697 −0.43265 0.665266
   Non-tumor 7.9643803±1.000720
   TP73 Tumor 7.9164301±2.009466 <0.0001 0.984879 −0.3943 0.69336
   Non-tumor 5.9146373±2.363484
   RGS19 Tumor 9.9259644±0.729683 <0.0001 1.143028 1.313883 0.188886
   Non-tumor 8.6949825±0.978522
   EIF4EBP1 Tumor 11.427944±1.081142 <0.0001 1.046433 0.666103 0.505345
   Non-tumor 9.8783883±0.927159
   CX3CL1 Tumor 9.7520038±1.782426 <0.0001 1.034169 0.77976 0.435532
   Non-tumor 12.084799±0.955857
   ITGB4 Tumor 14.253622±1.551662 <0.0001 0.936675 −1.33788 0.180934
   Non-tumor 13.188612±1.541228
   [82]Open in a new tab

   Abbreviations: ARG, autophagy-related genes; TCGA, The Cancer Genome
   Atlas.

Figure 5.

   [83]Figure 5
   [84]Open in a new tab

   The correlation between three genes included in prognostic signature
   and bladder cancer patients’ survival. Kaplan–Meier plots summarize
   results from analysis of correlation between (A) JUN expression level
   and patient survival, using best separation, (B) JUN expression level
   and patient survival, using median separation, (C) MYC expression level
   and patient survival, using best separation, (D) MYC expression level
   and patient survival, using median separation, (E) ITGA3 expression
   level and patient survival, using best separation, (F) ITGA3 expression
   level and patient survival, using median separation.

   The correlations between expression level of the three genes and
   clinicopathological parameters in BC are summarized in [85]Table 3. We
   observed significant correlations between JUN overexpression and tumor
   subtype of non-Papillary (P<0.001), high histological grade (P=0.033),
   advanced pathologic stage (P<0.001), and advanced pathological T stage
   (P=0.002). Elevated MYC was closely linked with tumor subtype of
   non-Papillary (P=0.001) and high histological grade (P=0.005). High
   expression of ITGA3 occurred in low histological grade (P<0.001), early
   pathological stage (P=0.033), early pathological T stage (P=0.008), and
   no lymph node metastasis (P=0.016).

Table 3.

   Relationships between three key prognostic genes expression and
   clinical parameters in bladder cancer by TCGA
   Parameters N JUN expression value MYC expression value ITGA3 expression
   value
   M±SD t value P-value M±SD t value P-value M±SD t value P-value
   Tissues Tumor 414 13.0028±1.1640 −7.410 <0.001 11.5414±1.5545 −4.822
   <0.001 13.8221±1.5157 2.405 0.017
   Non-tumor 19 15.0279±1.1851 13.2877±1.2655 12.9754±1.0911
   Age ≥60 321 13.0359±1.1327 1.076 0.150 11.6119±1.5304 1.571 0.117
   13.7787±1.5272 −1.665 0.097
   <60 87 12.8852±1.2530 11.3177±1.6171 14.0830±1.4561
   Gender Male 301 12.9868±1.1674 −0.495 0.621 11.4625±1.5666 −1.898 0.058
   13.7850±1.5661 −1.311 0.191
   Female 107 13.0515±1.1410 11.7930±1.4901 14.0084±1.3574
   Subtype Non-papillary 271 13.1979±1.0811 4.903 <0.001 11.7343±1.5708
   3.405 0.001 13.7762±1.5897 −1.492 0.136
   Papillary 132 12.6091±1.2284 11.1787±1.4654 14.0137±1.2921
   Tumor status With tumor 151 13.0248±1.0119 0.459 0.647 11.6697±1.5948
   1.351 0.177 13.7276±1.5787 1.763 0.079
   Tumor free 201 12.9720±1.2459 11.4571±1.4702 13.9996±1.4314
   Histological grade High grade 384 13.0359±1.1529 2.143 0.033
   11.5759±1.5844 3.002 0.005 13.7990±1.5338 −6.996 <0.001
   Low grade 21 12.4819±1.1653 11.0637±0.6886 14.8231±0.5669
   Pathologic stage Ⅲ～Ⅳ 274 13.1826±1.0448 4.206 <0.001 11.6576±1.4617
   1.944 0.053 13.7353±1.5386 −2.136 0.033
   I~II 132 12.6367±1.3030 11.3387±1.7148 14.0770±1.4496
   T stage T3–T4 252 13.1958±1.0249 3.064 0.002 11.6811±1.4423 1.643 0.101
   13.7290±1.5227 −2.684 0.008
   T1–T2 123 12.8322±1.1819 11.4088±1.6306 14.1726±1.4599
   N stage N1–3 129 13.1086±0.9329 1.265 0.207 11.4309±1.4607 −1.517 0.130
   13.6208±1.4944 −2.427 0.016
   N0 237 12.9643±1.2176 11.6798±1.5195 14.0205±1.5112
   M stage M1 11 12.8009±1.0016 −0.191 0.849 11.9098±1.8601 1.245 0.215
   13.5755±1.3148 −0.792 0.429
   M0 196 12.8728±1.2223 11.3199±1.5104 13.9395±1.4919
   [86]Open in a new tab

   Abbreviations: M±SD, mean±standard deviation; TCGA, The Cancer Genome
   Atlas.

Construction and definition of the PI

   The formula of PI is as follows: PI=(0.1643 × expression value of
   JUN)+(0.1555 × expression value of MYC)+(−0.1505 × expression value of
   ITGA3). It is noticed that the coefficient of ITGA3 is negative,
   indicating that the expression of JUN and MYC were negatively related
   with the survival time of BC patients, while the JUN was positively
   related with OS. Based on the median expression value of PI, the BC
   patients were stratified into high- and low-risk groups. We also
   calculated the expression levels of the three prognostic genes between
   high- and low-risk groups. Remarkably higher expression was noted for
   JUN and MYC in the high-risk groups, while lower expression was
   observed for ITGA3 in the high-risk groups ([87]Figure 6). These
   findings also hint that JUN and MYC were risk factors, while ITGA3 was
   a protective factor for the progression of BC patients.

Figure 6.

   [88]Figure 6
   [89]Open in a new tab

   Different expression of the three key genes between the high risk group
   and low risk group.

   In the meantime, the relationships between clinicopathological
   parameters and PI were also investigated. The results of independent
   sample t-tests showed that the PI values were higher in elder than in
   younger patients (P=0.009; [90]Figure 7A), higher in non-papillary than
   in papillary bladder cancer (P<0.001; [91]Figure 7C), higher in TIII–IV
   than in TI–II (P<0.001; [92]Figure 7D), higher in histological stage
   III–IV than in I–II (P<0.001; [93]Figure 7G), and higher in high grade
   than in low grade (P<0.001; [94]Figure 7H). No difference of PI value
   was observed between male and female (P=0.494; [95]Figure 7B), N1–3
   stage and N0 stage (P=0.250; [96]Figure 7E), or M1 and M0 stage
   (P=0.254; [97]Figure 7F).

Figure 7.

   [98]Figure 7
   [99]Open in a new tab

   The clinicopathological significance of prognostic index (PI) in
   bladder cancer. PI value in different (A) ages, (B) genders, (C) tumor
   subtypes, (D) pathological T stages, (E) pathological N stages, (F)
   pathological M stages, (G) pathological stages, (H) histological
   grades.

   To identify the performance of PI in predicting the clinical outcome of
   BC patients, the K-M plots were plotted to analyze the different
   survival time between the high- and low-risk groups. The results of K-M
   analysis indicated that the median OS for the high-expression group was
   734 days; the median OS for the low-expression group was 1,423 days.
   Patients in the high-risk group suffered significantly worse survival
   than those in the low-risk group (HR=1.610, 95% CI=1.200–2.160,
   P=0.002, [100]Figure 8A). [101]Figures 8B–[102]F show the PI
   distribution of patients in the training dataset, the number of
   patients in different risk groups, the OS of patients in the TCGA
   dataset, the number of censor patients, and the heatmap of the three
   genes expression profiles in the TCGA dataset. Furthermore, PI remained
   as an independent prognostic indicator for BC patients in multivariate
   analyses, after adjusting for clinicopathological features such as age,
   gender, tumor subtype, pathologic stage, and histological grade
   (HR=2.355, 95% CI=1.483–3.739, P<0.001, [103]Table 4).

Figure 8.

   [104]Figure 8
   [105]Open in a new tab

   Autophagy-related prognostic index (PI) of bladder cancer patients. (A)
   Kaplan–Meier plot represents that patients in the high-risk group had
   significantly shorter overall survival time than those in the low-risk
   group. (B) The PI distribution of patients in the training dataset. (C)
   The number of patients in different risk groups. (D) The overall
   survival of patients in the TCGA dataset. (E) The number of censor
   patients. (F) The heatmap of the three key genes expression profiles in
   the TCGA dataset.

   Abbreviations: TCGA, The Cancer Genome Atlas.

Table 4.

   Univariate and multivariate analyses of OS in bladder cancer patients
   of TCGA
   Variables Univariate analysis Multivariate analysis
   Hazard ratio (95%CI) P-value Hazard ratio (95% CI) P-value
   Age 1.033 (1.017–1.049) <0.001 1.029 (1.014–1.045) <0.001
   Gender 0.872 (0.631–1.203) 0.404 0.818 (0.588–1.139) 0.235
   Subtype 1.458 (1.030–2.065) 0.033 1.084 (0.756–1.553) 0.661
   Pathologic stage 1.707 (1.412–2.065) <0.001 1.617 (1.321–1.978) <0.001
   Histologic grade 2.968 (0.734–11.995) 0.127 0.931 (0.221–3.918) 0.922
   Prognostic index 2.717 (1.764–4.184) <0.001 2.355 (1.483–3.739) <0.001
   [106]Open in a new tab

   Notes: Age, stage, and prognostic index were coded as continuous
   variables. Specifically, stage was coded as I=1, II=2, III=3. IV=4. The
   risk factors of gender, subtype,and histologic grade are male,
   non-papillary, and high grade.

   Abbreviations: OS, overall survival; TCGA, The Cancer Genome Atlas.

Integrated prognostic signature by combining the PI with clinical parameters

   Based on the multivariate Cox regression analysis with TCGA dataset,
   age, stage, and PI were suggested as independent prognostic factors
   with complementary value. To further improve accuracy of PI in
   predicting OS of BC patients, we integrated age, pathological stage,
   and PI to derive an ACPI as (0.028 × age)+(0.467 × stage)+(0.834 × PI
   score). Similarly, patients were divided into high- and low-risk groups
   based on the median value of ACPI. As expected, ACPI stratified BC
   patients into two groups with a significantly different prognosis
   (HR=2.669, 95% CI=1.986–3.587, P<0.001; [107]Figure 9A). To evaluate
   how well the ACPI predicts the prognoses of BC patients, the
   time-dependent ROC curve analysis was carried out. The AUC for the ACPI
   was 0.689 ([108]Figure 9B), demonstrating the competitive performance
   of the ACPI for survival prediction in the TCGA dataset. The prognostic
   value of ACPI was also validated by [109]GSE13507 and [110]GSE31684.
   Consistent with the findings based on the TCGA dataset, patients in the
   high-risk group had significantly shorter overall survival than those
   in the low-risk group based on [111]GSE13507 (HR=7.389, 95%
   CI=3.645–14.980, P<0.001; [112]Figure 9C), and the AUC for the ACPI was
   0.864 ([113]Figure 9D). A similar trend was observed in [114]GSE31684
   (HR=1.665, 95% CI=0.872–3.179, P=0.122; [115]Figure 9E), and the AUC
   for the ACPI was 0.624 ([116]Figure 9F).

Figure 9.

   [117]Figure 9
   [118]Open in a new tab

   The prognostic value of autophagy-clinical prognostic index (ACPI) of
   bladder cancer patients. (A) Kaplan–Meier (K–M) survival curve showing
   overall survival (OS) outcomes according to relative high-risk and
   low-risk patients based on The Cancer Genome Atlas (TCGA) database. (B)
   Time-dependent ROC curve analysis for survival prediction by the ACPI
   based on the TCGA database. (C) K–M survival curve showing OS outcomes
   according to relative high-risk and low-risk patients based on the
   [119]GSE13507 dataset. (D) Time-dependent ROC curve analysis for
   survival prediction by the ACPI based on the [120]GSE13507 dataset. (E)
   K–M survival curve showing OS outcomes according to relative high-risk
   and low-risk patients based on the [121]GSE31684 dataset. (F)
   Time-dependent ROC curve analysis for survival prediction by the ACPI
   based on the [122]GSE31684 dataset.

Meta-analysis

   A total of 19 eligible studies were involved, including Blaveri Bladder
   2, Modlich Bladder, Sanchez Carbayo Bladder 2, TCGA, [123]GSE3167,
   [124]GSE13507, [125]GSE76211, [126]GSE2109, [127]GSE7476,
   [128]GSE30522, [129]GSE31189, [130]GSE37815, [131]GSE52519,
   [132]GSE65635, [133]GSE37817, [134]GSE100926, [135]GSE24152,
   [136]GSE19915 ([137]GPL3883 and [138]GPL5186), and [139]GSE40355. The
   results of meta-analysis and the diagnostic tests of meta-analysis were
   also updated. The expression of JUN and MYC in BC tissues were lower
   than that in non-BC tissues (I^2[JUN]=91.6%, P[JUN]<0.001;
   I^2[MYC]=90.3%, P[MYC]<0.001) ([140]Figures 10A and [141]C), while the
   expression of ITGA3 was opposite (I^2[ITGA3]=92.7%, P[ITGA3]<0.001)
   ([142]Figure 10E), the same expression tendence with that in TCGA. In
   addition, the diagnostic tests of meta-analysis showed that the AUC of
   the sROC of JUN, MYC, and ITGA3 were 0.91, 0.87, and 0.74, respectively
   ([143]Figures 10B, [144]D, and [145]F). Among the 19 studies involved
   in meta-analysis, only two microarrays ([146]GSE37137 and
   [147]GSE35824) and TCGA contained gene expression data from non-muscle
   invasive bladder cancer (NMIBC) and muscle invasive bladder cancer
   (MIBC) tissues. Then we performed a meta-analysis to evaluate the
   expression of JUN, MYC, and ITGA3 between NMIBC and MIBC tissues. The
   expression of JUN and ITGA3 in MIBC tissues were higher than that in
   NMIBC tissues (I^2[JUN]=72.1%, P[JUN]=0.028; I^2[ITGA3]=92.8%,
   P[ITGA3]<0.001), while the expression of MYC was opposite
   (I^2[MYC]=67.2%, P[ITGA3]=0.047) (data not shown). The heterogeneity of
   meta-analysis was significant due to the small size of cases involved
   in the microarrays and TCGA.

Figure 10.

   [148]Figure 10
   [149]Open in a new tab

   Meta-analysis. (A) Forest plot of JUN expression in bladder cancer with
   six datasets. (B) sROC curve for JUN expression in bladder cancer with
   six datasets. (C) Forest plot of MYC expression in bladder cancer with
   six datasets. (D) sROC curve for MYC expression in bladder cancer with
   six datasets. (E) Forest plot of ITGA3 expression in bladder cancer
   with six datasets. (F) SROC curve for ITGA3 expression in bladder
   cancer with six datasets.

Discussion

   BC is a major lethal malignancy worldwide. The stalled advance in
   molecular targeted therapy and no effective molecular biomarkers for BC
   prognosis monitoring warrants a better understanding of the molecular
   mechanisms that underlie this condition.[150]^19^,[151]^20 Exploration
   of autophagy mechanism opens new perspectives for BC.[152]^21^–[153]^25
   However, most research only focused on autophagy via studying a signal
   gene. To capture the genes necessary for BC from the perspective of
   autophagy, we screened ARGs and identified three key prognostic ARGs,
   all of which may offer additional potential therapeutic targets. We
   further leveraged the complementary value of molecular and clinical
   characteristics and showed that combining both could provide a more
   accurate estimation of overall survival in BC. This integrated study of
   multiple databases contributed to our novel understanding of BC biology
   and delineated potential therapeutic intervention possibilities.

   Given great advances in high-throughput sequencing recently, several
   large-scale databases emerged, such as TCGA and GEO, which have
   provided effective measures for selecting gene signatures. In the
   current study, we deeply mined the expression profiles of ARGs from
   TCGA and aimed to search molecular biomarkers for detecting the
   prognosis of BC patients. We first screened 27 differentially expressed
   ARGs between BC and non-tumor tissues. Considering these genes may be
   depth involved in the initiation of BC, we performed GO and KEGG
   analysis of these genes. Interestingly, functional analysis revealed
   that the most significant KEGG pathway (pathways in cancer) of these
   enriched genes was decreased. Based on the results, we hypothesized
   that autophagy may act as the tumor suppressor in the process of tumor
   initiation. Autophagy caused great concern; of particular interest was
   its multi-faceted character in cancers. Initially, the
   tumor-suppressive role of autophagy in cancers was proposed for
   autophagy inhibited by activation of mutations in oncogenes or
   inactivation of tumor suppressor genes.[154]^26 Furthermore, systemic
   mosaic deletion of autophagy genes in the setting of certain mouse
   models can result in the initiation of neoplasia.[155]^27
   Interestingly, autophagy turns to the guardian of malignant tumor cells
   after tumors are established.[156]^28 However, the role change of
   autophagy is not immutable and varies in different
   tumors.[157]^29^,[158]^30

   The result of univariate survival analysis revealed that four ARGs were
   associated with OS in the TCGA database. Further multivariate survival
   analysis helped us determine three key prognostic ARGs (JUN, MYC,
   ITGA3) to develop the PI, which could be an independent prognostic
   indicator for BC patients. JUN encodes c-Jun, which is the first
   discovered oncogenic transcription factor,[159]^31 involving diverse
   cellular processes, such as cell cycle progression,[160]^32
   anti-apoptotic, and tumorigenesis. Previous studies have suggested that
   up-regulation of c-Jun proteins was predictive of inferior OS for BC
   patients.[161]^33 However, well-informed insights of the functional
   mechanism of JUN in BC only has little coverage. The MYC protein is a
   multifunctional, nuclear phosphoprotein, and shows its evil face in the
   progress of a variety of tumors, including BC.[162]^34^–[163]^37
   Massari et al[164]^38 also found that c-Myc could exert excellent
   ability in stratifying patients with muscle invasive bladder urothelial
   carcinoma into high-risk and low-risk groups significantly for
   survival. In addition, several studies found that c-Myc knockdown could
   inhibit proliferation, migration, and invasion of bladder cancer
   cells.[165]^37 ITGA3 belongs to a family of the integrins, which
   triggers cell survival, proliferation, or migration events.[166]^39 The
   present study demonstrated that JUN and MYC overexpression were
   significantly associated with advanced pathological stage and high
   grade. Additionally, up-regulation of JUN and MYC indicated inferior
   OS. The opposite pattern was observed in the relationships between
   ITGA3 and clinical significance. Hence, we speculated that JUN and MYC
   may function as major driving forces of tumor progression, while ITGA3
   exerted its tumor suppressor role.

   To date, some prognostic signature of cancers based on expression
   profiles were proposed by the aid of advances in a large-scale public
   database. For example, Bao et al[167]^40 analyzed the RNA-Seq data of
   234 BC patients from TCGA and managed to obtain a four-lncRNA
   signature, which exerted a prognosis predicting value. Zhong et
   al[168]^41 also proposed a prognostic signature with six genes as a
   potential survival prediction marker for ER-positive breast cancer
   patients. However, these studies only focused on molecular biomarkers
   and overlooked the traditional clinical parameters. We attached much
   weight on molecular mechanisms and clinical perspective at once. Thus,
   the prognostic signature is promising to be converted into clinical
   application. However, a limitation of this study is its retrospective
   nature. Due to the lack of enough cases, we failed to evaluate the
   expression of JUN, MYC, and ITGA3 between NMIBC and MIBC tissues. In
   addition, other potential prognostic variables correlated to OS in BC,
   such as body mass index (BMI), residual tumor at tur,
   neutrophil-to-lymphocyte ratio (NLR), and lymphovascular invasion
   (LVI), should be investigated. Last, the changes of before and after
   the treatment, such as chemotherapy or Bacillus Calmette-Guérin (BCG)
   refractory, should also be considered to find the potential markers for
   predicting the treatment effect and prognosis.

Conclusion

   In conclusion, based on the comprehensive analyses with ARGs expression
   profiles and corresponding clinical features, three prognostic ARGs
   (JUN, MYC, and ITGA3) were identified. The genes identified in
   autophagy pathways also provide new possibilities for bladder cancer
   therapeutic intervention. By combining molecular signature and clinical
   characteristics, we constructed a novel risk score model ACPI which can
   robustly estimate BC patients’ survival. Also, the ACPI risk score
   model was validated by large sample size. However, further prospective
   experiments can be expected to test the clinical utility and aid in the
   search for optimal personalized targeted therapies.

Acknowledgments