Immune and Stromal Signature-Based Prognostic Gene of Bladder Cancer


 Background: Immune and stromal cells in the tumor microenvironment have exhibit critical relevance with tumorigenesis and progression. Therefore, our analysis was conducted to explore the potential prognostic factors and the therapeutic targets of bladder cancer (BC) that affiliate to immune and stromal infiltration signature.Methods: Immune and stromal score were calculated for BC patients from TCGA database using ESTIMATE algorithm to predict the level of infiltration. Kaplan-Meier curves were utilized to assess the correlation of immune/stromal infiltration with survival outcomes. Differential expression genes (DEGs) were identified. Enrichment analyses, Kaplan-Meier curves, protein-protein interaction (PPI) network construction were performed to describe and screen the core module genes, whose prognostic value was further validated by an independent GEO dataset. Transcriptional factors (TFs) and ncRNA. correlated with the core module were identified using RAID 2.0 and TRRUST 2.0 database, and drugs potentially regulative for BC were accordingly screened using DrugBank.Results: 393 and 93 BC patients were enrolled from TCGA and GEO respectively. Higher stromal infiltration indicated worse overall survival (P = 0.015), and higher immune infiltration indicated an improvement on overall survival (P = 0.042). 562 DEGs were identified and 123 of them associated with survival outcomes. PPI has identified the core module genes, in which EFNB2 was the only core prognostic gene that was validated by both TCGA (P = 0.042) and GEO dataset (P = 0.036). Four TFs and Five ncRNAs (e.g. HIF1A) were associated with the regulation of the core module genes, and six drugs were screened as potential candidates to regulate BC.Conclusion: Higher immune infiltration in bladder tumor was correlated with improved survival and higher stromal infiltration corelated with worse survival. Furthermore, a higher expression of EFNB2 was tied with poorer prognosis of BC, which was validated by two independent datasets.


Introduction
The tumor microenvironment of malignant solid tumor is believed to consist not only tumor cells but also normal cells. For instance, stromal cells and immune cells have been revealed to play vital roles in tumor growth, progression 1, 2 , metastasis and drug resistance 3,4 . The algorithm ESTIMATE (Estimation of STromal and Immune cells in MAlignant Tumor tissues using Expression data) was developed based on this concept 5 to predict the in ltration level of immune and stromal cells. The algorithm has been proved accurate across 11 different tumor types and applied to effectively describe immune/stromal signature and explore prognostic genes in glioblastoma 6 and gastric cancer 7 .
Bladder cancer is one of the most common urological malignant tumors, causing roughly 200,000 death in 2018, ranking 13 th among all tumors 8 . A past study has demonstrated the association of the in ltration of stromal and immune cells with the survival outcomes of muscle-invasive bladder cancer (MIBC) 9 . Therefore, one of our aims is to validate and expand the outcomes of the previous study in both non-muscle-invasive and muscle-invasive bladder cancer applying the ESTIMATE algorithm in TCGA and GEO database. Moreover, since several studies have indicated that the current TNM stage system provides limited prognostic value for malignant tumors 10,11 , it is expected that our study could identify the biomarkers of bladder cancer at a molecular level and introduce more information to the current TNM stage system and help physician in treatment decision making.
Method And Materials 2.1 Data extraction and calculation of immune / stromal score The study methodologies conformed to the standards set by the Declaration of Helsinki. Level 3 gene expression data of bladder cancer along with clinicopathological variables including gender, age, race, histology classi cation, tumor stage and survival outcomes were downloaded from The Cancer Genome Atlas (TCGA) Genomic Data Commons data portal (https://portal.gdc.cancer.gov/). Another independent dataset GSE31684 was also downloaded from Gene Expression Omnibus (GEO) for further external validation. Patients with survival less than one month were excluded. Immune and stromal score were calculated to present the level of in ltration of immune and stromal cells for each patient using the ESTIMATE algorithm 5 .

Correlation of immune and stromal score with tumor stage and prognosis
Correlation between immune/stromal score with tumor stage was assessed using Kruskal-Wallis test.
Optimal cut-off value of immune and stromal score was identi ed by "Maxstat" R package, a method previously described by Hothorn 12 , and patients were categorized into high or low immune/stromal score groups. Correlation of survival outcome with immune/stromal score was assessed with Kaplan-Meier curves.
2.3 Identi cation of differential expression gene (DEG) and GO and KEGG enrichment analyses of DEG Fold change >2 and P value < 0.05 were set to identify differential expression genes (DEG) with R package "limma". DEGs were visualized by heatmap and intersected with Venn diagram. Gene Ontology (GO) biological process including cellular component, molecular function and biological process, and KEGG pathway (Kyoto Encyclopedia of Genes and Genomes) enrichment analyses with false discovery rate (FDR) < 0.05 were performed towards the intersected DEGs in Venn diagram to screen the pathways that DEGs were correlated with.
2.4 Screening of prognostic DEGs, construction of protein-protein interaction (PPI) network and identi cation of the core module The prognostic value of DEGs was con rmed by Kaplan-Meier curve with a P value < 0.05 indicating that the expression of a particular gene was associated with the prognosis of bladder cancer. Prognostic DEGs were subsequently put into the STRING database (https://string-db.org/) and a protein-protein interaction (PPI) network was constructed. Densely connected modules with no less than 3 nodes in the network were de ned as the core modules.

Enrichment analyses and validation (GEO dataset) of the core module
GO and KEGG pathway enrichment analyses with false discovery rate < 0.05 were re-performed for the core modules. In addition, Kaplan-Meier curves were applying to the GEO dataset GSE31684 to validate the prognostic value of the core modules identi ed in TCGA database.
2.6 Transcription factors (TF) and ncRNAs associated with the core module and potential drugs exploration for bladder cancer Moreover, pivot nodes were analyzed to assess the association of the core module with ncRNA and with transcription factors (TF) using RAID 2.0 and TRRUST v.2 database, respectively 13,14 . On the basis of the core module and its associated ncRNA/TF, we further screened drugs that could potentially impact bladder cancer using DrugBank database (https://www.drugbank.ca/) 15 .

Calculation of immune and stromal score and its correlation with tumor stage and prognosis
A total of 393 and 93 bladder cancer patients were identi ed respectively from the TCGA and GSE31684 dataset. Immune and stromal score were calculated and the optimal cut-off value were identi ed as 1895.774 for immune score and -858.8528 for stromal score ( gure 1A). Higher stomal in ltration (46 patients) was correlated with worse overall survival but higher immune in ltration (226 patients) was correlated with improved overall survival ( gure 1B and 1C). Furthermore, stromal in ltration was associated with tumor stage (P < 0.001) while immune in ltration was not ( gure 1D and 1E).

Identi cation, function enrichment analysis and prognostic gene screening of DEGs
Volcano plot revealed differential expression genes between the high/low score groups. 136 DEGs were identi ed in the stromal groups and 468 DEGs were found in the immune groups ( gure 2A and 2B), which were then visualized in the heatmaps ( gure 2C and 2D). Venn diagram intersected all of those DEGs and found 42 of them were differentially expressed in both immune and stromal groups. Therefore, a total of 562 DEGs were eventually identi ed ( gure 2E). Subsequent enrichment analysis revealed the potential roles of DEGs in multiple pathways. Particularly, those DEGs were found to play a role in the urogenital system development, renin-angiotensin system and P53 signaling pathways ( gure 3). Kaplan-Meier survival curve screened 123 prognostic genes out of the 562 DEGs (P < 0.05). Part of the prognostic genes were displayed in gure 4.
3.3Construction of PPI network, identi cation and function enrichment analysis of the core module, validation of the prognostic meaning of the core module in GEO dataset A PPI network was constructed on the basis of the 123 prognostic genes to better understand the interactions between those genes. The core module with 34 nodes (MARK) of the PPI network was identi ed ( gure 5A). Function enrichment analysis was performed repeatedly and revealed that those 34 prognostic genes correlated with multiple pathways, for instance, HIF-1 signaling and TNF signaling pathways ( gure 5B and 5C). The validation of those 34 core prognostic genes in the GSE31684 dataset indicated that only the outcome of EFNB2 was consistent with TCGA (P = 0.036) ( gure 5D).
3.4 Regulation of the core module by ncRNA and TF and potential drugs for bladder cancer On the basis of the ncRNA-mRNA interaction and TF-mRNA interaction, pivot function of the core module was revealed using RAID 2.0 and TRRUST 2.0 database. Five ncRNAs (MALAT1, ABALON, TUG1, FENDRR and ANCR) and four TFs (HIF1A, IRF1, RUNX3 and STAT3) were associated with the regulation of the core module (table 1), which were also visualized in gure 5E. Moreover, we screened potential drugs for bladder cancer according to the pivot function of core module genes using DrugBank. Six drugs were found to be correlated with the core module and potentially regulable for bladder cancer, including HIF inhibitors like FG-2216, ENMD-1198, 2-Methoxyestradiol and PX-478, Carvedilol, a drug related with drug resistance of bladder tumor cell, and Emricasan, an investigative drug for non-alcoholic steatohepatitis (table 2).

Discussion
Several components of the tumor microenvironment were approved to create various compounds to regulate and participate in the progression of bladder tumor 16 . A recent study has revealed the signature of immune and stromal in ltration of MIBC could impact the treatment response and correlated signi cantly with the survival outcomes 9 . However, given that the majority of BC are diagnosed as nonmuscle-invasive and the patterns of immune and stromal in ltrates are varied among tumors, the necessity to investigate the correlation of immune and stromal in ltration with both non-muscle invasive and muscle-invasive bladder cancer as a whole would be elevated. Therefore, we employed the ESTIMATE algorithm, a method which has been veri ed to accurately predict the level of in ltration of immune and stromal cells in the tumor microenvironment across eleven types of tumors 5 and applied to effectively explore prognostic biomarkers in glioblastoma 6 and gastric cancer 7 , to characterize the immune / stromal in ltration signature and screen potential biomarker for bladder tumors with 393 and 93 patients respectively identi ed from TCGA and GSE31684.
It has been reported that higher in ltration of immune cell was associated with the response of immune checkpoint blockade in metastatic bladder cancer 17,18 and improved disease-speci c survival in muscleinvasive bladder cancer receiving trimodality therapy 9 . Despite the underlying mechanism still remains unclear, the studies imposed that treatment such as radiation could enhance immune-mediated tumor cells killing through promoting additional antigen release or chemokine secretion in the favorable microenvironment abundant of immune cells prior to treatment 19,20 . Similarly, regardless of treatment options, our study found that higher immune in ltration was correlated with improved overall survival of BC. In terms of the stromal in ltration, evidence has shown that the components of the stroma could contribute to the direct exclusion or inactivation of chemotherapy 21 and higher stromal in ltration signature was associated with worse survival after radical cystectomy 9 . In consistent, we found that higher stromal in ltration was associated worse overall survival of BC. Moreover, our outcomes also indicated that the in ltration of stromal cell was associated with bladder tumor stage, supporting the idea that stromal cells play critical roles in tumor progression 1,22 .
Through the signature of immune and stromal in ltration, we found that multiple DEGs were associated with the prognosis of bladder cancer but EFNB2 was the only core prognostic gene after external validation with GEO dataset. EFNB2 encodes the B class ephrin, a kind of membrane-bound ligands that binds to Eph receptors that comprises the largest protein-tyrosine kinase receptor family. Several studies have reported the association of ephrin with tumor progression, invasion and migration [23][24][25] . Particularly, FENB2 was reported as a poor prognostic indicator of ovarian and esophageal squamous cell carcinoma 26,27 . Recent studies have also demonstrated the role of EFNB2 in bladder cancer by showing that higher expression of EFNB2 was associated with poorer survival with samples from Johns Hopkins Hospital 28 .
Moreover, they also claimed that FENB2 signi cantly correlated with TP53 mutation 28 , a well-known tumor suppressor gene, and we performed GO enrichment analysis to nd P53 signaling pathway was signi cantly associated with the core module genes. Another study also showed that EFNB2 is a part of the PPI network centered by broblast growth factor receptor 3 which TP53 was involved with as a TF 29 .
Together with previous evidences, our study might promote the role of immune/stromal in ltration and EFNB2 as predictors or therapeutic targets for bladder cancer.
The hypoxia inducible factor 1 alpha (HIF-1A or HIF-1α) functions as an essential regulator of the transcription of many genes involving tumor angiogenesis. Multiple studies have demonstrated the correlation between HIF-1A with bladder cancer [30][31][32] . Earlier studies have reported that HIF-1A is associated with unfavorable prognosis of bladder cancer 30,31 , while a more recent research claims that the higher expression of HIF-1A can improve the local relapse-free survival of MIBC after radiation therapy plus carbogen and nicotinamide 32 . Another study also states that HIF-1A might play a critical role in gemcitabine resistance of bladder cancer 33 , supported by the evidence that the HIF-1A/MDR1 pathway confers chemoresistance to cisplatin in bladder cancer 34 . Similarly, we found that HIF-1A plus IFR1, RUNX3 and STATS were TFs regulating the core module of bladder cancer. We then screened six kinds of drugs that could potentially impact bladder cancer including Carvedilol, FG-2216, ENMD-1198, PX-478, 2-Methoxyestradiol and Emricasan. Carvedilol is a beta-adrenoceptor antagonist which has been reported to be associated with reversal of the abnormal regulation of HIF-1A and reverse of multidrug resistance to anticancer drugs(good job) 35,36 , while FG-2216, ENMD-1198, PX-478 and 2-Methoxyestradiol are all inhibitors of HIF-1A [37][38][39][40][41] . Emricasan, a caspase inhibitor and investigational drug for non-alcoholic steatohepatitis, was previously proven to be able to augment killing acute myeloid leukemia tumor cells by birinapant 42 . Our ndings enhanced the idea that HIF-1A could be a potential target for the management of bladder cancer and provided more thoughts of alternative drugs treating bladder cancer.
This was the rst study applying ESTIMATE algorithm in describing the immune and stromal in ltration signature BC. However, limitations should be noted before interpreting our ndings. The rst limitation is the retrospective nature of our study. With two independent public datasets employed for analysis, multiinstitution samples would improve the reliability. Furthermore, samples that were used in this study were from the core region of the tumor tissue. Those tissue samples might alienate themselves from the samples that extrapolated from different parts of the tumor in other analysis. Together, Those limitations should suggest the potential of uncharted territories for future studies.

Conclusion
The current study revealed that higher immune in ltration in bladder tumor was associated with improved survival outcomes while higher stromal in ltration was correlated with worse survival. Furthermore, higher expression of EFNB2 was validated by two independent datasets to be associated with poorer prognosis of BC.    Prognostic value of DEGs. 123 of those 562 DEGs were found to correlate with the overall survival of bladder cancer and a part of them were displayed here. DEG = differential expression gene.