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 prespecified analysis plan are detailed in S2 Text.

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, were treated with atezolizumab 1,200 mg IV every 3 weeks, and provided written consent according to Institutional Review Board-approved protocols permitting tissue and blood collection, sequencing, and correlative studies. Patient tumor samples were assessed prospectively and centrally (by HistoGeneX, Brussels, Belgium) for PD-L1 expression by immunohistochemistry with the SP142 assay (Ventana, AZ, USA) [1]. The PD-L1 tumor-infiltrating immune cell (IC) status was defined by the percentage of PD-L1–positive ICs in the tumor microenvironment: IC0 (<1%), IC1 (≥1% but <5%), and IC2/3 (≥5%), as defined in the original study. Six patients had multiple samples evaluated for PD-L1 IC status; for 4 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 2, the PD-L1 status of 1 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]. PBMCs were collected pretreatment 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.

DNA extraction and high-throughput 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 tumors 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. 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 approximately 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 tumor-infiltrating T lymphocytes (TIL) content was calculated as previously described [6–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 nonproductive rearrangements (a stop codon or out of frame CDR3 region was generated during VDJ recombination in 1 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β sequence 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 2 samples against the null hypothesis that the population abundance of the clone is identical in the 2 samples. We corrected for multiple testing to control the false discovery rate (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-exome-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 (Kapa Bio) and sequenced in separate lanes by paired-end 150 bp reads, using the Illumina HiSeq 4000. Whole exome sequencing results for 1 sample (for patient 4072) were excluded after failing to meet coverage requirements.

Somatic variant calling

DNA sequencing data for the tumor and normal samples were aligned to the GRCh37 reference using bwa-mem (v. 0.7.10) with default settings. The resulting BAMs were processed through Picard MarkDuplicates and the GATK (v. 3.5–0) pipeline, including Base Quality Score Recalibration and Indel Realignment. Single nucleotide variants (SNVs) were called from Mutect (v. 1.1.6) and Strelka (v. 1.0.14) with default settings. Variants from either call were included and the variants calls were further filtered to those with depth (in normal and tumor samples) ≥7 reads, >10% tumor variant allele frequency (VAF), and ≤3% normal VAF [10]. Mutations per megabase was computed by normalizing the number of mutations by the number of exonic loci with ≥7 reads in normal and tumor samples, calculated using Pageant (https://github.com/hammerlab/pageant).

Variants were annotated as missense variants by Varcode (v. 0.5.10, https://github.com/hammerlab/varcode) and PyEnsembl (v. 1.0.3, https://github.com/hammerlab/pyensembl) using Ensembl Release 75 and annotated as deleterious using PolyPhen (v. 2.2.2). DNA damage response (DDR) genes were gathered from [11,12].

Mutational signatures were inferred from the somatic mutation calls using deconstructSigs (v 1.6.0).

Human leukocyte antigen typing

Human leukocyte antigen (HLA) types for each patient were computed from the normal sequencing data using OptiType (v. 1.0.0).

RNA sequencing

RNA was extracted from 26 FFPE tumor samples and evaluated for quality and quantity using the Agilent RNA pico chip. Each sample was prepared for sequencing by constructing an Illumina Tru-Seq Stranded RNA kit, according to the manufacturer’s protocol. The resulting libraries were amplified by PCR, quantitated, pooled, and processed through a hybrid capture intermediate using the IDT xGen Exome reagent (as above). The captured fragments were quantitated, diluted, and were sequenced using 2 x 150 bp paired-end reads on the Illumina HiSeq 4000.

The RNA sequencing (RNA-seq) data were aligned to the GRCh37 reference in Ensembl Release 75 using STAR (v. 2.4.1d), and transcript quantification was performed using kallisto (v. 0.42.3). The STAR alignment was only used for identifying variant-supporting reads in the RNA. For gene-level analysis, the transcript quantifications were aggregated to the gene level using tximport (http://f1000research.com/articles/4-1521/v1).

Expressed mutations and neoantigens were computed using Isovar (v. 0.0.6, https://github.com/hammerlab/isovar), based on the RNA reads overlapping each mutation. Sleuth (v. 0.28.1) was used for differential expression analysis, and Gene Set Enrichment Analysis (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.3.1). Cohorts (v. 0.4.0, https://github.com/hammerlab/cohorts) and Biokepi (https://github.com/hammerlab/biokepi) were used to orchestrate the analysis. The Mann Whitney and Fisher’s Exact test were performed using the Python scientific computing library, SciPy (v. 0.18.1). Kaplan-Meier curves were computed with Lifelines (v. 0.9.1.0). Survival and logistic regression models were estimated using PyStan (v. 2.12.0.0), and the Stan statistical computing software (v. 2.12.0). Survival analyses utilized a proportional hazards piecewise exponential model with a random walk baseline hazard. The analysis for presence of a time-varying covariate effect was performed in R using survival (v. 2.39.5) to look for the association of scaled Schoenfeld residuals with log(time), whereas the estimation of the time-varying covariate effect was performed using Stan. This analysis estimated the covariate effect at each timepoint with a random-walk prior. In some cases, alternative specifications of models written in Stan were interrogated as sensitivity analyses; see the project’s GitHub repository (https://github.com/hammerlab/multi-omic-urothelial-anti-pdl1) and S3 Text for details.

All analysis code is available at https://github.com/hammerlab/multi-omic-urothelial-anti-pdl1 for open access by readers.

Mutation calls, TCR-seq, and RNA-seq data are available at http://doi.org/10.5281/zenodo.546110. Additional data are available at https://github.com/hammerlab/multi-omic-urothelial-anti-pdl1.

Patient characteristics

Twenty-nine patients with metastatic urothelial cancer from a single institution, treated with atezolizumab, as part of a single-arm phase II study (IMvigor 210, NCT 02108652), were included in the analyses. The patients displayed characteristics typical of the metastatic urothelial cancer population studied in IMvigor 210: 25 of 29 were males with urothelial cancers of bladder origin, and 18 of 29 had a reported prior smoking history (Table 1). Patients had an Eastern Cooperative Oncology Group (ECOG) performance status of 0 or 1 and had 0 to 3 prior regimens of chemotherapy. Of this group, 25 patients had sufficient tumor tissue for WES, 26 for RNA-seq, and 24 for TCR-seq. Twenty-nine had a pretreatment peripheral blood sample on which TCR-seq could be performed; 24 had 1 pretreatment and at least 1 posttreatment peripheral blood collection.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Baseline characteristics.

https://doi.org/10.1371/journal.pmed.1002309.t001

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 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-cytotoxic T-lymphocyte-associated protein 4 (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-seq of tumors and PBMCs at serial time points in our cohort. Due to limitations in sample availability (Methods), this analysis included tumors from 24 patients and peripheral blood from 29 patients, including pretreatment samples in all patients, and between 1 and 8 total time points. A median of 141,255 (range 43,052–335,089) T cells were analyzed per peripheral blood sample, including 82,636 (range 24,095–207,860) unique TCRs, with a median clonality of 0.080 (range 0.014–0.37) and a median T cell proportion of 0.31 (range 0.082–0.64). In the tumors, the corresponding values included 1,402 (range 63–133,167) T cells, 1,086 (range 67–56,273) unique TCRs, clonality of 0.096 (range 0.033–0.34), and T cell proportion of 0.097 (range 0.0098–0.33).

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 durable clinical benefit (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 progression-free survival (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 results 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), continuous PFS (split by median clonality, n = 24, log-rank p = 0.51), or continuous OS (split by median clonality, n = 24, log-rank p = 0.47). Tumors with less than the median TIL proportion or TIL clonality, considered jointly as 1 feature, were less likely to display DCB (Fig 1B, 25% of patients with DCB versus 81% of patients without DCB, n = 24, Fisher's Exact p = 0.021). Considering TIL proportion and TIL clonality separately, when using median as a threshold, did not result in a significant difference in terms of DCB in either case (S1A Fig). It remains unclear whether TIL clonality adds to TIL proportion in its association with DCB in this study (TIL proportion and TIL clonality versus TIL proportion alone, n = 24, log-likelihood p = 0.100).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. T cell receptor (TCR) clonality and treatment response.

(A) Tumor-infiltrating T lymphocytes (TIL) proportion alone was associated with durable clinical benefit (DCB), with a median of 0.21 (range 0.049–0.33) in tumors from patients who had DCB versus 0.069 (range 0.0098–0.24) in tumors from patients who did not (n = 24, Mann-Whitney p = 0.047). TIL clonality alone was not significantly associated with DCB, with a median of 0.12 (range 0.047–0.34) in tumors from patients with DCB and a median of 0.092 (range 0.033–0.22) in tumors from patients without DCB (n = 24, Mann-Whitney p = 0.10). (B) Tumors with less than the median TIL proportion or TIL clonality, considered jointly as 1 feature, were less likely to display DCB (Fig 1B, 25% of patients with DCB versus 81% of patients without DCB, n = 24, Fisher's Exact p = 0.021). (C) Patients with pretreatment peripheral TCR clonality less than the median exhibited improved progression-free survival (PFS, n = 29, log-rank p = 0.048). (D) Patients with pretreatment peripheral TCR clonality less than the median exhibited improved overall survival (OS, n = 29, log-rank p = 0.011). (E) There was a significant association between TCR clonality in the peripheral blood prior to initiating treatment and overall survival greater than 12 months (DCB-OS) (DCB-OS: TCR clonality 0.060 [range 0.022–0.21]; OS less than 12 months: 0.15 [range 0.031–0.35, n = 29, Mann-Whitney p = 0.0061]). (F) There was no significant association between pretreatment peripheral TCR clonality and DCB (DCB: TCR clonality 0.068 [range 0.027–0.21]; no DCB: 0.14 [range 0.022–0.35, n = 29, Mann-Whitney p = 0.25]). (G) Expansion of TCR clones found in TIL (orange bars) occurred in the peripheral blood 3 weeks after initiating treatment in all patients. (H) The number of TCR clones found in TIL that expanded in the peripheral blood 3 weeks after initiating treatment was 8.00 (range 4.00–12.00) in patients with DCB and 2.50 (range 1.00–18.00) in non-DCB patients (n = 22, Mann-Whitney p = 0.022).

https://doi.org/10.1371/journal.pmed.1002309.g001

We next examined pretreatment 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 TCR clonality would be inversely associated with response. We found that low pretreatment peripheral TCR clonality was associated with improved PFS (split by median clonality, n = 29, log-rank p = 0.048), overall survival (OS, split by median clonality, n = 29, log-rank p = 0.011), and OS greater than 12 months (DCB-OS, n = 29, Mann-Whitney p = 0.0061; Fig 1C, 1D and 1E), although not with DCB (Fig 1F, n = 29, Mann-Whitney p = 0.25).

Finally, we explored the relationship between intratumoral and peripheral TCR clonality. Individual T cell clones present in tumors can be tracked in the peripheral blood during treatment (examples in S1B Fig). Expansion of tumor-associated TCRs occurred in the peripheral blood in all patients (Fig 1G). However, a more pronounced expansion of intratumoral TCR clones was observed in DCB patients at 3 weeks after initiation of treatment (second dose of therapy) (n = 22, Mann-Whitney p = 0.022, Fig 1H) that was not significant at 6 weeks after therapy initiation (n = 20, Mann-Whitney p = 0.17, S1C Fig). Interestingly, all patients with low pretreatment peripheral TCR clonality and high TIL clonality survived greater than 12 months (DCB-OS, S1D Fig).

Association of tumor genetic features with PFS or OS

To further examine intratumoral factors associated with therapeutic efficacy, we performed WES on 25 FFPE archived tumor samples. Mean target coverage was 129 (range 44–194) in tumors and 73 (range 59–91) in normal tissue. SNVs 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 DCB-OS (n = 25, Mann-Whitney p = 0.37, S2A Fig). In a survival analysis for time to disease progression or mortality, the estimated hazard ratio associated with increase in missense SNV count per megabase was 0.92 (95% CI [0.78, 1.09]). These results are not surprising given that the present sample size (n = 25) is underpowered to detect an effect of magnitude similar to that observed by Rosenberg and colleagues [2] (power = 0.2, assuming median of 12.4 versus 6.4 mutations per megabase among patients with DCB versus non-DCB response).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Single nucleotide variants (SNVs) and treatment response.

(A) SNVs, premature stop codons, transversions, mutations in start or stop codons and splice site variants, as well as transitions and transversions were called for all samples. (B) Median mutations per megabase of 3.24 (range 0.038–11.46) in tumors from patients who progressed at or after 6 months as compared to 0.45 (range 0.019–9.90) in those who progressed in less than 6 months (n = 25, Mann-Whitney p = 0.22). (C) Median expressed neoantigens in tumors from patients who progressed at or after 6 months was 1.32 (range 0.00–6.06) versus 0.29 (range 0.00–5.70) in those who progressed before 6 months (n = 25, Mann-Whitney p = 0.29).

https://doi.org/10.1371/journal.pmed.1002309.g002

When filtering to expressed mutations, we found a median of 0.79 (range 0.00–3.36) expressed mutations per megabase for patients with DCB and a median of 0.16 (range 0.00–3.34) expressed mutations per megabase for patients without DCB (n = 25, Mann-Whitney p = 0.26, S2B Fig). Consistent with the known importance of specific variant calling pipelines to 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 postprocessing (via Base Quality Score Recalibration [BQSR] and depth/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 antitumor T cell response and are more common with increased mutation load. After performing in silico HLA typing (Methods), we examined predicted neoantigens that were 8 to 11 amino acids in length, resulting from the missense mutations of patients treated with atezolizumab. There was no significant association between predicted neoantigens per megabase and either DCB or DCB-OS. Patients with DCB had a median 4.58 (range 0.037–39.48) predicted neoantigens per megabase, while patients without DCB had 1.35 (range 0.00–20.22) (n = 25, Mann-Whitney p = 0.55, S2C Fig and S2D Fig). Filtering of predicted neoantigens to focus only on those expressed in RNA (Methods) also demonstrated no significant association between expressed predicted neoantigens and clinical benefit with atezolizumab (n = 25, Mann-Whitney p = 0.29, 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 PFS 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 reduced mortality or disease progression more than 3 months after treatment (n = 11, hazard ratio (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 those who survived less than or equal to 3 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 PFS 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 2 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 prespecified 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 (S3C Fig; see S3 Text and S3D Fig for further model interrogation). That said, we note that the 95% confidence intervals around the hazard ratios prior to 3 months include 1 (all p-values > 0.05; aggregate p(HR > 1) = 0.21), while those after 3 months are significantly less than 1 (all p-values < 0.05; aggregate p(HR > 1) = 0.020; p = 0.07 for interaction comparing aggregate HRs). This suggests that our threshold of 3 months, which was selected based on clinical experience, may be a convenient summary of these results.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Time-dependent relationship between mutation load and treatment response.

(A) There was a stronger association between somatic mutation load and progression-free survival (PFS) for events occurring more than 3 months following therapy (blue box: hazard ratio [HR] = 0.69, 95% CI [0.38, 0.99]), as compared to those in the first 3 months (red box: HR = 0.91, 95% CI [0.75, 1.07]). (B) There was a stronger association between somatic mutation load and overall survival (OS) for events occurring more than 3 months following therapy (blue box: HR = 0.80, 95% CI [0.60, 1.00]), as compared to those in the first 3 months (red box: HR = 1.02, 95% CI [0.79, 1.22]). (C) Kaplan-Meier estimates for PFS among patients with missense single nucleotide variants (SNVs) per megabase above the median observed value of 1.03 (red) and those with counts below the median value (blue). Time (Days) is plotted on a log-scale. For context, the frequency of observed events (progression and/or mortality) is plotted below the x-axis among patients with missense SNV per megabase above (red) and below (blue) the median value, with the approximate schedule of follow-up scans per the study protocol (see Methods) shown as vertical dotted lines.

https://doi.org/10.1371/journal.pmed.1002309.g003

These data suggest that in patients with rapidly progressive disease, factors other than mutation load likely determine their outcome. This observation is not surprising in that clinical factor analysis of this disease state has identified a heterogeneous population of patients, with 5 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 higher 5-factor score (3.00 (range 2.00–4.00) in contrast to 1.50 (range 0.00–4.00) in patients who survived longer than 3 months (n = 26, Mann-Whitney p = 0.018, S3A Fig). Patients surviving less than 3 months were much more likely to have liver metastases: 100% in patients surviving less than or equal to 3 months and 22% in patients surviving longer than 3 months (n = 29, Fisher's Exact p = 0.00097, S3B Fig). There were no significant differences in these patients with respect to Bacillus Calmette–Guérin (BCG) exposure (n = 29, Fisher's Exact p = 0.20), missense SNV load (n = 25, Mann-Whitney p = 0.26), and pretreatment peripheral TCR clonality (n = 29, Mann-Whitney p = 0.12). Keeping in mind the limited sample size of this cohort, these data suggest that there is a subset of nearly end-stage patients with cancer in whom clinical variables may negate immunological response, despite the presence of 1 or more favorable tumor-associated biomarkers. The inclusion of these clinical variables is warranted in future studies.

Examination of the tumor microenvironment shows evidence for adaptive immunity and suppression in responding tumors

Several studies have suggested that an “inflamed” tumor microenvironment, tumor, or IC PD-L1 expression increase the likelihood of response to checkpoint blockade. As seen in the published IMVigor 210 cohort, PD-L1 IC expression was significantly associated with DCB in this subset (n = 29, Spearman rho = 0.48 p = 0.0083, S4A Fig). We quantified immune infiltration from RNA-seq using ESTIMATE [19]. The immune score, while associated with the TIL proportion estimated through TCR-seq (S4B Fig), was estimated to be 764.37 (range −1195.08 to 1509.65) in patients with DCB and 263.49 (range −1100.78 to 1734.28) in patients without DCB but was not significantly different (n = 26, Mann-Whitney p = 0.33, S4C Fig). When we performed GSEA using the Hallmark Geneset [20], we did not observe any differentially expressed gene sets between tumors from patients with DCB versus no DCB. Furthermore, RNA expression of PD-L1 did not correlate with reported IC PD-L1 staining level (n = 26, Spearman rho = 0.045 p = 0.83, S4D Fig). We did not observe a difference in tumor MHC class I expression according to DCB (S4E Fig, HLA-A: n = 26, Mann-Whitney p = 0.26, HLA-B: n = 26, Mann-Whitney p = 0.36, HLA-C: n = 26, Mann-Whitney p = 0.24).

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 up-regulation 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), HAVCR2 exhibited significantly higher expression in tumors from patients who experienced DCB than in those who did not (S4F Fig). Interestingly, of the 4 IC2 tumors among HAVCR2-high patients (HAVCR2 expression greater than the median), 2 had missense SNV loads at or below the median (2 and 57); the other 2 had 180 and 412 SNVs. Additionally, although Rosenberg and colleagues [2] found that luminal cluster II showed a significantly higher response rate among the 4 subtypes of RNA expression from The Cancer Genome Atlas (TCGA), no significant association was found here between the 4 clusters and DCB (n = 20, Fisher's Exact p = 0.36) (S4G Fig) nor between the luminal/basal subcategorization and DCB (n = 20, Fisher's Exact p = 1.00), possibly due to sample size.

Relative importance of somatic, immune, and clinical factors in resistance and response to PD-L1 blockade

Unanswered questions that arise from the many studies of biomarker correlates of checkpoint blockade response are whether measures such as mutation load, PD-L1 staining, and others reflect the same “tumor state” or if each confers an independent effect on outcome?

When examined in conjunction with mutation load, the greater the expression of PD-L1, the more negative the association of mutation load with hazard (i.e., higher mutation load was associated with longer survival). Among patients with tumors showing little-to-no expression of PD-L1 (IC0 rated), each unit increase in missense SNV count per megabase was associated with a negligible change in hazard (n = 4, HR = 1.43, 95% CI [0.75, 2.98]). Among patients with tumors expressing PD-L1 at moderate or high levels (IC1 or IC2 staining), missense SNV count per megabase was associated with lower risk for disease progression or mortality (among IC1: n = 11, HR = 0.75, 95% CI [0.47, 1.14]; among IC2: n = 10, HR = 0.73, 95% CI [0.48, 1.06]). 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 the association of somatic mutation load (missense SNV count per megabase) and PFS, 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 Fig) than patients with liver metastasis (n = 9, HR = 0.96, 95% CI [0.66, 1.37], S5B Fig).

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, pretreatment 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, Fig 4C and 4D).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Associations between measured somatic, immune, and clinical variables.

(A) Although both programmed death-ligand 1 (PD-L1) staining and mutation load each were weakly associated with response, these variables were not correlated with each other (n = 25, Spearman rho = 0.14 p = 0.51). (B) Pretreatment peripheral T cell receptor (TCR) clonality did not correlate with mutation load (n = 25, Pearson r = 0.0017 p = 0.99). (C) Tumor-infiltrating T lymphocytes (TIL) proportion as estimated by TCR sequencing was associated with PD-L1 immune cell (IC) staining (n = 24, Spearman rho = 0.51 p = 0.010). (D) TIL clonality was associated with PD-L1 IC staining (n = 24, Spearman rho = 0.48 p = 0.017). (E) Hazard associated with log(pretreatment peripheral TCR clonality) by level of IC PD-L1 expression (IC0, IC1, or IC2). (F) Multivariate survival analysis of various clinical, peripheral, and intratumoral biomarkers for association with time to disease progression or mortality, utilizing a varying-coefficient model, which allows the hazard associated with a 1-unit increase in a biomarker’s value to vary according to the level of intratumoral PD-L1 expression (IC score).

https://doi.org/10.1371/journal.pmed.1002309.g004

In an analysis to see whether the association between pretreatment peripheral TCR clonality and PFS 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 pretreatment peripheral TCR clonality and PFS (among IC0: n = 4, HR [mean] = 1.86, HR [median] = 1.55, 95% CI [0.50, 4.99], p(HR > 1) = 0.21; among IC1: n = 11, HR [mean] = 0.69, HR [median] = 0.58, 95% CI [0.15, 1.84], p(HR < 1) = 0.19). Among patients with high levels of PD-L1 expression, by comparison, we observed almost complete separation of PFS according to pretreatment peripheral TCR clonality (among IC2: n = 10, HR [mean] = 89.88, HR [median] = 23.41, 95% CI [2.43, 506.94], p(HR>1) = 0.0014; Fig 4E). Similar results were seen in analyses with respect to OS and in a logistic regression analysis for DCB (S5C Fig, S5D Fig, S5E Fig).

To resolve the hypothesis that those patients with low peripheral TCR clonality simply were healthier, we examined the association between 5-factor score and pretreatment 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 IC PD-L1 expression (S2 Table).