Analysis of Autophagy-Related Signatures in The Tumor Immune Microenvironment and Identication of Clinical Prognostic Markers in Bladder Carcinoma

Objectives: Bladder carcinoma (BLCA) is one of the most common malignant diseases of urinary system. Our study aimed to investigate the autophagy-related signatures in the tumor immune microenvironment and construct effective prognosis prediction model. Methods: RNA sequencing data and corresponding clinical information of BLCA were downloaded from The Cancer Genome Atlas (TCGA) and Gene Expression Omnibus (GEO). Autophagy-related genes were extracted from TCGA dataset for consensus clustering analysis, and differences in survival rate were analyzed. STIMATE algorithm was used to analyze the tumor microenvironment (TME) and immune cell inltration was compared between different clusters. Differentially expressed genes (DEGs) between different clusters were identied, followed by function annotation. Independent prognostic signatures were further revealed to construct prognostic prediction model. Results: We identied 35 autophagy-related genes associated with prognosis. Survival rate of samples in cluster 1 was signicant lower than that in cluster 2. Cluster 2 had markedly lower tumor purity and signicantly higher estimate score and stromal score than cluster 1. The proportions of T cells CD8, macrophages M1, T cells CD4 memory activated, NK cells activated, and dendritic cells activated were higher in cluster 1. There were 1,275 DEGs which were mainly enriched in functions and pathways related to immune response and cancer. Seven genes (ATF6, CAPN2, NAMPT, NPC1, P4HB, PIK3C3, and RPTOR) were further identied as the independent prognostic signatures to construct risk score prediction model, which had good prediction performance. Conclusion: Prognosis prediction model based on 7 autophagy-related genes may have great value in BLCA prognosis prediction.


Introduction
Bladder carcinoma (BLCA) is one of the most common malignant diseases of urinary system and its incidence is increasing worldwide, with an estimated of 573,000 new cases diagnosed and approximately 213,000 deaths caused in 2020 [1,2]. Among the diagnosed BLCA patients, 70% of patients present with super cial tumors at the time of diagnosis and 30% of patients present as muscle-invasive disease associated with a high risk of death [3]. Despite great progress has been made in cancer biology and treatment for BLCA, such as surgery, immunotherapy, radiotherapy, and chemotherapy, the survival outcomes of the patients diagnosed with muscle-invasive and advanced BLCA remain poor [4,5].
Therefore, it is urgently to nd the prognostic markers to guide patient treatment selection and construct the effective prognosis prediction model for BLCA.
The tumor microenvironment (TME) has great importance in the tumor progression and the drug resistance of BLCA due to the heterogeneity and complexity of tumor-associated stromal cells and the extracellular milieu [6]. It has been shown that strati cation of TME immune types is bene cial for selection of immunotherapeutic strategies in tumor [7]. What's more, novel insights into the molecular mechanisms of immune evasion in BLCA are essential for developing immunotherapy that is targeting speci c components of TME [8]. According to an analysis of 542 patients with muscle-invasive bladder cancer, quantities and spatial organization of stromal tumor-in ltrating lymphocytes within TME could predict stages of tumor in ammation and patient survival [9].
Autophagy is an evolutionary catabolic mechanism to maintain cellular environmental homeostasis by recycling and degrading internal constituents, including macromolecules and organelles [10]. Recent ndings have revealed that autophagy plays an important role in the malignant transformation and is proposed as a signi cant factor in tumor development [11,12]. Furthermore, autophagy may affect TME by preventing cytotoxicity and providing sustainable source of biomolecules and cellular energy demands under stressful conditions [13,14]. In the study of Li et al., immune-related prognostic genes for BLCA were identi ed by analyzing the relationship between TME and gene expression pro les of one dataset from Gene Expression Omnibus (GEO), while the roles of autophagy in TME was not explored based on gene expression pro les [15].
In order to further reveal autophagy-related signatures in TME and uncover the potential molecular mechanisms, our current study included 394 BLCA samples as training dataset from The Cancer Genome Atlas (TCGA), and four datasets from GEO as validation dataset. Autophagy-related genes were rstly revealed and samples were clustered based on autophagy-related signatures, followed by differential expression and tumor immune microenvironment analysis. Moreover, prognostic prediction model was constructed based on autophagy-related signatures and compared with clinical prognostic prediction model ( Figure 1).  [17]. A total of 214 autophagy-related genes with non-zero expressions were extracted from TCGA dataset. Afterwards, autophagy-related genes with median absolute deviation (MAD) < 0.05 were excluded. The remained autophagy-related genes were submitted to univariate Cox regression analysis by Survival package (http://bioconductor.org/packages/survivalr/) [18]. Clustering analysis was performed for autophagyrelated genes with MAD > 0.05 and log-rank p value < 0.05 by ConsensuClusterPlus [19] in R based on Euclidean and Ward's linkage, with cycle computation 1,000 times to ensure stability and reliability.

Materials And Methods
Samples were classi ed into two different clusters by optimal k-means clustering and further con rmed using t-distributed stochastic neighbour embedding (t-SNE) algorithm [20]. The differences in survival rate between different clusters were analyzed using the Kaplan-Meier method [21].

Differentially expressed genes and function annotation
DEGs between the two clusters were identi ed by Limma package [17], with the cut-off criteria of |log2 fold change (FC) |> 0.5 and adjusted p < 0.05. Gene ontology (GO) functional enrichment in terms of biological process (BP), cellular component (CC), and molecular function (MF) as well as Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis were conducted for differentially expressed genes via the Database for Annotation, Visualization and Integrated Discovery (DAVID, https://david.ncifcrf.gov/) [22] with a threshold of p < 0.05.

Immune cell in ltration comparison between two clusters
We applied the ESTIMATE algorithm [23] and downloaded R script (https://sourceforge.net/projects/estimateproject/) to calculate the estimate scores, stromal scores, and immune scores, which were used to predict tumor purity and analyze the TME. In order to compare the differences in immune cell subtypes between the two clusters, the CIBERSORT package was used to calculate the proportions of 22 immune cell subtypes for TCGA samples. Samples with p < 0.05 in CIBERSORT analysis were used for further analysis, and the Mann-Whitney U test was used to compare the differences between the two clusters.

Construction of prognostic prediction model based on autophagy-related genes
Autophagy-related genes related to BLCA risk was identi ed by Lasso Cox regression analysis via the glmnet software package of R, and independent prognostic signatures were further revealed by multivariable Cox regression. The risk score was calculated based on the identi ed prognostic signature genes according to the following formula: Where β i referred to the multivariable Cox regression coe cients of signatures, and Exp i denoted the expression levels of signatures.
The receiver operating characteristic (ROC) curves were plotted by survivalROC package in R, and area under the curve (AUC; AUC = 0.5: no discriminatory power, AUC = 1: perfect discriminatory power) was used to evaluate the prognostic accuracy of the prediction model.

Statistical analysis
R software (Version 4.1.0; https://bioconductor.org/packages/release/bioc/html/limma.html) was used for statistical analysis. The differences in OS between groups were analyzed by Kaplan-Meier curves and log-rank test. Lasso Cox regression was applied to reveal the correlations of autophagy-related genes with prognosis, and multivariable Cox regression was used to identify the independent prognostic signatures. The performance of prediction model was evaluated by ROC curves and AUC. Mann-Whitney U test was used to compare the differences between the two clusters. All statistical p-values were bilateral, and p < 0.05 was considered statistically signi cant.

Results
Consensus clustering of autophagy-related genes in two clusters with different clinical outcomes of BLCA Among the 232 autophagy-related genes in HADb, 214 genes with MAD > 0.5 extracted from TCGA dataset were used to perform univariate Cox regression analysis and 35 candidate genes associated with prognosis with p < 0.05 were obtained, such as APOL1, CTSB, P4HB, ITGB1, and EIF4G1 (Table 1). By using the ConsensusClusterPlus in R, the TCGA samples were clustered into different groups based on the expression levels of these 35 autophagy-related genes. The crossover between the BLCA samples was the least when the consensus matrix k value was 2 (two subclasses: cluster 1 and cluster 2), the consensus matrix heat map showed a clear boundary (Figure 2A), and the crossover between the BLCA samples was the least (Figure 2B, 2C). The samples clustered into cluster 1 or cluster 2 were listed in supplementary table 3.
Meanwhile, t-SNE algorithm was used to reduce the dimensionality of features, and the two subclasses were largely consistent with the two dimensional pattern of t-SNE distribution ( Figure 2D). Heatmap of two clusters de ned by the autophagy-related genes showed that the expressions between the two clusters were different ( Figure 2E). Kaplan-Meier curves showed that the survival rate of samples in cluster 1 was signi cant lower than that in cluster 2 with a log-rank p < 0.01 ( Figure 2F).

Immune characteristics of the two clusters
After the estimate scores, immune scores, stromal scores, and tumor purity of the two clusters were achieved, we found the estimate score (p = 0.0325) and stromal score (p = 1.46e-04) of cluster 2 were obviously higher than those of the cluster 1, while signi cantly lower tumor purity (p = 0.0325) was observed in cluster 2 ( Figure 4A).
The proportions of 22 immune cell subtypes were calculated by CIBERSORT package, and cluster 2 had signi cantly higher proportions of the macrophages M2, macrophages M0, T cells CD4 memory resting, T cells regulatory (Tregs), and B cells naive ( Figure 4B). The proportions of T cells CD8, macrophages M1, T cells CD4 memory activated, NK cells activated, and dendritic cells activated were lower in cluster 2 by comparing with the cluster 1( Figure 4B) .

Prognostic model based on autophagy-related genes had good prediction performance
Lasso Cox regression analysis was performed for the 35 autophagy-related genes to minimize the risk of over tting and 24 genes were obtained ( Figure 5A, B). The regression coe cients of these 24 genes were listed in supplementary table 6. Multivariable Cox regression indicated 7 genes (ATF6, CAPN2, NAMPT, NPC1, P4HB, PIK3C3, and RPTOR) were the independent prognostic signatures ( Table 2; Figure 5C). The coe cients and expressions of these 7 genes were used to calculate the risk score. Samples were separated into low-risk or high-risk groups with the median risk score as cutoff value. Kaplan-Meier curves showed that high-risk group had poorer prognosis (p = 0.00019; Figure 5D). With the increasing of risk score, the number of dead patients elevated ( Figure 5E). According to the univariate and multivariate Cox regression analysis, age, stage and risk score were found to be independent prognostic factors ( Figure 5F). ROC curves of risk score, age, gender, stage and grade indicated that risk score had the best prediction performance in 1 year survival probability (AUC = 0.671) ( Figure 5G).

Autophagy-related genes in GEO dataset
Autophagy-related genes were extracted from GSE13507, GSE48075, GSE48276 and GSE69795. Univariate Cox regression analysis was used to reveal the autophagy-related genes associated with prognosis, and MTOR was the overlapping gene with the 35 autophagy-related genes in TCGA dataset (Supplementary gure 2). The ROC curves showed that MTOR expression levels were signi cantly related to tumor risk with p < 0.05 in GEO datasets ( Figure 6).

Discussion
BLCA is one of the top ten most common forms of cancer in the world, and the treatment burden of BLCA is high due to drug resistance, frequent surveillance strategies, high recurrence and metastasis rate [24,25]. Cancer stem cells factors, genetic heterogeneity, autophagy, and TME are involved in the drug resistance of cancer cells, and autophagy may affect the regulation of cancer stem cell homeostasis and TME, [26,27]. Therefore, it is of vital importance to understand the potential mechanism of autophagy and TME regulation in the BLCA.
In our study, 35 autophagy-related genes associated with prognosis were revealed for sample clustering, and the survival rate of samples in cluster 1 was signi cant lower than that in cluster 2 based on the Kaplan-Meier curves. According to ESTIMATE algorithm, the estimate score and stromal score were signi cantly higher in cluster 2 and tumor purity were markedly lower. In the study of Li et al, they found stromal scores were associated with clinical characteristics and prognosis of BLCA patients [15]. This result of our study is consistent with previous study, which indicated that the TME composition may in uence the prognosis of BLCA patients.
Previous researchers have indicated that macrophages in tumor tissues play an indispensable role in tumor immunity and progression depending on the cytokine microenvironment [28]. However, the speci c molecular mechanisms of tumor-associated macrophages in the modulation of TME remain poorly understood. In our current study, comparison of the immune cell proportions in the TME between the two clusters showed that the cluster 2 had signi cantly higher proportions of the macrophages M2, macrophages M0, T cells CD4 memory resting, T cells regulatory (Tregs), and B cells naïve. Meanwhile, the proportions of T cells CD8, macrophages M1, T cells CD4 memory activated, NK cells activated, and dendritic cells activated were lower in cluster 2. These results suggested that the immune cell in ltration was activated in patients with poor prognosis.
A total of 1,275 DEGs were identi ed between cluster 1 and cluster 2. Functional enrichment analysis showed that DEGs were enriched in the GO terms related to immune response and cancer, such as positive regulation of I-kappaB kinase/N-kappaB signaling, interferon-gamma-mediated signaling pathway, apoptotic process, cell proliferation, and angiogenesis. Bacterial immunotherapy for BLCA has the ability to induce tumor immunity causing anti-tumor effects that result from enhanced CD4 T cells, and ultimately requires tumor-intrinsic interferon-gamma signaling [29]. Moreover, DEGs were signi cantly involved in pathways related to tumor such as pathways in cancer, VEGF signaling pathway, cell cycle, proteoglycans in cancer, HTLV-I infection, MAPK signaling pathway, and PI3K-Akt signaling pathway. These results indicated that the DEGs between cluster 1 and cluster 2 classi ed by 35 autophagy-related genes participate in BLCA by affecting the immune microenvironment.
By using Lasso Cox regression and multivariable Cox regression analysis, 7 genes (ATF6, CAPN2, NAMPT, NPC1, P4HB, PIK3C3, and RPTOR) among the 35 autophagy-related genes were identi ed as the independent prognostic signatures to establish risk score prediction model. Kaplan-Meier curves showed that high-risk group had poorer prognosis. Activating transcription factor 6 (ATF6) is a member of the leucine zipper family and regulates endoplasmic reticulum stress [30]. It has shown that high expression of otubain 1 (OTUB1) promotes BLCA progression by stabilizing ATF6 under endoplasmic reticulum stress [31]. Cancer cell chemoresistance is associated with the ATF6-induced autophagy, resulting poor overall survival [32]. Nicotinamide phosphoribosyl transferase (NAMPT) is the rate limiting enzyme in procancer NAD+ synthesis, and it showed that inhibition of nicotinamide phosphoribosyl transferase (NAMPT) can suppress the growth of differentiated cancer cells through autophagy [33]. Genetic variants in NAMPT have been con rmed associated with BLCA risk and could be considered as genetic biomarker for the BLCA prognosis [34]. Previous study demonstrated that MPT0L145 which is a rst-in-class PIK3C3/FGFR inhibitor has potent anti-bladder cancer activity, which enhances autophagosome formation by FGFR inhibition and impairs autophagy ux via PIK3C3 inhibition [35]. What's more, blocking pro-survival autophagy by targeting PIK3C3 with MPT0L145 can signi cantly improve the sensitivity of cancer cells to targeted drugs or chemotherapeutic drugs, and promote the rational combination strategy of cancer treatment [36]. These results indicated that these 7 genes might be the potential targets in developing novel therapies for BLCA.
Univariate and multivariate Cox regression analysis indicated that age, stage, and risk score were independent prognostic factors. ROC curves of risk score, age, gender, stage and grade in 1, 3, 5 years survival probability indicated that risk score had the best prediction performance in 1 year survival probability and had good prediction performance in 3, and 5 years survival probability. Among the 35 autophagy-related genes related to prognosis in TCGA dataset, MTOR was the overlapping gene between GEO datasets and TCGA dataset. The ROC curves showed that MTOR expression levels were signi cantly related to tumor risk in GEO datasets.
In conclusion, 35 autophagy-related genes associated with prognosis were revealed and 1,275 DEGs were identi ed between cluster 1 and cluster 2 classi ed by 35 autophagy-related genes. Cluster 2 had markedly lower tumor purity and signi cantly higher estimate score and stromal score than cluster 1. The proportions of T cells CD8, macrophages M1, T cells CD4 memory activated, NK cells activated, and dendritic cells activated were higher in cluster 1. The identi ed DEGs were mainly enriched in functions and pathways related to immune response and cancer. Moreover, 7 genes (ATF6, CAPN2, NAMPT, NPC1, P4HB, PIK3C3, and RPTOR) were further identi ed as the independent prognostic signatures and risk score prediction model was constructed, which had good prediction performance.

Declarations
Ethics approval and consent to participate: Not applicable.
Consent for publication: Not applicable.  Figure 1 Flow chart of the study.    Survival curve of MTOR in GEO datasets GSE13507, GSE48075, GSE48276 and GSE69795.

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download.