Three Prognostic Biomarkers Correlate With Immunotherapy Response in Bladder Urothelial Carcinoma

ya guo (  gy8569851@xjtu.edu.cn ) Xi'an Jiaotong University Second A liated Hospital Department of Oncology https://orcid.org/00000002-0343-0976 Yin Bin Zhang Xi'an Jiaotong University Second A liated Hospital Yi Li Xi'an Jiaotong University Second A liated Hospital Wang Hui Su Xi'an Jiaotong University Second A liated Hospital Shan He xi'an jiaotong university second a liated hospital Shu Pei Pan Xi'an Jiaotong University Second A liated Hospital Kun Xu xi'an jiaotong university second a liated hospital Wei Hua Kou xi'an jiaotong university second a liated hospitla

reliable and repeatable biomarkers to assess the prognosis of BLCA patients (4)(5)(6)(7). Therefore, there is an urgent need to develop reliable molecular markers for predicting prognosis in BLCA.
In recent years, immunotherapy has attracted considerable attention for its in uence on the treatment of locally advanced and metastatic BLCA (8). However, only a proportion of patients respond to treatment with immune-checkpoint inhibitors (9,10). Therefore, identi cation of biomarkers to effectively predict prognosis and evaluate the bene t of immunotherapy is of great signi cance. Previous studies have reported that highly in ltrating lymphocytes are related to prognosis of BLCA (11). The evidence indicates that CD8+ T cells are involved in tumor adaptive immunity, and their in ltration is associated with prognosis (12,13). Programmed death ligand 1 (PD-L1), also known as B7-H1 or CD274, is involved in inhibiting T cell-mediated antitumor immunity through interaction with PD-1 (14) (15). Previous studies have reported that high expression of PD-L1 is associated with worse cancer outcomes in various malignancies (16). CD8+ T cell in ltration and tumor mutation burden (TMB) have been reported to be correlated with response to atezolizumab (anti-PD-L1) in metastatic urothelial cancer (mUC) (17)(18)(19).
Therefore, identi cation of biomarkers related to CD8+ T cell in ltration, PD-L1 expression, and TMB will help to predict prognosis and immunotherapy response in patients with BLCA.
In this study, we rst identi ed DEGs from the GEPIA (Gene Expression Profiling Interactive Analysis) and Oncomine databases, and then identi ed the overlapping DEGs between them using a Venn diagram. We further performed gene enrichment and protein-protein interaction (PPI) analyses to select hub genes. Moreover, prognosis-related hub genes were identi ed using the survival R package and con rmed in three online databases. We then explored the expression of the key prognosis-related genes with different clinical factors using the UALCAN database. Finally, we evaluated whether three key genes related to CD8+ T cell in ltration, TMB, and PD-L1 expression could be used to predict prognosis and monitor immunotherapy response in BLCA based on comprehensive bioinformatics analysis.

Methods
Identi cation of DEGs GEPIA (http://gepia.cancer-pku.cn/) is an online network tool based on data from The Cancer Genome Atlas (TCGA) and GTEx, which can be used to study interactions between DEGs, as well as for survival analysis, pro le plotting, and detection of similar genes (20). Oncomine is an online resource containing microarray data (https://www.Oncomine.org) (21). In this study, we rst used the GEPIA and Oncomine databases to identify DEGs by comparison of tumor samples with normal samples. In the GEPIA database, mRNAs with q<0.01 and |log2 fold change (FC)|≥2 were chosen as DEGs. In the Oncomine database, the selection criteria for DEGs are P<0.01, |log2FC|≥2, and gene rank ≤10%. Then, we identi ed the overlapping DEGs between them using a Venn diagram (http://bioinformatics.psb.ugent.be/webtools/Venn/) (22).

Functional Analysis and Pathway Enrichment Analysis
Metascape (http://metascape.org/) is an online resource for gene annotation and analysis (23). In the present study, Metascape was used to perform gene ontology (GO) and pathway analyses of 69 common hub genes. Pathway and process enrichment analyses were conducted based on several sources, including GO biological processes, The Kyoto Encyclopedia of Genes and Genomes pathways, reactome gene sets, and CORUM. Terms with a P-value less than 0.01, a minimum count of 3, and an enrichment factor greater than 1.5 were considered to represent signi cant processes or pathways.

Construction of PPI Network and Identi cation of Hub Genes
We evaluated PPI information of common genes using the STRING online database (https://stringdb.org/cgi/input.pl) and then visualized the resulting interaction network using Cytoscape software (http://www.cytoscape.org/) (24,25). A con dence score greater than 0.4 was de ned as signi cant. The Molecular Complex Detection (MCODE) plugin in Cytoscape was used to further screen key genes in the PPI network with degree cutoff = 5, K-score = 2, and node score cutoff = 0.2.

Prognostic Value of Hub Genes
We downloaded gene expression pro le and clinical data from TCGA (https://portal.gdc.cancer.gov/).
The prognostic value of 11 identi ed hub DEGs that were analyzed using the R survival package (26). PROGgenesV2 (http://genomics.jefferson.edu/proggene/ lter.php) is a web resource that allows researchers to study the correlations between genes and overall survival (OS) in multiple cancers based on TCGA and GEO data (27). PrognoScan (http://www.prognoscan.org/) was used to evaluate the associations between gene expression and patient prognosis, according to measures including OS and disease-free survival (DFS), across a large collection of publicly available cancer microarray datasets (28). The OSblca database (http://bioinfo.henu.edu.cn/BLCA/BLCAList.jsp) provides a useful tool to assess novel prognostic biomarkers in bladder cancer, based on data from 1,075 bladder cancer patients, including OS, disease-speci c survival (DSS), disease-free interval, and progression-free interval (29). In this study, we further con rmed the prognostic value of key genes based on the abovementioned three databases. The hub genes identi ed in this way were considered to be key prognosis-related genes.

Association between Prognosis-related Key Genes and Clinical Characteristics
The UALCAN database (http://ualcan.path.uab.edu/) is a comprehensive web resource for analyzing cancer OMICS data (TCGA and MET500) (30,31). A previous study reported that the expression and prognostic value of DEGs were associated with clinical characteristics, including TNM stage, smoking history, lymph invasion, and histological type (32). Therefore, we explored the relationship between the three key prognosis-related genes and clinical characteristics using UALCAN.
Immune Cell In ltration and TMB analysis TIMER is a user-friendly web portal for the systematic analysis of immune infiltrates across different types of cancer (https://cistrome. shinyapps.io/timer/) (33). In this study, we used TIMER to analyze the associations between the three key genes and immune cell in ltration. P<0.05 was considered to be statistically signi cant. In addition, we downloaded 33 cancer expression pro les with associated mutation data from TCGA, and assessed the association between the prognosis-related hub genes and TMB (26,34).
TISIDB Analysis TISIDB (http://cis.hku.hk/TISIDB) is a publicly available resource that allows the user to explore the function of a gene and its role in tumor-immune features (35). TISIDB consists of 10 modules: function, literature, screening, immunotherapy, lymphocyte, immunomodulator, chemokine, subtype, clinical, and drug. In this study, the lymphocyte module was used to evaluate the relationships between the expression of selected genes and tumor-in ltrating lymphocytes. The immunomodulator module was used to examine the associations between PDL1 and selected genes. The subtype module was used to assess the distribution of selected genes' expression across immune and molecular subtypes. We used the screening module to explore whether the expression of prognosis-related genes showed signi cant differences between responders and non-responders to immunotherapy. Using the drug module, we analyzed the drugs in the DrugBank database targeting these three genes.

Identi cation of Hub Genes
A total of 750 DEGs were identi ed from the GEPIA database, and 1,881 DEGs were identi ed from Oncomine. Sixty-nine common genes were screened out using the Venn diagram ( Fig. 1 A). We then performed GO and pathway enrichment analyses for the common genes using Metascape. The results showed that these common genes were involved in 20 main GO terms and pathways, of which the top ve were mitotic cell cycle, cGMP-PKG signaling pathway, extracellular matrix organization, muscle contraction, and response to hydrogen peroxide ( Fig. 1 B). Finally, a PPI network of DEGs was constructed using the STRING and Cytoscape software, containing 69 nodes and 83 edges ( Fig. 1 C). One signi cant module was identi ed using MCODE. This module contained 11 genes, Aurora-A kinase (AURKA), BIRC5, CENPA, CKS1B, ECT2, MYBL2, NUF2, RRM2, TK1, TPX2, and UBE2C, which were considered to be hub genes ( Fig. 1 D).

Prognosis-related Hub Genes
The association between the expression of 11 hub genes and overall survival was evaluated using the R survival package. High expression of three genes (AURKA, BIRC5, CKS1B) was related to an unfavorable prognosis ( Fig. 2 A-K). We further validated the prognostic value of the three genes in BLCA using the PROGgenesV2, PrognoScan, and OSblca databases. Our results showed that BLCA patients with higher expression levels of the three hub genes exhibited poorer OS, DFS, and DSS, indicating that the three hub genes may be associated with unfavorable prognosis (Table 1). In summary, our data suggested that the three key genes could serve as biomarkers of poor prognosis.

Association between Three Key Genes andClinical Parameters in BLCA
A previous study reported that the OS of BLCA patients was signi cantly associated with clinical characteristics, including TNM stage, smoking history, lymph invasion, and histological type (29). The three hub genes identi ed here were associated with various clinical characteristics including age, histological subtype, molecular subtype, nodal metastasis status, sample type, smoking, cancer stage, and TP53 mutation status ( Table 2). Most importantly, overexpression of the three hub genes was positively correlated with histological subtypes (Fig. 3 A-C). The expression of the three hub genes was higher in the basal squamous and neuronal subtypes than in the luminal subtype ( Fig. 3 D-F). BLCA patients with TP53 mutations also showed high expression of the three hub genes (Fig. 3 G-I). Taken together, these results suggest that increased expression of the three key genes might predict poor prognosis in patients with BLCA.

CD8+ T Cell In ltration Predicts Poor Prognosis in BLCA
The TIMER and TISIDB databases were used to explore the relationship between the three prognosisrelated genes and tumor-infiltrating immune cells. The three genes were positively associated with levels of in ltrating CD8+ T cells, neutrophils, and dendritic cells. Expression of BIRC5 was negatively correlated with in ltration of B cells (Table 3, Fig. 4 A-F). High levels of in ltration of CD8+ T cells were associated with poor prognosis (Fig. 4 G). These results suggest that these genes may affect prognosis via regulation of CD8+ T cells.

Three Key Genes Correlated with Immunotherapy Response
Previous studies have suggested that TMB and PD-L1 expression are correlated with response to atezolizumab in mUC (17). Our results indicated that the three key genes were associated with TMB in multiple cancers, including adrenocortical carcinoma, BLCA, and breast invasive carcinoma ( Fig. 5 A-C). We further discovered that the expression levels of the three hub genes were positively correlated with in ltration of PD-L1 (CD274) expression (Fig. 5 D-F). Finally, our results showed that the three hub genes exhibited a signi cant difference in expression between responders and non-responders to atezolizumab in urothelial cancer via TISIDB analysis (Table 4). Taken together, these results suggest that the impact of the three hub genes on response to immunotherapy in BLCA may be associated with TMB and PD-L1 expression. The three key genes may thus represent a promising immunotherapeutic target. Finally, we found that the three prognosis-related genes had the highest expression levels in the C2 (IFN-gamma dominant) subtype and the lowest in the C3 (in ammatory) subtype (Fig. 5 G-I).

Drug-Gene Interaction Network
Drugs targeting three key genes were collected from the DrugBank database. AURKA and 18 other targets were correlated with 15 drugs. BIRC5 and two targets were related to 15 drugs. CKS1B and another three targets interacted with three drugs (Fig. 6 A-C). Our results may contribute to the development of new targets for BLCA immunotherapy.

Discussion
BLCA is one of the most common urinary cancers in the world, with high recurrence and mortality rates that limit the e cacy of treatment (36). It is therefore essential to understand the molecular mechanism of BLCA. With the development of bioinformatics technology, increasing numbers of studies are using bioinformatics analysis to develop biomarkers and explore the molecular mechanism of BLCA (3, 7). However, there are still no reliable biomarkers associated with prognosis of BLCA patients. In recent years, immunotherapy has attracted considerable attention owing to its in uence on the treatment of locally advanced and metastatic BLCA (8), but limited information is available about biomarkers related to immunotherapy response. Overall, there are no satisfactory signatures to effectively predict prognosis and evaluate the bene t of immunotherapy in BLCA patients.
AURKA is a member of the serine/threonine kinase family and has a role in regulation of the cell cycle (37). Accumulating evidence indicates that AURKA is overexpressed in various cancers, including breast cancer, head and neck cancer, esophagus cancer, hematological malignancies, colorectal cancer, stomach cancer, pancreatic cancer, and ovarian and prostate cancers (38)(39)(40). Pathological overexpression of AURKA is correlated with shorter survival of cancer patients. According to previous reports, high expression of Aurora-A in tumor cells is closely related to poor prognosis (37,40,41). BIRC5 is a member of the inhibitor of apoptosis gene family, which has dual roles in promoting cell proliferation and preventing apoptosis (42). Overexpression of BIRC5 has been reported in several malignancies, and higher BIRC5 expression was also found to be associated with decreased survival (43). CKS1B is an oncogene that has been reported to show increased expression in various tumors (44). In accordance with previous studies, our results demonstrated that high expression of the three hub genes was associated with poor prognosis (Fig. 2 L-V). Further analysis revealed that these three key genes were overexpressed in non-papillary tumors, the basal squamous subtype, and TP53 mutation patients (Fig. 3 A-I). Moreover, our results indicated that the three hub genes were highly correlated with TMB and PD-L1 expression ( Fig. 5 A-F). Expression of the three unfavorable-prognosis-related genes was highest in the C2 (IFN-gamma dominant) subtype and lowest in the C3 (in ammatory) subtype (Fig. 5 G-I). Previous studies have reported that TP53 mutation is associated with poor prognosis (45). TMB has prognostic roles in various cancer types, including BLCA (46). A previous study reported that the C2 subtype had the highest levels of CD8+T cells, as well as having less favorable outcomes, whereas the C3 (in ammatory) subtype had the best prognosis (47). Taken together, these results indicate that the expression of these three key genes may represent a marker of poor prognosis in BLCA, as well as in non-papillary tumors, the basal squamous subtype, TP53 mutation patients, tumor with high TMB, and the C2 subtype.
Many types of human tumors are affected by tumor-in ltrating immune cells (TIICs) (48). A considerable number of studies have indicated thatTIICs affect the prognosis of patients and alter response to immunotherapy (49) (50). Higher in ltration levels of CD8+ T cells are associated with poor prognosis in gastric cancer (51). In line with previous studies, our results demonstrated that high in ltration of CD8+ T cells was associated with poor prognosis (Fig. 4 G). Moreover, we found that the three prognosis-related genes were positively related to CD8+ T cells (Fig. 4 A-F). These data suggest that CD8+ T cell in ltration has prognostic value in BLCA, and that the three genes comprising our prognosis-related signature may affect prognosis via CD8+ T cell in ltration. Increased CD8+T cell in ltration has been reported to be correlated with better immunotherapeutic effect (49). A previous study showed that biomarkers related to CD8+ T cell in ltration could facilitate the monitoring of immunotherapy response and the exploration of the immune in ltration mechanism in clear cell renal cell carcinoma (19). AURKA is overexpressed in cancer cells but not in normal tissues, making it a potential target for immunotherapy. AURKA-speci c CD8(+) T cells can selectively lyse leukemia cells (52). BIRC5 is correlated with T cell survival and proliferation, which can increase the accumulation and persistence of CD8(+) T cells following an encounter with Ag (53). A miR-197/CKS1B/STAT3-mediated network that promotes tumor PD-L1 expression and miR-197 replacement therapy may be a potential treatment for chemotherapy-resistant non-small-cell lung cancer (54). In this study, a positive association between the expression of three key genes and immune cells, particularly CD8+ T cells, was examined using the TIMER/TISIDB databases. Previous studies have reported that high expression of PD-L1 is associated with worse cancer outcomes (16). TMB has been investigated in various malignancies and found to be correlated with response to atezolizumab in mUC (17)(18)(19). We found that the expression levels of the three hub genes were positively correlated with PD-L1 expression and TMB (Fig. 5 A-F). Patients have previously been classi ed into immunotherapeutic responders and non-responders based on number and distribution of TIICs (55). We found that the three prognosis-related genes were overexpressed in responders to atezolizumab among patients with urothelial cancer (Table 4). Finally, our results demonstrated that AURKA and another 18 targets were correlated with 15 drugs. BIRC5 and two targets were related to 15 drugs. CKS1B and another three targets interacted with three drugs (Fig. 6 A-C). Collectively, these data suggest that the response of the three hub genes to PD-L1 in BLCA may be associated with CD8+ T cells, TMB, and PD-L1 expression. The three key genes may represent a promising immunotherapeutic target. Drug-gene interaction network analysis indicated that these genes and their related drugs could be used in the development of new targets for BLCA immunotherapy.

Conclusions
In summary, three key genes in BLCA were found to be correlated with poor prognosis and immune cell in ltration, especially that of CD8+ T cells. The responses of these prognosis-related genes to immunotherapy in BLCA may be associated with CD8+ T cells, TMB, and PD-L1 expression. These three key genes and their related drugs may help to develop new targets for BLCA immunotherapy. Further investigations and experiments using methods such as quantitative real-time PCR and western blotting, and clinical data analyses are required to validate the results of the present study. We will perform such further experiments to support our results.

Availability of data and materials
The data that support the ndings of this study are openly available from GEPIA (http://gepia.cancerpku.cn/) and Oncomine (https://www.Oncomine.org).  Tables  Table 1 Con rmation of the associations of three hub genes with prognosis in three different databases PROGgeneV2, PrognoScan, and OSblca databases were used to con rm the prognostic value of three hub genes in BLCA. HR, hazard ratio.            Evaluation of association between three prognosis-related key genes and clinical factors using UACLAN.

Ethics approval and consent to participate
Expression of three key prognosis-related genes based on different sample types, according to (A-C) histological subtypes, (D-F) molecular subtypes, and (G-I) TP53 mutation status.

Figure 3
Evaluation of association between three prognosis-related key genes and clinical factors using UACLAN.
Expression of three key prognosis-related genes based on different sample types, according to (A-C) histological subtypes, (D-F) molecular subtypes, and (G-I) TP53 mutation status.

Figure 3
Evaluation of association between three prognosis-related key genes and clinical factors using UACLAN.