Contribution of systemic and somatic factors to clinical response and resistance in urothelial cancer: an exploratory multi-omic analysis

Background: Inhibition of programmed death-ligand one (PD-L1) with atezolizumab can induce durable clinical benefit (DCB) in patients with metastatic urothelial cancers, including complete remissions in patients with chemotherapy refractory disease. Although mutation load and PD-L1 immune cell (IC) staining have been associated with response, they lack sufficient sensitivity and specificity for clinical use. Thus, there is a need to evaluate the peripheral blood immune environment and to conduct detailed analyses of mutation load, predicted neoantigens and immune cellular infiltration in tumors to enhance our understanding of the biologic underpinnings of response and resistance. Methods and Findings: The goals of this study were to (1) evaluate the association of mutation load and predicted neoantigen load with therapeutic benefit, and (2) determine whether intratumoral and peripheral blood T cell receptor (TCR) clonality inform clinical outcomes in urothelial carcinoma treated with atezolizumab. We hypothesized that an elevated mutation load in combination with T cell clonal dominance among intratumoral lymphocytes prior to treatment or among peripheral T cells after treatment would be associated with effective tumor control upon treatment with anti-PD-L1 therapy. We performed whole exome sequencing (WES), RNA sequencing (RNA-seq), and T cell receptor sequencing (TCR-seq) of pre-treatment tumor samples as well as TCR sequencing of matched, serially collected peripheral blood collected before and after treatment with atezolizumab. These parameters were assessed for correlation with DCB (defined as progression free survival (PFS) > 6 months), PFS, and overall survival (OS), both alone and in the context of clinical and intratumoral parameters known to be predictive of survival in this disease state. Patients with DCB displayed a higher proportion of tumor infiltrating T lymphocytes (TIL) (n=24, Mann-Whitney p=0.047). Pre-treatment peripheral blood TCR clonality below the median was associated with improved PFS (n=29, log-rank p=0.048) and OS (n=29, log-rank p=0.011). Patients with DCB also demonstrated more substantial expansion of tumor-associated TCR clones in the peripheral blood 3 weeks after starting treatment (n=22, Mann-Whitney p=0.022). The combination of high pre-treatment peripheral blood TCR clonality with elevated PD-L1 IC staining in tumor tissue was strongly associated with poor clinical outcomes (n=10, HR (mean)=89.88, HR (median)=23.41, 95% CI (2.43, 506.94), p(HR>1)=0.0014). Marked variations in mutation loads were seen with different somatic variant calling methodologies, which in turn impacted associations with clinical outcomes. Missense mutation load, predicted neoantigen load and expressed neoantigen load did not demonstrate significant association with DCB (n=25, Mann-Whitney p=0.22, n=25, Mann-Whitney p=0.55, and n=25, Mann-Whitney p=0.29 respectively). Instead, we found evidence of time-varying effects of somatic mutation load on progression-free survival in this cohort (n=25, p=0.044). A limitation of our study is its small sample size (n=29), a subset of the patients treated on IMvigor 210 (NCT02108652). Given the number of exploratory analyses performed, we intend for these results to be hypothesis-generating. Conclusions: These results demonstrate the complex nature of immune response to checkpoint blockade and the compelling need for greater interrogation and data integration of both host and tumor factors. Incorporating these variables in prospective studies will facilitate identification and treatment of resistant patients.


Introduction
Atezolizumab has demonstrated responses in 15-25% of patients with advanced urothelial carcinoma and improved survival compared to historical expectations [1,2] . Similar to predictive factor analyses in melanoma, colon cancer and non-small cell lung cancer studies with other checkpoint blockade agents, Rosenberg and colleagues reported a statistically significant association between mutation load and response to atezolizumab in urothelial cancer patients [2] . However, mutation load in the atezolizumab study was predicted based on an estimate using a targeted panel and not with WES. Similar to findings from prior studies, the association between this predicted mutation load and outcomes in patients with urothelial cancer was not dichotomous; there were tumors from patients with elevated mutation load that did not respond to therapy, and vice versa. Additionally, positive PD-L1 staining of infiltrating immune cells by immunohistochemistry was associated with, but poorly predicted, response. A statistical model suggested that both PD-L1 staining and mutation load impacted the likelihood of response.
However, the authors did not recommend their clinical use.
Collectively, studies to date imply that a combination of immune parameters are necessary to gain further precision in determining the likelihood of benefit from these immunotherapies and that a single biologic marker will be insufficient. There have been few attempts to integrate molecular and immunologic data from patients treated with checkpoint blockade and their tumors. Consequently we performed whole exome (WES), RNA, and T cell receptor (TCR) sequencing of tumor samples from patients treated with atezolizumab as well as TCR sequencing of matched, serially collected peripheral blood.

Ethics Statement
All research involving human participants was approved by the authors' Institutional Review Board (MSKCC IRB), and all clinical investigation was conducted according to the principles expressed in the Declaration of Helsinki . Written informed consent was obtained from the participants.

Analysis plan
The analyses that were and were not included in the pre-specified analysis plan are detailed in S1 Pre-Specified Analysis Plan.

Patients and clinical characteristics
All patients had locally advanced or metastatic urothelial carcinoma and were treated at Memorial Sloan Kettering Cancer Center (n= 29 ) on protocol NCT02108652 [2] . All patients initiated therapy in 2014 and were treated with atezolizumab 1200 mg IV every 3 weeks, and as defined in the original study. Six patients had multiple samples evaluated for PD-L1 IC status; for four of these patients, the sample used for PD-L1 IC status in this analysis was the same as the sample that was whole exome sequenced. For the remaining two, the PD-L1 status of one patient's tumor (1994) used a separate primary tumor sample that also agreed in status with a metastatic sample; the other patient's tumor (6229) used a metastatic sample site that agreed in status with another metastatic sample site. Smoking status was evaluated using previously completed self-reported smoking questionnaires or review of medical records. One patient was excluded because the patient did not consent to correlative studies beyond PD-L1 testing that was performed as part of the clinical trial.
Tumor and blood samples All tumor tissue used for sequencing was obtained prior to dosing with atezolizumab. Tumor samples used for whole exome sequencing were all formalin-fixed paraffin embedded (FFPE).
The presence of tumor tissue in the sequenced samples was confirmed by examination of a representative hematoxylin and eosin-stained slide by a genitourinary pathologist (H.A.).
Peripheral blood mononuclear cells (PBMCs) were isolated and stored as previously described [3] . PBMC were collected pre-treatment and during treatment.

Clinical efficacy analysis
Tumor responses to atezolizumab were evaluated by CT scan every 9 weeks for the first 12 months following day 1 of cycle 1. After 12 months, tumor assessments were performed every 12 weeks. The response evaluation criteria in solid tumors (RECIST) version 1.1 was used to define objective clinical responses by the institutional radiologist. TCRβ sequencing   Tumor samples from patients 0522 and 6800 were excluded from tumor TCR analyses after   failing quality control. Patients 0979, 7592, and 8214 did not have available tumor for TCRβ sequencing and were therefore excluded from tumor TCR analyses as well. Genomic DNA was purified from total PBMCs and tumor samples using the Qiagen DNeasy Blood extraction kit.

DNA extraction and high-throughput
The TCRβ CDR3 regions were amplified and sequenced using immunoSEQ® (Adaptive Biotechnologies, Seattle, WA) as previously described [4] . In brief, bias-controlled V and J gene primers were used to amplify rearranged V(D)J segments for high-throughput sequencing at ~20X coverage. After correcting sequencing errors via a clustering algorithm, CDR3 segments were annotated according to the International ImMunoGeneTics Collaboration [5] to identify the V, D, and J genes that contributed to each rearrangement. A mixture of synthetic TCR analogs in each PCR was used to estimate the absolute template abundance (i.e., the number of cells bearing each unique TCR sequence) from sequencing data, as previously described [6] . The estimated TIL content was calculated as previously described [6][7][8] . To determine TIL content in FFPE samples as a T cell fraction, we amplified several housekeeping genes and quantitated their template counts to determine the amount of DNA usable for TCRβ sequencing.
ImmunoSEQ then amplifies and sequences the molecules with rearranged TCRβ chains.
Because the immunoSEQ assay aligns sequences to the IMGT database, sequences are annotated as complete VDJ rearrangements or non-productive rearrangements (a stop codon or out of frame CDR3 region was generated during VDJ recombination in one of the alleles); all downstream analysis in this work proceeded with complete, productive sequences. To estimate the number of starting templates that were in the sample, the number of sequence reads for each TCR β sequences is measured. Synthetic control templates were also spiked into each sample, thereby enabling quantitation of input TCR β templates from the read counts. For each sample, Shannon entropy was also calculated on the clonal abundance of all productive TCR sequences in the data set. Shannon entropy was normalized to the range by dividing Shannon entropy by the logarithm of the number of unique productive TCR sequences in the data set.
This normalized entropy value was then inverted to produce the clonality metric. Those T cell clones whose frequencies differed between samples from a given subject taken at different time points, or between cell populations (e.g., between total PBMCs and tumor) were computationally identified as previously described [9] . The input data consisted of the absolute abundance for each TCR clone in each sample. Fisher's exact test was used to compute a p-value for each clone across the two samples, against the null hypothesis that the population abundance of the clone is identical in the two samples. We corrected for multiple testing to control FDR using the Benjamini-Hochberg procedure and employed a significance threshold of 0.01 on adjusted p-values.

Whole Exome Sequencing
Twenty-six FFPE-derived tumor and frozen PBMC-derived normal paired samples were sequenced by exome hybrid capture using the IDT xGen Whole Exome Panel ( https://www.idtdna.com/pages/products/nextgen/target-capture/xgen-lockdown-panels/xgen-ex ome-panel ) and standard protocols. Briefly, each sample was used to create a barcoded Illumina library, tumor samples were pooled at optimal multiplex to create an equimolar pool into the hybrid capture reaction, which was performed according to the manufacturer's suggested protocol. Similarly, normal samples were pooled and introduced to the hybrid capture reaction.
Following the recovery of captured library fragments, PCR amplification was performed, the resulting pools of fragments were quantitated using qPCR ( [11,12] . Mutational signatures were inferred from the somatic mutation calls using deconstructSigs (v 1.6.0).
HLA typing HLA types for each patient were computed from the normal sequencing data using OptiType (v. Sleuth (v. 0.28.1) was used for differential expression analysis and GSEA was used for pathway enrichment analysis. ESTIMATE was used to quantify immune and stromal scores from RNA-seq data.

Neoantigen calling
Neoantigens were computed from all nonsynonymous mutations using Topiary (v. 0.1.0 , https://github.com/hammerlab/topiary ) and NetMHCCons (v. 1.1) with HLA alleles calculated by OptiType. As with expressed mutations, expressed neoantigens were those supported in the RNA with at least 3 uniquely-mapped reads matching the cDNA sequence.

Statistical analysis
All statistical analysis was performed in Python and R (v. 3 pre-treatment peripheral blood sample on which TCR-seq could be performed; 24 had one pre-treatment and at least one post-treatment peripheral blood collection. ^Visceral metastasis defined as liver, lung, bone, or any non-lymph node or soft tissue metastasis.*Based on [13] †Based on [14] Intratumoral and peripheral TCR features associate with durable clinical benefit The importance of T cells to the anti-tumor response has long been known [15] ; the relevance of intratumoral and peripheral T cell receptor (TCR) clonality to the anti-tumor response is an area of active study. A single previous study of melanoma patients treated with anti-PD-1 therapy demonstrated that patients whose tumors featured both high levels of tumor infiltrating T lymphocytes (TIL) along with high TIL clonality were more likely to experience radiographic response to therapy [8] . A separate study examined the peripheral TCR repertoire in anti-CTLA-4-treated patients with prostate cancer or melanoma, and found that clonotype stability was associated with response [16] . To our knowledge, no prior study has reported both intratumoral and peripheral TCR clonality in a single population treated with checkpoint blockade therapy.
We performed TCR sequencing of tumors and peripheral blood mononuclear cells In our patient group, we first asked whether there was an association between outcome and either TIL clonality or TIL proportion, or with clonality in the peripheral blood. Consistent with the data from Tumeh and colleagues [8] , tumors from patients who experienced a DCB exhibited a higher TIL proportion than those patients who experienced progressive disease, with a median of 0.21 (range 0.049-0.33) in tumors from patients who had PFS greater than 6 months, versus 0.069 (range 0.0098-0.24) in tumors from patients who did not ( n=24, Mann-Whitney p=0.047 , Fig 1A ). The consistency of result is notable given that Tumeh and colleagues used a different methodology (IHC) to assess TIL proportion than was used in our study. However, TIL proportion was not associated with continuous PFS ( n=24, log-rank p=0.32 ) or OS ( n=24, log-rank p=0.26 ) when split by median proportion. TIL clonality alone was not significantly associated with DCB ( n=24, Mann-Whitney p=0.10 , Fig 1A)  We next examined pre-treatment peripheral blood clonality and its relationship to DCB.
Because a diverse TCR repertoire in circulation may increase the likelihood that a tumor specific T cell population is present, we hypothesized that T cell receptor clonality would be inversely associated with response. We found that low pre-treatment peripheral TCR clonality was associated with improved PFS (split by median clonality, n=29, log-rank p=0. variants were identified and annotated as silent, missense or nonsense mutations (Fig 2A ).
There was no significant association between median missense mutation load and DCB (median mutations per megabase 3.24 (range 0.038-11.46) in patients with DCB compared to 0.45 (range 0.019-9.90) in those without DCB, n=25, Mann-Whitney p=0.22 , Fig 2B ). There was also no significant association between missense mutation load and overall survival greater  output [17,18] , we found that different filtering techniques impacted the association with DCB (S1 Table). Missense mutation load, when counting only mutations that were removed after post-processing (via Base Quality Score Recalibration (BQSR) and depth/variant allele frequency (VAF) filtering), was predictive of response ( n=25, Mann-Whitney p=0.0078 ).
One hypothesis for explaining the association between mutation load and outcome to treatment with checkpoint blockade is the generation of neoantigens, altered peptides presented by the major histocompatibility complex that are capable of eliciting an anti-tumor T cell response and are more common with increased mutation load.  Fig 2C and S2B Fig). Here, too, we acknowledge the limitations in statistical power to detect associations due to the sample size of our study.
The association between mutation load and response likelihood strengthens over time Given that the mutation load and outcomes were weakly associated in the complete IMvigor210 dataset and not statistically significantly associated in this cohort, we embarked upon an exploration of additional factors, including tumor microenvironmental and systemic measures, which may modify the importance of this variable or independently affect outcomes.
To this end, we examined the time-varying association between mutation load and PFS to see if mutation load had a differential association with early hazards in contrast to late hazards. We found evidence of time-varying effects of somatic mutation load on progression-free survival in this cohort ( n=25, p=0.044 , see Methods). We estimated these effects to consist of a stronger association of somatic mutation load with mortality or disease progression more than 3 months after treatment ( n=11, HR=0.69, 95% CI (0.38, 0.99) ) as compared to that during the first 3 months ( n=25, HR=0.91, 95% CI (0.75, 1.07) ; Fig 3A ). This effect estimate yielded a p-value for interaction of 0.1, which does not contradict the test for presence of the effect ( n=25, p=0.044 ) because that test is better powered. When a similar analysis was performed for time-varying association with OS, the evidence in support of the existence of time-varying effects was similar ( n=25, p=0.082 ; notable, despite p>0.05, given the PFS results above). In terms of the estimate of these effects, patients who survived longer than 3 months exhibited a stronger association between the number of somatic mutations per megabase and a lower risk of subsequent mortality ( n=11, HR=0.80, 95% CI (0.60, 1.00) ) as compared to that during the first three months ( n=25, HR=1.02, 95% CI (0.79, 1.22) , Fig 3B ).
This suggests that the time-varying effect is not likely an artifact of differential association with survival versus progression. Looking at the Kaplan-Meier estimates of progression-free survival among patients with mutation load per megabase above and below the median value of 1.03, it is apparent that there is very little separation of these two populations until approximately 3 months after treatment, and that there is a high frequency of both progression and mortality events at this time ( Fig 3C ). While we report these results using a pre-specified threshold of 3 months, in a nonparametric analysis we found that the reduction in risk associated with somatic mutation load increased steadily over time without the emergence of a clear inflection point  clinical factors distinguishing those likely to experience a rapid and early death from those more likely to survive longer [13] . We hypothesized that such patients might simply be too clinically and systemically unwell to mount the necessary immune response, despite some of them harboring tumor biomarkers thought to confer a likelihood of DCB, including elevated mutation load. When we examined the 5-factor score in this subset relative to the rest of the dataset, we found that indeed patients who survived less than or equal to 3 months exhibited a significantly Given that such agnostic approaches did not reveal a clear association between tumor microenvironment factors and response, we pursued a hypothesis-driven approach examining the genes that show upregulation at the cell surface during T cell exhaustion. When categorized by DCB, there was no significant difference in expression of such genes, including CTLA-4, TIGIT, HAVCR2 (TIM-3) or LAG-3 [21] . When grouped by PD-L1 staining, we found low expression of all markers in the PD-L1 low group (IC0), as expected. However, in the PD-L1 high group (IC2), HAVCR exhibited significantly higher expression in tumors from patients who experienced DCB than in those who did not (S4F Fig). Interestingly, of the three IC2 tumors, two had missense SNV loads significantly below the median (17 and 57); the third had 412 SNVs.
Although our limited sample size precludes making an assertion that mutation load is associated with survival in any particular subgroup (e.g. when looking among IC1 and IC2 tumors alone); our data do support the presence of an interaction among these variables (p= 0.046 for interaction; S5A Fig). Given the plausibility of the finding that somatic mutation load may correlate better with survival among patients with an inflamed tumor microenvironment, the addition of somatic mutation load to PD-L1 IC staining warrants further study.
We found a similar albeit weaker interaction effect when looking at association of somatic mutation load (missense SNV count per megabase) and progression-free survival according to the presence/absence of liver metastasis prior to treatment administration (p= 0.14 for interaction). Among patients without liver metastasis, somatic mutation load was associated with a lower risk for disease progression or mortality ( n=16, HR=0.73, 95% CI (0.50, 1.02) , S5B To our surprise, although both PD-L1 staining and mutation load were each associated with response in the original study [2] these variables did not correlate with each other (Fig 4A ).
Furthermore, pre-treatment peripheral TCR clonality did not correlate with mutation load ( Fig   4B ). The lack of association between these variables suggests that each might have an independent or semi-independent role in determining the likelihood of response to therapy. TCR clonality and infiltration did, however, correlate with PD-L1 IC score: those tumors with higher clonality or higher infiltration also featured higher PD-L1 staining (p=0.02 and p=0.01, respectively, Figs 4C , 4D ). In an analysis to see whether the association between pre-treatment peripheral TCR clonality and progression-free survival varied by PD-L1 IC score, we found some evidence of an interaction (p= 0.015 for interaction; Fig 4E ). Among patients with low levels of PD-L1 expression, there was little association between pre-treatment peripheral TCR clonality and  Fig 4E ). Similar results were seen in analyses with respect to OS, and in a logistic regression analysis for DCB (S 5C F i g , S 5D F i g , S 5E F i g ).
To resolve the hypothesis that those patients with low peripheral TCR clonality simply were healthier, we examined the association between 5-factor score and pre-treatment peripheral TCR clonality and did not find such an association ( n=26, Spearman rho=0.25 p=0.22 , S5F Fig).
In a multivariate survival model for time to disease progression or mortality, which allows the effect of each biomarker to vary according to intratumoral PD-L1 IC score, we find that the correlation of each intratumoral, peripheral or clinical biomarker with disease progression or mortality is relatively independent of the others (Fig 4F , S5G Fig). Perhaps with the notable exception of the association between liver metastatic status and time to progression or survival, the correlation of each intratumoral or peripheral biomarker with outcome is strongest in the group with the highest levels of immune cell PD-L1 expression (S2 Table).

Discussion
The treatment of previously incurable metastatic solid tumors with checkpoint blockade agents has led to dramatic success in a minority of patients, a finding that has generated substantial excitement in the field, with associated correlative studies and drug development. Here, we undertook the in-depth characterization of tumors and peripheral blood from 29 patients treated on IMvigor 210, a Phase II study in which 310 patients were treated with the anti-PD-L1 agent atezolizumab. In this cohort of patients, we illustrate the importance of host immune factors, including intratumoral and peripheral T cell receptor clonality, infiltration and expansion, to clinical outcomes. We did not find a significant association between mutation or expressed neoantigen load and progression free survival or DCB (defined as progression free survival (PFS) > 6 months). However, we did observe a time-dependent relationship between mutation load and outcome, wherein a relationship between mutation burden and outcome could only be detected in those patients surviving greater than 3 months. This analysis implies that patients who experienced rapid progression may display systemic indicators of immune deficiency despite elevated mutation load in the tumors. Calculation of the hazard ratios for each measured biomarker and clinical factor underscores the concept that a complex interaction of both host and tumor variables determines whether a patient will experience clinical benefit from anti-PD-L1 therapy. Although the overall study found significant associations between mutation load as measured by the Foundation Medicine targeted sequencing panel and radiographic response [2] , there was no statistically significant association between mutation load and durable clinical benefit or survival in the patient subset studied here, despite the similarity of our study population to the parent study. This contrast may be due to a combination of factors. First, though statistically significant, the association in the overall study was not categorical: as in other studies of mutation load, this factor alone was not predictive of response. Second, we have less power to detect this association in our smaller subset compared with the larger studied cohort. Third, standardized definitions and calculations of mutation load do not exist as yet; each published study has used differing methodologies [2,10,22,23] . Indeed, in this study, depending on the method used, the association between mutation load and clinical outcomes varied from p<0.08 to p>0.4 (AUCs and p-values in S1 Table). In an attempt to deepen our understanding of the biology of response and resistance, we studied additional factors. We found that even in this small dataset, TCR clonality below the median in the peripheral blood prior to treatment, expansion of tumor-associated TCR in the periphery 3 weeks after initiating treatment, and higher TIL proportion are all associated with clinical benefit. These data suggest that TCR sequencing provided additional insights into response and resistance beyond mutation load and PD-L1 staining. With respect to biomarker development, our study implies that non-invasive metrics such as pre-treatment peripheral TCR clonality and known prognostic features such as the presence of liver metastases may be worthy of further study in urothelial cancer patients treated with PD-L1 blockade.
From a mechanistic perspective, these findings imply an important relationship between circulating and intratumoral immunity upon PD-L1 blockade. We hypothesize that low TCR clonality in the peripheral blood prior to treatment increases the likelihood that a patient harbors one or several clones capable of tumor recognition, whether of neoantigens or tumor-associated antigens. The expansion of tumor-associated TCRs in the peripheral blood underscores the continuity of the tumor and blood compartments, and suggests that the activity of PD-L1 blockade may involve circulating T cells more than was previously thought. Indeed, this raises the possibility that anti-tumor T cells may home from the periphery to the tumor before later recirculating.
Finally, though limited in power by the small sample size, we attempted to integrate the importance of the studied variables. This analysis demonstrated both hypothesized and unexpected interactions. For example, while mutation load seemed to be associated with outcome more significantly in PD-L1 IC1 and IC2 tumors, high PD-L1 IC staining in the setting of high peripheral TCR clonality was associated with a substantial hazard for poor outcome.
Given the significance of PD-L1 expression in mediating response to anti-PD-L1 therapy, the presence of these interactions may argue in their favor as predictive rather than prognostic biomarkers. Further analysis is required to elucidate the role of these biomarkers in mediating response to checkpoint blockade.
This study has several limitations. The patients under study were treated at a single institution and represent a small subset of the overall study, limiting statistical power. As a single arm Phase II study, there is no control arm for comparison. Tumor samples were FFPE and were not collected immediately prior to treatment initiation. Only one sample per patient was utilized, which does not necessarily capture the heterogeneity of each tumor. Finally, a number of analyses were performed on a small number of patients without an independent validation cohort; although most analyses were pre-specified, there is a risk of Type II error and no adjustments were made for multiple testing.
In conclusion, we have demonstrated the potential value of pursuing an integrated study of somatic, immune and clinical features in order to elucidate the biological mechanisms underlying response to checkpoint blockade and ultimately improve our ability to practice precision medicine. We hope this work will motivate further multi-omics studies of response to checkpoint blockade.

Supplementary Information, Captions
Supplementary Tables S1 Table. Choice of depth and VAF filtering, as well as whether or not to run Base Quality Score orange indicates TCRs present in the tumor and expanded in the blood with treatment. (C) There was no significant expansion of TIL-associated TCR clones between pre-treatment ( 3.00 (range 1.00-9.00) ) and 6 weeks post-treatment ( 2.00 (range 1.00-8.00) ), n=20, Mann-Whitney