Comprehensive Molecular Characterization of Urothelial Bladder Carcinoma

Demographic, clinical and pathological data

Samples, (from 19 tissue source sites), consisted of 131 chemotherapy-naïve, muscle-invasive, high-grade urothelial tumours (T2-T4a, Nx, Mx), as well as peripheral blood (n=118) and/or tumour-adjacent, histologically normal-appearing bladder tissue (n=23). Cases were retained only if they met the following criteria: tumour nuclei constituted ≥ 60% of all nuclei, tumour necrosis was ≤ 20% of the specimen, and variant histologies (squamous or small cell) were ≤50% (Supplementary Text S1). Clinical and demographic characteristics are described in Supplementary Data File S1.1. Five expert genitourinary pathologists re-reviewed all of the cases for multiple parameters, including the extent of variant histology (Supplementary Text S1 and Fig. S1.1a).

Somatic DNA alterations

The tumours displayed a large number of DNA alterations, slightly fewer than in lung cancer and melanoma, but more than in other adult malignancies studied by TCGA (Fig. 1)⁹. On average, there were 302 exonic mutations, 204 segmental alterations in genomic copy number and 22 genomic rearrangements per sample. We analyzed somatic copy number alterations (SCNAs) using both SNP 6.0 arrays and low-pass whole genome sequencing; the two were strongly concordant (Supplementary Methods S6.1 and Supplementary Fig. S6.1). There were 22 significant arm-level copy number changes (Supplementary Data File S6.1.1), and GISTIC (Supplementary Methods S6.2) identified 27 amplified and 30 deleted recurrent focal SCNAs (Supplementary Data Files S6.2.1 and S6.3.1). Focal amplifications involved genes previously reported to be altered in bladder cancer (Fig. 1c, Supplementary Fig. S6.2.1) and some not previously implicated. The latter included PVRL4, BCL2L1 and ZNF703. The most common recurrent focal deletion, seen in 47% of samples, contained CDKN2A (9p21.3) and correlated with reduced expression (Fig. 1 and Supplementary Fig. S2.7). Other focal deletions containing <10 genes appeared to target PDE4D, RB1, FHIT, CREBBP, IKZF2, FOXQ1, FAM190A, LRP1B and WWOX.

Whole-exome sequencing of 130 tumours and matched normal samples targeted 186,260 exons in 18,091 genes (mean coverage 100×, with 82% of target bases covered >30×). MuTect¹⁰ identified 39,312 somatic mutations (including 38,012 point mutations and 1,138 indels), yielding mean and median somatic mutation rates of 7.7 and 5.5 per megabase (Mb), respectively (Fig. 1a, Supplementary Table S2.1.1). Thirty-two genes showed statistically significant levels of recurrent somatic mutation (Fig. 1b, Supplementary Table S2.1.2) by analysis using MutSig 1.5¹¹ (Supplementary Methods S2.2). Three other genes identified by MutSig were not considered further because of low or undetectable expression (Supplementary Fig. S2.1.1). A similar analysis considering only mutations in the COSMIC database² identified three more significantly mutated genes (SMGs): ERBB2, ATM and CTNNB1 (Supplementary Table S2.1.3). We validated the mutation findings in three ways: targeted re-sequencing of all SMG mutations, comparison with RNA-Seq data for 123 samples and comparison with whole genome sequence data for 18 samples. Overall, the validation rate was > 99% by a combination of the methods (Supplemental Methods S2.4.)

Nearly half (49%) of the samples had TP53 mutations (Fig. 1b), which were mutually exclusive in their relationship with amplification (9%) and overexpression (29%) of MDM2; hence, TP53 function was inactivated in 76% of samples. Most RB1 mutations were inactivating, were associated with significantly reduced mRNA level (Supplementary Fig. S2.7) and were mutually exclusive with CDKN2A deletions (Supplementary Fig. S2.8 and Table S2.8.1). FGFR3 mutations (12%) typically affected known kinase-activating sites. PIK3CA mutations were relatively common (20%), clustering in the helical domain near E545 (Supplementary Fig. S2.4). Most TSC1 mutations (8%) were truncating, and six were homozygous (allele fraction > 0.5).

Many of the 32 genes identified in Fig. 1b have not previously been reported as statistically significantly mutated in bladder cancer: MLL2 (27%), CDKN1A* (14%), ERCC2* (12%), STAG2 (11%), RXRA* (9%), ELF3* (8%), NFE2L2 (8%), KLF5* (8%), TXNIP (7%), FOXQ1* (5%), RHOB* (5%), FOXA1 (5%), PAIP1* (5%), BTG2* (5%), ZFP36L1 (5%), RHOA (4%) and CCND3 (4%). The nine genes marked with asterisks have not been reported as SMGs in any other TCGA cancer type or reported in another study as mutated at >3% frequency². CDKN1A (p21^CIP1), a cyclin-dependent kinase inhibitor¹², had predominantly null or truncating mutations, implying loss of function. Fifteen of 16 mutations in ERCC2, a nucleotide excision repair gene¹³, were deleterious missense mutations, suggesting dominant negative effects. ERCC2-mutant tumours also had significantly fewer C>G mutations than did ERCC2-wild type tumours (Supplementary Figs. S2.3.1 and S2.3.2), and they trended toward higher overall mutation rate (Supplementary Figure S2.12). Seven of 12 mutations in RXRA (retinoid × nuclear receptor alpha)¹⁴ occurred at the same amino acid (five S427F; two S427Y) in the ligand-binding domain. Those seven tumours showed increased expression of genes involved in adipogenesis and lipid metabolism (Supplementary Fig. S2.6 and Data Files S2.6.1- S2.6.3), suggesting that the mutations cause constitutive activation.

Eleven tumours (8%) had deleterious missense mutations in the Neh2 domain of NFE2L2, a transcription factor that regulates the anti-oxidant program in response to oxidative stress¹⁵. Those tumours showed dramatically increased expression of genes involved in genotoxic metabolism and the reactive oxygen species (ROS) response (Supplementary Figs. S2.5.1-S2.5.3 and Date File S2.5.2). Furthermore, nine samples had mutations in redox regulator TXNIP¹⁶ (5 of them inactivating) and were mutually exclusive of samples with NFE2L2 mutations, providing another mechanism for dysregulation of redox metabolism. Predominant inactivating mutations were seen in STAG2, an X-linked cohesin complex component required for separation of sister chromatids during cell division¹⁷ (Supplementary Fig. S2.4).

Unsupervised clustering by non-negative matrix factorization of mutations and focal SCNAs in 125 samples identified three distinct groups (Fig. 1a, Supplementary Fig. S2.1.2). Group A (red), labeled as ‘focally-amplified’, is highly enriched in focal SCNAs in several genes, as well as mutations in MLL2 (Fig. 1; Supplementary Tables S2.1.4 and S2.1.5). Group B (blue), labeled as ‘papillary CDKN2A-deficient FGFR3-mutant‘, is enriched in papillary histology. Nearly all Group B samples show loss of CDKN2A, and the majority have one or more alterations in FGFR3. Group C (green), labeled as ‘TP53/cell-cycle-mutant’, shows TP53 mutations in nearly all samples, as well as enrichment with RB1 mutations and amplifications of E2F3 and CCNE1 (Fig. 1, Supplementary Table S2.1.4). Those differences in pattern of mutation suggest the possibility of different oncogenic mechanisms.

Seventy-two per cent of the cancers in this study were from current or past smokers, consistent with extensive epidemiological studies indicating an association between smoking and urothelial cancer risk. In contrast with lung cancer, however, there was no statistically significant association between smoking status and the mutational spectrum, frequency of mutation in any SMG, focal SCNAs or expression subtype (Supplementary Tables S2.9.1 and S2.9.2). Never-smokers did have a slightly higher fraction of C>G mutations than did current/former smokers (28.5% vs. 23.8%, p = 0.032; Supplementary Figs. 2.3.2 and 2.3.3). However, unsupervised clustering of promoter CpG island DNA methylation data revealed a major subgroup (34%) of tumours characterized by cancer-specific DNA hypermethylation (CIMP) (Supplementary Fig.S7.1). Multivariate regression analysis with age, sex and tumour stage as covariates identified smoking pack-years as the only significant predictor of CIMP phenotype, as has also been reported for colorectal cancer¹⁸.

Fifty-one per cent of mutations overall were Tp*C->(T/G) (Supplementary Table S2.1.1), a class of mutation recently reported to be mediated by one of the DNA cytosine deaminases, APOBEC^19,20. APOBEC3B was expressed at high levels in all of the tumours, suggesting a major role for APOBEC-mediated mutagenesis in bladder carcinogenesis (Supplementary Figs. S12.1 and S12.2).

Four genes involved in epigenetic regulation were SMGs: MLL2, ARID1A, KDM6A and EP300 (Fig. 1). Truncating mutations were significantly enriched in each of those genes (Supplementary Fig. S2.2 and Data Files S2.2.1-2). Three of them had previously been identified as mutated in urothelial cancers⁴, but mutation of MLL2, which encodes a histone H3 lysine 4 (H3K4) methyltransferase, is a novel finding. Several other chromatin-regulating genes had mutation rates ≥10% but were not statistically significant by MutSig analysis: MLL3, MLL, CREBBP, CHD7 and SRCAP. Many other epigenetic regulators were mutated at lower frequency but were also enriched with truncating mutations, suggesting functional significance (Supplementary Fig. S2.2 and Data Files S2.2.1 and S2.2.2). Non-silent mutations in chromatin regulatory genes overall were significantly enriched in bladder cancer in comparison with the entire exome, in contrast with all other epithelial cancers studied to date in the TCGA project (Supplementary Table S2.10). Mutations in MLL2 and KDM6A (the latter encoding a histone H3 lysine 27 (H3K27) demethylase) were mutually exclusive (Supplementary Fig. S2.8 and Table S2.8.1), suggesting that mutations in the two genes have redundant downstream effects on carcinogenesis or that the combined loss is synthetically lethal.

Chromosomal rearrangements and viral integration

To identify structural variations and pathogen sequences, we used low-pass, paired-end, whole-genome sequencing (WGS; 6-8× coverage) of 114 tumours and RNA sequencing of all tumours. We detected 2,529 structural aberrations, including 1,153 that involve gene-gene fusions. Among the translocations, 379 were inter-chromosomal, 237 were intra-chromosomal, 274 were the result of inversions and 263 resulted from deletions (Supplementary Table S3.1). We found several recurrent translocations of likely pathogenic significance, including an intra-chromosomal translocation on chromosome 4 involving FGFR3 and TACC3 (n=3). The breakpoints were in intron 16 (2 cases) or exon 17 (1 case) of FGFR3 and intron 10 of TACC3 (confirmed by DNA sequencing and RNA-seq). All three lead to fusion mRNA products whose predicted proteins include the N-terminal 758 amino acids of FGFR3 fused with the C-terminal 191 amino acids of TACC3 (Fig. 2a). Based on the structure of the FGFR3-TACC3 fusion protein, we predict that it can auto-dimerize, leading to constitutive activation of the kinase domain of FGFR3. FGFR3-TACC3 fusion, which was recently described in both glioblastoma²¹ and bladder cancer^7,8, represents a promising therapeutic target. The ERBB2 gene was also involved in translocations in four tumours, all with different fusion partners and all confirmed by DNA sequencing, RNA-Seq or both. In one case, exons 4 to 29 of ERBB2 were fused to the promoter plus exon 1 of DIP2B, and the fusion product was amplified (Fig. 2b). Two other fusion products resulted in novel mRNA products whose biological significance is not known.

We identified viral DNAs in 7 of 122 tumours (6%), and viral transcripts in 5 of 122 (4%). Three tumours expressed cytomegalovirus (CMV) transcripts (encoding RL5A, RNA2.7, RL9A, RNA1.2, UL5 and UL22A), one expressed BK polyoma virus and one expressed human papilloma virus 16 (HPV16). HPV16 and human herpes virus 6B DNA were each identified in one other sample but without expression. None of the tumours expressing CMV showed evidence of CMV integration into the host genome, suggesting the presence of a stable episome. In the BK-positive tumour, two BK genes were integrated into GRB14, a signaling adapter protein for receptor tyrosine kinases. In the HPV 16-expressing case, the virus integrated into BCL2L1, an apoptosis-regulating gene (Fig. 2c). In that tumour, BCL2L1 was amplified (∼6×) and overexpressed (∼10× median; > 2× any of the other samples). Overall, those findings suggest that viral infection may play a role in the development of a small percentage of urothelial carcinomas.

mRNA, microRNA and protein expression

Analysis of RNA-seq data from 129 tumours identified four clusters (clusters I-IV) (Fig. 3 and Supplementary Fig. S4.1). Cluster I (‘papillary-like’) is enriched in tumours with papillary morphology (p=0.0002), FGFR3 mutations (p=0.0007, q=0.02), FGFR3 copy number gain (p=0.04, q=0.1) and elevated FGFR3 expression (p<0.0001) (Fig. 3a). It includes all three samples with FGFR3-TACC3 fusions. Cluster I samples also show significantly lower expression of miR-99a and miR-100, microRNAs that down-regulate FGFR3 expression (p=0.0002, Figs. 3aS5.3)²². They also show lower expression of miR-145 and miR-125b, which have been reported as frequently downregulated in bladder cancer²³. Tumours with FGFR3 alterations, and perhaps other tumours that share the Cluster I expression profile, may respond to inhibitors of FGFR or its downstream targets.

Reverse-phase protein array (RPPA) data indicate that clusters I and II express high HER2 (ERBB2) levels and an elevated estrogen receptor beta (ESR1) signaling signature, indicating potential targets for hormone therapies such as tamoxifen or raloxifene (Fig. 3d). In fact, the HER2 protein levels in a subset of the tumours are comparable to those found in TCGA HER2-positive breast cancers²³.

For comparison, we asked whether any of the four clusters show gene signatures similar to those identified in any other tumour type(s) among the first 11 analyzed by TCGA. We found that the signature of bladder cancer cluster III (‘basal/squamous-like’) is similar to that of basal-like breast cancers, as well as squamous cell cancers of the head and neck and lung (Supplementary Fig. S4.2)^24,25. All four of those cancer types express characteristic epithelial lineage genes, including KRT14, KRT5, KRT6 and EGFR. Basal-like subtype²⁶ and squamous cell subtype²⁷ of urothelial carcinoma, have been independently reported. Many of the samples in bladder cluster III express cytokeratins (i.e. KRT14 and KRT5) that were recently reported to mark stem/progenitor cells²⁶. Some of those samples also show a level of variant squamous histology (Fig. 3b). Bladder clusters I and II show features similar to those of Luminal A breast cancer, with high mRNA and protein expression of luminal breast differentiation markers, including GATA3 and FOXA1 (Fig. 3c). Markers of urothelial differentiation such as the uroplakins are also highly expressed in clusters I and II, as are the epithelial marker E-cadherin and members of the miR-200 family of microRNAs (which target multiple regulators of epithelial-mesenchymal transition)²⁸ (Fig. 3c). Taken together, those observations suggest that, despite their diverse tissue origins, some bladder, breast, head and neck and lung cancers share common pathways of tumour development.

To determine if the expression-based clusters could be seen in other datasets, we used the muscle-invasive bladder cancer samples from Sjodahl, et al.²⁷, hierarchically clustering them with the genes used in our analysis. From the sample dendrogram, we identified four groups (Supplemental Figure S4.3a). The four groups identified in the Sjodahl data set correlated well with the four clusters identified in our TCGA data. (Supplemental Figure S4.3b).

When we analyzed the RNA-seq data for transcript splice variation using SpliceSeq²⁹ (Supplemental Material S11), one finding of interest was an average of 3% PKM1 and 97% PKM2 transcripts in the tumour samples. The PKM2 isoform of pyruvate kinase is the principal driver of a shift to aerobic glycolysis in tumours (the Warburg effect)³⁰. Therefore, urothelial bladder cancers (and other cancer types) may prove sensitive to inhibition of glycolysis or related metabolic pathways.

Pathway analysis and therapeutic targeting

Integrated analysis of the mutation and copy-number data revealed three main pathways as frequently dysregulated in bladder cancer: cell cycle regulation (altered in 93% of cases); kinase and PI3K signaling (72%); chromatin remodeling, including mutations/SCNAs in histone modifying genes (89%); and components of the SWI/SNF nucleosome remodeling complex (64%) (Fig. 4a). To complement those results for well-defined pathways, we applied network analysis methods to examine other possible interactions between genes and pathways (Fig. 4b). In particular, we used the TieDIE algorithm to search for causal regulatory interactions within the PARADIGM network, which connects mutated genes to active transcriptional hubs^31,32. The analysis identified a sub-network linking mutated histone-modifying genes to a large array of activated transcription factors, suggesting potential far-reaching effects of histone modification on other pathways (Supplementary Fig. S8.2.1), converging on MYC/MAX regulation. Both MYC and MAX showed similar levels of pathway activity, independent of mutations in chromatin genes, suggesting that mutations in histone-modifying genes provide just one mechanism for disruption of the MYC/MAX hub. In contrast, tumours with chromatin-related mutations showed differential activity of transcription factors FOXA2 and SP1, implicating de-differentiation processes as a result of the mutations. Our network analysis also identified HSP90AA1 as a critical signaling hub, suggesting that inhibitors of HSP90 may have therapeutic value in urothelial carcinoma. Although the linkages between mutations and transcriptional changes were statistically significant in terms of their proximity in the network (as determined by permutation tests; see Supplementary Fig. S8.2), further studies will be needed to assess the biological relevance of the findings.

Integrated analysis also identified mutations, copy number alterations or RNA expression changes affecting the PI3-kinase/AKT/mTOR pathway in 42% of the tumours (Fig.5a). Included were activating point mutations in PIK3CA (17%; potentially responsive to PI3K inhibitors), mutation or deletion of TSC1 or TSC2 (9%; potentially responsive to mTOR inhibitors) and overexpression of AKT3 (10%; potentially responsive to AKT inhibitors). We also observed mutations, genomic amplifications or gene fusions that affect the RTK/RAS pathway in 44% of the tumours (Fig. 5b). Included were events that can activate FGFR3 (17%; potentially responsive to FGFR inhibitors or antibodies), amplification of EGFR (9%; potentially responsive to EGFR antibodies or inhibitors), mutations of ERBB3 (6%; potentially sensitive to ERBB kinase inhibitors) and mutation or amplification of ERBB2 (9%; potentially sensitive to ERBB2 kinase inhibitors or antibodies). ERBB3 mutations in bladder cancer have been noted previously⁴, but statistically significant mutation of ERBB2 in bladder cancer has not been reported. Both genes are potential therapeutic targets in other diseases^33-35 (Fig. 5c). Interestingly, ERBB2 alterations were approximately as frequent in this study as in TCGA breast cancers, but with fewer amplifications and more mutations (Fig. 5d)²⁴.