Single-Cell Transcriptome Analysis Reveals Gene Signatures Associated with T-cell Persistence Following Adoptive Cell Therapy

Adoptive cell therapy (ACT) using tumor-infiltrating lymphocytes (TIL) has shown significant efficacy in patients with metastatic melanoma (1). Striking responses have been observed in individual patients with metastatic cholangiocarcinoma, cervical cancer, colorectal cancer, or breast cancer after ACT (2–5). However, clinical responses to ACT have been limited in the majority of patients with solid epithelial cancers (3, 4, 6). To increase the efficacy of the current treatments, multiple aspects of ACT are being explored, including the selection of effective targets and the avidity of antigen-reactive T cells (7, 8). Our previous studies have shown that neoantigen recognition by T cells is a major factor associated with effective TIL-ACT (2, 9–11). Other studies have also demonstrated that effective immune checkpoint blockade immunotherapies correlate with the abundance of neoantigens expressed by tumors (12–15). Despite efforts to expand and transfer selective TIL populations recognizing neoantigen targets, many of the TIL treatments have been ineffective (6, 16).

The persistence of specific T-cell clonotypes in the peripheral blood after T-cell infusion has become an important parameter to monitor during the course of ACT (17). In addition to the antigen avidity of T cells, several factors can also influence T-cell persistence in the peripheral blood. Examples include the levels of homeostatic cytokines, endogenous T-cell reconstitution, tumor antigen expression, and the tumor microenvironment (17–20). Therefore, it is difficult to study persistence data because of many potential variables, unless the analysis is from a large number of patients treated with nearly the same cell product (17).

In a previous study, we observed that patient-4095 with metastatic colorectal cancer, who was treated with mutated KRAS(G12D)-reactive TILs, experienced a partial response for 9 months after ACT (4). A single nonresponding pulmonary lesion was resected, and the patient has remained disease-free for 3 years after treatment. We found that the infused TILs comprised three dominant clonotypes (49.5%, 19.1%, and 6.9%, respectively). The first clonotype, 9.1-NP, recognized a KRAS(G12D) 9-mer epitope (GADGVGKSA); the second clonotype, 10-P, recognized a KRAS(G12D) 10-mer epitope (GADGVGKSAL); and the third clonotype, 9.2-P, recognized the same 9-mer epitope as 9.1-NP. The most dominant clonotype, 9.1-NP, did not persist (NP) at all on day 40, but the 10-P and 9.2-P clonotypes continued to persist (P) up to day 266 after ACT. The 9.1-NP and 9.2-P clonotypes shared the same specificity and similar T-cell receptor α/β (TCRα/β) sequences, except one amino acid at the CDR3α region and six amino acids at the CDR3β region, and no significant difference in TCR avidity was observed (4).

Because these three clonotypes shared the same in vitro culture environment and in vivo host environment, we sought to investigate whether different gene expression profiles were associated with different durations of persistence after ACT. Because the majority of biological assays cannot distinguish the difference between 9.1-NP and 9.2-P T cells, we performed single-cell TCR and transcriptome analysis for this single-patient study. The results suggested that TILs that could persist in patient-4095 had distinguishable gene expression profiles, including key genes encoding surface markers and transcription factors.

Patients were enrolled in a clinical trial of TIL therapy (ClinicalTrial.gov ID: NCT01174121). This trial was approved by the Institutional Review Board of the NCI, and the written-informed consent was obtained from the patients, following NIH guidelines and Declaration of Helsinki.

The characteristics, treatment, and clinical response for patient-4095 with metastatic colorectal cancer have been published previously (4). We have also reported the summary of characteristics for patient-4007, -4071, and -4081 with metastatic colorectal cancer and patient-4069 with pancreatic cancer (6, 16). Briefly, TILs were generated from the metastatic tumors of patients. TIL cultures were selected based on the reactivities against tumor-specific mutations, and selected TIL cultures were expanded for treatment. The patients were treated with a lymphodepleting chemotherapy regimen, a single infusion of TILs, followed by several doses of IL2. The peripheral blood lymphocyte (PBL) samples were obtained from the patients every 2 to 3 days during hospitalization and during follow-up visits. PBL samples were cryopreserved and stored in a liquid nitrogen container before use. The percentage of individual T-cell clonotypes in samples was obtained by an Immunoseq Assay service, provided by Adaptive Biotechnologies.

TILs used for this study were generated by methods described previously (21). Briefly, metastatic tumors were resected from patients, and tumor fragments were excised and cultured in RPMI medium supplemented with 10% in-house human AB serum, 2 mmol/L l-glutamine, 25 mmol/L HEPES, gentamicin (10 μg/mL), and IL2 (6,000 IU/mL; Clinigen). TIL cultures were grown for 2 to 4 weeks and then screened for recognition of tumor-specific mutations (22). The screening results for patient-4095, -4007, -4071, -4081, and -4069 have been published (4, 6, 16). The mutation-reactive TIL cultures were selected and expanded using a rapid expansion protocol (REP) to large numbers for patient infusion (23). The REP culture contained 5 × 10⁶ TILs and 5 × 10⁸ irradiated PBMC feeder cells in RPMI/AIM V medium (50%/50% mixed), supplemented with 7.5% in-house human AB serum, 2 mmol/L l-glutamine, IL2 (3,000 IU/mL), and OKT3 antibody (30 ng/mL; Miltenyi Biotec) in a G-Rex100 flask (Wilson Wolf). A small portion of TILs were cryopreserved and stored in a liquid nitrogen container for the experiments shown in this report.

The samples of infused TILs were thawed and recovered overnight in RPMI/AIM V medium (50%/50%) supplemented with 5% human AB serum (Valley Biomedical). TILs were resuspended in PBS and then subjected to a 10X Chromium instrument (10X Genomics) for the single-cell analysis, as described below. Note that 1 × 10⁷ PBLs from patient-4095 were stained with KRAS-9mer tetramer or KRAS-10mer tetramer (NIH tetramer core facility), together with CD8 antibody (clone RPA-T8, BD Biosciences) for 40 minutes. Stained cells were washed twice with PBS containing 5% FBS (SAFC). After washes, tetramer-positive cells were sorted by BD FACSAria II and subjected to a 10X Chromium instrument for the single-cell analysis.

We used the standard protocol and reagent kit for single-cell V(D)J analysis, provided by 10X Genomics. Up to 8 samples/reactions in 8 channels can be processed simultaneously by a 10X Chromium instrument. Briefly, 10,000 cells per reaction/channel were loaded, with the targeted cell recovery of 6,000 cells. For TIL4095, a total of 4 channels were loaded to obtain the desired number of cells. For patient-4095′s PBLs on day 12 after ACT, KRAS-9mer tetramer⁺ cells were loaded on 3 channels and KRAS-10mer tetramer⁺ cells were loaded on 1 channel. For patient-4095′s PBLs on day 40 after ACT, KRAS-9mer tetramer⁺ cells and KRAS-10mer tetramer⁺ cells were loaded on 2 channels each. Single cells were captured and lysed, and mRNA was reverse transcribed to cDNA using the provided reagents (10X Genomics). Single-cell cDNA samples were barcoded and amplified by 12 PCR cycles. The amplified cDNA samples were divided into two fractions for use in two analyses. The first fraction was target-enriched for TCRs (10X Genomics proprietary), and barcoded samples were pooled and sequenced by an Illumina NextSeq 550 sequencer (High Output v2 kit, Read1: 150 b.p., Read2: 150 b.p.). The second fraction was processed for 5′ gene expression libraries by following the manufacturer's instruction (10X Genomics). Each channel was processed and sequenced separately by an Illumina NextSeq 550 sequencer (High Output v2 kit, Read1: 26 b.p., Read2: 98 b.p.).

TIL4095 (TILs from patient-4095) single cells were processed by two identical Fluidigm C1 systems (Fluidigm), following the protocol published previously (24). Briefly, the infused TIL sample from patient-4095 was thawed and recovered overnight in RPMI/AIM V medium (50%/50%) supplemented with 5% human AB serum. TILs were stained with KRAS-9mer tetramer or KRAS-10mer tetramer, together with CD8 antibody for 40 minutes. After washes, tetramer-positive cells were sorted by BD FACSAria II and subjected to Fluidigm C1 systems for the single-cell analysis. Twelve thousand KRAS-9mer tetramer⁺ cells or KRAS-10mer tetramer⁺ cells were loaded separately on an IFC plate (Fluidigm). Because of the low abundance of 9.2-P clonotype, KRAS-9mer tetramer⁺ cells were loaded on 6 additional IFC plates, and 9.2-P single cells were selected based on the CDR3 sequence information. Single-cell samples were processed by the Fluidigm C1 system and barcoded by Nextera XT DNA Library Preparation Kit (Illumina). Single-cell samples were first sequenced by an Illumina MiSeq sequencer (High Output v3 kit, Read1: 250 b.p., Read2: 250 b.p.) to identify CDR3 sequences. Single-cell samples with correct CDR3 sequences (9.1-NP, 9.2-P, or 10-P) were sequenced again by an Illumina HiSeq2500 sequencer (Read1: 100 b.p., Read2: 100 b.p.), in order to obtain sufficient sequence information for data analysis.

Single-cell sequencing data were first processed by Cell Ranger pipelines (v2.1.1; 10X Genomics), in order to obtain T-cell clonotypes and gene expression profiles associated with individual single-cell barcodes. T-cell clonotypes were defined based on TCR Vβ CDR3 nucleotide sequences. Single cells with two different TCR Vβ CDR3 nucleotide sequences were likely doublets and were excluded. Low-quality single cells, with less than 200 detectable genes, were also excluded. Genes with less than 10 read counts, as well as mitochondrial genes, were removed from the subsequent analysis. Noncoding RNA genes were not investigated in this study.

To construct the gene expression matrix for the datasets, we merged all uniquely aligned reads with the same unique molecular identifiers (UMI), cell barcode, and gene annotation in the HDF5 format. The raw single-cell data were then normalized with convolution methods (25). The normalized data were decomposed using randomized principal component analysis (26). The principal components were used to generate t-distributed stochastic neighbor embedding (t-SNE) projections, which were implemented with R package Seurat v.2.3 (https://github.com/satijalab/seurat). The differential expression comparisons were generated with edgeR package with selected genes (FDR < 0.05).

The Persistence_Up and Persistence_Down gene sets were constructed based on the shared, differentially expressed genes (FC > ±1.5, FDR < 0.05), obtained from the comparisons of 9.1-NP versus 9.2-P and 9.1-NP versus 10-P (Supplementary Table S1). Gene set enrichment analysis (GSEA; ref. 27) software/scripts and instruction are publicly available from Broad Institute (http://software.broadinstitute.org/gsea/index.jsp). We made necessary in-house modifications for the scripts, in order to execute the GSEA scripts on the NIH Biowulf server ((http://hpc.nih.gov).

Single-cell whole-transcriptome sequences obtained by the Fluidigm C1 were first analyzed by an in-house bioinformatics pipeline to obtain TCR-CDR3 nucleotide sequences (24). The sequencing data of single cells with correct CDR3 were further analyzed by a bioinformatics pipeline, including the alignment using STAR (v2.4.2) and raw counts by subread (v1.4.6). Good-quality single cells were defined as cells in which more than 1,500 genes were expressed and detected. Individual genes expressed in more than two cells were selected for the subsequent analysis. The gene expression matrix was generated, and the P values were calculated by the Wilcoxon signed-rank test.

The data discussed in this publication have been deposited in NCBI's Gene Expression Omnibus (28) and are accessible through GEO Series accession number GSE136394 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE136394).

Statistics for single-cell data and GSEA are described in the above sections.

To study the gene expression profiles of persistent and nonpersistent clonotypes, the infused TILs from patient-4095 were subjected to single-cell TCR/transcriptome sample preparation by a 10X Genomics single-cell instrument, followed by next-generation sequencing. Each single cell contained a unique barcode, which could be used to link the gene expression profiles to one of the three dominant clonotypes (9.1-NP, 9.2-P, or 10-P), based on their unique TCR-CDR3β sequences. After filtering out single cells with low-quality data, we obtained 9,707 single cells from the 9.1-NP clonotype, 1,737 single cells from the 9.2-P clonotype, and 2,771 single cells from the 10-P clonotype for the downstream bioinformatics analysis (Fig. 1A).

Because 9.1-NP and 9.2-P clonotypes recognized exactly the same 9-mer epitope with similar avidity, we first compared the gene expression profiles of single cells from nonpersistent 9.1-NP and persistent 9.2-P clonotypes. We identified a total of 472 differentially expressed genes with FDR less than 0.05 (Fig. 1B and C; Supplementary Table S2). By comparing the fold changes (FC) of gene expression between 9.1-NP and 9.2-P, we identified the top 10 upregulated genes and the top 10 downregulated genes. These genes are listed in Fig. 1D, and these 20 genes are also highlighted in the heatmap and volcano plot (Fig. 1B and C). The brief descriptions of their known functions are listed in Supplementary Table S3, based on the summaries obtained from online databases Genecards (www.genecards.org) and OMIM (www.omim.org). Gene functions in 18 out of the 20 genes have been reported in previous studies, and 9 of them have been shown to function significantly in T cells. However, the biological roles of C1orf162 and TPRG1 have been largely unknown (Fig. 1D).

We then compared the single-cell gene expression profiles from nonpersistent 9.1-NP and persistent 10-P single cells, with the caveat that 10-P clonotype recognized a longer epitope with one additional leucine at the C-terminus. Similarly, a total of 528 differentially expressed genes (FDR < 0.05) were identified in this comparison (Fig. 1E and F; Supplementary Table S2). The top 10 upregulated genes and the top 10 downregulated genes based on the FC of gene expression are listed in Fig. 1G. These 20 genes are also highlighted in the heatmap and volcano plot (Fig. 1E and F). The brief descriptions of their known functions are listed in Supplementary Table S4. Gene functions in all 20 genes have been reported in previous studies, and 9 of them have been shown to function significantly in T cells (Fig. 1G).

By assessing these two sets of differentially expressed genes, a total of 411 genes were shared between these two comparisons (Fig. 2A). A total of 11 genes in the top upregulated and downregulated gene lists were shared between both groups (bold text, Fig. 1D and G). Because these gene signatures were also found in another persistent clonotype 10-P, these genes might be more significant in the biology of persistence than the 61 unique genes in the 9.2-P versus 9.1-NP comparisons. Using the shared 411 genes identified from these two comparisons, a t-SNE plot was generated. The majority of 9.2-P and 10-P single cells were separated from the 9.1-NP single cells (Fig. 2B).

Among these shared gene signatures, we focused on cell surface markers and transcription factors in the top 10 upregulated and downregulated gene lists (Fig. 1D and G). The cell surface markers could be potentially utilized to monitor or enrich T-cell populations associated with the persistence phenotype based on antibody staining, whereas the transcription factors could be the potential key factors defining the persistence phenotype. As shown in Fig. 2C, the gene expression of granzyme K (GZMK) was much higher in the 9.1-NP nonpersistent single cells, compared with the 9.2-P persistent single cells (log₂FC = −3.48, FDR = 0). Similarly, the GZMK expression was much higher in the 9.1-NP single cells, compared with the 10-P persistent single cells (log₂FC = −3.55, FDR = 0; Fig. 2C). The expressions of other granzymes, including GZMB, GZMH, and GZMM, were also significantly higher in the 9.1-NP, compared with the 9.2-P or 10-P single cells (Supplementary Table S2). These results suggested that 9.1-NP single cells showed a more effector-like phenotype, compared with 9.2-P or 10-P single cells. Conversely, natural killer cell marker KLRB1 (killer cell lectin-like receptor subfamily B, member 1, CD161) had high expression in 9.2-P single cells, compared with 9.1-NP single cells (log₂FC = 3.43, FDR = 0; Fig. 2D). The 10-P single cells were also expressed higher KLRB1 than 9.1-NP single cells, but the difference was less than the 9.1-NP versus 9.2-P comparison (log₂FC: 3.43 vs. 1.69). The integrin beta 1 (ITGB1/CD29) gene also had high expression in both 9.2-P and 10-P single cells, compared with 9.1-NP single cells (log₂FC = 1.97 and 1.83, respectively; Fig. 2E). In these comparisons, we also found other differentially expressed genes related to T-cell mobility, including integrin ITGA1 (integrin alpha-1, CD49a), ITGAX (integrin alpha X, CD11c), chemokine CCL3, XCL2, and chemokine receptor CXCR4 and CCR5 (Supplementary Table S2).

IL7 has been known as one of the master regulators for T-cell homeostasis (20). Higher expression of IL7 receptor (IL7R/CD127) was observed in 9.2-P single cells, compared with 9.1-NP single cells (log₂FC = 1.31, FDR = 5.19 × 10⁻⁷³; Fig. 2F). Similarly, higher IL7R expression was observed in 10-P single cells, compared with 9.1-NP single cells (log₂FC = 1.51, FDR = 2.03 × 10⁻¹⁵²). These results suggested that these 9.2-P T cells, as well as 10-P T cells, might be better equipped to compete for the limited resource of IL7, leading to better persistence (20). Sphingosine-1-phosphate receptor 1 (S1PR1) has been described as the predominant receptor for S1P, which is the central mediator of T-cell egress from the thymus or lymphoid organs into the blood (29). The expression of S1PR1 was higher in 9.2-P single cells, compared with 9.1-NP single cells (log₂FC = 1.22, FDR = 5.38 × 10⁻⁷²; Fig. 2G). Similarly, higher S1PR1 expression was observed in 10-P single cells, compared with 9.1-NP single cells (log₂FC = 1.60, FDR = 1.47 × 10⁻²⁰⁵). C10orf54, also known as B7H5 or VISTA (the V domain-‐containing immunoglobulin suppressor of T-cell activation), is one of the B7 family members. Previous studies have suggested that VISTA can negatively regulate T-cell functions in antitumor responses, similar to other checkpoint inhibitors (30). The C10orf54/VISTA expression was lower in 9.1-NP single cells, compared with 9.2-P single cells or 10-P single cells (log₂FC = 1.19 and 1.41, respectively; Fig. 2H). The expression levels of other negative costimulatory molecules were low (Supplementary Fig. S1). Cell surface markers CD5 and CD6 have been suggested to negatively regulate TCR signaling, but the evidence for this is not conclusive. CD5 and CD6 were both downregulated in both 9.2-P and 10-P single cells, compared with 9.1-NP single cells (Fig. 2I and J).

Several transcription factors were differentially expressed in this analysis, including Zinc finger protein 683 (ZNF683), sex-determining region Y-box 4 (SOX4), CCAAT/enhancer binding protein δ (CEBPD), Kruppel like factor 2 (KLF2), and eomesodermin (EOMES; Fig. 3A–E; Supplementary Table S2). ZNF683, also known as HOBIT (homolog of BLIMP1 in T cells), has been suggested to serve as a key transcription factor expressed in tissue-resident memory T cells (31). The expression of ZNF683 was upregulated in 9.2-P single cells, compared with 9.1-NP single cells (log₂FC = 2.71, FDR = 1.92 × 10⁻¹¹⁹; Fig. 3A). Similarly, the ZNF683 expression was upregulated in 10-P single cells, compared with 9.1-NP single cells (log₂FC = 3.97, FDR = 0). These results suggested that ZNF683 might also play an important role in T-cell persistence. Sex-determining region Y-box 4 has been shown to play a key role for CXCL13 production in human CD4⁺ T cells under inflammatory conditions (32). CEBPD has been known to work synergistically with CEBPB to control the induction of proinflammatory cytokines IL6 and TNFα in macrophages (33). However, its role in T cells is largely unclear. Both SOX4 and CEBPD were downregulated in 9.2-P single cells, compared with 9.1-NP single cells (log₂FC = −1.89 and −1.79, respectively; Fig. 3B and C). Both SOX4 and CEBPD were also significantly downregulated in 10-P single cells, compared with 9.1-NP single cells (log₂FC = −1.74 and −2.48, respectively).

It has been reported that the transcription factor Klf2 can promote the expression of S1pr1 and restrain the differentiation of CD4⁺ T follicular helper cells in mice (34). The expression of KLF2 was higher in 9.2-P single cells, compared with 9.1-NP single cells (log₂FC = 1.15, FDR = 3.94 × 10⁻⁶⁹; Fig. 3D). The KLF2 expression in 10-P single cells was also higher than 9.1-NP single cells (log₂FC = 1.68, FDR = 7.62 × 10⁻²⁵⁰). Lastly, the transcription factor EOMES has been identified as one of the key markers among exhausted T-cell populations (35). Significantly less EOMES expression was observed in 9.2-P single cells, compared with the 9.1-NP single cells (log₂FC = -0.89, FDR = 3.59 × 10⁻⁹⁴; Fig. 3E). The EOMES expression was also significantly downregulated in 10-P single cells, compared with the 9.1-NP single cells (log₂FC = −1.57, FDR = 0).

To test whether we could observe the same gene expression patterns using a different approach, we utilized another single-cell system developed by Fluidigm. This system uses a different molecular approach for the single-cell transcriptome analysis. Because Fluidigm C1 is a low-throughput system, TILs were first enriched using the binding of KRAS-9mer or KRAS-10mer tetramer and then subjected to the Fluidigm C1 single-cell system. The full-length, whole-transcriptome sequences from each single cell were obtained, including TCR-CDR3 information for clonotype identification (24). We obtained a total of 123 good-quality single cells from 9.1-NP, 42 single cells from 9.2-P, and 72 single cells from 10-P. As illustrated in Figs. 2 and 3, significant differences (P < 0.05) were observed in 11 of the 13 genes when comparing 9.1-NP versus 9.2-P, except IL7R and KLF2, possibly due to the low number of cells obtained from 9.2-P single cells from this low-throughput approach (Supplementary Fig. S2). Significant differences were observed in all 13 genes comparing 9.1-NP versus 10-P (Supplementary Fig. S2).

To further observe the gene expression profiles of T cells after ACT, the patient's PBLs were stained with KRAS-9mer tetramer or KRAS-10mer tetramer in the attempt to enrich the 9.1-NP, 9.2-P, or 10-P populations. The nonpersistent clonotype 9-1-NP was detectable on day 12, but not detectable on day 40 after ACT. Using the methods previously described (Fig. 1A), single-cell samples were prepared from the enriched populations by a 10X Genomics single-cell instrument (Fig. 4A). Using the peripheral blood sample from day 12 after ACT, we obtained a total of 7,365, 525, and 3,487 good-quality single cells associated with 9.1-NP, 9.2-P, and 10-P clonotypes, respectively. From the peripheral blood sample on day 40 after ACT, we acquired a total of 7,919 and 5,806 good-quality single cells associated with 9.2-P and 10-P clonotypes, respectively. In this analysis, TILs were plotted together with PBLs to observe the changes of gene expression following ACT (Supplementary Table S5).

As previously demonstrated, the 9.2-P and 10-P TILs expressed minimum amounts of GZMK, compared with 9.1-NP TILs. On day 12 after cell transfer, significant expression of GZMK was observed in the peripheral blood 9.2-P single cells, although the expression was slightly lower than 9.1-NP single cells in PBLs (log₂FC = −0.16, FDR = 9.69 × 10⁻³; Fig. 4B). On day 40, 9.2-P single cells continued to express a high GZMK. Similar observations were found in 10-P single cells, compared with 9.1-NP single cells in the PBLs on day 12 (log₂FC = −1.70, FDR = 0). KLRB1 had high expression in 9.2-P TILs and 10-P TILs. After ACT, 9.2-P T cells in PBLs continued to express high KLRB1 on day 12 and day 40 (Fig. 4C). However, 10-P T cells lost the expression of KLRB1 on day 12 after ACT, and 10-P T cells continued to express low KLRB1 on day 40 (Fig. 4C). Similar to the observation in TILs, ITGB1 expression was higher in both 9.2-P and 10-P, compared with 9.1-NP on day 12 (log₂FC = 1.12 and 0.96, respectively). For both 9.2-P and 10-P, the ITGB1 expression was slightly higher on day 40, compared with day 12 (log₂FC = 0.21 and 0.31, respectively; Fig. 4D). As one of the key regulators for T-cell homeostasis (20), IL7R expression remained higher in 9.2-P single cells, compared with 9.1-NP single cells on day 12 (log₂FC = 0.62, FDR = 2.62 × 10⁻¹²). The 9.2-P T cells significantly upregulated IL7R on day 40, compared with 9.2-P T cells on day 12 (log₂FC = 1.66, FDR = 6.21 × 10⁻¹²⁷; Fig. 4E). Similar results were observed in peripheral blood 10-P T cells, compared with 9.1-NP T cells (log₂FC = 0.86, FDR = 6.4 × 10⁻¹²²). These results suggested that 9.2P and 10-P were more competitive in terms of obtaining limited homeostatic cytokine IL7 in vivo. The S1PR1 expression was not significantly different between peripheral blood 9.1-NP and 9.2-P T cells on day 12, but the S1PR1 expression was significantly higher in 9.2-P T cells on day 40, compared with 9.2-P T cells on day 12 after ACT (log₂FC = 0.95, FDR = 3.43 × 10⁻¹⁵; Fig. 4F). On the other hand, the S1PR1 expression was significantly higher in peripheral blood 10-P T cells, compared with 9.1-NP T cells on day 12 (log₂FC = 0.42, FDR = 8.06 × 10⁻²³). The S1PR1 expression was significantly higher in 10-P T cells on day 40, compared with 10-P T cells on day 12 after ACT (log₂FC = 0.83, FDR = 4.55 × 10⁻⁶⁵). The expression of C10orf54 (VISTA) was slightly higher in 9.2-P T cells, compared with 9.1-NP T cells on day 12 after ACT (log₂FC = 0.27, FDR = 1.48 × 10⁻²; Fig. 4G). The expression of other negative costimulatory molecules was low (Supplementary Fig. S1). Both CD5 and CD6 expression was significantly lower in 9.2-P T cells, compared with 9.1-NP T cells on day 12 (log₂FC = −1.23 and −1.24, respectively; Fig. 4H and I).

Similar to the observation in TILs, the transcription factor ZNF683 also had high expression in 9.2-P T cells on day 12, compared with 9.1-NP T cells on day 12 (log₂FC = 1.83, FDR = 2.14 × 10⁻⁷²; Fig. 5A). Similarly, 10-P T cells expressed higher ZNF683, compared with 9.1-NP T cells on day 12 (log₂FC = 2.51, FDR = 0). Both 9.2-P and 10-P T cells on day 40 expressed higher ZNF683, compared with 9.2-P and 10-P T cells on day 12 (log₂FC = 0.54 and 0.61, respectively). The expression of SOX4 and CEBPD was minimum among the peripheral blood 9.1-NP, 9.2-P, and 10-P T cells (Fig. 5B and C). The expression of transcription factor KLF2 was higher in 9.2-P T cells, compared with 9.1-NP T cells in the peripheral blood on day 12 (log₂FC = 0.34, FDR = 4.97 × 10⁻⁴; Fig. 5D). Similarly, 10-P T cells expressed higher KLF2 than 9.1-NP T cells on day 12 (log₂FC = 0.81, FDR = 2.17 × 10⁻¹⁰⁹). Lower levels of EOMES expression were observed in 9.2-P on day 12, compared with 9.1-NP (log₂FC = −0.49, FDR = 7.38 × 10⁻⁶). However, the expression of EOMES in 9.2-P increased on day 40, compared with day 12 (log₂FC = 0.26, FDR = 2.13 × 10⁻²; Fig. 5E). Similarly, lower EOMES expression was observed in 10-P on day 12, compared with 9.1-NP (log₂FC = −0.58, FDR = 3.05 × 10⁻⁴³). However, the expression of EOMES in 10-P also increased on day 40, compared to day 12 (log₂FC = 0.18, FDR = 1.3 × 10⁻⁵; Fig. 5E). Taken together, we observed the continuing high expression of IL7R, ITGB1, ZNF683, and KLF2 in peripheral blood 9.2-P and 10-P T cells following ACT. On the other hand, 9.1-NP T cells had relatively lower expression of these genes in PBLs, compared with 9.2-P or 10-P T cells.

The persistence of T cells can be influenced by the host condition, including the levels of homeostatic cytokines, endogenous T-cell reconstitution, the tumor microenvironment, as well as the binding affinity between TCR and antigen/HLA (17–20). With these caveats in mind, we attempted to test whether or not other infused TILs expressed a set of genes associated with persistence. From the comparisons of 9.1-NP versus 9.2-P and 9.1-NP versus 10-P, 38 shared genes were significantly upregulated, and 41 shared genes were significantly downregulated (FC > ±1.5, FDR <0.05), which were defined as the Persistence_Up gene set and Persistence_Down gene set, respectively (Supplementary Table S1). To perform the GSEA analysis (27), we selected four autologous TILs, 4007, 4069, 4071, and 4081, obtained from patients with metastatic gastrointestinal cancers (6, 16). These autologous TILs were used to treat these 4 patients within the 5-month window prior to the treatment of patient-4095. The same as the majority of infused TILs in this clinical trial, these four TILs contained predominantly (32.3%–91.9%) neoantigen-reactive CD8⁺ T cells but persisted poorly 1 month after ACT (Fig. 6A; Supplementary Table S6; refs. 4, 16). Following the same process as the previous experiments, the single-cell transcriptome analysis was performed to study gene expression profiles linked to neoantigen-reactive clonotypes. These single cells associated with neoantigen-reactive clonotypes were first compared with the persistent clonotypes (9.2-P and 10-P), and the differentially expressed genes were subjected to the GSEA analysis. As shown in Fig. 6B–E, strong differences were observed in all eight GSEA tests for Persistence_Up and Persistence_Down gene sets. Therefore, these TILs obtained from 4 patients did not express the gene sets associated with the persistence, suggesting that these TILs might lack the “fitness” to persist in vivo.

In this study, we discovered a gene expression profile associated with T-cell persistence in a patient treated for KRAS(G12D)-expressing tumors. It has been known that IL7 is one of the key cytokines involved in T-cell homeostasis (20). It is thus likely that T cells with high IL7R could have a survival advantage over the low IL7R cells. The 9.2-P and 10-P persistent T cells expressed high IL7R prior and after ACT, compared with the 9.1-NP nonpersistent T cells. Similarly, a previous study suggested that in vitro–expanded MART1- or gp100-specific T-cell clones that expressed the IL7R were more likely to persist longer after adoptive transfer to patients (36). Next, we studied S1PR1 and KLF2, the key regulators to promote T-cell egress from the thymus or lymphoid organs into the blood (29). A previous mouse study showed that Klf2 could induce the expression of S1pr1 and limit the differentiation of CD4⁺ T follicular helper cells (34). Another study suggests that CD8⁺ resident memory T cells express low Klf2 and S1pr1 (37). Consistent with these findings, we found that 9.2-P and 10-P persistent TILs expressed high KLF2 and S1PR1, whereas 9.1-NP nonpersistent TILs expressed low KLF2 and S1PR1. After ACT, 9.2-P and 10-P T cells continued to express KLF2 and S1PR1. Because of the limitation in the human study, it is unclear whether or not T cells can downregulate KLF2 and S1PR1, and subsequently enter the tumor.

One of the most significant changes observed in this study was the expression of ZNF683 in T cells that were infused and persisted. Znf683 (Hobit) was identified as a homolog of Blimp1 and a key transcription factor for CD8⁺ tissue-resident memory T cells in mice (31). Zundler and colleagues demonstrate that Znf683 and Blimp-1 double-knockout mice are protected in CD4⁺ tissue-resident, memory T-cell–mediated experimental colitis models (38). Conversely, a human study observes that ZNF683 (HOBIT) has high expression in viral-specific CD8⁺ effector T cells in peripheral blood, and it was proposed that ZNF683 might be a key transcription factor for long-lived effector-type T cells based on this observation (39, 40). However, no further evidence has been published to support this hypothesis. In this study, we found that the ZNF683 expression in the infused CD8⁺ TILs was associated with their persistence in vivo. ZNF683 was highly expressed in persistent 9.2-P and 10-P T cells in peripheral blood 12 and 40 days after cell infusion. Therefore, ZNF683 might serve as a key transcription factor for T-cell persistence, potentially working together with other transcription factor(s), such as KLF2. Conversely, it is possible that not only a few key transcription factors, but the entire gene expression profile is responsible for the fitness of TILs.

It has been well-established that “younger” T cells, including naïve T cells, stem cell memory T cells, and central memory T cells, can persist better in vivo and induce stronger antitumor responses, compared with the “older,” terminally differentiated, effector memory T cells (41). However, it is puzzling that terminally differentiated T cells still can induce durable tumor regressions following TIL-ACT or checkpoint immunotherapy. To explain this, it has been proposed that mouse exhausted T cells can be divided into two subsets: (i) “Early” exhaustion subset, defined as Eomes^loPD-1^intTcf1⁺ T cells; (ii) “Terminal” exhaustion subset, defined as Eomes^hiPD-1^hiTbet^loTcf1⁻ T cells (42). The Eomes^hi subset has poor proliferative potential but retains partial cytotoxicity, whereas Eomes^lo subset retains moderate proliferative capacity but has limited cytotoxicity (35). The Eomes^lo early exhaustion subset can convert to the Eomes^hi terminal exhaustion subset after antigen stimulation and proliferation, and this process is not reversible (35). More importantly, PD-1 blockade can reverse and reinvigorate the exhausted phenotype of Eomes^lo subset, but not the Eomes^hi subset (43). However, the majority of evidence supporting this hypothesis was obtained from the LCMV chronic infection model, and there is little evidence in human cancer studies to support this hypothesis directly (42). In the present study, the persistent 9.2-P and 10-P TILs had lower EOMES expression than the nonpersistent 9.1-NP TILs. Similarly, 9.2-P and 10-P TILs also expressed low GZMK, known to be associated with cytotoxicity. Following ACT, the persistent 9.2-P and 10-P T cells in the peripheral blood had lower EOMES than the nonpersistent 9.1-NP T cells, and on day 40, 9.2-P and 10-P T cells upregulated EOMES, compared with 9.2-P and 10-P T cells on day 12 after ACT. Therefore, our observation in this study was in line with this hypothesis. However, further evidence is required, such as tracking the T-cell phenotypes in the peripheral blood and tumors through the course of immunotherapy, which can be difficult in human studies. It has also been hypothesized that Tcf1⁺ “memory-like” exhausted T cells are proliferation-competent and responsible for the antivirus or likely antitumor responses after PD-1 blockade immunotherapy (44–46). However, the expression of TCF1 (TCF7) was minimum in our single-cell analysis, and no significant difference of TCF1 (TCF7) expression was observed between 9.1-NP versus 9.2-P or 9.1-NP versus 10-P.

It remains an open question whether or not the persistence of T cells in the peripheral blood is required for tumor regression and clinical responses or is merely a byproduct of vigorous immune responses. Our previous ACT studies using melanoma TILs demonstrate a significant correlation between the persistence of TIL clonotypes and clinical responses (47). Conversely, in another ACT trial using peripheral blood–derived, MART-1–specific cytotoxic T cells, followed by a standard course of anti–CTLA-4, there was no relationship between the persistence of the infused T cells and the tumor responses (48). In ACT studies using CD19 chimeric antigen receptor–modified T cells for patients with B-cell leukemias or lymphomas, several studies show the correlation between persistence and responses (17, 18, 49–52). However, no significant correlation was found in one of the studies (53). In one ACT study using NY-ESO-1 TCR-transduced T cells for patients with synovial sarcoma, a correlation between persistence and clinical responses was observed (54). However, the correlation was not observed in another study for patients with melanoma or synovial sarcoma (55). Taken together, the outcome of ACT is determined by many factors, which may or may not include the persistence of T cells. More direct evidence is needed to further understand these complicated observations.

In summary, we have identified two different TIL populations with distinctive gene expression profiles, which were associated with the persistence of TILs in a patient. Although we could obtain significant knowledge in a relatively defined environment in a single patient, the observation might not be generalized because of the diversity of patients. We plan to do this single-cell analysis for a large number of patients as more patients are accrued in this clinical trial. Given these caveats, the information obtained in this study has opened the possibility to improve TIL therapies against epithelial cancers by enriching TIL populations based on one or the combinations of the cell surface markers, such as KLRB1, ITGB1, IL7R, S1PR1, or C10orf54. Manipulating the expression of key transcription factors, such as ZNF683 and KLF2, may lead to better T-cell products for immunotherapy. These hypotheses need to be further tested in human studies.

E. Tran is a consultant/advisory board member for PACT Pharma, Genocea Biosciences, and Tailored Therapeutics. No potential conflicts of interest were disclosed by the other authors.

Conception and design: Y.-C. Lu, Z. Zheng, E. Tran, P.F. Robbins, S.A. Rosenberg

Development of methodology: Y.-C. Lu, Z. Zheng

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): Y.-C. Lu, Z. Zheng, E. Tran

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): Y.-C. Lu, L. Jia, Z. Zheng

Writing, review, and/or revision of the manuscript: Y.-C. Lu, Z. Zheng, E. Tran, P.F. Robbins, S.A. Rosenberg

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): Z. Zheng

Study supervision: Y.-C. Lu, S.A. Rosenberg

The authors would like to thank Nicholas Restifo, Steven Shema, Valery Bliskovsky, Michael Kelly, and Jinguo Chen for suggestions and technical support. The authors also thank Doug Joubert, NIH Library Writing Center, for article editing assistance. Next-generation sequencing was conducted at the CCR Genomics Core at the NCI. This work utilized the computational resources of the NIH HPC Biowulf cluster (http://hpc.nih.gov).

This work was supported by the Intramural Research Program of NCI.

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