Transcriptional programs of neoantigen-specific TIL in anti-PD-1-treated lung cancers

1The efficacy of PD-1- and PD-L1-blocking agents is predicated upon CD8 T cell-mediated anti-tumour immunity. Early studies focused on tumour-associated antigens, whereas recent work has shifted attention to T cell recognition of mutation-associated neoantigens (MANA), owing to the large numbers of somatic mutations acquired by many cancers during their development. The association of improved anti-PD-1 and anti-PD-L1 clinical responses with high tumour mutational burden strongly suggests that MANA are important targets of anti-tumour immunity induced by PD-1 blockade.

2Despite the success of immune checkpoint blockade (ICB) in improving clinical outcomes, most cancers still do not respond. Improving response rates to ICB will require an understanding of the functional state of tumour-specific T cells, particularly in the tumour microenvironment. However, a fundamental limitation in the current understanding of the T cell functional programs that underpin the response to ICB has been the absence of transcriptional profiling of true MANA-specific TIL. A related problem is the paucity of information regarding the differences between MANA-specific TIL in ICB-responsive versus resistant tumours. Indeed, MANA-specific T cells represent a small fraction of total TIL,, particularly in lung cancer, in which they have been shown to selectively upregulate CD39. This highlights the challenges confronting characterization of the cells responsible for the activity of T cell-targeting immunotherapies.

3For this study, we used peripheral blood and tissue biospecimens obtained from the first-in-human clinical trial of neoadjuvant anti-PD-1 (nivolumab) in resectable non-small cell lung cancer (NSCLC; ClinicalTrials.gov identifier: NCT02259621; Fig. 1a, top) to study the transcriptional programs of MANA-specific TIL. Nine out of 20 patients with NSCLC (45%) treated in this trial had a major pathologic response (MPR) at the time of resection, defined as no more than 10% viable tumour at the time of surgery; previous studies have established an association between MPR and improved overall survival–. A schematic of the study design and experimental approach is shown in Fig. 1a, bottom. Combined single-cell RNA sequencing (scRNA-seq) and T cell receptor sequencing (TCR-seq) was performed on TIL (n = 15), paired adjacent normal lung (n = 12), tumour-draining lymph nodes (TDLN, n = 3) and a distant metastasis (Extended Data Fig. 1a, Supplementary Tables 1–3). In total, 560,916 T cells passed quality control (Fig. 1b, Supplementary Table 3) and were carried forward for analyses.

4Uniform manifold approximation and projection (UMAP) analysis of cells from all samples on the basis of filtered and normalized transcript counts defined 15 T cell clusters (Fig. 1b, c, Extended Data Fig. 1b–e, Supplementary Data 1.1). Expression of subset-defining markers and T cell checkpoints was visualized in red scale on the UMAP (Fig. 1d). The two clusters designated as TRM had the highest expression of the canonical TRM genes, ZNF683 (also known as HOBIT) and ITGAE (also known as CD103), and the highest expression of a TRM gene set (Extended Data Fig. 1f, Supplementary Data 2.1). Principal component analysis (PCA) of samples based on concatenated cell-cluster-level pseudobulk profiles distinguished adjacent normal-lung T cells from TIL (Fig. 1e), but did not distinguish MPR from non-MPR TIL (Fig. 1f). We did not observe notable differentially expressed gene programs between MPR and non-MPR TIL (Supplementary Data 3), indicating that gene expression profiling of total TIL has limited sensitivity in distinguishing the pathologic response to PD-1 blockade.

5We next performed the MANA functional expansion of specific T cells assay (MANAFEST) on 9 of the 16 individuals on whom scRNA-seq–TCR-seq was conducted. This assay detects in vivo antigen-experienced T cell responses and identifies the clonal identity of the T cell receptor (TCR) corresponding to these cells. Of these nine, four were classed as MPR and five were non-MPR (results from one individual have been previously described). Putative MANA (Supplementary Tables 4–6), peptide pools representing influenza matrix and nucleoproteins, and a pool of major histocompatibility complex (MHC) class I-restricted cytomegalovirus (CMV), Epstein–Barr virus (EBV) and influenza virus epitopes (CEF) were queried for CD8+ T cell reactivity in parallel (Supplementary Tables 6, 7). From 7 (3 with MPR and 4 without MPR) of the 9 individuals, 72 total unique MANA-specific TCRs, 33 unique CEF-specific TCRs, and 52 unique influenza-specific TCRs were identified (Extended Data Fig. 2, Supplementary Tables 8, 9, Supplementary Data 4, 5). Out of 33 CEF-specific TCRs, 6 matched known public EBV-specific TCRs and 3 matched known public influenza-specific TCRs. No CMV-reactive TCRs were mapped from our CEF-specific TCRs. Notably, 4 of the 41 MANA-specific TCRVβ complementarity-determining region 3 (CDR3) clonotypes identified in a patient without MPR (patient ID MD01-004) (Extended Data Fig. 2) were specific for a MANA (MD01-004-MANA12) derived from a p53 R248L hotspot mutation, and were found at appreciable frequency in the pre- and post-treatment tumour (Extended Data Fig. 3), despite the tumour not attaining MPR. Most MANA-specific clones were detected at very low frequency (median: 0.001%) in the peripheral blood across all available time points (Fig. 2a, Extended Data Fig. 3). Overall, pathologic response was not associated with the prevalence, frequency or intratumoral representation of MANA-specific T cells (Extended Data Figs. 2, 3 Supplementary Table 9). In fact, more MANA-specific TIL were observed among non-MPR TIL than among MPR TIL. No consistent pattern was observed for the frequency of viral-specific T cells in the tissue or peripheral blood (Extended Data Figs. 2, 3).

6Ten MANA-specific clonotypes, for which the TCRα could be confidently identified from the single-cell analysis, were selected for validation of MANA recognition via TCR cloning and introduction into a Jurkat–NFAT luciferase reporter system. Seventy per cent of tested clonotypes (representing 95.2% of total cells bearing TCRs identified by MANAFEST) were validated as MANA-specific (Extended Data Fig. 4a–c). Peptide–human leukocyte antigen (HLA) binding assays demonstrated that two MANA peptides—MD01-005-MANA7 and MD01-004-MANA12—displayed comparably high MHC class I affinity (measured dissociation constants (Kd) = 5.1 nM and 17.5 nM, respectively) and stability (Extended Data Fig. 4d, e).

7We next evaluated the transcriptional programming of MANA- and viral-specific CD8+ T cells. Refined clustering of all CD8+ T cells (n = 235,851) identified 15 unique clusters (Fig. 2b, Extended Data Fig. 5a, Supplementary Data 1.2). Clusters were named on the basis of previously defined T cell states from single-cell transcriptomic studies. Six clusters had gene expression programs consistent with TRM T cells, characterized by high expression of HOBIT, LINC02446, CD103 and a previously published TRM gene set (Extended Data Fig. 5b). Selective genes and linkage to the global CD3 T cell clusters shown in Fig. 1 were visualized (Extended Data Figs. 5c, d). The six TRM subsets were heterogenous in their expression of an exhaustion gene set described previously in NSCLC (Extended Data Fig. 5e, Supplementary Data 2.2). None of the most frequent tumour-infiltrating clonotypes were restricted to a single cluster (Extended Data Fig. 5f). Among all tested individuals, a total of 28 MANA-specific CD8 clonotypes (1,350 total cells from 3 patients with MPR and 3 patients without MPR) as identified by MANAFEST were detected in the single-cell data, of which 20 clonotypes (890 cells) were in the tumour (Fig. 2c, Supplementary Table 8). Of the viral-specific T cell clonotypes, 23 influenza-specific (866 cells) and 2 EBV-specific (281 cells) clones were found in the CD8 single-cell analysis.

8Overlay of these clonotypes onto the CD8+ T cell UMAP demonstrated a marked distinction between the clonotypes with different antigen specificities (Fig. 2c, Extended Data Fig. 6a–c). EBV-reactive T cells primarily resided in effector T (Teff) cell clusters, whereas influenza- and MANA-specific T cells largely occupied distinct TRM clusters. Notably, because influenza is a respiratory virus, influenza-specific T cells may be considered the archetypal lung-resident memory T cells. None of the patients in our study were symptomatic for influenza in the six weeks preceding surgery. It is therefore not surprising that influenza-specific CD8 cells were TRM rather than Teff cells. By contrast, EBV-specific T cells exclusively occupied Teff clusters, consistent with periodic acute stimulation upon latent EBV reactivation. Whereas influenza-specific cells were the most abundant in normal lung, MANA-specific CD8 cells were more numerous in the tumour (Extended Data Fig. 6d, e).

9There was considerable shared expression of selective cytotoxic T lymphocyte (CTL) activation genes between MANA- and EBV-specific T cells, in particular genes encoding T cell activation and CTL activity, such as HLA-DRA, GZMH, IFNG and NKG7 (Fig. 2d, Supplementary Data 1.3). However, genes encoding certain canonical cytolytic molecules, such as GZMK, were expressed at low levels in MANA-specific TIL. Most notably, EOMES, which encodes a transcription factor that is critical for CTL activity, was expressed in EBV-specific CD8 cells but was minimally expressed in most MANA-specific cells. Multiple checkpoints were significantly upregulated in MANA-specific TIL compared with EBV-specific TIL. Notably, MANA-specific cells expressed higher levels of PRDM1, which encodes BLIMP-1 and has been reported to participate in coordinated transcriptional activation of multiple checkpoint genes, including PD-1 (also known as PDCD1), LAG3, TIGIT and HAVCR2. TOX, which encodes a chromatin modifier important for exhaustion programs of chronic virus-specific and tumour-specific T cells in mouse models,, was only marginally increased in MANA-specific cells, whereas its homologue, TOX2, which has also been reported to drive T cell exhaustion, showed much higher upregulation in MANA-specific versus EBV-specific CD8 TIL. HOBIT, which is selectively upregulated in TRM T cells, was also upregulated in MANA-specific TIL, even relative to influenza-specific TRM (Fig. 2e). Indeed, MANA-specific T cells demonstrated the highest immune checkpoint and exhaustion signatures (Extended Data Fig. 6f). These findings demonstrate that MANA-specific CD8 T cells in the tumour have an unconventional hybrid transcriptional program characterized by incomplete activation of effector programs and significant upregulation of checkpoint molecules such as PD-1, CTLA-4, TIM3, TIGIT and CD39. Genes encoding each of these checkpoint molecules were more highly expressed among MANA-specific CD8 cells than either influenza- or EBV-specific CD8 cells, with CD39 being the most highly differentially expressed (Fig. 2d, e), congruent with previous flow cytometry findings on MANA-specific lung cancer TIL.

10Influenza-specific TRM were distinguished from MANA-specific TRM by low levels of both activation and effector CTL programs and had lower expression of multiple checkpoint molecules, but had the highest levels of genes associated with T memory stem cells, such as TCF7 and IL-7R (Fig. 2e, f). Indeed, IL-7R expression was 4.6-fold higher on influenza-specific TIL relative to MANA-specific TIL. In TIL obtained from patient without MPR MD01-004, culture with titrating concentrations of IL-7 in vitro induced much higher levels of IL-7R-regulated genes (Supplementary Data 2.3) in influenza-specific TIL than in MANA-specific TIL (Fig. 2g, Extended Data Fig. 7). Nonetheless, supraphysiological levels of IL-7 induced appreciable upregulation of IL-7R-induced genes in MANA-specific TIL. Given the distinct transcriptional programs of the identified MANA-specific CD8 cells, we hypothesized that other CD8 T cells in the same TRM cluster showing differential expression relative to influenza-specific T cells (Fig. 2g) may also recognize MANA that were not detected by the MANAFEST assay. We cloned seven TCRs corresponding to CD8+ T cells with highly differential gene expression relative to influenza-specific T cells. We screened each TCR with a library of candidate MANA (Supplementary Table 6) and confirmed MANA recognition in three of these TCRs, one TCR each from patients MD01-004, MD01-005 and MD043-011 (Extended Data Fig. 8a–d).

11To next investigate the ligand-dependent TCR signalling capacity of antigen-specific T cells, we performed a dose–response curve with cognate peptides matched to the ten total Jurkat-validated MANA-specific TCRα–TCRβ pairs (Supplementary Table 10). Peptide dose–response curves of MPR-derived TCRs were comparable to those of EBV- and influenza-specific TCRs, suggesting that these TCRs were capable of strong ligand-dependent signalling (sometimes referred to as functional avidity). However, the peptide dose–response curves of TCRs derived from patients without MPR were markedly lower (approximately 2 log10 leftward shift in peptide dose–response curve) (Fig. 2h, Extended Data Fig. 8e). Together, our data show that despite similar measured MANA–HLA binding affinities (Extended Data Fig. 4c, d), TCR from expandable MANA-specific clones from the patient with MPR had significantly higher functional avidity than MANA-specific clones from patients without MPR.

12To explore determinants of ICB sensitivity, we examined differences in gene expression patterns between MPR and non-MPR MANA-specific TIL. The neoadjuvant clinical trial format enabled us to make this distinction through pathological analysis of surgically resected tissue. In total, we compared 45 MPR TIL transcriptomes (39 from MD01-005, 2 from MD043-003 and 4 from NY016-025) with 885 non-MPR TIL transcriptomes (782 from MD043-011, 62 from MD01-004 and 22 from NY016-014; Extended Data Fig. 9, Supplementary Table 8). We observed highly significant differences between pathologic MPR and non-MPR tumours (Fig. 3a, Supplementary Data 1.4). Significantly higher levels of genes associated with T cell dysfunction such as TOX2, CTLA4, HAVCR2 and ENTPD1 were observed for non-MPR MANA-specific T cells, whereas MPR MANA-specific T cells had higher expression of genes associated with memory (IL7R and TCF7) and effector function (GZMK) (Fig. 3a–c). Both the checkpoint score and exhaustion score were higher in MANA-specific TIL from patients without MPR (Fig. 3d, Extended Data Fig. 10a, b). Of note, CXCL13 is one of the genes most highly correlated with checkpoint-associated genes in non-MPR MANA-specific TIL, and was also found to be highly expressed in MANA-specific cells relative to virus-specific cells among CD8 TIL (Fig. 2d–f).

13A number of genes encoding T cell inhibitory molecules were more highly correlated with a composite immune checkpoint score of MANA-specific TIL from patients without MPR than those from patients with MPR (Fig. 3e, Extended Data Fig. 10c). In two patients without MPR (MD01-004 and MD043-011) and one patient with MPR (MD01-005), we also detected MANA-specific cells upon single-cell profiling of CD8 T cells from TDLN (Extended Data Fig. 10d, e). Tracking the MANA-specific CD8 clonotypes from the primary tumour, we detected those clones among TIL from a brain metastasis resected from patient MD043-011 24 months after primary tumour resection (Extended Data Fig. 10f). Relative to the primary tumour, even-higher levels of three checkpoints—LAG3, TIGIT and HAVCR2—were expressed on MANA-specific TIL in the metastasis (Extended Data Fig. 10g, Supplementary Data 1.5).

14Going back to overall TIL transcriptomic patterns, we hypothesized that MANA-specific T cells and/or a MANA-specific T cell-like signature might correlate with response to ICB, even though total TIL single-cell transcriptomic patterns did not (Fig. 1e). Among CD8 TIL from six MPR tumours and nine non-MPR tumours, the greatest correlation with pathologic response status was observed by combining four TIL clusters most highly enriched in MANA-specific cells, whereas the expression profile of total CD8 TIL did not distinguish MPR from non-MPR (Extended Data Fig. 11). These data suggest that additional T cells with this profile may contribute to the anti-tumour response.

15We next performed scRNA-seq–TCR-seq of serial peripheral blood T cells from patient with MPR MD01-005 after enriching for expression of the TCR-Vβ genes corresponding to this patient’s MANA-specific TCRs using fluorescence-activated cell sorting (FACS) (Fig. 4a–c, Extended Data Fig. 12a). Nine out of ten MANA-specific clones mapped to a TRM-like cluster (Tmem(3); Tmem, memory T cell), with some transcriptional features of TRM, such as expression of HOBIT) two weeks after the initiation of anti-PD-1 treatment (Fig. 4d). By four weeks (time of tumour resection), a significant diversification of phenotype was observed (P ≤ 0.021; Methods). Half of the MANA-specific cells were in Teff clusters (Fig. 4e). By 11 weeks (7 weeks after tumour resection), the MANA-specific cells were below the limit of detection in the blood, consistent with known TRM patterns in the peripheral blood. Using RNA velocity, we observed a clear bidirectional flow of TRM-like memory MANA-specific T cells in the Tmem(3) cluster towards either an activated effector (Teff(3)) or a Tmem(2) transcriptional profile (Fig. 4f). Genes associated with Teff cell function and activation, T cell homing and migration, and tissue retention were upregulated along the pseudotime from Tmem(3) to Teff(3), whereas there was a decrease in genes associated with resting memory T cells (Fig. 4g, h). Gene Ontology (GO) analysis revealed significant enrichment of an IFNγ-mediated signalling pathway along the differentiation trajectory (Extended Data Fig. 12b–f). Although all these tissue compartments were only available for one MPR, these findings are consistent with our hypothesis that, upon activation, functional effector MANA-specific T cells enter the blood and traffic into tissues, including normal lung, in search of micro-metastatic tumour, and are compatible with a previous study showing that TRM cell plasticity can influence systemic memory T cell responses.

16Here we describe the transcriptional programming of MANA-specific TIL after ICB in lung cancer, and further, differential gene programs between patients whose tumours show MPR versus those that do not. Using the MANAFEST platform, MANA-specific CD8 T cells in peripheral blood were detected in the majority of patients who were treated with anti-PD-1; these were also found among TIL in roughly a third of these individuals. Detection of these T cells was independent of tumour response, suggesting that factors in the tumour microenvironment affecting T cell function probably contribute to anti-tumour responsiveness. Indeed, the most frequent MANA-specific clonotype, representing 782 TIL, was observed in a patient with no MPR. This tumour had dual KRAS and STK11 oncogenic mutations, which are known to be highly associated with non-response to PD-1 blockade. Consistent with an earlier study, CD39 expression was a key difference between MANA-specific and viral-specific T cells. Among MANA-specific CD8 TIL, roughly 90% were TRM with high expression of HOBIT that also displayed a partial but incompletely activated Teff program, along with upregulation of several targetable checkpoints in non-MPR tumours. MANA-specific T cells also express far less IL-7R relative to influenza TRM, translating functionally into poor IL-7 responsiveness. These features may all contribute to their limited tumour-specific responsiveness in contrast to anti-viral responses. Future studies are warranted to assess the diminished functional capacity of MANA-specific T cells that was suggested by the transcriptomic profiles observed in our study.

17One hypothesis for the lack of ICB response in some patients is that tumour-specific T cells exhibit low activity owing to poor avidity or affinity of their TCR for its cognate peptide MHC. Our finding comparing the ligand-induced TCR signalling of three MANA-specific TCRs from MPR TIL with seven from patients without MPR supports this notion, although additional studies of this type are necessary to definitively test the hypothesis. An overall limitation of these studies is the modest number of MANA-specific cells among TIL that we were able to detect, representing three responders and three non-responders. Indeed, identification of MANA-specific cells is experimentally challenging, and only a few studies have successfully identified these cells in NSCLC,,,,, yet none of these profiled the transcriptome of MANA-specific T cells at single-cell resolution. Among the 930 MANA-specific transcriptomes that we identified in TIL, there was high consistency among cells from each response group in highly differential expression of key genes known to regulate T cell function. These findings inform on potential ICB combination therapies to overcome anti-PD-1 resistance that occurs even in the presence of potent MANA-specific T cells. For example, our data demonstrated reduced activation of transcriptional programs downstream of IL-7 ligation in MANA-specific TIL relative to influenza-specific TIL, but the MANA-specific TIL retain their ability to respond to supraphysiological levels of IL-7. Because IL-7 signalling is a requisite for maintenance of T cell homeostasis and long-lived memory, it is conceivable that targeting the IL-7 pathway could enhance ICB response. Our findings thus provide a platform for follow-up studies to more rigorously test the generalizability of our conclusions in the setting of resectable and metastatic NSCLC.

18No statistical methods were used to predetermine sample size. The experiments were not randomized. The investigators were not blinded to allocation during experiments and outcome assessment.

19This study was approved by the Institutional Review Boards (IRB) at Johns Hopkins University (JHU) and Memorial Sloan Kettering Cancer Center and was conducted in accordance with the Declaration of Helsinki and the International Conference on Harmonization Good Clinical Practice guidelines. The patients described in this study provided written informed consent. All biospecimens were obtained from patients with stage I-IIIA NSCLC who were enrolled to a phase II clinical trial evaluating the safety and feasibility of administering two doses of anti-PD-1 (nivolumab) before surgical resection. Pathological response assessments of primary tumours were reported previously,. Tumours with no more than 10% residual viable tumour cells were considered to have a MPR.

20Cryobanked T cells were thawed and washed twice with pre-warmed RPMI with 20% FBS and gentamicin. Cells were resuspended in PBS and stained with a viability marker (LIVE/DEAD Fixable Near-IR; ThermoFisher) for 15 min at room temperature in the dark. Cells were the incubated with Fc block for 15 min on ice and stained with antibody against CD3 (BV421, clone SK7) for 30 min on ice. After staining, highly viable CD3+ T cells were sorted into 0.04% BSA in PBS using a BD FACSAria II Cell Sorter. Sorted cells were manually counted using a hemocytometer and prepared at the desired cell concentration (1,000 cells per μl), when possible. The Single Cell 5′ V(D)J and 5′ DGE kits (10X Genomics) were used to capture immune repertoire information and gene expression from the same cell in an emulsion-based protocol at the single-cell level. Cells and barcoded gel beads were partitioned into nanolitre-scale droplets using the 10X Genomics Chromium platform to partition up to 10,000 cells per sample followed by RNA capture and cell-barcoded cDNA synthesis using the manufacturer’s standard protocols. Libraries were generated and sequenced on an Illumina NovaSeq instrument using 2 × 150-bp paired end sequencing. 5′ VDJ libraries were sequenced to a depth of ~5,000 reads per cell, for a total of 5 million to 25 million reads. The 5′ DGE libraries were sequenced to a target depth of ~50,000 reads per cell.

21Genomic data for most individuals in our study were reported previously, and whole-exome sequencing, variant calling and neoantigen predictions for individuals MD043-003 and NY016-025 were performed prospectively for the present study. Whole-exome sequencing was performed on pre-treatment tumours unless otherwise noted (Supplementary Table 4) and matched normal samples. DNA was extracted from tumours and matched peripheral blood using the Qiagen DNA kit (Qiagen). Fragmented genomic DNA from tumour and normal samples was used for Illumina TruSeq library construction (Illumina) and exonic regions were captured in solution using the Agilent SureSelect v.4 kit (Agilent,) according to the manufacturers’ instructions as previously described. Paired-end sequencing, resulting in 100 bases from each end of the fragments for the exome libraries was performed using Illumina HiSeq 2000/2500 instrumentation (Illumina). The depth of total and distinct coverage is shown in Supplementary Table 4. Somatic mutations, consisting of point mutations, insertions, and deletions across the whole exome were identified using the VariantDx custom software for identifying mutations in matched tumour and normal samples as previously described,. Somatic mutations, consisting of nonsynonymous single base substitutions, insertions and deletions, were evaluated for putative MHC class I neoantigens using the ImmunoSelect-R pipeline (Personal Genome Diagnostics) as previously described. Somatic sequence alterations are listed in Supplementary Table 5.

22We used the MANAFEST assay to evaluate T cell responsiveness to MANA and viral antigens. In brief, pools of MHC class I-restricted CMV, EBV and influenza peptide epitopes (CEFX, jpt Peptide Technologies), pools representing the matrix protein and nucleoprotein from H1N1 and H3N2 (jpt Peptide Technologies), and putative neoantigenic peptides defined by the ImmunoSelect-R pipeline (jpt Peptide Technologies; Supplementary Table 6) were each used to stimulate 250,000 T cells in vitro for 10 days as previously described. The time point of peripheral blood collection used for each MANAFEST assay is described in Supplementary Tables 2, 7. In brief, on day 0, T cells were isolated from PBMC by negative selection (EasySep; STEMCELL Technologies). The T cell-negative fraction was co-cultured with an equal number of selected T cells in culture medium (IMDM/5% human serum with 50 μg ml−1 gentamicin) with 1 μg ml-1 relevant neoantigenic peptide, 1 μg ml−1 of an MHC class I-restricted CMV, EBV, and influenza peptide epitope pool (CEFX, jpt Peptide Technologies), 1 μg ml−1 of pools representing the matrix protein and nucleoprotein from H1N1 and H3N2 (jpt Peptide Technologies), or no peptide (to use as a reference for non-specific or background clonotypic expansion). On day 3, half the medium was replaced with fresh medium containing cytokines for a final concentration of 50 IU ml−1 IL-2 (Chiron), 25 ng ml−1 IL-7 (Miltenyi) and 25 ng ml−1 IL-15 (PeproTech). On day 7, half the medium was replaced with fresh culture medium containing cytokines for a final concentration of 100 IU ml−1 IL-2 and 25 ng ml−1 IL-7 and IL-15. On day 10, cells were harvested, washed twice with PBS, and the CD8+ fraction was isolated using a CD8+ negative enrichment kit (EasySep; STEMCELL Technologies). DNA was extracted from each CD8-enriched culture condition using the Qiamp micro-DNA kit according to the manufacturer’s instructions. TCR sequencing was performed on each individual peptide-stimulated T cell culture using survey-level sequencing (max depth ~60,000 reads) by Adaptive Biotechnologies using their established platform or by the Sidney Kimmel Comprehensive Cancer Center FEST and TCR Immunogenomics Core (FTIC) facility using the Oncomine TCR Beta short-read assay (Illumina) and sequenced on an Illumina iSeq 100 using unique dual indexes, for a maximum of ~40,000 reads per sample.

23Data pre-processing was performed to eliminate non-productive TCR sequences and to align and trim the nucleotide sequences to obtain only the CDR3 region. Sequences not beginning with C or ending with F or W and having less than seven amino acids in the CDR3 were eliminated. TCR sequencing samples with less than 1,000 productive reads were excluded from downstream analysis. MD043-011-MANA_22 was the only such sample in the present study (see Supplementary Table 7). Resultant processed data files were uploaded to our publicly available MANAFEST analysis web app (http://www.stat-apps.onc.jhmi.edu/FEST) to bioinformatically identify antigen-specific T cell clonotypes.

24Bioinformatic analysis of productive clones was performed to identify antigen-specific T cell clonotypes meeting the following criteria: (1) significant expansion (Fisher’s exact test with Benjamini–Hochberg correction for false discovery rate (FDR), P < 0.05) compared to T cells cultured without peptide, (2) significant expansion compared to every other peptide-stimulated culture (FDR <0.05) except for conditions stimulated with similar neoantigens derived from the same mutation, (3) an odds ratio >5 compared to the no peptide control, and (4) present in at least 10% of the cultured wells to ensure adequate distribution among culture wells. A lower read threshold of 300 was used for assays sequenced by the FTIC and a lower threshold of 30 was used for samples sequenced by Adaptive Biotechnologies. In MANAFEST assays testing less than 10 peptides or peptide pools, cultures were performed in triplicate and reactive clonotypes were defined as being significantly expanded relative to T cells cultured without peptide (FDR <0.05) in two out of three triplicates, and not significantly expanded in any other well tested. When available, TCRseq was also performed on DNA extracted from tumour, normal lung, and lymph node tissue obtained before treatment and at the time of surgical resection, as well as serial peripheral blood samples. The assays performed on each biospecimen are outlined in Supplementary Table 2.

25Peptide affinity for cognate HLA molecules was assessed using a luminescent oxygen channeling immunoassay (LOCI; AlphaScreen, Perkin Elmer) as previously described. This is a proximity-based system using a donor and acceptor bead, each conjugated with an epitope tag. When the donor bead is excited with light at 650 nm and can activate an acceptor bead, resulting in a signal at 520–620 nm, which can be quantified per second as a surrogate of affinity. A higher number of counts per second indicates higher affinity of the peptide:HLA pair. The stability of peptide loaded complexes was measured by refolding MHC with peptide and subsequently challenging complexes with a titration of urea. The denaturation of MHC was monitored by ELISA as described previously.

26Ten MANAFEST+ TCR sequences for which the TCRα chain could be enumerated (>3 cells in single-cell data with the same TCRα–TCRβ pair) were selected for cloning. In addition, seven clones (from three individuals: MD01-004, MD01-005 and MD043-011) that have high composite signature (using the AddModuleScore function) consisting of differential gene programs of MANA-specific T cell relative to influenza-specific T cells in the TRM were selected for cloning. Relevant TCRs were analysed with the IMGT/V-Quest database (http://www.imgt.org). The database allows us to identify the TRAV and TRBV families with the highest likelihood to contain the identified segments which match the sequencing data. To generate the TCRs, the identified TCRA V-J region sequences were fused to the human TRA constant chain, and the TCRB V-D-J regions to the human TRB constant chain. The full-length TCRA and TCRB chains were then synthesized as individual gene blocks (IDT) and cloned into the pCI mammalian expression vector, containing a CMV promoter, and transformed into competent Escherichia coli cells according to the manufacturer’s instructions (NEBuilder HiFi DNA Assembly, NEB). Post transformation and plasmid miniprep, the plasmids were sent for Sanger sequencing to ensure no mutations were introduced (Genewiz).

27To generate a Jurkat reporter cell in which we could transfer our TCRs of interest, the endogenous TCR α- and β-chains were knocked out of a specific Jurkat line that contains a luciferase reporter driven by an NFAT response element (Promega) using the Alt-R CRISPR system (Integrated DNA Technologies, IDT). Two sequential rounds of CRISPR knockout were performed using crDNA targeting the TCRα constant region (AGAGTCTCTCAGCTGGTACA) and the TCRβ constant region (AGAAGGTGGCCGAGACCCTC). Limiting dilution was then used to acquire single cell clones and clones with both TCRα and TCRβ knocked out, as confirmed by Sanger sequencing and restoration of CD3 expression only by the co-transfection of TCRα or TCRβ chains, were chose. CD8α and CD8β chains were then transduced into the TCRα−TCRβ− Jurkat reporter cells using the MSCV retroviral expression system (Clontech). Jurkat reporter cells were then co-electroporated with the pCI vector encoding the TCRB and TCRA gene blocks, respectively, using ECM830 Square wave electroporation system (BTX) at 275 V for 10 ms in OptiMem media in a 4-mm cuvette. Post electroporation, cells were rested overnight by incubating in in RPMI 10% FBS at 37 °C, 5% CO2. TCR expression was confirmed by flow cytometric staining for CD3 on a BD FACSCelesta and 50,000 CD3+ T cells were plated in each well of a 96-well plate. Reactivity of the TCR-transduced Jurkat cells was assessed by co-culturing with 1 × 105 autologous EBV-transformed B cells, loaded with titrating concentrations of MANA peptides, viral peptide pools or negative controls. After overnight incubation, activation of the NFAT reporter gene was measured by the Bio-Glo Luciferase Assay per manufacturer’s instructions (Promega). Jurkat cells were routinely tested for mycoplasma contamination. No cell line authentication was performed.

28gBlocks (IDT) encoding HLA A*68:01, p53(R248L) and wild-type p53 were cloned into pcDNA3.4 vector (Thermo Fisher Scientific, A14697). COS-7 cells were transfected with plasmids at 70–80% confluency using Lipofectamine 3000 (Thermo Fisher Scientific, L3000015) and incubated at 37 °C overnight in T75 flasks. A total of 30 μg plasmid (1:1 ratio of HLA plasmid per target protein plasmid in co-transfections) was used. Post transfection, COS-7 cells were plated with TCRαβ-transfected JurkaT cells containing NFAT reporter gene at a 1:1 ratio. After overnight incubation, activation of the NFAT reporter gene was measured by the Bio-Glo Luciferase Assay per manufacturer’s instructions (Promega).

29Cell Ranger v3.1.0 was used to demultiplex the FASTQ reads, align them to the GRCh38 human transcriptome, and extract their cell and unique molecular identifier (UMI) barcodes. The output of this pipeline is a digital gene expression (DGE) matrix for each sample, which records the number of UMIs for each gene that are associated with each cell barcode. The quality of cells was then assessed based on (1) the number of genes detected per cell and (2) the proportion of mitochondrial gene/ribosomal gene counts. Low-quality cells were filtered if the number of detected genes was below 250 or above 3× the median absolute deviation away from the median gene number of all cells. Cells were filtered out if the proportion of mitochondrial gene counts was higher than 10% or the proportion of ribosomal genes was less than 10%. For single-cell VDJ sequencing, only cells with full-length sequences were retained. Dissociation/stress associated genes,, mitochondrial genes (annotated with the prefix ‘MT-’), high abundance lincRNA genes, genes linked with poorly supported transcriptional models (annotated with the prefix ‘RP-’) and TCR (TR) genes (TRA/TRB/TRD/TRG, to avoid clonotype bias) were removed from further analysis. In addition, genes that were expressed in less than five cells were excluded.

30Seurat (3.1.5) was used to normalize the raw count data, identify highly variable features, scale features, and integrate samples. PCA was performed based on the 3,000 most variable features identified using the vst method implemented in Seurat. Gene features associated with type I Interferon (IFN) response, immunoglobulin genes and specific mitochondrial related genes were excluded from clustering to avoid cell subsets driven by the above genes. Dimension reduction was done using the RunUMAP function. Cell markers were identified by using a two-sided Wilcoxon rank sum test. Genes with adjusted P <0.05 were retained. Clusters were labelled based on the expression of the top differential gene in each cluster as well as canonical immune cell markers. Global clustering on all CD3 T cells and refined clustering on CD8 T cells were performed using same procedure. To select for CD8+ T cells, SAVER was used to impute dropouts by borrowing information across similar genes and cells. A density curve was fitted to the log2-transformed SAVER imputed CD8A expression values (using the ‘density’ function in R) of all cells from all samples. A cut-off is determined as the trough of the bimodal density curve (that is, the first location where the first derivative is zero and the second derivative is positive). All cells with log2-transformed SAVER imputed CD8A expression larger than the cut-off are defined as CD8+ T cells. TRB amino acid sequences were used as a biological barcode to match MANA, EBV or influenza A-specific T cell clonotypes identified from the FEST assay with single-cell VDJ profile and were projected onto CD8+ T cell refined UMAP.

31PCA was performed on a standardized pseudobulk gene expression profile, where each feature was standardized to have a mean of zero and unit variance. In global CD3 and CD8 TIL PCA, for each cell cluster we first aggregated read counts across cells within the cluster to produce a pseudobulk expression profile for each sample and normalized these pseudobulk expression profiles across samples by library size. Combat function in the sva R package, was applied to address potential batch effects on the normalized pseudobulk profile. Highly variable genes (HVGs) were selected for each cell cluster by fitting a locally weighted scatterplot smoothing (LOESS) regression of standard deviation against the mean for each gene and identifying genes with positive residuals. For each sample, all cell clusters were then concatenated by retaining each cluster’s HVGs to construct a concatenated gene expression vector consisting of all highly variable features identified from different cell clusters. Each element in this vector represents the pseudobulk expression of a HVG in a cell cluster. Samples were embedded into the PCA space based on these concatenated gene expression vectors. Canonical correlation, between the first two PCs (that is, PC1 and PC2) and a covariate of interest (that is, tissue type or response status) was calculated. Permutation test was used to assess the significance by randomly permuting the sample labels 10,000 times. In the MANA-specific PCA (Extended Data Fig. 11), MANA-enriched cell clusters, defined by clusters with MANA-specific T cell frequency at least two fold higher than randomly expected, were aggregated as one combined cell cluster. Then, a similar procedure by first identifying HVGs, computing the first 2 PCs and then calculating the canonical correlation was repeated for the combined MANA-enriched cell cluster and each of the other CD8 clusters.

32The gene expression read counts were adjusted by library size. SAVER was used to impute the dropouts, and further log2-transformed the imputed values after adding a pseudocount of 1. A linear mixed-effect model was constructed to identify genes that are significantly differential between MPR and non-MPR among total CD8/CD4 TIL and by each cell cluster, respectively. The B-H procedure was used to adjust the P values for multiple testing, and the statistical significance is determined using a cut-off of FDR <0.05.

33Differential-expression tests for antigen-specific T cells were performed using FindAllMarkers functions in Seurat with Wilcoxon rank-sum test on SAVER imputed expression values. Genes with >0.25 log2-fold changes, at least 25% expressed in tested groups, and Bonferroni-corrected p values <0.05 were regarded as significantly differentially expressed genes (DEGs). Antigen-specific (MANA versus influenza versus EBV) T cell marker genes were identified by applying the differential expression tests for upregulated genes between cells of one antigen specificity to all other antigen-specific T cells in the dataset. MANA-specific T cell genes associated with response to ICB were identified by applying the differential expression tests comparing MANA-specific T cells from MPR versus those from non-MPR. Top ranked DEGs (by log-fold changes) with a log2-fold changes >0.8 and DEGs relating to T cell function were extracted for further visualization in a heat map using pheatmap package. SAVER-imputed expression values of selective marker genes (transcriptional regulators, memory markers, tissue-resident markers, T cell checkpoints, effector and activation markers) were plotted using the RidgePlot function in Seurat.

34Cryopreserved TIL from patient MD01-004 were thawed, counted, and stained with the LIVE/DEAD Fixable Aqua (ThermoFisher) viability marker and antibodies specific for CD3 (PE, clone SK1) and CD8 (BV786, clone RPA-T8). Thirty-thousand CD8+ T cells per condition were sorted on a BD FACSAria II Cell Sorter into a 96-well plate. Autologous T cell-depleted PBMC were added as antigen presenting cells (APC) at 1:1 ratio. The cells were stimulated with either influenza A or MD01-004-MANA 12 peptide and titrating concentrations of recombinant human IL-7 (Miltenyi) for 12 h in a round-bottomed 96-well plate.

35Following 12 h of antigen and IL-7 stimulation, cells were spun down, counted and re-suspended in 1% BSA at desired concentration. scRNA-seq and VDJ libraries were prepared using 10X Chromium single cell platform using 5′ DGE library preparation reagents and kits according to manufacturer’s protocols (10X Genomics) and as described above. MANA- or influenza-specific T cell clonotypes from the single-cell dataset were identified by using the TRB amino acid sequences as a biological barcode. SAVER imputed gene expression was scaled and centred using the ScaleData function in Seurat. A composite score for the IL-7-upregulated gene set (Supplementary Data 2.3) expression was computed using the AddModuleScore function and subsequently visualized using ridgeplot. Mean ± standard error was used to show the dose–response curve of the IL-7-upregulated gene-set score by antigen-specific T cells and peptide-stimulation groups.

36To characterize dysfunctional CD8 MANA TIL, six best-characterized (and clinically targeted) checkpoints: CTLA4, PDCD1, LAG3, HAVCR2, TIGIT and ENTPD1, were used to compute the T cell checkpoint score, and a published gene list from exhausted T cells was used to compute the T cell exhaustion score, using AddModuleScore function in Seurat. Applying T cell checkpoint score as an anchor, genes that were maximally correlated to the score were identified using linear correlation in MANA-specific TIL from MPR and non-MPR, respectively. Top-30 genes (from HVG selected using FindVariableGenes function in Seurat and excluded the 6 genes included in immune checkpoint score generation) with the highest correlation coefficients were plotted as a bar plot. The difference of correlation coefficients of the above genes was additionally computed between MPR and non-MPR and visualized using waterfall plot.

37Peripheral blood T cells from patient MD01-005 were sorted based on expression of CD8 and TCRVβ2, followed by scRNA-seq–TCR-seq and clustering on conventional CD8+ T cells (MAIT cells excluded). To evaluate whether there was a statistically significant change in the cell types of MANA cells between week 2 (W2) and week 4 (W4) samples in Fig. 4d, e, we first conducted a Fisher’s exact test, which yields a P = 0.021, indicating a statistically significant phenotype change in MANA-specific cells (Fig. 4e). We also conducted a more sophisticated test that adjust for potential background differences in cell type abundance between W2 and W4 samples. In this test, we let mc,t denote the probability that a MANA-specific T cell collected at time point t (W2 or W4) comes from cell type c, and let pc,t denote the proportion of all cells in time point t that come from cell type c. We evaluated the ratio Rc,t = mc,t/pc,t, which characterizes the relative abundance of MANA-specific T cells in each cell type. We compared the null model where this ratio does not change over time (H0: Rc,W2 = Rc,W4 for all cell type c) versus the alternative model where W2 and W4 T cells have different ratios (H1: Rc,W2 ≠ Rc,W4). To do this, we computed the test statistic using the observed data and compared it to its null distribution obtained using Monte Carlo simulations. To construct the null distribution for, we pooled cells from W2 and W4 together and treated them as one sample to estimate the common ratio Rc,W2 = Rc,W4 = Rc shared by W2 and W4, and then derived the probability that a MANA-specific T cell collected at time point t comes from cell type c under the null model H0, which is proportional to pc,t Rc (that is, the product of the sample-specific background cell type proportion pc,t and the common MANA-abundance ratio Rc shared between samples). The MANA-specific T cells at time point t were then redistributed to different cell types randomly based on a multinomial distribution with this expected MANA-specific T cell type proportion (that is, the expected probability that a MANA-specific T cell at time point t comes from cell type c under H0 is ), while keeping the total number of MANA-specific T cells at each time point the same as the observed MANA-specific T cell number at that time point. The test statistic S was then computed using this simulated sample. We repeated this simulation 10,000 times to derive the null distribution of S. Comparing the observed S to its null distribution yields a P < 10−4.

38The RNA velocity analysis was performed by first recounting the spliced reads and unspliced reads based on aligned bam files of scRNA-seq data using the velocyto Python package. The calculation of RNA velocity values for each gene in each cell and embedding RNA velocity vector to low-dimension space were done using the SeuratWapper workflow for estimating RNA velocity using Seurat (https://github.com/satijalab/seurat-wrappers/blob/master/docs/velocity.md). The first two diffusion components from Diffusion map were used to construct the coordinates along with velocity. TSCAN (v.1.7.0) was used to reconstruct the cellular pseudotime on diffusion maps space for the PBMC T cells from three time points (samples) of one patient (MD01-005). Based on velocity analysis, the Tmem(3) cluster was specified as the starting cluster for the pseudotemporal trajectory which has branches. For each branch, log2-transformed and library size-normalized SAVER-imputed gene expression values were used for analysing gene expression dynamics along the pseudotime. 10,325 genes with normalized expression ≥ 0.01 in at least 1% of cells were retained. For each gene g, the gene expression along pseudotime t in each sample s was described as a function which was obtained by fitting B-spline regression to the gene’s normalized expression values in single cells. The red curves in Fig. 4h are the mean of the function of the three samples. In order to test whether the gene expression shows a significant change along pseudotime, we compared the above model with a null model in which is assumed to be a constant over time. The likelihood ratio statistic between the two models was computed. To determine the P value, the null distribution of the likelihood ratio statistic was constructed by permuting the pseudotime of cells in each sample, refitting the models and recomputing the likelihood ratio statistic. The P value was calculated as the number of permutations out of a total of 1,000 permutations that produce a likelihood ratio statistic larger than the observed one. The P values from all genes were converted to FDR by Benjamini–Hochberg procedure to adjust for multiple testing. Genes with FDR < 0.05 were considered as dynamic genes with statistical significance. k-Means clustering was applied to group genes with similar dynamic expression patterns into clusters. topGO (v.2.42.0) was used to identify the enriched Gene Ontology terms by comparing the genes in each cluster to all 10,325 genes as background.

39Further information on research design is available in the Nature Research Reporting Summary linked to this paper.

40Any methods, additional references, Nature Research reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at 10.1038/s41586-021-03752-4.