. Author manuscript; available in PMC: 2020 Jan 29.

Published in final edited form as: Nat Med. 2019 Jul 29;25(8):1251–1259. doi: 10.1038/s41591-019-0522-3

Clonal replacement of tumor-specific T cells following PD-1 blockade

Kathryn E Yost ^1,^#, Ansuman T Satpathy ^1,^2,^3,^#,^*, Daniel K Wells ³, Yanyan Qi ¹, Chunlin Wang ⁴, Robin Kageyama ³, Katherine McNamara ^5,^6,⁷, Jeffrey M Granja ^1,^6,⁸, Kavita Y Sarin ⁹, Ryanne A Brown ^2,⁹, Rohit K Gupta ¹⁰, Christina Curtis ^5,^6,⁷, Samantha L Bucktrout ³, Mark M Davis ^3,^11,^12,¹³, Anne Lynn S Chang ^9,^*, Howard Y Chang ^1,^3,^6,^9,^13,^*

PMCID: PMC6689255 NIHMSID: NIHMS1531727 PMID: 31359002

Abstract

Immunotherapies that block inhibitory checkpoint receptors on T cells have transformed the clinical care of cancer patients¹. However, whether the T cell response to checkpoint blockade relies on reinvigoration of pre-existing tumor infiltrating T cells (TILs) or on recruitment of novel T cells remains unclear^2–4. Here, we performed paired single-cell RNA (scRNA) and T cell receptor (TCR)- sequencing on 79,046 cells from site-matched tumors from patients with basal cell carcinoma (BCC) or squamous cell carcinoma (SCC) pre- and post-anti-PD-1 therapy. Tracking TCR clones and transcriptional phenotypes revealed a coupling of tumor-recognition, clonal expansion, and T cell dysfunction marked by clonal expansions of CD8⁺CD39⁺ T cells, which co-expressed markers of chronic T cell activation and exhaustion. However, this expansion did not derive from pre-existing TIL clones; rather, it was comprised of novel clonotypes not previously observed in the same tumor. Clonal replacement of T cells was preferentially observed in exhausted CD8⁺ T cells and evident in BCC and SCC patients. These results demonstrate that pre-existing tumor-specific T cells may have limited reinvigoration capacity, and that the T cell response to checkpoint blockade derives from a distinct repertoire of T cell clones that may have just recently entered the tumor.

Main

We generated droplet-based 5’ scRNA- and TCR-seq libraries from 11 patients with advanced BCC before and after anti-PD-1 treatment in site-matched primary tumors (Fig. 1a, Supplementary Table 1, Methods). CD3 immunohistochemistry (IHC) and whole exome sequencing (WES) supported an immunological response to checkpoint blockade, including increased CD3⁺ T cell infiltration (Fig. 1b) and mutational loss following treatment affecting both clonal and sub-clonal mutations and neoepitopes, suggestive of tumor immunoediting (Fig. 1c, Extended Data Fig. 1a–c, Supplementary Tables 1–3)⁵.

Figure 1. — **(a)** Workflow for sample processing and scRNA-seq analysis of advanced BCC samples collected pre- and post-PD-1 blockade. Graphics courtesy of the Parker Institute for Cancer Immunotherapy. **(b)** Immunohistochemistry staining for CD3⁺ cells in representative BCC tumors before and after PD-1 blockade. Tumor boundaries denoted with dashed lines. All scale bars represent 100 μm. IHC staining was performed once for each sample (n = 16 samples). **(c)** Bar plot of neoepitope burden pre- and post-treatment based on exome sequencing. Variants were classified as predicted neoepitopes if the peptide was found to bind to the MHC allele with less than 500 nM binding strength and its wildtype cognate bound to the same allele with greater than 500 nM binding strength. **(d)** UMAP of all tumor-resident cells pre- and post-therapy for all 11 BCC patients. Clusters denoted by color are labeled with inferred cell types, which include 2 malignant clusters, 2 CD4⁺ T cell clusters, 3 CD8⁺ T cell clusters, and proliferating T cells, endothelial cells, melanocytes, myofibroblasts, and cancer-associated fibroblasts (CAFs), dendritic cells (DCs), macrophages, and plasmacytoid dendritic cells (pDCs), 3 B cell clusters, and 1 NK cell cluster. **(e)** UMAP of tumor-resident cells colored by patient identity (top left), FACS sort markers (top right), anti-PD1 treatment status (bottom left), and TCR detection (bottom right). **(f)** Inferred CNV profiles based on scRNA-seq data. Non-immune, non-malignant cells (fibroblasts and endothelial cells, n = 2,122) were used as normal reference for malignant cell CNV inference (n = 3,548). **(g)** Representative examples of hematoxylin and eosin (H&E) staining of different BCC subtypes. All scale bars represent 100 μm. H&E staining was performed once for each sample (n = 9 samples). **(h)** UMAP of malignant cells colored by patient (left) and clinical subtype (right). **(i)** UMAP of malignant cells colored by enrichment of basal and squamous cell carcinoma gene signatures (from Atwood *et al*., 2015 and Hoang *et al*., 2017) (top). Malignant cells ordered based on the difference between basal and squamous signatures (bottom). The clinical diagnosis associated with each cell and expression of signature associated genes are indicated below.

In total, we obtained scRNA-seq profiles from 53,030 malignant, immune, and stromal cells with paired TCR sequences in 28,371/33,106 T cells (86%; Fig. 1d–e, Extended Data Fig. 1d, Methods). We identified 19 cell clusters based on scRNA-seq profiles, including 2 malignant clusters, 6 T cell clusters, 4 stromal clusters, 3 myeloid clusters, 3 B cell clusters, and 1 NK cell cluster. Immune cell classifications were consistent with surface markers used to isolate cells and bulk RNA-seq from reference populations⁶ (Fig. 1e, Extended Data Fig. 2a–d). Notably, immune cells from different patients clustered together, indicating consistency in TME immune cell types across patients. Single-cell copy number variation (CNV) estimation revealed patient-specific CNVs only in malignant cells, which were consistent with CNVs detected by WES and previously described CNVs in BCC⁷ (Fig. 1f, Extended Data Fig. 3a, Methods).

Re-clustering of 3,548 malignant cells revealed clustering by patient and BCC subtype (Fig. 1g–h), indicating significant intertumoral heterogeneity as observed in other cancers^8,9. We identified a shared malignant gene expression program that included EPCAM, BCAM and TP63, known markers of BCC (Extended Data Fig. 3b)^10–12. We identified 577 genes differentially expressed between patient tumors, which included Ras signaling genes, suggesting aberrant squamous cell pathway activation in BCC, as previously reported¹³. Scoring of malignant cells for enrichment of BCC and SCC expression signatures from bulk RNA-seq^14,15 revealed a differentiation continuum with basal signature enrichment in nodular BCCs and squamous signature enrichment in infiltrative and metatypical BCCs (Fig. 1i). Altogether, these results demonstrate that gene expression in BCC is driven by patient-specific malignant pathways, but largely does not influence immune cell phenotypes in the TME.

We next focused on TILs to understand the clonal T cell response to checkpoint blockade. First, we re-clustered 33,106 TILs and identified 9 distinct T cell clusters containing cells from multiple patients and pre- and post-treatment timepoints (Methods, Fig. 2a–b, Extended Data Fig. 4a–d). CD4⁺ clusters included regulatory T cells (Tregs), follicular helper T cells (Tfh), and T helper 17 cells (Th17). CD8⁺ clusters included naïve cells, memory cells, effector memory cells, activated cells, chronically activated/exhausted cells, hereafter referred to as exhausted cells, and intermediate ‘exhausted/activated’ cells, which co-expressed activation and exhaustion associated genes³. Notably, we observed an increased frequency of Tfh cells and activated, exhausted, and exhausted/activated CD8⁺ T cells post-treatment, supporting reports that PD-1 blockade primarily impacts CD8⁺ T cells (Fig. 2b)^16,17.

Figure 2. — **(a)** UMAP of tumor-infiltrating T cells present in BCC samples pre- and post-PD-1 blockade. Clusters denoted by color labeled with inferred cell types (left). UMAP also colored by patient (top right) and anti-PD-1 treatment status (bottom right). **(b)** Heatmap of differentially expressed genes (rows) between cells belonging to different T cell subsets (columns). Specific genes associated with different T cell clusters are highlighted. Bars at top of heatmap indicate the number of cells, post-therapy enrichment, and number of patients in each cluster. **(c)** Diffusion map of naïve, memory, activated and exhausted CD8⁺ T cells using the first two diffusion components (left). Cells colored based on cluster identities from Fig. 2a. Cells are also colored by diffusion pseudotime and treatment status (top right). Average expression of selected core activation and exhaustion genes is quantified along diffusion components 1 and 2 (bottom right). **(d)** Co-expression analysis of differentially expressed genes (n = 146 genes) between activated, exhausted and activated/exhausted CD8⁺ T cells (n = 5454 cells). Inset indicates core exhaustion module identified by hierarchical clustering, with canonical exhaustion genes highlighted. **(e)** Diffusion map of CD8⁺ T cell subsets colored by clone size (left) and boxplot of Gini indices for each CD8⁺ T cell cluster calculated for each patient (right), showing significant clonal expansion within exhausted CD8⁺ T cells (n = number of patients, unpaired t-test, one-tailed, relative to basemean; box center line, median; box limits, upper and lower quartiles; box whiskers, 1.58× interquartile range, here and throughout). Exhausted refers to both exhausted and exhausted/activated clusters. **(f)** Activation score (based on expression of top 50 genes most correlated with *IFNG* expression) versus exhaustion score (based on expression of top 50 genes most correlated with *HAVCR2* expression) for all CD8⁺ T cells (n = 17,561), colored by expression levels of indicated genes. **(g)** Activation score versus exhaustion score enrichment for TCR clones with >1 cell (n = 6,422) based on average activation and exhaustion scores of individual cells belonging to that clone, colored by the most frequent assigned phenotype for cells belonging to that clone, and size based on clone size (top right) or cell cycle score (bottom right).

We used diffusion maps to visualize the relationship between CD8⁺ T cell clusters and order cells in pseudotime (Methods, Fig. 2c)¹⁸. The first diffusion component separated activated and exhausted cells and was highly correlated with T cell exhaustion genes, including PDCD1 and HAVCR2, while the second diffusion component separated naïve and memory cells from activated and exhausted cells and was highly correlated with T cell activation genes, including IFNG and TNF (Fig. 2c, Extended Data Fig. 5a,b). We used co-expression analysis to identify a core T cell exhaustion signature in the context of checkpoint blockade, which included known exhaustion markers (HAVCR2, TIGIT), tissue resident memory T cell (Trm) markers (ITGAE, CXCR6)^19,20, and markers of tumor-reactive CD8⁺ TILs including CD39 (ENTPD1)^21–24 (Fig. 2d). These results suggest that exhausted CD8⁺ TILs increase after PD-1 blockade and express gene signatures of chronic activation, T cell dysfunction, and tumor reactivity.

Since tumor antigen-specific CD8⁺ T cells expand clonally during a productive immune response, we analyzed single-cell TCR sequences to identify clonally-expanded cells as an indicator of tumor-specificity. We grouped cells by TCRα and TCRβ sequence and noted large clone sizes in exhausted T cells compared to other CD8⁺ clusters (Fig. 2e). We measured clonality using Gini index and observed significantly higher clonality in exhausted T cells compared to all other CD8⁺ T cells (Fig. 2e). Analysis of CD4⁺ T cells demonstrated increased clonality in Tfh cells post-treatment, which in one patient was accompanied by an increase in B cells expressing germinal center markers (Extended Data Fig. 5c–g).

To examine the link between clonal expansion and exhaustion, we scored all CD8⁺ cells for activation and exhaustion signatures based on the top 50 genes correlated with IFNG and HAVCR2 expression, respectively (Methods, Fig. 2f). We found that T cells with a high exhaustion signature exhibited gene expression patterns associated with tumor reactivity, including CD39 (ENTPD1) and CD103 (ITGAE) expression, and absence of KLRG1 expression^21–24. To characterize individual clones, we averaged exhaustion and activation scores for all cells in a clone and observed high exhaustion gene signatures in the largest clones (Fig. 2g). Exhausted clones also exhibited a high proliferation signature, as reported in melanoma (Methods, Fig. 2g)²³.

We next analyzed lineage relationships between T cell phenotypes and clonotypes. Globally, we found that cells grouped by clonotype were more likely to share a common phenotype and were more correlated in gene expression than randomly grouped cells, in line with prior studies (Fig. 3a–c, Extended Data Fig. 6a–c)^25–27. We used GLIPH (grouping of lymphocyte interactions by paratope hotspots) to identify ‘TCR specificity groups,’ clusters of distinct TCR sequences that likely recognize common antigens via shared motifs in CDR3 sequence²⁸. T cells expressing distinct TCRs within a specificity group were more likely to share a common phenotype and exhibit highly correlated gene expression, compared to randomly grouped clones (Fig. 3a–c, Extended Data Fig. 6d). These results suggest that clonally-expanded TILs are highly correlated in cellular phenotype and that PD-1 blockade does not promote phenotype instability within a clone. Moreover, specificity group analysis suggests that antigen specificity also contributes to T cell fate.

Figure 3. — **(a)** UMAP of tumor-infiltrating T cells colored by selected TCR clones (left). UMAP of T cells colored by TCRβ clones belonging to the same TCR specificity (GLIPH) group (right). **(b)** Phenotypes of single cells belonging to the same TCR clone or TCR specificity group. Shown are the top five most abundant clones (top and middle) larger than 10 cells for each patient. Each bar is colored by individual phenotypes of single cells within the clone. The bottom plots show phenotypes of distinct TCR clones within a TCR specificity group. Both analyses show substantial phenotypic similarity among single cells belonging to a clone or group. **(c)** Distribution of the proportion of cells within each clone, or TCRβ clones within each TCR specificity group, (>=3 cells) that share a common cluster identity compared to randomly selected and size matched groups of T cells from the same sample (left, n = number of TCRβ clones or TCR specificity groups, unpaired t-test, two-tailed). Distribution of cell-cell correlations between cells that belong to the same TCR clone or cells within the same TCR specificity group but different clonotypes, compared to randomly selected and size matched groups of T cells from the same sample (bottom, n = number of cells or clones, unpaired t-test, two-tailed). Cell-cell correlations were calculated using log-transformed expression of differentially expressed genes. **(d)** Heatmap showing the fraction of clonotypes belonging to a primary phenotype cluster (rows) that are shared with other secondary phenotype clusters (columns). **(e)** Heatmap of all observed phenotype transitions for matched clones during PD-1 blockade for clones with at least 3 cells for each timepoint. **(f)** *TCF7*⁺/stem-like score (from Im *et al*., 2016) versus exhaustion score for all CD8⁺ T cells (n = 17,561), colored by expression of indicated genes (left). *TCF7*⁺/stem-like score versus exhaustion score for cells belonging to primarily exhausted clones, colored by phenotype (top right). Violin plot of *TCF7*⁺/stem-like score for memory and exhausted clones separated by change in clone abundance following treatment (bottom right, n = number of clones, unpaired t-test, two-tailed). Clones were defined as expanded or contracted if they significantly changed in abundance by a Fisher exact test (P < 0.05 and fold change > 0.5), and persistent if they did not significantly change in abundance and at least one cell was detected at each timepoint. Exhausted refers to both exhausted and exhausted/activated clusters. **(g)** Pie charts showing clone size and distribution of phenotypes for matched clones pre- and post-therapy. Selected clones had a primarily exhausted phenotype pre-therapy and increased in abundance post-therapy, and are separated by the presence of a high *TCF7*⁺ signature prior to treatment. **(h)** Pie chart of all clones with an exhausted phenotype post-treatment, colored by whether the clone contained *TCF7* high cells pre-therapy.

We asked whether divergent phenotypes within clones could inform lineage transitions between T cell phenotypes. We aggregated all clonotypes in a given cluster (‘primary phenotype’) and measured the fraction shared with another cluster (‘secondary phenotype’) (Fig. 3d). Broadly, we noted significant overlaps between CD8⁺ T cell phenotypes, including memory and activated T cells, suggesting common transitions between activation states. We detected minimal clonotype sharing between exhausted and effector cells, supporting a strict bifurcation between these phenotypes²³. CD4⁺ T cell clones were largely restricted to single phenotypes, suggesting limited plasticity between CD4⁺ cell states. We also observed a non-random distribution of phenotypes of individual clones within specificity groups, suggesting that specific T cell phenotypes may result from different TCR signal strength thresholds (Extended Data Fig. 6e). We noted overlaps between CD8⁺ and CD4⁺ phenotypes within specificity groups, such as specificity groups containing CD8⁺ exhausted, and CD4⁺ Treg and Tfh clones, suggesting that CD4⁺ and CD8⁺ TILs responding to the same antigen may arise from distinct clonotypes.

To track clonal cell fates following PD-1 blockade, we matched clonotypes between treatment timepoints based on TCRβ sequences and compared the primary phenotypes at each timepoint for all matched clones (Fig. 3e, Extended Data Fig. 7b). We observed stability among CD4⁺ clusters and frequent transitions among CD8⁺ clusters, similar to clonotype sharing in individual timepoints. While we observed frequent transitions between memory and effector to activated states, pre-treatment exhausted clones did not transition to non-exhausted phenotypes post-treatment, suggesting that exhausted TILs have limited capacity for phenotype transition after PD-1 blockade.

Prior studies have identified stem-like T cells expressing the transcription factor TCF7 that proliferate following PD-1 blockade^17,29–32. To ask whether similar cells existed in our dataset, we scored CD8⁺ cells for exhaustion or TCF7⁺ stem-like signatures (Fig. 3f). We observed a small population of cells (28% of exhausted cells, 1.5% of CD8⁺ T cells) with both TCF7⁺ and exhaustion signatures. We found that for both memory and exhausted phenotypes, persistent clones had a significantly higher TCF7⁺ signature pre-treatment compared to clones that contracted (Fig. 3f). However, this analysis was limited to only two exhausted clones that significantly expanded, prompting us to identify 10 exhausted clones that increased in frequency post-therapy but were previously excluded due to low clone size and limited expansion (Fig. 3g). The majority of these clones remained exhausted, but those with a high TCF7⁺ signature pre-treatment expanded more substantially. Nevertheless, only a small fraction (10.3%) of post-therapy exhausted clones were derived from clones containing TCF7⁺ cells, suggesting that post-treatment exhausted clones may be derived from additional sources (Fig. 3h).

Since few pre-existing exhausted T cells showed expansion post-therapy, we asked how clone abundance changed globally following treatment by comparing pre- and post-treatment frequencies of each clone (Fig. 4a). We identified significantly expanded clones post-therapy, many of which were not detected prior to treatment (68% of significantly expanded clones). Integration of scRNA-seq data revealed strikingly different persistence patterns for each phenotype. Namely, post-therapy exhausted clones were significantly enriched for novel clonotypes; on average 84% of exhausted clones were derived from novel clonotypes for each patient, compared to only 40% of naïve, activated, memory, or effector memory clones (Fig. 4a–b). Across all patients, 55% of post-treatment exhausted clones containing at least 5 cells were derived from novel clonotypes, compared to 20% of memory clones (Fig. 4c). We next asked how expansion of novel clones, a phenomenon we termed ‘clonal replacement,’ contributed to exhausted T cell frequency in each patient (Fig. 4d, Extended Data Fig. 8a). We found that 7/11 patients had an increased exhausted CD8⁺ T cell frequency following treatment, and in 6/7 patients, the majority were derived from novel clonotypes (Fig. 4d). Interestingly, post-treatment exhausted clones were enriched for novel TCR specificity groups, suggesting that novel clones may represent new antigen specificities (Extended Data Fig. 8b).

Figure 4. — **(a)** Scatterplots comparing TCRβ clone frequencies pre- and post-treatment measured by scRNA+TCR-seq for all BCC patients (n = 11 patients). Clones that were significantly expanded or contracted post-treatment based on a Fisher exact test (P < 0.05) are highlighted on the left. Clones where the majority of cells exhibit an exhausted CD8⁺ phenotype (middle, red) or a memory CD8⁺ phenotype (right, blue) highlighted. In this and subsequent panels, exhausted refers to both exhausted and exhausted/activated clusters. **(b)** Boxplot of the fraction of novel clones detected by scRNA+TCR-seq within each cluster following treatment (n = number of patients, unpaired t-test, two-tailed). Clones with only one cell detected and cells from su003 with no clonotype overlap between timepoints were excluded. **(c)** Lorenz curve of TCRβ clone frequencies based on scRNA+TCR-seq for exhausted CD8⁺ T cell clones (left) and memory CD8⁺ T cell clones (middle) greater than or equal to 5 cells, colored by presence of each clone prior to treatment. Proportion of novel clones in each phenotype quantified on the right. **(d)** Fraction of exhausted cells out of total T cells detected by single cell RNA+TCR-seq for each patient, separated by treatment status. Cells belonging to novel clones detected post-treatment are highlighted. **(e)** Scatterplots comparing TCRβ clone frequencies pre- and post-treatment measured by bulk TCR-seq (n = 8 patients). Clones that were significantly expanded or contracted post-treatment based on a binomial test (two-sided, Bonferroni corrected P value < 0.01) are highlighted on the left, with expanded clones further separated based on their detection pre-treatment. Clones where the majority of cells share an exhausted CD8⁺ phenotype based on scRNA-seq (middle, red) or a memory CD8⁺ phenotype (right, blue) highlighted. **(f)** Bar plot of fraction of clones with significant expansion post-treatment based on bulk TCR-seq, separated by phenotype and colored by replacement status. **(g)** Overlap between TCRβ clones in peripheral blood and tumor infiltrating T cells detected by bulk TCR-seq (n = 5 patients, left). Fraction of TIL clones detected in peripheral blood, separated by sample (top right). Fraction of novel exhausted TIL clones detected in PBMCs, separated by treatment status (bottom right). **(h)** Violin plot of clone enrichment (tumor frequency / PBMC frequency) detected by bulk TCR-seq, separated by phenotype and treatment status (data from 5 patients, n = number of clones, unpaired t-test, one-tailed). **(i)** Characteristics of squamous cell carcinoma (SCC) samples treated with anti-PD-1 (left) and UMAP of tumor-infiltrating T cells present in SCC samples pre- and post-PD-1 blockade (right). Clusters denoted by color are labeled with inferred cell types. Graphics courtesy of the Parker Institute for Cancer Immunotherapy. **(j)** Fraction of exhausted cells out of total T cells detected by single-cell RNA+TCR-seq for each patient, separated by treatment status. Novel clones detected post-treatment are highlighted (bottom left). Sample su010-S derived from an SCC lesion from patient su010 who presented with both BCC and SCC lesions. **(k)** Bar plot of fraction of clones with significant expansion based on bulk TCR sequencing post-treatment, separated by phenotype and colored by replacement status (bottom right).

To increase sensitivity for detecting rare clonotypes, we performed bulk TCR-seq on remaining biopsy material from 8/11 patients (Methods, Fig. 4e). Similar to scTCR-seq analysis, we observed a substantial number of novel expanded clones that were either not expanded (<5 cells) or undetected pre-treatment. Compared to all other CD8⁺ phenotypes, exhausted cells had a higher proportion of significantly expanded clones following treatment, and the majority of expanded clones were derived from novel clonotypes, which was consistently observed across patients with expanded exhausted cells post-therapy (Fig. 4f, Extended Data Fig. 8c). To address whether this resulted from sampling bias, we performed bulk TCR-seq on site-matched samples from one patient twice pre-therapy and twice post-therapy at approximately two-month intervals (Extended Data Fig. 8d). We only observed clonal replacement of exhausted clones when comparing pre- to post-treatment samples, suggesting that TCR dynamics of exhausted cells were mainly influenced by PD-1 blockade, not tumor biopsy timing or location.

To determine whether novel clonally-expanded TILs could be detected in peripheral blood, we performed bulk TCR-seq on 10 blood samples from 5 patients. Overall, 41% of TIL TCR clonotypes could also be detected in blood, but only represented 6% of blood clonotype diversity (Fig. 4g, Extended Data Fig. 9a). Importantly, blood clonotypes represented all TIL phenotypes, including novel exhausted TILs (Extended Data Fig. 9a,b). Overall, 35.5% of novel exhausted TIL clonotypes were detected in the peripheral blood post-treatment, and surprisingly, 11.8% of novel exhausted TIL clonotypes were detected in peripheral blood pre-treatment, despite being undetectable by deep TCR sequencing in the tumor pre-treatment (Fig. 4g). We compared clonotype enrichment in the tumor over peripheral blood by comparing clonotype frequency in each location. We noted a significant increase in the enrichment of exhausted clones and specificity groups relative to other phenotypes post-treatment, suggesting preferential expansion and retention in the tumor, supporting their tumor-specificity (Fig. 4h, Extended Data Fig. 9c). These results suggest that it may be feasible to monitor the clonal tumor-specific T cell response to checkpoint blockade in the blood and that novel TIL clones may be recruited from peripheral sources^33,34.

Finally, we asked whether clonal replacement of exhausted cells could be observed in a different cancer type. We generated scRNA+TCR-seq profiles of 26,016 TILs from serial tumor biopsies in 4 patients with SCC treated with anti-PD-1 (Fig. 4i, Extended Data Fig. 10a). SCC samples were obtained on average 31 days post-treatment, enabling analysis of TIL dynamics relatively early after treatment. We first confirmed our prior findings regarding: 1) TIL phenotypes, including clonal expansion of exhausted CD8⁺ T cells which expressed tumor-specificity markers, including CD39^21–24 (Extended Data Fig. 10b–d) and 2) TIL clonotype dynamics including stability of clone phenotypes (Extended Data Fig. 10e–i). Next, we observed that a considerable proportion of exhausted cells detected post-treatment were derived from novel clonotypes (Fig. 4j). Integration of scRNA-seq data with bulk TCR-seq confirmed clonal replacement preferentially in exhausted T cells compared to other phenotypes, with overall 50% of expanded exhausted clones derived from novel clonotypes compared to only 29% of other expanded CD8⁺ clones (Extended Data Fig. 10j, Fig. 4k). Notably, we observed a similar degree of expansion of pre-existing clones at early timepoints in SCC as we did at later timepoints in BCC.

Here, we performed single-cell profiling of clinical tumor biopsies including integration of TCR clonotype and scRNA-seq phenotype which revealed that: 1) clonally-expanded cells were highly enriched in exhausted CD8⁺ T cells and expressed markers of tumor specificity, including CD39 and CD103^21–24, and 2) the clonal repertoire of exhausted CD8⁺ T cells was largely replaced by novel clones after therapy compared to other phenotypes. These results suggest that chronic activation and exhaustion of pre-existing TILs limits their re-invigoration following checkpoint blockade⁴, and that the T cell response to immunotherapy derives from a distinct repertoire of tumor-specific T cell clones.

Clonal replacement of tumor-specific T cells is consistent with several prior findings in the context of PD-1 blockade, including limited reinvigoration of exhausted T cells due to broad epigenetic remodeling^4,35, proliferation of CXCR5⁺CD8⁺ T cells in lymphoid organs but not other tissues²⁹, and loss of anti-tumor T cell responses following chemical inhibition of T cell migration³⁶. Importantly, our study did not identify the source of novel T cell clones, and clones could derive from tumor-extrinsic sources, including lymphoid organs, or rare unexpanded clones within the TME or tumor periphery. Although our bulk TCR sequencing supports the first possibility, further work will be required to identify the source of novel T cell clones and their impact on clinical response. Importantly, both possibilities are compatible with the potential derivation of these cells from a TCF7⁺ precursor population^17,29–32.

Furthermore, expansion of novel TCR clones and specificity groups following PD-1 blockade, coupled with neoepitope loss, suggests that novel T cell clones may initiate a distinct immuno-editing wave³⁷. The antigen identities recognized by each wave require further investigation, perhaps using high-throughput tumor specificity assays^38,39. Finally, our results suggest that improved checkpoint blockade activity in immune-infiltrated (‘hot’) vs immune-desert (‘cold’) tumors may result from an intrinsic ability to constantly attract new T cells⁴⁰, rather than reactivation of pre-existing TILs. In summary, this study reveals insights into the clonal T cell response to checkpoint blockade in human cancer, which has important implications for the design of checkpoint blockade immunotherapies.

Methods

Human subjects

This study was approved by the Stanford University Administrative Panels on Human Subjects in Medical Research, and written informed consent was obtained from all participants. All patients had histologically-proven advanced or metastatic BCC or SCC and were not good candidates for surgical resection⁴¹. Exclusion criteria included 1) prior exposure to checkpoint blockade agents, 2) systemic immunosuppressant use, treatment with radiotherapy or other anti-cancer agents within 4 weeks of first biopsy⁴¹. Patients were treated with 200 mg pembrolizumab every 3 weeks or 350 mg cemiplimab every 2 weeks. A subset of patients received ongoing treatment with 150 mg vismodegib daily (Supplementary Table 1). Response was assessed by RECIST version 1.1⁴².

Sample Collection and Processing

Fresh biopsies were collected from the primary tumor site. A portion of the tumor was stored in RNAlater for whole exome sequencing and bulk TCR sequencing. The remaining tissue was processed for single cell RNA sequencing.

H&E and Immunohistochemical (IHC) staining

For H&E staining, formalin-fixed, paraffin-embedded tissue cut at 4 microns and stained using the Tissue-Tek Prisma automated slide stainer. Immunostaining was performed on the Ventana Benchmark Ultra platform (CD3) and the Leica Bond platform (CD8, PD-L1). Antibodies used include anti-human CD3 (cat. no. 103A-76, Cell Marque), anti-human CD8 (cat. no. M7103, Dako) and anti-human PD-L1 (cat. no. 13684S, Cell Signaling Technology).

Exome Sequencing

Exome capture was performed by Accuracy and Content Enhanced (ACE) augmented exome strategy (Personalis) and sequenced on an Illumina HiSeq 2500 with paired-end 100-bp sequencing, with an average of 110-fold coverage (range 82–138).

HLA Typing

All samples were genotyped for HLA-A, -B, -C, -DPA1, -DPB1, -DQA1, -DQB1, -DRB1 and -DRB3/4/5 loci using the MIA FORA NGS FLEX HLA Typing 11 Kit 96 Tests (Immucor, Inc. Norcross, GA, USA), following manufacturer’s semi-automated protocol and as described previously^43,44. Briefly, paired-end sequence reads were generated using an NGS-based HLA genotyping method targeting 11 HLA genes with extensive coverage of the HLA genomic region by long-range polymerase chain reaction (PCR). Coverage for HLA class I loci is >200 bp 5’UTR to 3’UTR ~200–400 bp. In the case of HLA-DQA1 locus is ~200 bp of the 5’UTR to ~200 bp of the 3’UTR and for HLA-DQB1 locus is ~70bp of the 5’UTR to ~100 bp of the 3’UTR. For the remaining class II loci specific key regions of the gene were amplified. For HLA-DPA1 locus, this coverage is from exon 1 to exon 4 and for HLA-DPB1 locus from exon 2 to exon 4. All HLA-DRB1/3/4/5 genes were co-amplified in two separate reactions. The coverage for HLA-DRB1/3/4 loci included ~300–500 bp of the 5’UTR to the first ~270 bp of intron-1 and the 3’ end of intron-1 (~250 bp) to exon-6. For the HLA-DRB5 gene exon 2 to exon 6 were amplified. All libraries were sequenced on an Illumina MiSeq. For assignment of HLA genotypes, NGS paired-end reads were analyzed using the MIA FORA FLEX v3.0 software (Immucor, Inc.). Final HLA genotyping calls were confirmed by manual review.

Bulk TCR Sequencing

Deep sequencing of the TCRβ gene was performed using the immunoSEQ platform (Adaptive Biotechnologies) on genomic DNA extracted from tumor biopsies or peripheral blood with input amounts ranging from 616 ng to 6,004 ng. Only data from productive rearrangements was exported from the immunoSEQ Analyzer for further analysis. On average, 26,066 TCR templates were detected from tumor samples (range 554–99,264), representing an average of 6,041 unique clonotypes (range 237–17,181), ~20-fold increase in sampling depth compared to scTCR-seq. For peripheral blood samples, on average 113,528 TCR templates were detected (range 24,679–257,772), representing an average of 36,536 unique clonotypes (range 7,041–71,462).

Tumor Dissociation

Fresh tumor biopsies were minced and digested in 5 mL digestion media (DMEM/F12 media + 250 μg/mL Liberase TL + 200 U/mL DNAse I) in a C-tube using the gentleMACS Octo system at 37°C for 3 hours at 20 rpm. Following digestion, 50 μL of 500 mM EDTA was added and sample collected by centrifugation at 300xg for 5 minutes. The cell suspension was then passed through a 70 μm filter and pelleted by centrifugation at 300xg at 4°C for 10 minutes. Cells were then resuspended in 1 mL of RPMI media and cryopreserved in FBS supplemented with 10% DMSO until further processing.

Cell Sorting

Cells were categorized as peri-tumoral T cells (CD45⁺CD3⁺), other peri-tumoral lymphocytes (CD45⁺CD3⁻) and malignant/stromal cells (CD45⁻CD3⁻). For patients su009, su010, su011, su012, su013, and su014, only peri-tumoral T cells were isolated and used for scRNA-seq. For sample su010-S, we additionally isolated CD8⁺CD39⁺ peri-tumoral T cells (CD45⁺CD3⁺CD8⁺CD39⁺). For bulk RNA-seq datasets, CD4⁺ T helper cells were sorted as naive T cells (CD4⁺CD25⁻CD45RA⁺), Treg (CD4⁺CD25⁺IL7R^lo), Th1 (CD4⁺CD25⁻IL7R^hiCD45RA⁻CXCR3⁺CCR6⁻), Th2 (CD4⁺CD25⁻IL7R^hiCD45RA⁻CXCR3⁻CCR6⁻), Th17 (CD4⁺CD25⁻IL7R^hiCD45RA⁻CXCR3⁻CCR6⁺), Th1–17 (CD4⁺CD25⁻IL7R^hiCD45RA⁻CXCR3⁺CCR6⁺), and Tfh subsets (CXCR5⁺ counterparts of each). Antibodies used included anti-human-CD45 conjugated to V500 (clone HI30, cat. no. 560779, lot no. 7172744, BD Biosciences), anti-human-CD3 conjugated to FITC (clone OKT3, cat. no. 11–0037-41, lot no. 2007722, Invitrogen), anti-human-CD8 conjugated to Pacific Blue (clone 3B5, cat. no. MHCD0828, lot no. 1964935, Invitrogen), anti-human-CD39 conjugated to APC (clone A1, cat. no. 328210, lot no. B268898, BioLegend), anti-human-PD-1 conjugated to APC/Cy7 (clone EH12.2H7, cat. no. 329921, lot no. B245235, BioLegend) and anti-human-HLA-DR conjugated to eVolve 605 (clone LN3, cat. no. 83–9956-41, lot no. 1949784, Affymetrix-Ebioscience), anti-human-CD45RA conjugated to PERCP-Cy5.5 (clone HI100, cat. no. 304107, lot no. B213966, BioLegend), anti-human-CD127 conjugated to Brilliant Violet 510 (clone A019D5, cat. no. 351331, lot no. B197159, BioLegend), anti-human-CD4 conjugated to APC/Cy7 (clone OKT4, cat. no. 317417, lot no. B207751, BioLegend), anti-human-CCR6 conjugated to PE (clone G034E3, cat. no. 353409, lot no. B203239, BioLegend), anti-human-CD25 conjugated to FITC (clone BC96, cat. no. 302603, lot no. B168869, BioLegend), anti-human-CXCR3 conjugated to Brilliant Violet 421 (clone G025H7, cat. no. 353715, lot no. B206003, BioLegend), anti-human-CXCR5 conjugated to Alexa-Fluor-647 (clone RF8B2, cat. no. 558113, lot no. 5302868, BD Pharmingen), and anti-human-CD3E conjugated to Pacific Blue (clone UCHT1, cat. no. 558117, lot no. 4341657, BD Biosciences). All antibodies were used at a 1:200 dilution, with the exception of anti-CD45 and anti-HLA-DR antibodies which were used at a 1:100 dilution. Propidium iodine (cat. no. P3566, Invitrogen) was used for live/dead staining at a final concentration of 2.5 μg/mL.

Single-cell RNA-seq Library Preparation

Single-cell RNA-seq and TCR-seq libraries were prepared using the 10X Single Cell Immune Profiling Solution Kit, according to the manufacturer’s instructions. Briefly, FACS sorted cells were washed once with PBS + 0.04% BSA and resuspended in PBS + 0.04% BSA to a final cell concentration of 100–800 cells/μL as determined by hemacytometer. Cells were captured in droplets at a targeted cell recovery of 500–7000 cells, resulting in estimated multiplet rates of 0.4–5.4%. Following reverse transcription and cell barcoding in droplets, emulsions were broken and cDNA purified using Dynabeads MyOne SILANE followed by PCR amplification (98°C for 45 sec; 13–18 cycles of 98°C for 20 sec, 67°C for 30 sec, 72°C for 1 min; 72°C for 1 min). Amplified cDNA was then used for both 5’ gene expression library construction and TCR enrichment. For gene expression library construction, 2.4–50 ng of amplified cDNA was fragmented and end-repaired, double-sided size selected with SPRIselect beads, PCR amplified with sample indexing primers (98°C for 45 sec; 14–16 cycles of 98°C for 20 sec, 54°C for 30 sec, 72°C for 20 sec; 72°C for 1 min), and double-sided size selected with SPRIselect beads. For TCR library construction, TCR transcripts were enriched from 2μL of amplified cDNA by PCR (primer sets 1 and 2: 98°C for 45 sec; 10 cycles of 98°C for 20 sec, 67°C for 30 sec, 72°C for 1 min; 72°C for 1 min). Following TCR enrichment, 5–50 ng of enriched PCR product was fragmented and end-repaired, size selected with SPRIselect beads, PCR amplified with sample indexing primers (98°C for 45 sec; 9 cycles of 98°C for 20 sec, 54°C for 30 sec, 72°C for 20 sec; 72°C for 1 min), and size selected with SPRIselect beads.

Sequencing

Single-cell RNA libraries were sequenced on an Illumina NextSeq or HiSeq 4000 to a minimum sequencing depth of 25,000 reads/cell using the read lengths 26bp Read1, 8bp i7 Index, 98bp Read2. Single-cell TCR libraries were sequenced on an Illumina MiSeq or HiSeq 4000 to a minimum sequencing depth of 5,000 reads/cell using the read lengths 150bp Read1, 8bp i7 Index, 150bp Read2.

Data Processing of exome libraries

Whole exome sequencing was preprocessed using a standard GATK approach⁴⁵. Briefly, both tumor and normal samples were aligned to GRCh37 using bwa-mem⁴⁶ and further processed to remove duplicates and recalibrate base quality scores. All processing was performed in FireCloud⁴⁷.

Mutation calling and neoepitope prediction

Small somatic variants were identified using Mutect2⁴⁸ and further annotated with the GATK. Somatic copy number variants were identified using the GATK best practices pipeline. HLA typing was performed on the germline whole exome sample using xHLA⁴⁹. Neoepitopes were identified using pVAC-seq⁵⁰, where a peptide-MHC pair was considered a neoepitope if the peptide was found to bind to the MHC allele with less than 500 nM binding strength and its wildtype cognate bound to the same allele with greater than 500 nM binding strength. We also examined additional neoepitope filters in Extended Data Fig. 1b including mutant peptide binding strength less than 500 or 50 nM and observed similar trends across patients and treatment conditions.

Tumor Clonal Composition Analysis

For the clonal evolution analysis, somatic single-nucleotide variants (SNVs) were called using Mutect 1.1.7⁴⁸ and the variant assurance pipeline (VAP)⁵¹ for filtering and rescuing. The VAP filters for FFPE and other artifacts and also leverages sequencing data from related samples to salvage false-negatives that would otherwise occur due to limits of the variant caller. VAFs (variant allele frequencies) were calculated for the detected and rescued variants by dividing the number of reads carrying the variant by the total number of reads spanning that position. For each case, mutations covered by less than 20 reads in any sample were removed, as were mutations where the alternate allele was not supported by at least four reads in at least one sample. TitanCNA⁵² was utilized to define local copy number and purity of the tumor samples. Observed VAFs were adjusted for local copy number and purity using the CHAT⁵³ framework in order to generate CCF (cancer cell fraction) estimates for each mutation in each sample. Case su002 was excluded from further clonal evolution analysis because it had a purity of < 15% (as inferred by TitanCNA) in both the pre-treatment and post-treatment samples, reducing the accuracy of imputed CCF values. Next, we used PyClone⁵⁴ to define mutational clusters and assess changes in cluster frequencies across treatment. PyClone’s Dirichlet process clustering was carried out on the functional mutations identified in each case. For case su006, since fewer functional mutations were identified compared to the other cases, all filtered mutations (i.e. including synonymous SNVs) that passed the depth of coverage thresholds described above were used to define mutational clusters. The PyClone beta binomial model was run using default parameters for 10,000 iterations with a burn-in of 1000. For visualization of each case, we plotted PyClone clusters comprising at least 1% of the total number of utilized mutations.

Data Processing of single-cell RNA-seq libraries

Single-cell RNA-seq reads were aligned to the GRCh38 reference genome and quantified using cellranger count (10X Genomics, version 2.1.0). Filtered gene-barcodes matrices containing only barcodes with UMI counts passing threshold for cell detection were used for further analysis. On average, we obtained reads from 1,862 genes per cell (median: 1,716) and 6,304 unique transcripts per cell (median: 4,777), comparable to prior droplet based scRNA-seq studies of human cancers.

Principal component analysis (PCA) and UMAP clustering

All additional analysis was performed using Seurat (version 2.3.4)⁵⁵. Cells with less than 200 genes detected or greater than 10% mitochondrial RNA content were excluded from analysis, with 79,046/83,583 cells passing filter (95%).

For clustering of all cell types in BCC TME, raw UMI counts were log normalized and variable genes called on each dataset independently based on average expression greater than 0.1 and average dispersion greater than 1. Variable T cell receptor and immunoglobulin genes were removed from the list of variable genes to prevent clustering based on variable V(D)J transcripts. To remove batch effects between samples associated with a heat shock gene expression signature, we assigned a heat shock score using the AddModuleScore function based on genes annotated with the GO biological process ontology term ‘cellular response to heat’. Additionally, we assigned scores for S and G2/M cell cycle phase based on previously defined gene sets⁸ using the CellCycleScoring function. Scaled z-scores for each gene were calculated using the ScaleData function and regressed against number of UMIs per cell, mitochondrial RNA content, S phase score, G2/M phase score and heat shock score. Scaled data was used an input into PCA based on variable genes. Clusters were identified using shared nearest neighbor (SNN) based clustering based on the first 20 PCs with k = 30 and resolution = 0.4. The same principal components were used to generate the UMAP projections^56,57, which were generated with a minimum distance of 1 and 20 neighbors.

For malignant cell and T cell specific clustering in BCC samples, we isolated subsets of cells from the complete data set identified as either malignant or T cells based on broad clustering. Cells were then re-clustered as described above, with the following modifications: For malignant cells, we did not observe cell-cycle associated effects and did not regress out cell cycle scores. Variable genes were called on the merged dataset based on average expression greater than 0.1 and average dispersion greater than 1.8. For UMAP visualization, we used the first 10 PCs, a minimum distance of 0.15 and 15 neighbors. For T cell clustering, we called variable genes on each dataset independently based on average expression greater than 0.15 and average dispersion greater than 3, and used the merged variable gene set for PCA. T cell clusters were identified using SNN-based clustering based on the first 16 PCs with k = 30 and resolution = 0.3. For UMAP visualization, we used the same PCs, a minimum distance of 0.05 and 20 neighbors.

T cell clustering of SCC samples was performed as described above with the following modifications: Variable genes were called on each dataset independently based on average expression greater than 0.15 and average dispersion greater than 2 and the merged variable gene set used for PCA. We observed three small outlier clusters based on initial clustering which expressed B cell and macrophage marker genes which were removed from further analysis. T cell clusters were identified using SNN-based clustering based on the first 16 PCs with k = 30 and resolution = 0.3. For UMAP visualization, we used the first 16 PCs, a minimum distance of 0.05 and 20 neighbors.

Cell Cluster Annotation

Clusters were annotated based on expression of known marker genes, including CD3G, CD3D, CD3E, CD2 (T cells), CD8A, GZMA (CD8⁺ T cells), CD4, FOXP3 (CD4⁺ T cells/Tregs), KLRC1, KLRC3 (NK cells), CD19, CD79A (B cells), SLAMF7, IGKC (Plasma cells), FCGR2A, CSF1R (Macrophages), FLT3 (Dendritic cells), CLEC4C (Plasmacytoid Dendritic cells), COL1A2 (Fibroblasts), MCAM, MYLK (Myofibroblasts), FAP, PDPN (CAFs), EPCAM, TP63 (Malignant cells), PECAM1, VWF (Endothelial cells), PMEL, MLANA (Melanocytes). Clusters were also confirmed by identifying differentially expressed marker genes for each cluster and comparing to known cell type marker genes. Finally, we downloaded bulk RNA-seq count data from sorted immune cell populations from Calderon et al., 2018⁶ and compared bulk gene expression to pseudo-bulk expression profiles from single cell clusters. UMI counts were summed for all cells in each cluster to generate pseudo-bulk profiles. Gene counts from aggregated single-cell and bulk data were then normalized and depth corrected using variance stabilizing transformation in DESeq2 (version 1.18.1). Genes with a coefficient of variation greater than 20% across bulk RNA-seq datasets were used to calculate the Pearson correlation between bulk datasets and pseudo-bulk profiles.

Data Processing of single-cell TCR-seq libraries

TCR reads were aligned to the GRCh38 reference genome and consensus TCR annotation was performed using cellranger vdj (10X Genomics, version 2.1.0). TCR libraries were sequenced to a minimum depth of 5,000 reads per cell, with a final average of 15,341 reads per cell. On average, 12,335 reads mapped to either the TRA or TRB loci in each cell. TCR annotation was performed using the 10X cellranger vdj pipeline as described. 85% of annotated T cells were assigned a TCR and only 0.18% of cells not annotated as T cells were assigned a TCR. For cells with two confident TCRs, both were considered in the analysis. Overall, 5.5% of T cells with TCR reads were assigned two productive TRB sequences and 5.7% of T cells with TCR reads were assigned two productive TRA sequences and both sequences were assigned to each cell and used for clonotype grouping. Only 1.6% of all cells were assigned two TRB sequences and two TRA sequences. We detected an average of 1,863 unique clonotypes on average in each patient (range 151 – 4,081). Of 27,956 total clonotypes detected, an average of 1.84 cells were assigned to each clonotype, 5,291 clonotypes comprised of greater than one cell, and clonotype sizes ranging from 1 cell to 564 cells.

GLIPH Analysis

To identify TCR specificity groups, GLIPH analysis was carried out as described previously²⁸. GLIPH clusters TCRs based on two similarity indexes: 1) global similarity, meaning that CDR3 sequences differ by up to one amino acid, and 2) local similarity, meaning that two TCRs contain a common CDR3 motif of 2, 3, or 4 amino acids (enriched over random sub-sampling of unselected repertoires). We performed GLIPH with the following modifications: 1) for clusters based on global similarity, CDR3b fragments within the same cluster are required to be at most one amino acid different, and this difference must be at the same amino acid location in all fragments within the cluster, and 2) for clusters based on local motifs, the starting positions of motifs of the same cluster within CDR3b fragments must be within 3 amino acids to be considered.

Single-cell CNV detection

Single-cell CNVs were detected using HoneyBADGER⁵⁸. Log-transformed UMI counts were used as input, after removing genes with mean expression lower than 0.1 normalized counts (7,189 genes passing filter, 75–753 genes per chromosome). Non-immune, non-malignant cells were used as a normal reference, including fibroblasts, endothelial cells and melanocytes (n = 2,122). CNVs were detected based on the average gene expression in sliding windows across each chromosome (n = 101 genes/window) relative to average expression in normal reference cells. CNV profiles of malignant and reference cells were visualized with z-score limits of −0.6 and 0.6.

Generation and Data Processing of bulk RNA-seq libraries

For bulk CD4⁺ T cell subset RNA-seq, cDNA library construction was performed using the SMART-Seq v4 Ultra Low Input RNA Kit (Clontech) with 2 ng of input RNA. Sequencing libraries were prepared using the Nextera XT DNA Library Prep Kit (Illumina), quantitated using the Qubit dsDNA HS Kit (Thermo Fisher Scientific), and pooled in equimolar ratios. Final pooled libraries were sequenced on an Illumina HiSeq 2500 with paired-end 50-bp read lengths. Paired end RNA-seq libraries from basal cell carcinoma tumors (Atwood et al., 2015¹⁴), cutaneous squamous cell carcinoma tumors (Hoang, et al., 2017¹⁵) and T cell subsets (Simoni, et al. 2018²¹ and this study) were aligned to the GRCh38 reference genome using STAR (version 2.6.1a) following adapter trimming by cutadapt (version 1.17). Uniquely-mapped reads were counted using featureCounts (version 1.6.2) using Ensembl GRCh38 GTF transcript annotations. Differential expression analysis was performed using DESeq2 to identify cell-type specific expression programs (version 1.18.1).

Gene expression signature scoring

Individual cells were scored for bulk RNA-seq expression programs derived from bulk RNA-seq data as follows. Raw UMI counts were used as input into the AUCell package to score each cell for gene set enrichment based on AUC scores to correct for gene dropout and library size differences⁵⁹. After building a gene expression ranking for each cell, the gene set enrichment was calculated for each cell using the area under the recovery curve using default parameters.

Activation and exhaustion signatures were derived by identifying variable genes across all CD8⁺ T cells using the FindVariableGenes function in Seurat with an average expression cutoff of 0.05 and dispersion cutoff of 0.5. The Pearson correlation between reference genes IFNG (activation signature) and HAVCR2 (exhaustion signature) and all variable genes across all CD8⁺ T cells was computed using scaled expression values. Exhaustion and activation signature genes were comprised of the top 50 genes with the highest correlation with reference genes IFNG and HAVCR2. The TCF7⁺/stem-like signature was obtained from processed data from Im et al., 2016²⁹. Individual cells were scored for enrichment of gene signatures using the function AddModuleScore in Seurat. Cell cycle scoring was performed as previously described⁸. Briefly, cells were scored for enrichment of cell cycle associated genes using the CellCycleScoring function in Seurat.

Diffusion map and pseudotime analysis

Single cells from BCC samples assigned to CD8⁺ T cell clusters were used for diffusion map and pseudotime analysis. Differentially expressed genes were used to recalculate principal components. Data was then exported to Scanpy (version 1.2.2)⁶⁰ for diffusion map and pseudotime analysis. Data was preprocessed by computing a neighborhood graph using 40 neighbors, the first 20 PCs. The first three components of the diffusion map were then computed. A randomly selected naïve T cell was used as the root cell for diffusion psuedotime computation using the first 3 diffusion components and a minimum group size of 10.

Sources for bulk RNA sequencing data

Reference bulk RNA-seq from sorted immune populations were obtained from GEO (GSE118165). Reference bulk RNA-seq data from CD8⁺ T cells were obtained from GEO (GSE113590). Reference bulk RNA-seq data from basal cell carcinomas were obtained from GEO (GSE58377). Reference bulk RNA-seq data from squamous cell carcinomas were obtained from the ArrayExpress database (E-MTAB-5678).

Statistical analysis

The statistical methods used for each analysis are described within the figure legends and in the Life Science Reporting Summary linked to this article.

Code Availability

All custom code used in this work is available upon request.

Reporting Summary

Further information on research design is available in the Life Sciences Reporting Summary linked to this article.

Data Availability

All ensemble and single-cell RNA sequencing data have been deposited in the Gene Expression Omnibus (GEO) and are available under accession GSE123814. Exome sequencing data has been deposited in the Sequence Read Archive (SRA) and are available under accession PRJNA533341. Bulk TCR sequencing data can be accessed through Adaptive Biotechnologies’ ImmuneACCESS database (doi:10.21417/KY2019NM; https://clients.adaptivebiotech.com/pub/yost-2019-natmed). All other relevant data are available from the corresponding author upon reasonable request.

Extended Data

Extended Data Fig. 3 — **(a)** Inferred CNV profiles for malignant cells separated by patient based on scRNA-seq (scCNV) and WES. Dashed line indicates a potential subclone identified by scCNV highlighted for su005. For all patients, pre and post-treatment malignant cells were analyzed together and exhibited similar CNV profiles, with the exception of su006. For su006, differences between timepoints were apparent in CNV profiles obtained from both scRNA-seq as well as exome, analogous to the changes in mutation composition identified in Extended Data Fig. 1a. **(b)** Heatmap of differentially expressed genes (rows, n = 577) across malignant BCC cells (n = 3,548) aggregated by patient (columns, n = 8). Cutoffs for differential expression were less than 0.01 adjusted P-value (Wilcoxon rank sum test, two-tailed, Bonferroni corrected), greater than 0.3 average log fold change and greater than 0.3 difference in fraction of positive cells. Core BCC genes that are differentially expressed between all malignant cells and other TME cells are shown in top cluster. Genes differentially expressed between patients are shown in the bottom clusters. Specific genes associated with cancer-associated pathways are highlighted.

Extended Data Fig. 4 — **(a)** Enrichment of bulk T cell subtype signatures for each T cell cluster identified in the BCC TME. T cell subtype signatures were derived from bulk datasets (from this study and Simoni *et al.*, 2018, n = 3–7 biologically independent samples from different donors) and single T cells from BCC dataset were scored for signature enrichment. Heatmaps represents the z-scored average signature enrichment for each cluster. **(b)** Heatmap of Pearson correlation between T cell clusters based on first 20 PCs used for clustering (n = 33,106 cells). **(c)** UMAP of all T cells colored by marker gene expression. **(d)** UMAP of all T cells separated by patient and colored by anti-PD-1 treatment status.

Extended Data Fig. 5 — **(a)** Violin plots of cell coordinates in diffusion components 1 and 2 separated by cluster identity (left, middle). Violin plot of pseudotime values separated by cluster identity (right). N = number of cells. **(b)** Heatmap of expression of genes with highest correlation with diffusion components 1 and 2 (rows) across cells belonging to each cell type cluster (columns). **(c)** Boxplot of Gini indices for each CD4⁺ T cell cluster separated by timepoint, showing clonal expansion of Tfh cells following treatment. Each point represents a patient with more than 10 cells belonging to a cluster at that timepoint, with the size proportional to the number of cells. **(d)** UMAP of all cells detected for patient su001 colored by treatment timepoint (left) and relative proportions of each immune cell type (right), showing increased frequency of B cells posttreatment. **(e)** UMAP of T cells detected for patient su001 colored by treatment timepoint (left) and relative proportions of CD4⁺ phenotype (right), showing increased frequency of Tfh cells post-treatment. **(f)** Bar plot of percent *AICDA* positive B cells, separated by patient. **(g)** H&E staining of post-treatment BCC tumor from patient su001 post-treatment demonstrating islands of BCC in sclerotic stroma with a peripheral cuff of dense lymphoid tissue. Scale bar for top image represents 400 μm and scale bar for bottom image represents 100 μm. H&E staining was performed once for each sample.

Extended Data Fig. 6 — **(a)** Distributions of the proportion of cells within each clone (>=3 cells) that share a common cluster identity, separated by patient (for patients with >3 clones with >=3 cells), compared to randomly selected and size matched groups of T cells (n = number of clones, unpaired t-test, two-tailed). **(b)** Distribution of the proportion of CD4⁺ cells (left) and CD8⁺ cells (right) within each clone or TCRβ clones within each TCR specificity group (>=3 cells) that share a common cluster identity, separated by treatment timepoint, compared to randomly selected and size matched groups of T cells from the same sample (left, n = number of clones, unpaired t-test, two-tailed). **(c)** Bar plot of T cell cluster assignments for all clones with greater than 10 cells, separated by patient and treatment status. **(d)** Bar plot of T cell cluster assignments for the largest 10 TCR specificity groups, separated by TCRβ clone. Conserved motifs between TCRβ clones identified by GLIPH highlighted in red. Representative TCRβ sequences shown for TCR specificity groups with more than four unique clonotypes. **(e)** Heatmap of the fraction of TCR specificity groups with clones belonging to a given primary phenotype (rows) that also contain clones belonging to a secondary phenotype (columns).

Extended Data Fig. 7 — **(a)** Heatmap of TCRβ clonotype overlap between all samples, indicating correct pairing of samples and a significant number of overlapping clones between timepoint within individual patients with the exception of one pair with limited cell numbers and no clonotype overlap (su003) **(b)** Bar plot of T cell cluster assignments for matched TCRβ clones between timepoints for top 60 clones with at least 3 cells per timepoint. Related to Figure 3e.

Extended Data Fig. 8 — **(a)** Scatterplots comparing TCRβ clone frequencies pre- and post-treatment measured by single-cell RNA+TCR-seq, separated by patient. Clones where the majority of cells share an exhausted CD8⁺ phenotype (red) or a memory CD8⁺ phenotype (blue) are highlighted. Patient su003 without no clonotype overlap between timepoints excluded. In this and subsequent panels, exhausted refers to both exhausted and exhausted/activated clusters. **(b)** Boxplot of the fraction of novel TCR specificity groups within each cluster following treatment for TCR specificity groups containing at least two distinct TCRβ sequences and at least 3 cells, separated by patient (n = number of patients). **(c)** Bar plot of fraction of clones with significant expansion post-treatment based on bulk TCRseq, separated by patient and phenotype and colored by replacement status. **(d)** Scatterplots comparing TCRβ clone frequencies between timepoints measured by bulk TCR-seq for sequential timepoints in patient su001, with clones where the majority of cells share an exhausted CD8⁺ phenotype (red) or a memory CD8⁺ phenotype (blue) highlighted. Novel clones emerging between timepoints are highlighted in dark red and are detected only in pre- and post-treatment comparisons, but not in comparisons between pre-treatment timepoints, suggesting that replacement is primarily a result of PD-1 blockade rather than time between sampling.

Extended Data Fig. 9 — **(a)** Pie chart of percentage of TCRβ clones detected in peripheral blood that were also detected by scRNA-seq, expanded to show distribution of phenotypes in tumor, as well as fraction of exhausted clones detected in peripheral blood, colored by replacement status in tumor. In this and subsequent panels, the exhausted category includes both exhausted and exhausted/activated clusters. **(b)** Bar plot of percentage peripheral T cells matching tumor-infiltrating TCRβ clones with exhausted phenotypes post-treatment as detected by scRNA-seq. **(c)** Violin plot of TCR specificity group enrichment (tumor frequency / PBMC frequency) detected by bulk TCRseq, separated by phenotype and treatment status (n = number of TCR specificity groups, unpaired t-test, one-tailed).

Extended Data Fig. 10 — **(a)** UMAP of tumor-infiltrating T cells present in SCC samples pre- and post-PD-1 blockade colored by patient (top right) and anti-PD-1 treatment status (bottom right). **(b)** Heatmap of correlation between averaged RNA expression between BCC and SCC T cell clusters. **(c)** Boxplot of Gini indices for each CD8⁺ T cell cluster calculated for each patient (n = number of patients). In this and subsequent panels, exhausted refers to both exhausted and exhausted/activated clusters, unless otherwise noted. **(d)** Abundance of the top 12 exhausted clones in sample su010-S identified by unsupervised clustering compared to the abundance of the same clones in sorted CD8⁺ CD39⁺ T cells, colored by assigned phenotype. **(e)** Distribution of the proportion of cells within each clone or TCRβ clones within each TCR specificity group (>=3 cells) that share a common cluster identity, separated by treatment timepoint, compared to randomly selected and size matched groups of T cells from the same sample (left, n = number of TCRβ clones or TCR specificity groups, unpaired t-test, two-tailed). **(f)** Heatmap of the fraction of clonotypes belonging to a given primary phenotype cluster (rows) that are shared with other secondary phenotype clusters (columns). **(g)** Heatmap of all observed phenotype transitions for matched clones during PD-1 blockade for clones with at least 3 cells for each timepoint. **(h)** *TCF7*⁺/stem-like score (from Im *et al*. 2016) versus exhaustion score for all CD8⁺ T cells, colored by gene expression (left). *TCF7*⁺/stem-like score versus exhaustion score for exhausted cells and cells of other phenotypes belonging to primarily exhausted clones, colored by phenotype (top right). Violin plot of *TCF7*⁺/stem-like score for exhausted cells and cells of other phenotypes belonging to primarily exhausted clones, demonstrating that the highest *TCF7*⁺/stem-like score is observed in cells with an exhausted phenotype (bottom right, n = number of cells). **(i)** Violin plot of *TCF7*⁺/stemlike score for memory and exhausted cells separated by change in clone abundance following treatment (left, n = number of cells, unpaired t-test, two-tailed). Clones were defined as expanded or contracted if they significantly changed in abundance by a Fisher exact test (P < 0.05 and fold change > 0.5), and persistent if they did not significantly change in abundance and at least one cell was detected at each timepoint. **(j)** Scatterplots comparing TCRβ clone frequencies pre- and post-treatment measured by single-cell TCR sequencing for all BCC patients. Clones that were significantly expanded or contracted post-treatment based on a binomial test (two-sided, Bonferroni corrected P value < 0.01) are highlighted on the left. Clones where the majority of cells share an exhausted CD8⁺ phenotype (middle, red) or a memory CD8⁺ phenotype (right, blue) are also highlighted.</Figure_Caption>

Supplementary Material

NIHMS1531727-supplement-1.pdf^{(76KB, pdf)}

NIHMS1531727-supplement-2.xlsx^{(2.8MB, xlsx)}

Acknowledgments

We thank members of the Chang laboratory for helpful discussions. We thank X. Ji, D. Wagh, and J. Coller at the Stanford Functional Genomics Facility, J. Sung and S. Fitch at the Stanford Clinical and Translational Research Unit Biobank, and D.D. Hiraki at the Stanford Blood Center Histocompatibility, Immunogenetics, and Disease Profiling Laboratory. We thank A.D. Colevas, S. Reddy, L. Van Der Bokke, R. Patel for clinical collaboration. We thank A. Pague and A. Valencia for assistance with clinical specimens. We thank Khriszha Quema-Yee at the Parker Institute for Cancer Immunotherapy for assistance with illustrations. This work was supported by the Parker Institute for Cancer Immunotherapy (A.T.S., D.K.W., R.K., S.L.B., M.M.D., and H.Y.C.), the Michelson Foundation (A.T.S.), and the National Institutes of Health (NIH) P50HG007735, R35CA209919 (H.Y.C.), K08CA23188-01 (A.T.S.), K23CA211793 (K.Y.S.) and 5U19AI057229 (M.M.D.). K.E.Y. was supported by the National Science Foundation Graduate Research Fellowship Program (NSF DGE-1656518) and a Stanford Graduate Fellowship. A.T.S. was supported by a Bridge Scholar Award from the Parker Institute for Cancer Immunotherapy, a Career Award for Medical Scientists from the Burroughs Wellcome Fund, and the Human Vaccines Project Michelson Prize for Human Immunology and Vaccine Research. Cell sorting for this project was done on instruments in the Stanford Shared FACS Facility. Sequencing was performed by the Stanford Functional Genomics Facility (supported by NIH grant S10OD018220). M.M.D. and H.Y.C. are investigators of the Howard Hughes Medical Institute.

Footnotes

Competing Interests

H.Y.C. is a co-founder of Accent Therapeutics and an advisor for 10x Genomics and Spring Discovery. A. L. S. C. was an advisory board member and clinical investigator for studies sponsored by Merck, Regeneron, Novartis, Galderma, Genentech Roche. A.T.S. and D.K.W. are advisors for Immunai.

References

1.Sharma P & Allison JP The future of immune checkpoint therapy. Science 348, 56–61 (2015). [DOI] [PubMed] [Google Scholar]
2.Sakuishi K et al. Targeting Tim-3 and PD-1 pathways to reverse T cell exhaustion and restore anti-tumor immunity. J. Exp. Med 207, 2187–2194 (2010). [DOI] [PMC free article] [PubMed] [Google Scholar]
3.Wherry EJ & Kurachi M Molecular and cellular insights into T cell exhaustion. Nat. Rev. Immunol 15, 486–499 (2015). [DOI] [PMC free article] [PubMed] [Google Scholar]
4.Pauken KE et al. Epigenetic stability of exhausted T cells limits durability of reinvigoration by PD-1 blockade. Science 354, 1160–1165 (2016). [DOI] [PMC free article] [PubMed] [Google Scholar]
5.Riaz N et al. Tumor and Microenvironment Evolution during Immunotherapy with Nivolumab. Cell 171, 934–949.e16 (2017). [DOI] [PMC free article] [PubMed] [Google Scholar]
6.Calderon D et al. Landscape of stimulation-responsive chromatin across diverse human immune cells. bioRxiv 409722 (2018). doi: 10.1101/409722 [DOI] [PMC free article] [PubMed] [Google Scholar]
7.Bonilla X et al. Genomic analysis identifies new drivers and progression pathways in skin basal cell carcinoma. Nat. Genet. 47, 398 (2015). [DOI] [PubMed] [Google Scholar]
8.Tirosh I et al. Dissecting the multicellular ecosystem of metastatic melanoma by single-cell RNA-seq. Science 352, 189–196 (2016). [DOI] [PMC free article] [PubMed] [Google Scholar]
9.Puram SV et al. Single-Cell Transcriptomic Analysis of Primary and Metastatic Tumor Ecosystems in Head and Neck Cancer. Cell 171, 1611–1624.e24 (2017). [DOI] [PMC free article] [PubMed] [Google Scholar]
10.Tellechea O, Reis JP, Domingues JC & Baptista AP Monoclonal antibody Ber EP4 distinguishes basal-cell carcinoma from squamous-cell carcinoma of the skin. Am. J. Dermatopathol 15, 452–455 (1993). [DOI] [PubMed] [Google Scholar]
11.Bircan S, Candir O, Kapucoglu N & Baspinar S The expression of p63 in basal cell carcinomas and association with histological differentiation. J. Cutan. Pathol 33, 293–298 (2006). [DOI] [PubMed] [Google Scholar]
12.Bernemann TM, Podda M, Wolter M & Boehncke WH Expression of the basal cell adhesion molecule (B-CAM) in normal and diseased human skin. J. Cutan. Pathol 27, 108–111 (2000). [DOI] [PubMed] [Google Scholar]
13.Ransohoff KJ, Tang JY & Sarin KY Squamous Change in Basal-Cell Carcinoma with Drug Resistance. (2015). doi: 10.1056/NEJMc1504261 [DOI] [PubMed] [Google Scholar]
14.Atwood SX et al. Smoothened variants explain the majority of drug resistance in basal cell carcinoma. Cancer Cell 27, 342–353 (2015). [DOI] [PMC free article] [PubMed] [Google Scholar]
15.Hoang VLT et al. RNA-seq reveals more consistent reference genes for gene expression studies in human non-melanoma skin cancers. PeerJ 5, e3631 (2017). [DOI] [PMC free article] [PubMed] [Google Scholar]
16.Wei SC et al. Distinct Cellular Mechanisms Underlie Anti-CTLA-4 and Anti-PD-1 Checkpoint Blockade. Cell 170, 1120–1133.e17 (2017). [DOI] [PMC free article] [PubMed] [Google Scholar]
17.Sade-Feldman M et al. Defining T Cell States Associated with Response to Checkpoint Immunotherapy in Melanoma. Cell 175, 998–1013.e20 (2018). [DOI] [PMC free article] [PubMed] [Google Scholar]
18.Haghverdi L, Buettner F & Theis FJ Diffusion maps for high-dimensional single-cell analysis of differentiation data. Bioinforma. Oxf. Engl 31, 2989–2998 (2015). [DOI] [PubMed] [Google Scholar]
19.Kumar BV et al. Human Tissue-Resident Memory T Cells Are Defined by Core Transcriptional and Functional Signatures in Lymphoid and Mucosal Sites. Cell Rep. 20, 2921–2934 (2017). [DOI] [PMC free article] [PubMed] [Google Scholar]
20.Savas P et al. Single-cell profiling of breast cancer T cells reveals a tissue-resident memory subset associated with improved prognosis. Nat. Med 24, 986 (2018). [DOI] [PubMed] [Google Scholar]
21.Simoni Y et al. Bystander CD8 + T cells are abundant and phenotypically distinct in human tumour infiltrates. Nature 557, 575 (2018). [DOI] [PubMed] [Google Scholar]
22.Duhen T et al. Co-expression of CD39 and CD103 identifies tumor-reactive CD8 T cells in human solid tumors. Nat. Commun 9, 2724 (2018). [DOI] [PMC free article] [PubMed] [Google Scholar]
23.Li H et al. Dysfunctional CD8 T Cells Form a Proliferative, Dynamically Regulated Compartment within Human Melanoma. Cell 176, 775–789.e18 (2019). [DOI] [PMC free article] [PubMed] [Google Scholar]
24.Thommen DS et al. A transcriptionally and functionally distinct PD-1 + CD8 + T cell pool with predictive potential in non-small-cell lung cancer treated with PD-1 blockade. Nat. Med 24, 994 (2018). [DOI] [PMC free article] [PubMed] [Google Scholar]
25.Azizi E et al. Single-Cell Map of Diverse Immune Phenotypes in the Breast Tumor Microenvironment. Cell 174, 1293–1308.e36 (2018). [DOI] [PMC free article] [PubMed] [Google Scholar]
26.Zhang L et al. Lineage tracking reveals dynamic relationships of T cells in colorectal cancer. Nature 1 (2018). doi: 10.1038/s41586-018-0694-x [DOI] [PubMed] [Google Scholar]
27.Zemmour D et al. Single-cell gene expression reveals a landscape of regulatory T cell phenotypes shaped by the TCR. Nat. Immunol 19, 291–301 (2018). [DOI] [PMC free article] [PubMed] [Google Scholar]
28.Glanville J et al. Identifying specificity groups in the T cell receptor repertoire. Nature 547, 94–98 (2017). [DOI] [PMC free article] [PubMed] [Google Scholar]
29.Im SJ et al. Defining CD8+ T cells that provide the proliferative burst after PD-1 therapy. Nature 537, 417–421 (2016). [DOI] [PMC free article] [PubMed] [Google Scholar]
30.Kurtulus S et al. Checkpoint Blockade Immunotherapy Induces Dynamic Changes in PD-1−CD8+ Tumor-Infiltrating T Cells. Immunity 50, 181–194.e6 (2019). [DOI] [PMC free article] [PubMed] [Google Scholar]
31.Siddiqui I et al. Intratumoral Tcf1+PD-1+CD8+ T Cells with Stem-like Properties Promote Tumor Control in Response to Vaccination and Checkpoint Blockade Immunotherapy. Immunity 50, 195–211.e10 (2019). [DOI] [PubMed] [Google Scholar]
32.Miller BC et al. Subsets of exhausted CD8 + T cells differentially mediate tumor control and respond to checkpoint blockade. Nat. Immunol 20, 326 (2019). [DOI] [PMC free article] [PubMed] [Google Scholar]
33.Huang AC et al. A single dose of neoadjuvant PD-1 blockade predicts clinical outcomes in resectable melanoma. Nat. Med 25, 454 (2019). [DOI] [PMC free article] [PubMed] [Google Scholar]
34.Kamphorst AO et al. Proliferation of PD-1+ CD8 T cells in peripheral blood after PD-1-targeted therapy in lung cancer patients. Proc. Natl. Acad. Sci. U. S. A 114, 4993–4998 (2017). [DOI] [PMC free article] [PubMed] [Google Scholar]
35.Ghoneim HE et al. De Novo Epigenetic Programs Inhibit PD-1 Blockade-Mediated T Cell Rejuvenation. Cell 170, 142–157.e19 (2017). [DOI] [PMC free article] [PubMed] [Google Scholar]
36.Spitzer MH et al. Systemic Immunity Is Required for Effective Cancer Immunotherapy. Cell 168, 487–502.e15 (2017). [DOI] [PMC free article] [PubMed] [Google Scholar]
37.Matsushita H et al. Cancer exome analysis reveals a T-cell-dependent mechanism of cancer immunoediting. Nature 482, 400–404 (2012). [DOI] [PMC free article] [PubMed] [Google Scholar]
38.Scheper W et al. Low and variable tumor reactivity of the intratumoral TCR repertoire in human cancers. Nat. Med 25, 89–94 (2019). [DOI] [PubMed] [Google Scholar]
39.Gee MH et al. Antigen Identification for Orphan T Cell Receptors Expressed on Tumor-Infiltrating Lymphocytes. Cell (2017). doi: 10.1016/j.cell.2017.11.043 [DOI] [PMC free article] [PubMed] [Google Scholar]
40.Li J et al. Tumor Cell-Intrinsic Factors Underlie Heterogeneity of Immune Cell Infiltration and Response to Immunotherapy. Immunity 49, 178–193.e7 (2018). [DOI] [PMC free article] [PubMed] [Google Scholar]
41.Chang ALS et al. Pembrolizumab for advanced basal cell carcinoma: an investigator-initiated, proof-of-concept study. J. Am. Acad. Dermatol (2018). doi: 10.1016/j.jaad.2018.08.017 [DOI] [PMC free article] [PubMed] [Google Scholar]
42.Eisenhauer EA et al. New response evaluation criteria in solid tumours: revised RECIST guideline (version 1.1). Eur. J. Cancer Oxf. Engl 1990 45, 228–247 (2009). [DOI] [PubMed] [Google Scholar]
43.Wang C et al. High-throughput, high-fidelity HLA genotyping with deep sequencing. Proc. Natl. Acad. Sci. U. S. A 109, 8676–8681 (2012). [DOI] [PMC free article] [PubMed] [Google Scholar]
44.Thorstenson YR et al. Allelic resolution NGS HLA typing of Class I and Class II loci and haplotypes in Cape Town, South Africa. Hum. Immunol (2018). doi: 10.1016/j.humimm.2018.09.004 [DOI] [PMC free article] [PubMed] [Google Scholar]
45.McKenna A et al. The Genome Analysis Toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data. Genome Res. 20, 1297–1303 (2010). [DOI] [PMC free article] [PubMed] [Google Scholar]
46.Li H & Durbin R Fast and accurate short read alignment with Burrows–Wheeler transform. Bioinformatics 25, 1754–1760 (2009). [DOI] [PMC free article] [PubMed] [Google Scholar]
47.Birger C et al. FireCloud, a scalable cloud-based platform for collaborative genome analysis: Strategies for reducing and controlling costs. bioRxiv 209494 (2017). doi: 10.1101/209494 [DOI] [Google Scholar]
48.Cibulskis K et al. Sensitive detection of somatic point mutations in impure and heterogeneous cancer samples. Nat. Biotechnol 31, 213–219 (2013). [DOI] [PMC free article] [PubMed] [Google Scholar]
49.Xie C et al. Fast and accurate HLA typing from short-read next-generation sequence data with xHLA. Proc. Natl. Acad. Sci 114, 8059–8064 (2017). [DOI] [PMC free article] [PubMed] [Google Scholar]
50.Hundal J et al. pVAC-Seq: A genome-guided in silico approach to identifying tumor neoantigens. Genome Med. 8, 11 (2016). [DOI] [PMC free article] [PubMed] [Google Scholar]
51.Between-region genetic divergence reflects the mode and tempo of tumor evolution. - PubMed - NCBI. Available at: https://www.ncbi.nlm.nih.gov/pubmed/28581503. (Accessed: 20th March 2019) [DOI] [PMC free article] [PubMed]
52.Ha G et al. TITAN: inference of copy number architectures in clonal cell populations from tumor whole-genome sequence data. Genome Res. 24, 1881–1893 (2014). [DOI] [PMC free article] [PubMed] [Google Scholar]
53.Li B & Li JZ A general framework for analyzing tumor subclonality using SNP array and DNA sequencing data. Genome Biol. 15, 473 (2014). [DOI] [PMC free article] [PubMed] [Google Scholar]
54.Roth A et al. PyClone: statistical inference of clonal population structure in cancer. Nat. Methods 11, 396–398 (2014). [DOI] [PMC free article] [PubMed] [Google Scholar]
55.Butler A, Hoffman P, Smibert P, Papalexi E & Satija R Integrating single-cell transcriptomic data across different conditions, technologies, and species. Nat. Biotechnol 36, 411–420 (2018). [DOI] [PMC free article] [PubMed] [Google Scholar]
56.McInnes L, Healy J & Melville J UMAP: Uniform Manifold Approximation and Projection for Dimension Reduction. ArXiv180203426 Cs Stat (2018). [Google Scholar]
57.Becht E et al. Dimensionality reduction for visualizing single-cell data using UMAP. Nat. Biotechnol. 37, 38–44 (2019). [DOI] [PubMed] [Google Scholar]
58.Fan J et al. Linking transcriptional and genetic tumor heterogeneity through allele analysis of single-cell RNA-seq data. Genome Res gr.228080.117 (2018). doi: 10.1101/gr.228080.117 [DOI] [PMC free article] [PubMed] [Google Scholar]
59.Aibar S et al. SCENIC: single-cell regulatory network inference and clustering. Nat. Methods 14, 1083–1086 (2017). [DOI] [PMC free article] [PubMed] [Google Scholar]
60.Wolf FA, Angerer P & Theis FJ SCANPY: large-scale single-cell gene expression data analysis. Genome Biol. 19, 15 (2018). [DOI] [PMC free article] [PubMed] [Google Scholar]

Associated Data

This section collects any data citations, data availability statements, or supplementary materials included in this article.

Supplementary Materials

NIHMS1531727-supplement-1.pdf^{(76KB, pdf)}

NIHMS1531727-supplement-2.xlsx^{(2.8MB, xlsx)}

PERMALINK

Clonal replacement of tumor-specific T cells following PD-1 blockade

Kathryn E Yost

Ansuman T Satpathy

Daniel K Wells

Yanyan Qi

Chunlin Wang

Robin Kageyama

Katherine McNamara

Jeffrey M Granja

Kavita Y Sarin

Ryanne A Brown

Rohit K Gupta

Christina Curtis

Samantha L Bucktrout

Mark M Davis

Anne Lynn S Chang

Howard Y Chang

Abstract

Main

Figure 1. Characterization of the BCC TME pre- and post-PD-1 blockade by single-cell RNA-seq.

Figure 2. Exhausted CD8+ T cells are clonally expanded and express markers of tumor-specificity.

Figure 3. Clonal dynamics and phenotype transitions of tumor-infiltrating T cells.

Figure 4. Clonal replacement of exhausted CD8+ T cells following PD-1 blockade.

Methods

Human subjects

Sample Collection and Processing

H&E and Immunohistochemical (IHC) staining

Exome Sequencing

HLA Typing

Bulk TCR Sequencing

Tumor Dissociation

Cell Sorting

Single-cell RNA-seq Library Preparation

Sequencing

Data Processing of exome libraries

Mutation calling and neoepitope prediction

Tumor Clonal Composition Analysis

Data Processing of single-cell RNA-seq libraries

Principal component analysis (PCA) and UMAP clustering

Cell Cluster Annotation

Data Processing of single-cell TCR-seq libraries

GLIPH Analysis

Single-cell CNV detection

Generation and Data Processing of bulk RNA-seq libraries

Gene expression signature scoring

Diffusion map and pseudotime analysis

Sources for bulk RNA sequencing data

Statistical analysis

Code Availability

Reporting Summary

Data Availability

Extended Data

Extended Data Fig. 1. Mutational landscape of BCC tumors following PD-1 blockade, Related to Figure 1.

Extended Data Fig. 2. Characterization of cell types present in BCC TME, Related to Figure 1.

Extended Data Fig. 3. Copy number alterations and gene expression of individual BCC tumors, Related to Figure 1.

Extended Data Fig. 4. Characterization of T cell subtypes present in BCC TME, Related to Figure 2.

Extended Data Fig. 5. Characterization of activation/exhaustion trajectories and increase in Tfh cell clonality accompanied by B cell expansion, Related to Figure 2.

Extended Data Fig. 6. Correlations between T cell clones or TCR specificity groups and scRNA-seq phenotype, Related to Figure 3.

Extended Data Fig. 7. Details of clone transitions, Related to Figure 3.

Extended Data Fig. 8. Clonal expansion in tumor and peripheral blood detected by bulk TCR sequencing, Related to Figure 4.

Extended Data Fig. 9. TCR overlap between peripheral blood and tumor detected by bulk TCR sequencing, Related to Figure 4.

Extended Data Fig. 10. Clonal replacement analysis in SCC TILs following PD-1 blockade, Related to Figure 4.

Supplementary Material

Acknowledgments

Footnotes

References

Associated Data

Supplementary Materials

Data Availability Statement

ACTIONS

PERMALINK

RESOURCES

Similar articles

Cited by other articles

Links to NCBI Databases

Figure 2. Exhausted CD8⁺ T cells are clonally expanded and express markers of tumor-specificity.

Figure 4. Clonal replacement of exhausted CD8⁺ T cells following PD-1 blockade.