Circulating Tumor DNA Reveals Clinically Actionable Somatic Genome of Metastatic Bladder Cancer

Urinary bladder cancer is the fifth most common malignancy in North America, and was diagnosed in approximately 430,000 patients worldwide in 2012 (1). A quarter of patients present with muscle-invasive bladder cancer, which is generally high-grade and progresses to metastatic disease in 50% of patients (2). Until recently, cisplatin-based chemotherapy was the only treatment with a proven survival benefit for metastatic bladder cancer, and the 5-year survival rate with this regimen is only 5% to 15% (3). Recent data suggest that immune checkpoint blockade elicits responses in 15% to 30% of patients with advanced bladder cancer that have progressed on cisplatin-based therapy (4–7). Furthermore, targeted agents directed against ERBB2 (HER2), FGFR3, or the PI3K/mTOR pathway may hold promise in patient subsets (8–10).

Patient molecular stratification will be critical for optimal use of second-line therapies in the metastatic setting. Although few metastatic tumors have been studied, primary bladder cancer is molecularly heterogeneous and characterized by high somatic mutation rates (the third highest of all studied cancers) and recurrent copy-number alterations. Common tumor suppressors (e.g., TP53, RB1, CDKN2A) are frequently lost, and the majority of patients have alterations to chromatin modifiers and/or the PI3K/mTOR and MAPK pathways (11–15). Furthermore, there are distinct molecular subtypes based primarily on RNA expression (13, 16–18). The influence of specific mutations or molecular subtypes on clinical outcomes is not yet clear, but there is evidence to suggest that ERCC2 and TSC1 mutations correlate with cisplatin and everolimus sensitivity respectively, and that RNA subtypes are linked to response to atezolizumab and cisplatin (4, 16, 19, 20). High somatic mutation rates also appear to be a biomarker for response to immune checkpoint inhibition, consistent with similar data from melanoma, colorectal, and NSCLC (21–23). The majority of current biomarker approaches rely upon archival tissue obtained at transurethral resection or cystectomy, which may have reduced applicability in the setting of second and subsequent lines of therapy. A method for real-time genetic analysis of metastatic tumors is urgently needed.

The analysis of circulating tumor DNA (ctDNA) in plasma has shed light on the somatic landscape of metastatic disease in several solid malignancies. Prognostic and predictive ctDNA-based biomarkers are beginning to emerge in prostate, colon, and pancreatic cancer (24–26). In NSCLC, a “companion diagnostic” ctDNA test for the EGFR T790M mutation recently received FDA approval. However, the abundance and utility of ctDNA in advanced bladder cancer remains largely unexplored, with prior studies limited by small cohort sizes and dependency on single-gene sequencing and/or digital droplet PCR (ddPCR) assays to detect specific missense mutations (27–29). While ddPCR provides the sensitivity required to detect extremely low levels of ctDNA (e.g., for predicting disease relapse after curative therapy (29), only broad next-generation sequencing assays can deliver a comprehensive analysis of clinically relevant alterations such as PI3K/mTOR pathway deregulation or somatic hypermutation. In this study, we applied a custom sequencing approach to plasma cell-free DNA (cfDNA) collected from a cohort of patients with aggressive bladder cancer. We reveal the landscape of somatic alterations detected across 50 clinically relevant driver genes in metastatic bladder cancer, demonstrating the promise of liquid biopsies as both a discovery tool and a biomarker for therapy selection.

We accrued 51 consecutive patients with muscle-invasive bladder cancer to our liquid biopsy program between December 2014 and December 2016 (Supplementary Table S1). Approval for this study was granted by the University of British Columbia Ethics Board, or the local ethics review board for participating sites in Edmonton, Canada and Tübingen, Germany. The study was conducted in accordance with the Declaration of Helsinki, and written informed consent was obtained from all participants prior to collection. Methods for plasma cfDNA and buffy coat fraction germline DNA (gDNA) extraction from whole blood are listed in the Supplementary Methods (30). Control DNA for testing the targeted sequencing panel was extracted from bladder cancer cell lines RT4v6 and RT112. Cell lines were obtained in 2008 from MD Anderson Cancer Center (Houston, TX) and genotyped using a Promega StemElite kit to PCR amplify eight short tandem repeat loci plus a gender determining marker, Amelogenin. Mycoplasma testing was performed in July 2016 prior to DNA extraction using PlasmoTest (InvivoGen).

We employed a targeted sequencing strategy capturing the exonic regions of 50 bladder cancer relevant genes, as well as the TERT promoter and FGFR3 introns (Supplementary Table S2). Genes were selected on the basis of the established landscape of recurrent mutations and copy-number alterations in bladder cancer (12–15). In total, 0.7 Mb of genomic sequence was captured (including off- and near-target regions). For each cfDNA sample, 10–25 ng of DNA was used for library preparation. A-tailing, end repair, Illumina-compatible adapter ligation, and PCR amplification (between 15 and 17 cycles) was performed. Library quantification was carried out with the NanoDrop spectrophotometer, and each library was run on an ethidium bromide gel to confirm success. Up to 25 purified sample libraries at a time were multiplexed to obtain single pools with a combined mass of 1 μg, allowing a minimum 40 ng input for each sample library. These pools were hybridized to the capture panel for a minimum of 16 hours at 47°C. The subsequent wash, recovery, and amplification of the captured regions was performed according to the NimbleGen SeqCap EZ system protocols. Final libraries were purified with Agencourt AMPure beads and quantitated using either the KAPA qPCR kit, or the Qubit 2.0 Fluorometer (Life Technologies) and Qubit dsDNA HS Assay Kit. Pools were diluted to 20 pmol/L, and were sequenced on Illumina MiSeq (V3 600 cycle kit) or HiSeq 2500 (V4 250 cycle kit) machines. Preparation and sequencing of gDNA libraries was performed as described for cfDNA, with the following exceptions: all gDNA samples were sheared into 180-bp fragments with a Covaris focused-ultrasonicator prior to library preparation; 100 ng of DNA were used per library; PCR amplification was limited to 12 cycles. Exome capture and sequencing was performed using libraries previously prepared for targeted sequencing, and following the identical protocols described above, utilizing the Roche Nimblegen SeqCap EZ MedExome kit.

Paired-end reads were aligned against the hg38 reference genome using Bowtie-2.2.4 (31). Optical and PCR duplicates were removed using samblaster-0.1.20 (32). Adapters were trimmed in paired mode using cutadapt-1.11 (33). Low-quality read tails (smoothed baseq < 30) were trimmed using an in-house algorithm. Per-base read coverages in target regions were counted using bedtools-2.25.0 (34). To verify sample identities and pairings, sample genotypes were compared at over 500 high-confidence genotyped SNPs.

Somatic mutations were called in cfDNA samples by searching for variants with a mutant allele fraction (MAF) of ≥ 1%, and ≥ 10 supporting reads. In addition, the allele fraction was required to be 25 times higher than the background error rate (i.e., the average allele fraction across all gDNA samples), and three times higher than the allele fraction in the paired gDNA sample. Protein-level consequences of variants were predicted using ANNOVAR (35), and all variants were annotated with their frequency in the COSMIC database. All variants were further annotated with their population frequency in the ExAC and KAVIAR databases, and somatic mutation candidates with a population frequency > 0.5% were discarded if paired gDNA exhibited low coverage at the site. Finally, all somatic mutation candidates were inspected in Integrative Genomics Viewer (IGV). Variants in low mappability regions and variants near indels were flagged for detailed scrutiny. Because of low coverage at TERT promoter mutation hotspots, these mutations were called by manual inspection of read pileups. PCR and Sanger sequencing were performed to confirm high allelic frequency somatic mutations in cases where flanking regions were amenable to primer design, and ample cfDNA was available.

The somatic mutation rate of a given sample (in mutations per megabase) was estimated as |${\sum _k} = {\,\,\frac{{1000000}}{{Bk}}$| where B_k was the number of genomic locations sequenced to sufficient depth for mutation k to be detected. The threshold for sufficient depth was R/F_k where R was the minimum number of alternate allele reads required for a mutation to be called, and F_k was the allele fraction of mutation k.

ctDNA fractions were estimated on the basis of the allele fractions of autosomal somatic mutations. We first observed that the MAF in diploid chromosomes is highest when the mutation is combined with loss-of-heterozygosity (somatic mutations in genes with a detectable amplification were omitted from ctDNA fraction estimation). In this situation, the MAF and ctDNA fraction are related as MAF = (ctDNA × 1)/[(1 – ctDNA) × 2 + ctDNA × 1], and so ctDNA = 2/((1 / MAF) + 1). As we could not detect deletions or loss-of-heterozygosity from low ctDNA fractions, we conservatively assumed that all mutations could have associated LOH. To deal with stochastic variation in observed mutant allele read counts, we modeled the mutant read count as arising from a binomial distribution, and conservatively calculated what the true MAF would be if the highest observed MAF was a 95% quantile outlier.

Reads were counted in all capture regions using bedtools-2.25.0 (34). Coverage logratios were calculated against a low-noise median reference derived from all gDNA samples. GC fraction was calculated for all target regions, and Loess regression was applied to bait logratios to correct for GC content bias. After GC bias correction, read counts were corrected for residual differences in overall coverage using median-of-ratios normalization.

As hemizygous copy-number changes are associated with a loss of one allele at heterozygous bases, we used gDNA samples to find all heterozygous germline SNPs in captured regions, and then calculated the allele fractions of those SNPs in the corresponding cfDNA samples. We then calculated a median divergence from heterozygosity, defined as median(abs(AF - 0.5)), for each gene and cfDNA sample.

A deletion was called for a gene when coverage logratio ≤ –0.1 and MAF ≥ 0.6, or regardless of MAF if coverage logratio ≤ –0.3. A gain was called for a gene when coverage logratio ≥ 0.1 and MAF ≥ 0.6, or regardless of MAF if logratio ≥ 0.3. These conservative thresholds were determined empirically by studying a plot of coverage logratios and heterozygous SNP allele fractions in samples with and without detectable ctDNA (Supplementary Fig. S1). Amplifications and homozygous deletions were called using very conservative logratio thresholds of 0.7 and –1.0, respectively.

ctDNA fraction corrected absolute gene copy-numbers were estimated using the formula C_abs = (2 × 2L – 2) / F + 2, where C_abs is the absolute copy-number, L is the median coverage logratio of targeted regions inside the gene, and F is the ctDNA fraction.

Although most large deleterious genomic rearrangements fall within introns and are likely undetectable by an exon-limited panel, we searched for those in or near exonic regions by splitting unaligned reads into two 30-bp anchors that were then aligned to the hg38 genome. Discordant anchor pairs were grouped into clusters by position, and anchor pair clusters with eight or more supporting reads were flagged as rearrangement candidates and manually curated using IGV and BLAT. For the FGFR3-ADD1 fusion gene, putative coiled coil oligomerization domains in the ADD1 protein were identified using the COILS software (36).

To span the heterogeneous clinical spectrum of aggressive bladder cancer, we collected 85 plasma samples from 51 patients, including 14 patients with localized muscle-invasive bladder cancer (at initial sample collection) and 37 with nodal and/or distant metastatic disease (Fig. 1). Two patients provided matched metastatic tissue samples (peritoneal carcinomatosis and lymph node). Our cohort included plasma collected from patients prior to, during, and after platinum-based chemotherapy in both the neoadjuvant/adjuvant localized and metastatic setting. We obtained significantly less cfDNA from chemotherapy-naïve localized patient samples (n = 13) compared with those with metastatic disease not currently receiving systemic therapy (n = 42; median 6.4 vs. 20.3 ng/mL plasma; mean, 9.6 vs. 43.2; P = 0.00038, rank-sum test; Fig. 2A). Samples from patients receiving platinum-based chemotherapy at time of collection (n = 26) had high cfDNA yields (median 20.8 ng/mL plasma; mean, 30.9; range, 5.5–110) regardless of disease state (Fig. 2A).

For cost-effective next-generation sequencing of cfDNA, it is not possible to retain high sensitivity for ultra-low ctDNA fractions while simultaneously covering a large genomic territory. However, while the proportion of cfDNA that is tumor-derived (the ctDNA fraction) is largely unknown in metastatic bladder cancer, it is often greater than 1% in advanced solid malignancies studied to date (24, 37, 38). In data from The Cancer Genome Atlas (TCGA), 98% of primary muscle-invasive bladder cancer had a nonsynonymous somatic mutation in at least one of 50 bladder cancer genes in our panel (13). As such, we reasoned that with a sequencing depth of 500 to 1,000×, somatic mutations would theoretically be detected in almost all patients should their ctDNA fraction exceed 5%. To test this, we spiked genomic DNA from bladder cancer cell line RT4v6 into gDNA from patient 10044. Read support for five tested variants was observed at the lowest level of dilution (1% cell line DNA; Supplementary Fig. S2). Using our conservative thresholds for somatic variant calling, all five variants were independently detected at 10% dilution, four of five at 5% dilution, and two of five at 2% dilution. Strong copy-number alterations such as CDKN2A biallelic loss and PIK3CA amplification (both present in RT4v6 cells) passed our logratio thresholds for detection at dilution factors of 20% to 30%, but not at lower dilutions, although we could not use germline SNPs to augment copy-number calling as we would for paired patient specimens (see Methods; Supplementary Fig. S3).

We applied our targeted sequencing methodology to 62 of 85 cfDNA samples (median depth 753×, IQR 640×–922×; Supplementary Figs. S4 and S5), and matched leukocyte gDNAs for all patients subjected to cfDNA sequencing (median depth 427×, IQR 357×–496×). This included all cfDNA samples (n = 52/52) from patients with metastasis together with representative samples (n = 10/33) from patients with localized muscle-invasive bladder cancer. As expected, given their lower tumor burden, only one of seven localized patients (3/10 samples, all originating from the same patient) had evidence of ctDNA above our detection thresholds (Fig. 2B–D). Among patients diagnosed with metastatic disease, 58% of cfDNA samples (30/52; 24/37 patients) had evidence of ctDNA above 2% of total cfDNA; however, 12 of these samples were collected from patients currently undergoing chemotherapy and/or without evidence of disease at time of collection. When considering samples obtained from patients with distant metastatic disease not currently receiving systemic therapy, 73% (16/22) had ctDNA fractions above our detection threshold (Fig. 2C and D).

Among the 33 samples with evidence of ctDNA, the ctDNA fraction (see Methods) ranged from 3.9% to 72.6% (mean, 26.6%; median, 19.8%; Fig. 2B; Supplementary Table S3). Eleven samples had a ctDNA fraction over 30%, and all were from patients with distant metastatic disease at the time of collection. Conversely, 12 of 23 samples with ctDNA fractions between 3.9% and 30% were from patients with distant metastatic disease (P = 0.0058 for comparison to 11/11; Fisher exact test). In patients with distant metastatic disease, regardless of current treatment modality, there was no apparent difference in ctDNA fraction between samples from patients who had received prior cystectomy (n = 21), and those where the primary tumor remained intact (n = 12; median 13.1% vs. 6.8%, P = 0.29, rank-sum test; Fig. 2D). Although these numbers are small, it suggests that the majority of ctDNA is being shed from metastatic sites in these patients.

Four cfDNA samples carried an estimated ctDNA fraction higher than 70% (Fig. 2B). In each of these samples, the high ctDNA fraction was supported by multiple high allele fraction somatic mutations (Fig. 2B) and strong copy-number logratios (Supplementary Fig. S6). All four patients also had distant metastases, and high cfDNA yields (between 60 and 369 ng/mL plasma). Sanger sequencing confirmed the presence of high allele fraction somatic mutations (Supplementary Fig. S7).

Across the 26 patients with quantifiable ctDNA in at least one cfDNA sample (or in two cases, metastatic tissue), we detected 281 somatic mutations including 121 protein-altering mutations (Fig. 3A; Supplementary Table S4). We identified 22 missense or truncating mutations in TP53 in 17 of 26 patients (65%; Fig. 3B). Incorporating copy-number results, 24 of 26 patients carried either a TP53 inactivating change (n = 19), RB1 inactivating change (n = 8), MDM2 gain (n = 2), or CDKN2A loss (n = 6), likely resulting in disrupted cell-cycle regulation (Fig. 3C). For the two patients without alterations in these four genes, E-006 had nonmetastatic disease, and BC-008 had a low ctDNA fraction not amenable to copy-number analysis. Four patients exhibited definitive evidence for biallelic inactivation of RB1 or CDKN2A (Fig. 3C). Over half the cohort (19/26 patients) had mutations or disrupting rearrangements in chromatin modifier genes (Fig. 3C), including eight truncating mutations within ARID1A that remove the DNA-binding domain and/or glucocorticoid receptor binding domain, and eight truncating mutations in KMT2D (MLL2). TERT promoter mutations were identified in 12 of 26 patients: 10 patients harbored the chr5:1295113:G>A mutation (reported in 65% of bladder cancer), while two carried the chr5:1295135:G>A mutation (reported in 10% of bladder cancer; ref. 39). The majority of patients had alterations to the PI3K/mTOR pathway, including six with hotspot missense mutations in PIK3CA [K111E, E542K, E545K (n = 2), E674Q or E726K], four with truncating mutations in TSC1 (mutually exclusive with PIK3CA hotspot mutations), one with PIK3R1 stopgain, one with PTEN stopgain, and one with a TSC2 truncating rearrangement (Fig. 3C). Six patients showed evidence for PTEN or PIK3R1 deletion.

Nine patients (35%) carried ERBB2 activating somatic changes, including amplification in five patients (Supplementary Fig. S8), and hotspot mutations (S310F, L755S, I767M, V777L) in four patients. Four patients carried ERBB3 activating alterations, including two amplifications and two hotspot mutations (M91I, V104L; Fig. 3C). Remarkably, one patient had an average of 71 copies ERBB2 in his ctDNA. Three more patients carried activating RAS mutations [KRAS G12D (n = 2), HRAS Q61R], and one patient carried a KRAS amplification. In total, 15 of 26 patients carried activating somatic alterations in the MAPK pathway.

Overall, we observed no significant differences in mutation frequency (across the 50 genes in our panel) between samples in our metastatic cohort and in the TCGA-localized muscle-invasive bladder cancer cohort (Supplementary Fig. S9; ref. 13). Even the distribution (within each gene) and nature of mutations (e.g., predominantly missense versus truncating) were highly consistent between the two cohorts, providing a strong measure of validation for mutations described here.

Two patients (T-015 and T-019) had both metastatic tissue and matched plasma cfDNA sequenced, but only T-015 had quantifiable ctDNA (estimated 5% fraction). Reassuringly, the solid and liquid biopsies for T-015 revealed the same five somatic mutations, although the CDKN2A missense mutation was detected at a variant allele frequency of only 0.9% in the liquid biopsy (Supplementary Fig. S10). Four patients in our cohort had evidence of ctDNA in more than one plasma collection obtained over the course of chemotherapy. Temporal samples from the three patients with distant metastatic disease were highly consistent, with the same patient-unique mutations detected over multiple timepoints (Supplementary Fig. S11).

Projecting our targeted sequencing data across the genome suggested a mean and median somatic mutation rate of 17.2 and 14.3 per Mb, respectively (Fig. 3A; Supplementary Table S4). Three of 26 (12%) of patients exhibited projected rates ≥ 30 per Mb: close to that detected in mismatch repair deficient colorectal cancers (21). Our findings are higher than the mean and median somatic mutation rates of 7.7 and 5.5 per Mb, respectively, observed in primary muscle-invasive bladder cancer (13).

To assess the degree to which our targeted sequencing–derived mutation rates recapitulate those observed over the entire exome, we subjected 11 cfDNA libraries from eight patients (the majority of cfDNAs harboring a ctDNA fraction over 25%, and their matched gDNA libraries) to exome sequencing (median depth 176×, IQR 171×–187×). Somatic mutation allele fractions and copy-number coverage logratios observed in exome sequencing data were consistent with the ctDNA fractions derived from targeted sequencing data (Fig. 4). When considering mutations robustly detectable by both panels (i.e., those with allele fraction > 10%) we observed a high correlation for mutation counts between the 50-gene targeted panel and whole-exome sequencing data (Pearson correlation = 0.83; Supplementary Fig. S12). Estimated somatic mutation rates (see Methods) were also correlated between the two panels (Pearson correlation = 0.71; Supplementary Fig. S12). Nevertheless, the targeted panel consistently overestimated somatic mutation rates relative to whole-exome sequencing (median 66% increase in mutation rate). This increase may be due to recurrent driver mutations in the targeted genes, although we observed a high proportion of synonymous or noncoding variants (unlikely to be under selective pressure) in the samples with the highest mutation rates (Supplementary Fig. S12B).

At the individual sample level, the exome data provided strong technical validation for mutation calls from the initial targeted sequencing data (Supplementary Fig. S12). The higher depth of sequencing in the targeted panel data inevitably detected more subclonal mutations than exome sequencing. For example, BC-005 (predicted ctDNA fraction = 26%) had 25 somatic mutations detected in targeted sequencing data (between AF 2.0%–17.3%). Only four of these mutations were detected in whole-exome sequencing with an allele fraction of 10% or higher (the detection threshold in exome data; Supplementary Fig. S12C). None of the 25 mutations showed a significantly different allele fraction between targeted and whole-exome sequencing (P > 0.01 for all, Fisher exact test). This suggests that the subclonal (< 10% allele fraction) mutations detected by targeted sequencing are real, but only detectable at high sequencing depth, highlighting the impact of sequencing depth on estimated somatic mutation rates.

Although the heterogeneity of our cohort prevents broad analysis of clinical outcomes, we note that six patients with metastasis in our cohort had received immunotherapy during the course of their disease (T-017, T-014, T-009: nivolumab; BC-017: atezolizumab; BC-014: durvalumab; 10044: pembrolizumab). BC-017 experienced a complete response to atezolizumab 2 years prior to plasma sampling (treatment was ceased due to immune-related colitis). T-014 had a partial response to nivolumab. Both patients exhibited high projected somatic mutation rates (22 and 39 Mb⁻¹, respectively). The remainder (all nonresponders) had either no/low ctDNA (T-017, T-009) or a low/median mutation rate (BC-014, 10044).

We did not detect any FGFR3 mutations, consistent with their predominance in non-muscle–invasive papillary bladder cancer (13). However, we identified a novel FGFR3-ADD1 fusion event in two cfDNA samples (collected 8 months apart) from patient T-002 (Fig. 5; Supplementary Table S5). The rearrangement juxtaposes the 3′-UTR of FGFR3 with the second intron of ADD1, likely resulting in an in-frame chimeric protein where the last FGFR3 exon is spliced out due to its lack of a splice donor site. As ADD1 contains a dimerization domain, it plausibly enables constitutive activation of FGFR3, functioning similarly to the well-described FGFR3-TACC3 fusion gene (40). We did not detect FGFR3-TACC3 in any patients here, although it was robustly identified in two cell lines (RT112 and RT4v6; Supplementary Fig. S13). In addition, we report two different ERBB2 fusions in patients with ERBB2 gene amplification. The fusion genes, ERBB2-CDK12 and MED1-ERBB2, likely represent secondary passenger events accompanying the formation of high-level gene amplification. Also detected (in four separate patients) were a 27-kb deletion in ARID1A, 34-bp deletion overlapping a TSC2 splice donor site, and truncating rearrangements in RB1 and RUNX3 (Supplementary Table S5). Repeat samples from the same patient, when available, supported the rearrangements.

Practical molecular stratification of patient tumors will be critical to guide the use of emerging targeted therapies and immunotherapy in metastatic bladder cancer. However, the limited availability of tissue from metastases in routine clinical practice hampers attempts to implement real-time precision oncology in this setting. In this context, we describe clinically relevant somatic alterations detected through minimally invasive plasma ctDNA profiling in the majority of patients with metastatic bladder cancer.

Confidence in our results can be drawn from several sources of evidence. First, the distribution and frequency of mutations, copy-number alterations, and chromosomal rearrangements was highly consistent with localized muscle-invasive bladder cancer (13, 14). Second, the only recurrent mutations detected in our study were known hotspots in PIK3CA, ERBB2, and the TERT promoter: all other mutations were unique to individual patients. Third, when multiple temporal samples from individual patients were analyzed, the same mutations were consistently observed. Fourth, somatic mutations with high allelic frequency from our targeted sequencing approach were independently verified by Sanger sequencing and matched exome sequencing. Copy-number calls and ctDNA fractions derived from matched exome sequencing data were also highly consistent with results from our targeted panel. Finally, similar to observations from other malignancies (28), we observed higher ctDNA fractions in patients with distant metastatic disease than those with local lymph node involvement or localized cancer.

Prior studies have suggested that bladder cancer releases abundant ctDNA, but were based on very small patient numbers and assays limited to specific missense mutations or single genes (27–29). In our study, the majority of patients with metastasis had ctDNA fractions above 2%, and one-third had fractions greater than 20%. This implies that samples from most patients with progressive metastatic disease are suitable for cost-effective targeted sequencing approaches, and a significant proportion for standard exome sequencing. However, even when excluding clinically metastatic patients likely to have a very low or absent tumor burden (e.g., those with pathologic lymph node disease), there remained some plasma cfDNA samples from patients with metastasis that yielded no ctDNA evidence. Our assay was also relatively uninformative in the localized setting, although the availability of tumor tissue collected at transurethral resection (or cystectomy) and urine-based biomarkers reduce the rationale for pursuing plasma ctDNA in this setting (27, 41). Nevertheless, future studies seeking to capture somatic information from all patients with metastasis could augment targeted panel sequencing with ddPCR for patient-specific mutations identified through archival tissue analyses or for recurrent hotspot mutations (e.g., in PI3KCA or TERT), as has been applied successfully to patients with clinically localized bladder cancer (27, 29). Alternatively, the incorporation of unique molecular identifiers during library preparation and subsequent ultra-deep sequencing has proved successful in other solid malignancies for detecting mutant alleles as rare as 0.02% (28, 42), but this can become expensive if applied beyond a few genes.

We observed high total cfDNA yields in patients receiving chemotherapy, in both the localized and metastatic setting. As this was not accompanied by increased (or any) ctDNA in most instances, it suggests elevated nonmalignant cell death in the wake of systemic cytotoxic therapy. In pregnant women, fetal cfDNA is cleared rapidly from circulation after childbirth (43), and ctDNA levels reduce after surgery in patients with NSCLC and colorectal cancer (44, 45). The same is likely true in patients with bladder cancer responding to chemotherapy. Future ctDNA studies in bladder cancer should carefully consider the timing of plasma collections, and be cautious extrapolating data from other cancers where overall cfDNA yield correlates with tumor burden (44).

In clinical trials to date, approximately 15% to 30% of patients with bladder cancer appear to respond to checkpoint blockade (2, 4–7). In the atezolizumab phase II trial, responding patients had a significantly higher mutation rate than nonresponders (4). In our study, approximately 12% of patients with metastasis had projected mutation rates (from targeted sequencing) above 30 per Mb: levels of hypermutation associated with response to immunotherapy in other cancers (21–23). It is likely that our targeted approach, encompassing frequently mutated genes, overestimates global mutation rates. However, overall trends in mutation counts from targeted sequencing were consistent with those derived from matched exome sequencing, suggesting that it is at least possible to distinguish between relatively high and low mutational burdens. Certainly, a minimally invasive assay for approximate mutational burden, as described herein, could prove valuable to the dissection of patient response to immune checkpoint blockade (4). It is plausible that ctDNA could be used in combination with IHC of PD-1/L1 and molecular subtype to help select patients most suitable for treatment.

The practical detection of mutations, copy-number changes and even novel fusion genes from plasma in metastatic bladder cancer offers exciting possibilities for precision oncology. We identified 11 patients with amplification and/or activating hotspot mutation to ERBB2/ERBB3. A recent phase II trial evaluated afatinib (an irreversible inhibitor of the ERBB family) in bladder cancer and found that activity was confined to patients with alterations in ERBB genes (46). The future of other HER2-directed agents such as trastuzumab and lapatinib in bladder cancer also likely depends upon biomarker-driven patient selection (10). Clinical trials of FGFR3 inhibitors are underway in bladder cancer, but dovitinib (for example) has shown no activity to date in patients with FGFR3 mutations or positive IHC (47). FGFR3 inhibitors are not yet tested in patients known to harbor activating gene fusions, even though bladder cancer cell lines expressing FGFR3 fusions are extremely sensitive to FGFR-targeted agents, more so than those carrying mutations (40). Our data demonstrate that ctDNA profiling of FGFR3 introns alone is sufficient to identify novel FGFR3 fusions. We identified several patients with TSC1/TSC2 truncating alterations, which are reported to confer sensitivity to the mTOR inhibitor everolimus (19). Indeed, the PI3K/AKT/mTOR pathway was aberrant in a large proportion of patients, suggesting opportunities for testing of AKT inhibitors such as MK2206 that has previously demonstrated activity in bladder cancer cell lines (14). We detected several mutations in DNA repair genes, including a missense mutation in ERCC2 and two truncating ATM mutations, raising the possibility that ctDNA profiling may also help guide implementation of cisplatin-based chemotherapy (48, 49).

The profiling of plasma ctDNA also shows promise for monitoring therapy response or disease progression in patients with bladder cancer. Here, we demonstrate proof-of-principle for ctDNA profiling to capture temporal ctDNA fraction and the prevalence of driver mutations in three patients with distant bladder cancer metastases undergoing chemotherapy. In the future, this approach could be leveraged to dissect associations between changes in ctDNA fractions during treatment and clinical outcomes. Indeed, a recent study suggested that presence of ctDNA in plasma samples collected postcystectomy in clinically localized bladder cancer was associated with disease recurrence (29), which is consistent with breast cancer where ctDNA detection after curative therapy portends rapid metastatic relapse (50). An added advantage of ctDNA profiling is the theoretical potential to survey intrapatient tumor heterogeneity, and the effect of the therapy on clonal dynamics. Advanced bladder cancer was recently shown to be characterized by high levels of intrapatient heterogeneity (51), but the comprehensive dissection of longitudinal heterogeneity in ctDNA samples will likely require broader genome-wide approaches and more matched metastatic tissue, than used in the current study.

Our study demonstrates that the liquid biopsy can provide a practical and cost-effective snapshot of bladder cancer driver genes status in patients with advanced disease. Given the relative ease of plasma collection, we propose that ctDNA analysis be incorporated in future clinical studies of metastatic bladder cancer, where possible. With the proviso of comprehensive somatic data, the liquid biopsy can inform valuable context to better dissect therapy response and accelerate rational implementation of precision oncology.

A. Stenzl is a consultant/advisory board member for Alere, Bristol-Myers Squibb, Ipsen, and Roche. No potential conflicts of interest were disclosed by the other authors.

Conception and design: T. Todenhöfer, B.J. Eigl, P.C. Black, A.W. Wyatt

Development of methodology: M. Annala, A.W. Wyatt

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): G. Vandekerkhove, T. Todenhöfer, W.J. Struss, A. Wong, K. Beja, S. Brahmbhatt, J. Hennenlotter, K.N. Chi, S. North, B.J. Eigl, P.C. Black

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): G. Vandekerkhove, M. Annala, W.J. Struss, K. Beja, E. Ritch, S. Volik, C.C. Collins, P.C. Black, A.W. Wyatt

Writing, review, and/or revision of the manuscript: G. Vandekerkhove, M. Annala, T. Todenhöfer, W.J. Struss, S. Volik, M. Nykter, K.N. Chi, A. Stenzl, B.J. Eigl, P.C. Black, A.W. Wyatt

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): G. Vandekerkhove, T. Todenhöfer, J. Hennenlotter, S. North, A. Stenzl, P.C. Black

Study supervision: M. Nykter, P.C. Black, A.W. Wyatt

This study was primarily supported by Bladder Cancer Canada (to A.W. Wyatt, B.J. Eigl, P.C. Black). Additional support was provided by the Canadian Urologic Oncology Group (to P.C. Black), the Emil Aaltonen Foundation (to M. Annala), the Academy of Finland#269474 (to M. Nykter), and the German Cancer Foundation (Deutsche Krebshilfe; to T. Todenhöfer).

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.