Racial differences in genomic features and new trans-ancestry prognosis subtypes in bladder cancer

The incidence and survival of bladder cancer vary greatly among different populations but 28 the influence of the associated molecular features and evolutionary processes on its clinical 29 treatment and prognostication remains unknown. Here, we analyze the genomic architectures 30 of over 500 bladder cancer patients from Asian/Black/White populations. We identify novel 31 association between AHNAK mutations and APOBEC-a mutational signature whose activities 32 vary substantially across populations. All significantly mutated genes but only half of arm- 33 level somatic copy number alterations (SCNAs) are enriched with clonal events, indicating 34 large-scale SCNAs as rich sources fostering bladder cancer clonal diversities. The prevalence 35 of TP53 and ATM clonal mutations as well as the associated burden of SCNAs is 36 significantly higher in Whites/Blacks than in Asians. We identify a trans-ancestry prognostic 37 subtype of bladder cancer characterized by: enrichment of non-muscle-invasive patients and 38 muscle-invasive patients with good prognosis, increased CREBBP/FGFR3/HRAS/NFE2L2 39 mutations, decreased intra-tumor heterogeneity and genome instability and activated tumor 40 microenvironment. 41 in higher levels of intra-tumor and higher rate of


INTRODUCTION 43 44
Bladder cancer is the sixth most common type of adult cancer with almost 430, 000 new 45 cases diagnosed per year 1 . Bladder cancer is classified into non-muscle invasive bladder 46 cancer (NMIBC) and muscle-invasive bladder cancer (MIBC) according to their 47 histopathological features and clinical behaviors. NMIBC, the major subtype of bladder 48 cancer that accounts for 70% to 80% of all newly diagnosed cases, has better overall survival 49 but high recurrence rate of 31-78% while the rest of patients with muscle-invasive or 50 metastatic diseases exhibit quite low survival rate 2,3 . The genetic or molecular basis 51 underlying such great disparities in clinical outcomes between NMIBC and MIBC remains 52 largely unclear. prognostic biomarkers for bladder cancer had been developed utilizing their clonal or sub-75 two signatures in Black and four signatures in White ( Fig. 1a and Supplementary Fig. 1). 126 Through clustering analysis of the cosine similarities and Pearson correlation values between 127 these signatures and the 30 Sanger signatures with reported roles in mutagenesis 14-16 , the five 128 mutational signatures identified in our cohorts could be grouped into categories resembling  Fig. 1b and Supplementary Fig. 3). These observations might suggest the 144 presence of population variations in mutational processes operating at different evolution 145 stages during bladder cancer development. To further assess the potential clinical association 146 of these mutational signatures, we performed unsupervised clustering of the patients in 147 Cohorts 1 and 2 according to their mutational signature activities and identified four 148 mutational signature clusters (MSig1 to MSig4, Fig. 1c). Different clusters of patients 149 showed varied activities of different mutational signatures (Fig. 1d) and were associated with 150 distinct clinical outcomes (Fig. 1e, P = 0.002). Cluster MSig2 showed enrichment in Asian 151 patients (Fig. 1c) and had high MMR-like signature activity (P < 0.001 & FDR< 0.1, Fig.  152 1d). Prognostic analysis of Cohort 1 revealed that MIBC samples in cluster MSig2 tended to 153 have the shortest overall survival (Fig. 1e) and this observation had to be further investigated 154 in future studies as the number of samples with survival information available in cluster 155 MSig2 was quite small in our study. Patients in the MSig3 and MSig4 clusters had high 156 levels of APOBEC and ERCC2 & AGE signature mutagenesis (Fig. 1d), respectively. 157 Multivariate Cox regression analysis of Cohort 1 patients showed that those patients in the 158 MSig3 and MSig4 clusters had better survival than those belonging to MSig1 cluster (Fig.  159   1f). 160

161
The Asian cohort in our study was consisted of patients from both Europe/USA areas (Cohort 162 1: MIBC) and Asia areas (Cohort 2: MIBC and NMIBC) and we compared the relative 163 contribution of the mutational signatures among these three Asian patient groups and found 164 that the activities of APOBEC-a and -b signatures were relatively higher in MIBC patients 165 from Cohort 2 than the other two groups of Asian patients. Both MIBC and NMIBC patients 166 from Cohort 2 showed significantly high levels of MMR-like signature activities compared 167 with MIBC samples from Cohort1 (Supplementary Fig. 3). Nevertheless, more studies and 168 samples were needed to validate the observed racial differences in the activity of MMR-like 169 signature and we could not rule out the potential influence of study batch effect on this 170 observation due to the following reasons: (i) the MMR-like signature with minor exposures 171 could only be observed in Cohort 2 generally; (ii) most of the mutations contributed to the 172 MMR-like signature were subclonal; (iii) the MMR-like signature had a relatively low level 173 of similarity to any known COSMIC signatures; and (iv) little previous evidence suggested 174 the presence of this signature in bladder cancer (Supplementary Fig. 4). 175

Gene mutations potentially associated with signature activities 177
To explore whether there were any genes whose mutations were associated with the 178 mutational signature activities, we performed signature enrichment analysis for those genes 179 mutated at elevated frequencies with permutation tests which could control the inflation that 180 an increase in the overall mutation burden in each sample or in the total number of mutations 181 in each gene usually correlated with increased signature activity. To determine the 182 significance level of associations between the elevated mutation frequency of each gene and 183 the increased activity of each signature, the mutation burden attributed to a given signature 184 from tumors that harbored non-silent mutations in the gene was compared with tumors that 185 did not by combining the two study cohorts (Online Methods). 186 187 Consistent with previous reports 17 , we showed that ERCC2 mutation state and smoking were 188 significantly correlated with ERCC2 & AGE signature activity (FDR < 0.1, P < 0.001, Fig  189   2a, b). Interestingly, we also found that ERCC2 mutation state was correlated with the 190 activity of APOBEC-b signature (FDR < 0.1), with an increase of 75 APOBEC-b mutations 191 in ERCC2-mutant patients (Fig. 2b). Notably, we identified that the APOBEC-a signature 192 activity, an indicator of good clinical outcome in bladder cancer 4 , was associated with 193 AHNAK mutations which appeared to act as a tumor suppressor through potentiating the 194 TGFβ signalling in a previous study at the protein level 18 (Fig. 2a). The median numbers of 195 APOBEC-a mutations in patients with and without non-silent AHNAK mutations were 235 196 and 62, respectively ( Fig. 2b, P < 0.001). Correlation between HRAS mutation states and the 197 MMR-like signature activity was also noticed in our study (Fig. 2a, b). 198

199
We further validated that whether racial differences in mutation frequencies of several key 200 genes or pathways may partially account for the racial differences in mutational signature 201 activities among Asian/Black/White populations by combining Cohorts 1 and 2. The 202 frequency of AHNAK mutations was the highest in the White patients who also showed the 203 highest APOBEC-a signature activity (9.6% in White vs 4.3% in Black & 6.4% in Asian) 204 while the frequency of HRAS mutations was the highest in the Asian patients (14.3% in Asian 205 vs 0.0% in Black & 3.4% in White) (Fig. 2c). At the pathway level, genes involved in the 206 DNA polymerase pathway and other conserved DNA damage response genes also had the 207 highest mutation frequencies in White population (Supplementary Fig. 5a). 208

209
To determine whether low-frequency mutations in genes involved in the DNA damage 210 responses were linked to mutational signature activities, we performed permutation tests for 211 the DNA damage response pathways from curated databases (Online Methods). We observed 212 that mutations of the nucleotide excision repair (NER) pathway which includes the ERCC2 213 gene were significantly associated with the activities of both ERCC2 & AGE and APOBEC-b 214 signatures (Supplementary Fig. 5b). Mutations in the DNA polymerase pathway were 215 identified to be associated with the APOBEC-a signature activity. Additionally, mutations of 216 a group of genes defined as other conserved DNA damage response genes were significantly 217 associated with the mutagenesis of three signatures, APOBEC-a, APOBEC-b and ERCC2 & 218 AGE (Supplementary Fig. 5b, c). 219 220 We next compared the expression levels of AHNAK, ERCC2 and HRAS among the three 221 racial populations in Cohort 1. We observed that the expression levels of both AHNAK and 222 ERCC2 were significantly lower in Asian than those in White population (P < 0.05 and P < 223 0.001, respectively) ( Supplementary Fig. 6), which provided further evidence supporting 224 that racial differences in the activities of AHNAK and ERCC2 might contribute to the 225 observed racial differences in the prevalence of APOBEC-a and ERCC2 & AGE signatures 226 in bladder cancer. However, HRAS expression showed no significant difference among 227 different race groups (Supplementary Fig. 6). 228 229

Racial disparities in cancer gene mutations and SCNAs 230
To depict the racial disparities in molecular features among different racial groups, we first 231 identified the significantly mutated genes (SMGs) by combining all samples from the two 232 study cohorts and then compared the prevalence of mutations in SMGs among different race 233 groups. A total of 57 SMGs were identified by MutSig2CV (q < 0.1) (Supplementary Data 234 2) and 26 of them were cancer census genes that were expressed in over 75% of TCGA 235 patients (supported by at least three read counts) and were kept for downstream analysis 19 . 236 Asian patients from Cohorts 1and 2 were consistently higher than Whites. Compared with 245 Asian patients, clonal mutations in ATM and ARID1A showed higher prevalence in Blacks 246 and Whites, respectively. The frequency of NRAS mutations was higher in Blacks than 247 Whites and three genes, including FGFR3, ELF3 and TSC1, showed higher frequencies of 248 clonal mutations in Whites than in Asian patients from Cohort 2. Interestingly, FGFR3 clonal 249 mutations were observed to occur more frequently not only in Asian patients in Cohort 1 than 250 in Cohort 2. FGFR3 mutations are the early driver events during NMIBC evolution and we 251 also observed that FGFR3 clonal mutations occurred more frequently in NMIBC than in 252 MIBC in Cohort 2 patients 23 (Fig. 3a, b). 253

254
We next compared the overall burden of SCNAs at the chromosome/arm/focal levels across 255 different racial groups (Online Methods). NMIBC patients generally had a lower burden of 256 SCNAs at the chromosome/arm levels than the MIBC patients. The overall arm-or 257 chromosome-level SCNA scores in Asian MIBC patients were lower than those in White or 258 Black MIBC patients (Fig. 3c). This observation was consistent with the finding of low 259 prevalence of clonal mutations in the TP53 and ATM gatekeeper genes governing DNA repair 260 in Asian MIBC patients. We further categorized each chromosome/arm level SCNA into 261 clonal or subclonal event and found that the majority of arm-level deletions either in clonal or 262 subclonal states showed substantial inter-population variations in their prevalence. A large 263 number of chromosomes/arms showed higher prevalence of clonal deletions in White/Black 264 MIBC patients than in the Asian MIBC patients. We also observed a large number of 265 chromosomes/arms showing higher prevalence of subclonal SCNAs in Whites than in Asian 266 patients from Cohort 2 (Fig. 3d). Moreover, we observed that only quite a few number of 267 clonal arm-level SCNAs, including del(3p) and amp(20), had higher prevalence but a number 268 of subclonal SCNAs showed significantly lower prevalence in Asian MIBC patients from 269 Cohort 2 than in Asian patients from Cohort 1 (Fig. 3e). 270 271 It has been reported that bladder cancer patients with different genetic backgrounds may be 272 presented with different disease stages at the time of diagnosis 10 . In our study cohorts, we 273 also observed that MIBC patients of Asian origins had higher prevalence of low-stage disease 274 while Blacks and Whites had higher prevalence of high-stage disease (Supplementary Fig.  275   7). To correct for the potential influence of heterogeneities in disease stages on the racial 276 disparities in genomic features of MIBC, we further analyzed the racial differences in the 277 mutation landscapes and SCNA profiles in low-staged (staged I & II) MIBC patients due to 278 the lack of sufficient numbers of high-staged patients in the Asian cohorts. Again, the 279 prevalence of clonal mutations in TP53 and ATM was significantly higher in Whites and 280 Blacks than in Asians, respectively, among MIBC patients with low-stage disease 281 ( Supplementary Fig. 8). The prevalence of clonal mutations in NFE2L2 was the highest in 282 low-staged Black patients and the fractions of White patients harboring clonal mutations in 283 FGFR3, TSC1 and CDKN2A were significantly higher than those of Asian patients from 284 Cohort 2. Clonal mutations in HRAS also showed significantly higher prevalence in Asians 285 than in Whites among the low-staged patients. We also confirmed that the prevalence of 286 FGFR3 clonal mutations was significantly higher in Asians from Cohort 1 than in Asians 287 from Cohort 2 among the low-stage MIBC patients. The trends of higher prevalence of clonal 288 deletions and subclonal SCNAs in Whites or Blacks were also observed among patients with 289 low stage MIBC. Differences in the frequencies of clonal deletions of several chromosome 290 arms including del(3p) and subclonal SCNAs affecting a number of other arms were also 291 observed between Asians from Cohorts 1 and 2 among the low-staged patients. 292 293

Genomic architectures and subtypes of bladder cancer 294
Our study included a large number of bladder cancer patients with different genetic 295 backgrounds and we then tried to perform molecular subtyping of bladder cancer by 296 integrating multiple somatic events according to their CCFs in all samples by combining 297 Cohorts 1 and 2. We inferred the genomic architectures of bladder cancer by focusing on the 298 26 SMGs with reported roles in the Cancer Gene Census database and the 34 frequent arm-299 level SCNAs altered in at least 30% of patients as defined above. To explore the contribution 300 of these SMGs and frequent arm-level SCNAs during bladder cancer evolution, we estimated 301 the CCF of each SMG and SCNA in each sample. Overall, the median CCFs of SMGs were 302 higher than those of SCNAs by combining the two study cohorts (Fig. 4a). All the 26 SMGs 303 were found to be enriched with clonal non-silent mutations (FDR < 0.1) while only about 304 50% (17/34) and 15% (5/34) of the frequent arm-level SCNAs were significantly enriched 305 with clonal and sub-clonal events (FDR < 0.1), respectively (Fig. 4a). These data suggested 306 that heterogeneities in the genomic architectures of bladder cancer were more likely to arise 307 from diversity in the clonity of arm-level SCNAs.  Fig. 9). Thus, we further evaluated the 317 prognostic value of the SMGs and frequent SCNA events according to their clonal or sub-318 clonal states by ignoring the racial origins of the patients to increase the sample sizes. We 319 combined the different racial groups in Cohort 1 together in Kaplan-Meier analysis and 320 multivariate Cox regression analysis. In total, we identified only four individual events with 321 potential prognostic value in Cohort 1 ( We next tried to stratify whether there were any subtypes of bladder cancer whose overall 327 genetic architectures had great influence on their clinical outcomes. We performed NMF 328 analysis of the CCFs of the 26 SMGs and 34 arm-level SCNAs by combining the two study 329 cohorts. In total, we identified two molecular subgroups (termed as clusters A and B), whose 330 clinical outcomes were divergent. Because long-term survival information was not available 331 for the patients from Cohort 2, we performed survival analysis for all the MIBC samples in 332 the TCGA dataset. TCGA MIBC patients in cluster A had worse survival rate than those in 333 cluster B (Fig. 4b) (HR = 1.64, 95% CI: 1.14 ~ 2.36; P = 0.007). Multivariate analysis further 334 proved that the genomic architectures could predict the clinical outcomes of bladder cancer 335 patients independently by adjusting for age, gender, TNM stage and mutational signature 336 (Fig. 4c). Further comparison of these two genomic subtypes with the molecular subtypes 337 defined by TCGA at different expression levels showed that genomic cluster B tended to be 338 clustered with the molecular subtypes mRNA_Luminal_papillary, lncRNA_3, miRNA_3 and 339 RPPA_1 in Cohort 1 patients 20 (Supplementary Fig. 10). It had been shown previously that To understand which factors contributed to the distinct prognosis of bladder cancer subtypes 349 defined by genomic architectures, we compared the prevalence of 60 potential driver events 350 (26 SMGs and 34 arm-level SCNAs) between clusters A and B. We showed that four SMGs 351 (CREBBP, FGFR3, HRAS and NFE2L2) but none of the frequent arm level SCNAs had 352 significantly higher prevalence of genetic aberrations in cluster B patients who had better 353 prognosis than in cluster A patients (FDR < 0.1, Fig. 4a and Supplementary Data 4). On the 354 other hand, three SMGs (TP53, RB1 and PRDM2) and about 91% of the frequent arm level 355 SCNAs were enriched with alterations in cluster A patients who showed poor prognosis 356 Fig. 4a and Supplementary Data 4). Notably, of the SMGs mutated more 357 frequently in cluster B than in cluster A, three showed race-specific enrichment of clonal 358 mutations among the low-staged patients with FGFR3, HRAS and NFE2L2 mutations being 359 enriched in Whites/Asians from Cohort 2, Asians from the two study cohorts and Blacks, 360 respectively (Fig. 4a) intermediate/late evolution stages, respectively (Fig. 4a). These data further supported that 370 the clonal diversities of both clusters A and B patients were largely contributed by the large-371 scale SCNAs acquired at late evolution stages. Several SCNAs, including del(9), 372 del(17p),del(8p) and del(11p), were early clonal events shared by patients in both clusters A 373 and B while arm-level SCNAs, including del(4q), del(15q), del(10q),del(11q) and del(12), 374 were acquired at the early and late stages in clusters A and B, respectively (Fig. 4a). 375

376
To explore whether there were any other differences in the overall genomic landscapes 377 between clusters A and B patients, we firstly compared their non-silent somatic mutation 378 burdens and found no significance difference. The extent of intra-tumor heterogeneities as 379 measured by MATH values and Shannon entropy of CCFs in cluster A patients was 380 significantly higher than that of cluster B patients (P< 0.001, Fig. 4a and Supplementary 381 Fig. 11). We next estimated the overall burden of SCNAs at the chromosome/arm/focal 382 levels for each patient and showed that cluster A patients had significant higher burden of 383 large-scale SCNAs (arm & chromosome levels) but not focal SCNAs (P < 0.001, Fig. 4d). 384 Whole-genome doubling (WGD) events were prevalent in cluster A and patients with Asian 385 origin were enriched in cluster B (both P < 0.001). To be more specific, although our cohorts 386 included only a small number of NMIBC patients from Cohort 2, cluster B were enriched 387 with Cohort 2 NMIBC patients (51.0% in cluster B vs. 22.9% in cluster A, P < 0.005) and 388 Cohort 1 MIBC patients who were less likely to develop metastatic disease (52.6% in cluster 389 B vs. 33.3% in cluster A, P < 0.001, Fig. 4a). Further analysis of the overall survival for 390 Cohort 1 MIBC patients showed that patients in cluster A had worse survival rate than those 391 in cluster B after adjusting for WGD, racial origins as well as other clinical and genomic 392 parameters (Fig. 4c). These observations suggested that it was feasible to distinguish the 393 heterogeneous MIBC patients into two subtypes with distinct clinical outcomes based on 394 their inherited genomic architectures. 395

396
To further study whether high intra-tumor heterogeneity and chromosomal instability of 397 cluster A patients had influence on the tumor microenvironment, we compared the relative 398 expression levels of the immune signature genes among two genomic prognostic clusters 399 using ssGSEA (Fig. 4e). In Cohort 1, the infiltrating levels of Th17 cells and CD8+ T cells 400 were observed to be significantly higher in cluster B patients (P < 0.001 and FDR < 0.1, Fig.  401   4e). Previous studies proved that Th17 cells as subsets of T helper lymphocytes can mediate 402 the antitumor immune responses through interacting with effector CD8+ T cells 24 . Cluster A 403 patients had significantly higher levels of infiltrating Th2 cells and Macrophages (Both P < 404 0.001 and FDR < 0.1, Fig. 4e), which had been proved to facilitate tumor growth and 405 metastasis leading to the reduced survival in the cancer 25-27 . These observations on 406 differences in tumor microenvironment between clusters A and B were also validated 407 independently in the Chinese cohort (Fig. 4e). and White populations based on their racial origins (Supplementary Table 2 and 3). 522 523

Identification of SNVs and SCNAs 524
Bam files for the tumor and matched normal samples were processed as input. Somatic 525 mutations were called by using Mutect2 software 12 and further filtered with a read depth of at 526 least 10x in the germline and tumor samples, a maximum of two variant supporting reads in 527 the germline, a minimum tumor variant allele frequency of 10% and a maximum germline 528 variant allele frequency of 2%. Allelic SCNAs were called with the ReCapSeg software 529 following the Broad's GATK documentation 530 (http://gatkforums.broadinstitute.org/categories/recapseg-documentation). 531 532

Estimation of intratumor genetic heterogeneity 533
To estimate the level of intratumor genetic heterogeneity in bladder cancer, we calculated 534 each tumor's MATH value from the median absolute deviation (MAD) and the median of its 535 mutant-allele fractions at somatically mutated loci: MATH=100*MAD/median 36 . 536 537

RNA-seq analysis and estimation of immune cell score 538
Bam files for paired samples were processed to quantify the transcript abundance levels by 539 using kallisto 37 . We then used tximport 38 R package to convert transcript level value into 540 gene level. We used signature genes from previous study 39 to infer the infiltration levels of and White) and 30 COSMIC signatures was performed using the standard hierarchial 567 clustering R package with a distance measurement of 'cosine' similarity or 'pearson' 568 correlation. The MSig clustering analysis was performed using a standard hierarchical 569 clustering in R, with a 'euclidian' distance for the signature activity matrix and a 'ward.D' 570 linkage tree. By using the Elbow method to look at the total within-cluster sum of square 571 (wss) as a function of the number of clusters, the number of MSig clusters (K = 4) was 572 chosen so that adding another cluster does not improve much better the total wss. 573

Signature enrichment analysis 575
To search for genes whose mutation statuses were associated with the activity of a specific 576 signature, we applied a permutation statistical test 17 to compare, for each gene, the signature 577 activities between samples with and without mutations in the gene. Firstly, to remove the 578 inflation in the number of mutations in each gene associated with the elevated background 579 mutation rates, we randomly produced a total of 10 4 permutated mutation matrix in which the 580 total counts of gene-specific and sample-specific mutations were kept the same as those 581 observed in our patient cohorts, following an approach described by Strona et al 49 . We used 582 the one-tailed Wilcoxon rank-sum P value to compare the signature activity between mutant 583 and wild type samples of a specific gene. To get the final P value for a given gene, we 584 calculated the fraction of permutated matrix with a test P value equal to or more extreme than 585 the observed matrix (Sum of Pobserved ≥ Prandom times/the total number of permutations). To 586 increase the computational efficiency, we analyzed only 174 genes with a non-silent mutation 587 frequency of 5% or above, 26 SMGs, the APOBEC family member gene sets and the human 588 DNA repair gene sets that had been reported to be associated with some mutational signature 589 activities (downloaded from https://www.mdanderson.org/documents/Labs/Wood-590 Laboratory/human-dna-repair-genes.html). We corrected for multiple-hypothesis testing 591 using the Benjamini-Hochberg procedure and used FDR q < 0.1 as the significance threshold. We calculated the SCNA scores according to a previous study with minor modifications 50 . 596 Our calculation of SCNA was based on the integer allelic copy number calls generated by 597 ABSOLUTE 51 which took the tumor purities into consideration. The absolute copy number 598 of each segment was determined and a chromosomal arm event would be called out if the 599 cumulative length of the deletions (genomic segments with minor alleles of 0 copy) or 600 amplifications (genomic segments with major alleles of 2 or more copies) on this arm was 601 greater than 50% of the chromosome arm. Chromosome level events were distinguished from 602 arm level events when both arms of a chromosome had the same copy number changes (in 603 sign). After excluding the arm-and chromosome-level copy number changes, the remaining 604 segments with SCNAs were defined as focal changes. Based on the copy number change sign 605 associated with each segment x of a chromosome, arm and focal event t, the assigned score S 606 were defined as: 607 Wilcox.test and chisq.test to generate the empirical P values, respectively. P values were 634 adjusted for multiple hypothesis tests using the R function p.adjust with the "fdr" option 635 636

Survival analysis 637
Asian bladder cancer patients from Asia cohort were excluded in survival analysis due to the 638 lack of some key clinical information. Chi-squre test statistics in Kaplan-Meier curves were 639 computed using log-rank tests. P values were also calculated from multivariate Cox 640 proportional-harzards regression models using the R package "survival".