Identication of RB1 Gene as a Immune-Related Prognostic Signature Based on Tumor Mutation Burdens in Bladder Cancer

Background: Bladder cancer (BC) is known as the eleventh most common malignant tumor all over the world, for either males or females. Developing effective regimens targeting more promising biomarkers aiming for better prognosis are required. Immune checkpoint inhibitors (ICI) have been demonstrated as a prospective and practical means to resist cancers. Theoretically, adequate inltration of immune cells indicates more immunotherapy targets and may promise better prognosis. Methods: Full transcriptome data (n=433), clinical information (n=581) and mutation sequencing (n=412) were obtained freely from The Cancer Genome Atlas and independent mutation sequencing data of 101 samples were acquired from International Cancer Genome Consortium. Statistical processing was conducted using R packages with R x64 4.0.2. Gene biologically functional research was performed with gene set enrichment analysis (GSEA) based on Kyoto Encyclopedia of Genes and Genomes (KEGG) database. 22 types of immune cell inltration were assessed and calculated in 398 samples of BC tumors. Results: Tumor mutation burdens (TMB) of mutant type groups were higher than wild type groups for 19 genes, except for FGFR3 and CREBBP verifying that genomic mutation associates positively with TMB in BC tumor. Kaplan-Meier analysis showed high mutation frequency on RB1 had a negative effect on prognosis of BC patients and RB1 was an independent prognostic factor (p=0.004, HR=1.776) in BC. It was also demonstrated that RB1 mainly participate in singling pathways of cell proliferation and cell cycle. Proportions and correlation of 22 types of immune cells in 433 samples were determined. Immune cells with similar function are inclined to co-exist in tumor microenvironment of BC. Among them, regulatory T cells (Tregs) were detected as a negatively correlated type immune cell to mutation of RB1 that probably increases the incidence of tumor immune escaping in BC. Conclusion: RB1 identied as independent prognostic predictor, and there is a chance for contribution to poor overall survival as the mutation occurs. more, RB1 a biomarker


Introduction
According to the latest research, bladder cancer (BC) is known as the eleventh most common malignant tumor all over the world [1]. However, in Europe and America, the incidence of BC is the fourth highest among males, ranking after prostate cancer, lung cancer and colon cancer [2]. BC can be classi ed into non-muscular invasive bladder cancer (NMIBC) and muscular invasive bladder cancer (MIBC), yet the latter contributes to worse prognosis and poorer response to radiotherapy or chemotherapy.
Unfortunately, almost 25% of patients are identi ed as MIBC at rst diagnosis, suffering from relatively poorer 5-year survival of 50% [3,4]. Considering that the existing therapy strategies, including chemotherapy, radiation therapy, and radical cystectomy depending on particular situations, are unable to bring about satisfying improvements of long-term outcomes, new regimens targeting promising biomarkers are urgently needed.
Currently, a novel therapy aimed at the discovery of immune checkpoint inhibitors (ICIs) has been demonstrated as a prospective and practical means to resist cancers [5]. For instance, inhibition of programmed cell death-1 (PD-1) can temporize the programmed death of immune cells, thus increasing immunotherapy targets against various tumors appear. However, the theory in BC remains to be proved [5]. As a consequence, exploring feasible immunosuppressive targets would provide prescient guidance when selecting treatment for patients with BC.
Notably, the associations among TMB, cancer immunotherapy and tumor-related prognosis have been persuasively demonstrated [6,7]. Rizvi NA, Hellmann MD et al argued that either in various oncogenes or tumor suppressor genes, the higher nonsynonymous mutation burden accumulated, the more potential neoantigens could emerge, which increases the immune response against advanced and metastatic tumor cells [8]. In terms of one certain carcinoma, the higher TMB is measured, the stronger response to immunotherapy could happen [9,10]. The Cancer Genome Atlas (TCGA) and International Cancer Genome Consortium (ICGC) were employed in order to probe into the relations among degrees of TMB, prognosis and immune response.

Materials And Methods
Data retrieval and TMB-statistical processing TCGA (https://cancergenome.nih.gov/) was applied to retrieve full transcriptome data of 433 BC samples, clinical information of 581 patients and somatic mutation data of 412 individuals. Somatic mutation data of extra 101 samples were outrightly acquired from ICGC (https://dcc.icgc.org/) databases. The landscape genome mutations of different genes from two databases were demonstrated as waterfall plots separately. Statistics were processed and analyzed with "oncoprint" R package.
Notably, only nonsynonymous missense mutation count for the calculation. According to the standard of mutation occurrence, the whole samples were divided into wild and mutation groups. The difference of TMB between two groups was performed as box plots, using two independent samples t-test, p<0.05 was supposed to be statistically signi cant. Statistics were processed and analyzed with "maftools" R package.

Prognostic factors analysis
Clinical data of 581 samples quali ed with available prognostic information and characteristic features were downloaded from TCGA database. For different genes, Kaplan-Meier curves (K-M curve) with logrank test were applied to characterize the difference of overall survival rates between wild cohort and mutation cohort. Hazard Ratios (HRs) were used as effective value, p<0.05 was supposed to be statistically signi cant. In order to identify speci c genes as potential risk factors of prognosis independent from other features of BC patients, random effects model was utilized to perform univariate Cox regression and multivariate Cox regression analysis. The clinical baseline variables included age, gender, tumor grade, TNM stage and TMB occurrence.

Gene functional exploration
To predict the pathways relevant genes most possibly participate in between the wild and the mutation groups of TCGA-BC clinical samples, Gene set enrichment analysis (GSEA) was applied to respectively assess the enrichment score of the wild and the mutation groups of gene sets, based on Kyoto Encyclopedia of Genes and Genomes (KEGG) database [11,12]. The function of analyzed genes thus could be re ected by mutative enrichment between the wild and the mutation on KEGG pathways. Statistics were processed and visualized with R package.

Immune cell in ltration assessment
The proportions of 22 types of immune cells in 433 samples were visualized as a bar plot in BC tumors, establishing those that play major roles in microenvironment of BC [13,14]. The 22 immune cell types include B cells, T cells, plasma cells, monocytes, macrophages, natural killer cells, eosinophils, neutrophils, etc. Correlations degree between each type of immune cell was visually presented with pairwise correlation analyses and evaluated by Spearman's rank correlation analyses, p<0.05 was considered to be statistically signi cant. The vioplot demonstrated the in ltrations of 22 types of immune cells in wild and mutation groups by using "vioplot" R package, of which Wilcoxon rank-sum test's p<0.05 was statistically signi cant criterion.

Description of TMB in BC
The research process was showed in Figure 1. Full transcriptome data belonging 433 individual BC tumor samples were retrieved from TCGA database with outrightly open acquisition access, among which 398 were included into study after processing. Then we analyzed the whole-exome sequencing data using the "maftools" package to visualize an exhaustive presentation of mutational information. The waterfall plots of mutational information from TCGA and ICGC were separately established in Figure 2A and Figure  2B. The color category beside the plot represents different genomic mutation types. Selectively, we focused on the missense mutations marked with green.
An intersection of 19 mutated genes were identi ed for research after selected simultaneously from the top 31 mutated genes of two databases, including TP53, PIK3CA, KDM6A, RB1 and so on ( Figure 2C).
Then the TMB calculation of 19 included genes was caught out to analyze the difference between wild type groups and mutant type groups. Statistically, TMBs of mutant type groups were higher than wild type groups for almost every set, except for FGFR3 and CREBBP ( Figure 2D), verifying that genomic mutation associates positively with TMB in BC.
Overall survival analysis and prognostic factors identi cation In order to gure out the in uence of mutated genes on prognosis, we performed survival analyses of the 19 genes mentioned above. 19 K-M curves of overall survival were drafted to compare survival rates and times showing as Figure 3A-S. We found that, among all of the 19 genes, wild group presents signi cant survival advantage to mutation group only in RB1 (log rank p=0.039), which meant the high mutation frequency on RB1 had a negative effect on prognosis of BC patients.
For further investigation of whether RB1 could be identi ed as an independent survival factor, randomsize effects models combining with clinical characteristics were separately performed to comprehensively analyze the HRs of multiple potential prognostic factors. As is shown in Figure 4A (Table 1). More evidently, according to the regression analyses, either TMB or mutation of RB1 associates negatively with prognosis. Moreover, this discovery is prognosticly consistent with the conclusion mentioned in Figure 2D that genomic mutation associates positively with TMB in BC.

Biological roles of RB1 in KEGG pathways
To further predict the potential biological roles of RB1 it plays in BC, the GSEA analysis was employed for KEGG pathways research. In accordance to the GSEA score, KEGG pathways like DNA replication, mismatch repair, cell cycle, oocyte meiosis and regulation of autophagy ( Figure 5A), which mostly relate to cell proliferation and cell differentiation, enrich on the mutation region of RB1 band. This nding signi es that mutation types of RB1 might be more inclined to take part in cell proliferation and development of BC.

Patterns of immune cell in ltration in BC tumor
The population proportions of 22 types of immune cells in 398 samples were visualized as a bar plot in BC tumors, and these types include B cells, T cells, plasma cells, monocytes, macrophages, natural killer cells, eosinophils, neutrophils, etc. As is presented in the bar plot, 398 samples from TCGA database display the respective distributions of 22 immune cell in ltration( Figure 6A). Generally, the proportions of every type of immune cell are neither even nor stable in different BC samples.
As for the correlation patterns, relevance between each of the 22 immune cells is multifarious, or even loose in general. While it is worth noting that CD4 memory activated T cell shows signi cantly positive correlations with CD8 T cell, NK resting cell, Macrophage 1 (M1) and follicular helper T cell. Whereas, both of CD8 T cell and CD4 memory activated T cell associate respectively negatively with CD4 memory resting T cell and Macrophage 0 (M0) ( Figure 6B). This nding indicates that immune cells cooperating together are inclined to in ltrate in tumors. In contrast, those with conversed function tend to exist mutually opposite to each other in the tumor microenvironment of BC.
In terms of RB1, fractions of 22 immune cell types were established in the vioplot. Wilcoxon rank-sum test reveals that the amount of regulatory T cells (Tregs) in ltrating in wild type group was signi cantly higher than in mutant type group (Wilcoxon rank-sum test, p = 0.031, Figure 6C). One probable inference is that the occurrence of genomic mutation on RB1 decreases the antigens Tregs used to be able to recognize, which suppresses its speci c immune response to the BC tumor while strengthen the tumor immune escaping.

Discussion
BC is the most prevalent malignant tumor of urinary system, taking major part in yearly mortality related to cancers. While cancer immunotherapy has become increasingly notable for its effectiveness and accuracy in treating advanced and metastatic carcinomas, immunotherapeutic strategies aimed at BC remain insu cient [14][15][16]. In 2017, Goodman AM, et al pointed out that tumor mutational burden can function as an independent predictor of response to immunotherapy in various cancers [10]. Subsequently, Fan CN, et al identi ed three speci c long-noncoding RNAs as correlated factors with TMB and clinical outcomes of triple-negative breast cancer in 2019, revealing that long-noncoding RNAs can be the predictor of TMB and prognosis [17]. Till 2020, Wang CL, et al complementary of microRNA as a marker of TMB and immune cell in ltration in lung adenocarcinoma (LUAD), which led to further re nement of the association between non-coding genes and TMB-mediated immune cell in ltration [18]. Moreover, the communication among coding RNAs expression, TMB level, prognosis and immune cell in ltration was currently described by Wang FR, et al in breast cancer [19]. However, researches of gene mutations as a TMB-related factor are rarely reported. Therefore, coming across this imperfection when constructing prognostic or immunotherapeutic models, we carried out this comprehensive analysis to explore relations among mutation, TMB and immune cell in ltration in BC tumor.
We found 19 genes equipped with high level of missense mutation in BC by intersecting TCGA (n = 412) and ICGC (n = 101) databases, of which 17 showed signi cant difference between wild type group and mutation type group. Based on the presence and absence of mutation, levels of TMB respectively appear high and low. The 17 signi cant genes included TP53, PIC3CA, KDM6A, TTN, SYNE1, SYNE2, ERCC2, ARID1A, MUC16, ERBB2, XIRP2, EP300, RB1, AHNAK2, CSMD3, HMCN1 and MACF1. Noticeably, TP53 and RB1 were widely known as important anti-oncogenes in diverse range of cancers, while KDM2A was acknowledged as a core gene that promotes tumor development by participating in histone methylation [20][21][22][23]. In a word, the variety of functions of the relevant local TMB-regulated genes suggests that changes in these TMB levels interfere with complex biological signal pathways.
In the further study of the in uence of mutation on prognosis, mutation occurring to RB1 distinctively contributed to poor clinical outcome. As is referred to before, RB1 represses the proliferation and malignant progression through negatively regulating the cell cycle in BC tumor. Consequently, it can be generally inferred that the high level of TMB occurring to RB1 sequence likely depresses its original function of cell cycle delaying. In other words, the increasing TMB in BC accelerates the proliferation and development of tumor cells and directly leads to poor survival rates of patients, which is in agreement with the argument of Luo C, et al [24]. The subsequent analysis of the KEGG pathways reinforces the ndings above. RB1 gene was more likely to participate in signal pathways related to cell proliferation and cell differentiation like DNA replication, oocyte meiosis and regulation of autophagy. Thus, mutation might alter its original roles in pathways.
In the tumor microenvironment, in ltration of in ammatory cells is one of typical features, which affects tumor proliferation and development to a large extent [25,26]. And mutual interaction of in ammatory cytokines and immune cells is a major research eld. Previous clinical studies revealed that the immune cell in ltration gives rise to different denouements of plenty of carcinomas [27][28][29]. Thus, the relationship between TMB and immune cell in ltration was taken into our consideration.  [30]. Interestingly, Tregs interact apparently negatively with pro-in ammatory cells, such as CD4 memory activated T cell, CD8 T cell, NK resting cell and (M1). Furthermore, Treg cell was found as a particular negatively regulatory target of RB1. Regulatory T cell has been acknowledged to have regulatory capability by suppressing the in ltration of other T cells. The speci c trait of Tregs allows it to depress the chronic immune responses against viruses, tumors and self-antigens [31]. According to our present research, mutations occurring in RB1 signi cantly reduce the immune in ltration of Tregs, implying that high level of TMB is capable of stimulating the malignant progression. From the perspective of the mechanism, chances are RB1 can inhibit the Tregs-associated in ammation to play its role in cancer promotion by mutations. This approach to undermining Treg-mediated chronic immune responses against tumor can be further deemed as a partial explanation to tumor immune escape of BC.
In the present research, according to bioinformatics analyses, we put forward that the high TMB of RB1 in BC can completely turnover the original antineoplastic in uence the wild type had before. Moreover, quanti cation of the TMB in BC was suggested to be a potential tool for clinical outcome prediction due to our evidence. Inevitably, there are some de ciencies in our study, mostly because of the retrospective data extraction from databases. Admittedly, prospective validations based on experimental results will make our opinions more persuasive. Biases would also generate during the progress of genes selection. For example, signi cance was observed on SYNE 2 in TMB calculation and survival analysis, but its association with cancers is still inconclusive.
To conclude, the high level of TMB in RB1 could be apparently measured in BC, and it could serve as an independent predictor of poor prognosis in BC patients. Meanwhile the more occurrences of mutations in RB1 are inclined to protect BCelimination, by chance of delaying the apoptosis in cell cycles and weaken the antitumor effects of Treg cell in ltration in tumor microenvironment, in turn resulting in further malignant development.

Conclusion
In conclusion, RB1 can be identi ed as an independent prognostic predictor, and there is a chance for contribution to poor overall survival as the mutation occurs. What's more, mutation of RB1 also functions as a biomarker that represses the in ltration of Tregs, increasing the incidence of tumor immune escaping in BC.   Figure 1 The ow chart of the whole research process.    GSEA of KEGG pathways RB1 takes part in after mutation.