Single-Cell Sequencing Reveals Trajectory of Tumor-Infiltrating Lymphocyte States in Pancreatic Cancer

Pancreatic ductal adenocarcinoma (PDAC) accounts for ∼85% of all diagnosed pancreatic cancers and, with a 5-year survival rate of only 9%, is the fourth-highest cause of cancer-related death in the United States (1, 2). Chemotherapy and surgery remain the chief therapy options, providing curative potential. However, only 10% to 20% of patients are eligible for surgery, and 80% of those will experience relapse (3). The unmet clinical need of this patient population has spurred efforts to expand therapeutic options.

Immunotherapy represents an exciting new therapeutic strategy, in part due to the success of immune-checkpoint blockade in melanoma, lung cancers, and renal cell carcinoma (4–7). Cellular immunotherapies also have found success. These include the adoptive transfer of T cells engineered to express chimeric antigen receptors (CAR-T) in several hematologic cancers (8, 9) and tumor-infiltrating lymphocytes (TIL) in metastatic melanoma (10, 11). However, immune-checkpoint blockade and CAR-T have provided limited benefit in PDAC so far (12–14).

Due to the variable response to immunotherapy across patients with cancer, it is important to understand the underlying immune landscape of the tumor microenvironment (TME) relevant to clinical response. Although the presence of TIL has been correlated with better prognosis in a variety of cancers, including PDAC (15–20), suppression by regulatory T cells (Treg) and attenuation of cytotoxic function by entering a dysfunctional state have been identified as mechanisms of tumor immune escape (21–23). Conversely, TIL populations such as stem cell–like memory T cells (T_SCM) and tissue-resident memory-like T cells (T_RM) have been identified as important for response or improved survival (24–26). These populations are often defined by a complex combination of surface markers and transcriptional states, and no consensus signature relevant for outcome has emerged yet. T-cell repertoire analysis can also provide insight into the antitumor response in the TME, given that expansion of T-cell clones at the tumor site can indicate an immune response and high-frequency clones have been found to be the tumor-reactive clones (27).

Single-cell RNA sequencing (scRNA-seq) has emerged as a powerful tool to interrogate the heterogeneity of complex cell populations such as TIL by combining the transcriptomic profile and T-cell receptor (TCR) repertoire in the same cells. Indeed, several scRNA-seq studies have revealed TIL diversity across multiple solid tumor types, identifying various states of T-cell dysfunction as well as potentially novel TIL populations associated with response and prognosis (24, 28–35). The lack of response of PDAC to immunotherapy thus far, despite the presence of tumor-reactive TIL (36–38), and the dire need for improved therapy options make PDAC an ideal candidate for deeper analysis using scRNA-seq.

Previous work by this group has shown the feasibility of ex vivo expansion of tumor-reactive PDAC TIL for adoptive cell therapy (ACT; ref. 38). Therefore, the purpose of this study was to characterize the PDAC TIL landscape with an eye toward informing TIL use in cellular immunotherapy of PDAC. To perform this study, freshly resected tumor and uninvolved tissue from patients with PDAC were analyzed using combined scRNA-seq and TCR sequencing to interrogate the heterogeneous TIL populations present to better understand the tumor-infiltrating T cells. Additionally, cultured PDAC TIL from a subset of tumor samples were analyzed using similar methods to investigate the contribution of each initial TIL population to the composition of the final TIL product with its therapeutic potential in mind. Overall, seven CD8⁺ TIL cell states and five CD4⁺ TIL cell states in the fresh tissue were elucidated along with their inferred relationship based on transcriptomic and TCR data. Furthermore, transcriptomic and TCR analysis of the ex vivo expanded TIL product showed preferential expansion of select effector T-cell states originally identified in situ.

To delineate TIL substates in PDAC, we performed scRNA-seq on seven primary, human PDAC samples and one uninvolved pancreas sample (hereafter referred to as “uninvolved”), focusing on the CD3⁺ TIL (Fig. 1A, MDA1 data set). To further increase the data sample size, CD3⁺ TIL scRNA-seq data from another cohort at the institution (MDA2 data set; 26 primary PDAC and 10 uninvolved) and from a previous study by Peng and colleagues (PUMCH data set; 24 PDAC and 11 normal pancreas) were combined (39). In total, 39,694 T cells were analyzed from 57 PDAC samples and 22 uninvolved/normal samples. Clustering of T cells from all samples identified 13 populations overall, with 5 CD4⁺ TIL cell states (CD4-FOXP3, CD4-CXCR4, CD4-CCR7, CD4-CXCL13, and CD4-MX1), 7 CD8⁺ TIL cell states (CD8-GZMK, CD8-CXCR6/IL7R, CD8-ZNF683, CD8-GZMB/PRF1, CD8-CXCL13, CD8-CCR7/IL7R, and CD8-MX1), and a cycling subpopulation (Fig. 1B). Cohort-specific and patient-specific batch effects were limited as evidenced by the even mixing between the cohorts and the patients as shown in uniform manifold approximation and projection (UMAP; Supplementary Fig. S1A and S1B, respectively). These cell states were annotated after the identification of defining T-cell genes that are expressed within that cell state (Supplementary Fig. S1C). The presence of all CD4⁺ and CD8⁺ TIL cell states, with the exception of CD8-CXCR6/IL7R, was validated in formalin-fixed paraffin-embedded primary PDAC tissues using RNA in situ hybridization (RNAScope), which enables visualization of single RNA molecules in individual cells (Supplementary Fig. S1D).

To understand the biological function of each T-cell cluster, the unbiased top differentially expressed genes were investigated as well as the expression of specific known transcription factors. This analysis revealed cluster-specific expression of key genes, indicating cell-state programming differences between the clusters (Fig. 1C; Supplementary Fig. S2A). Although the TIL cell states were primarily composed of αβ T cells, marked by expression of the TCR alpha-chain gene TRAC, γδ T cells (expressing the TCR gamma-chain gene TRDC) were found to be a small component of the CD8-GZMB/PRF1 cell state (Supplementary Fig. S2B). In addition to the transcriptional states, there was a distinct cycling T-cell state with the expression of MKI67 and other cell cycle–related genes (Fig. 1B; Supplementary Fig. S2C).

Regarding the potential functional identity of the TIL states, the CD4-FOXP3 cluster was classified as Tregs due to the expression of FOXP3, IL2RA (CD25), CTLA4 (Supplementary Fig. S2A), and the Treg-associated transcription factors IKZF2 (Helios) and IKZF4 (Ikaros; Fig. 1C). An activated CD4⁺ helper T-cell population, CD4-CXCR4, was also identified based on CD69, GZMA, and CCL5 expression (Supplementary Fig. S2A). However, low expression of the transcription factors TBX21 (TBET) and GATA3 made the classification as a Type 1 or Type 2 helper T-cell unclear (Fig. 1C). Within the CD8⁺ TIL population, the CD8-ZNF683 cluster appeared to be a T_RM-like population as suggested by expression of the T_RM transcription factor ZNF683 and the _TRM-associated genes XCL1 and GNLY, and low expression of KLF2, as previously reported (Fig. 1C; Supplementary Fig. S2A; ref. 29). More­over, a cytotoxic CD8⁺ cell state was identified (CD8-GZMB/PRF1) by high expression of the well-characterized effector cell markers GZMB and PRF1 (Supplementary Fig. S2A) and the transcription factor TBX21 (Fig. 1C).

Beyond these unique cell states, the single-cell data showed CD8⁺ and CD4⁺ clusters that shared expression of several defining genes. One of these sets, CD4-CXCL13 and CD8-CXCL13, showed expression of the transcription factor RBPJ and chemokine CXCL13, which are both identified as part of a reputed dysfunctional T-cell phenotype (Fig. 1C; Supplementary Fig. S1A; refs. 35, 40). These two cell states also showed expression of the well-known immune-checkpoint marker genes CTLA4, PDCD1 (PD1), HAVCR2 (TIM3), TNFRSF18 (GITR), LAG3, and TIGIT (Supplementary Fig. S2D). Another set of cell states, CD4-CCR7 and CD8-CCR7/IL7R, shared expression of transcription factors LEF1 and TCF7 (Fig. 1C) and surface markers CCR7, IL7R, and SELL (CD62L; Supplementary Fig. S2A), consistent with a central-memory T-cell phenotype (T_CM).

Across the cohorts, T-cell states were detected at varying frequencies (Fig. 1D; Supplementary Fig. S2E). For example, in tumor samples with ≥50 T cells detected, the CD8-CXCL13 population was observed in 43/53 PDAC samples, but only 7 patients had frequencies of this cluster ≥5% of their total TIL (Fig. 1D). This contrasts with the CD8-GZMK state that was identified in all tumors, and 50/53 patients had frequencies ≥5% of their total T cells. Additionally, uninvolved/normal samples predominantly clustered away from PDAC samples, suggesting global differences in their T-cell cluster composition (Fig. 1D). The CD8-CXCR6/IL7R cell state was significantly enriched in the uninvolved/normal tissue compared with PDAC samples (P < 0.001), whereas the reverse was found in the CD8-ZNF683 cell state (P < 0.001, permutation test, 1,000 permutations). Furthermore, to establish if the cell states are tumor tissue–specific or reflecting the organ from which the tumor originated, the cell-state composition was compared in a subset of 10 patients who had matched tumor and uninvolved tissue samples (Fig. 1E). In these patients, it was found that the CD4-FOXP3 and CD4-CXCR4 cell states were enriched in the tumor over the uninvolved tissue (P = 0.021 and P = 0.035, respectively, two-sided, paired Wilcox test), whereas the CD8-CXCR6/IL7R and CD8-GZMB/PRF1 cell states were enriched in the uninvolved versus the tumor (P = 0.041 and P = 0.011, respectively, two-sided, paired Wilcox test). However, these changes were not statistically significant after P-value correction for multiple testing.

To further verify if the cell states were tumor-specific or overlapping with cell states found in the periphery, their transcriptional profiles were compared with that of peripheral blood T cells (Supplementary Fig. S3). Publicly available sequencing data from healthy donors’ and PDAC patients’ peripheral blood mononuclear cells (PBMC) were combined to make the comparison (Supplementary Fig. S3A; ref. 41). Shared nearest-neighbors clustering revealed five CD4⁺ T-cell states and four CD8⁺ T-cell states as defined by their differentially expressed genes (Supplementary Fig. S3B; Supplementary Table S1). Comparing the top 20 differentially expressed genes (by P value) of the tissue T-cell states and the PBMC cell states showed top-marker similarity only between three states, the CD8⁺ T_CM and the CD4⁺ T_CM cell states (both defined by CCR7) and the CD8⁺ cytotoxic cell states (defined by GZMB and PRF1; Supplementary Fig. S3C). As CCR7 can identify less-differentiated T cells and naïve T cells, it is not surprising that the two are found in both the tumor and the blood. Nonetheless, this highlights that the majority of the defined cell states are unique to the tissue.

Because T cells are known for their capacity to transition between different functional and differentiation states, we applied pseudotemporal inference to understand the relationship between the identified cell states. Because CD4⁺ and CD8⁺ T cells are not developmentally related beyond thymic selection, trajectory analysis on these major subsets was performed separately (Fig. 2A and B). Additionally, due to the high complexity of the T-cell development landscape, the lines displayed on the graph are weighted by the density of cells around each node. This density is proportional to the amount of data supporting the trajectory at that point, and therefore the confidence of that location in the graph.

For the CD8⁺ TIL, the CD8-CCR7/IL7R T_CM state and the CD8-CXCR6/IL7R state were inferred to be related in the trajectory. However, that portion of the graph had low cellular density support, suggesting the other CD8⁺ TIL states did not entirely arise from the T_CM population detected in the tumor (Fig. 2A). In contrast, the transition of TIL through the other cell states was strongly supported by the trajectory analysis, based on the closeness on the path, and the thickness and length of their branches. For example, the CD8-CXCR6/IL7R and CD8-ZNF683 cell states occupied the same long branch, indicating a high degree of similarity (Fig. 2A). This contrasts with the latter part of the pseudotime trajectory beyond CD8-ZNF683, where the shorter branch lengths and multiple forks indicate that the CD8⁺ TIL can potentially enter distinct states. The trajectory also suggests that CD8-GZMK may be an intermediate cell state between the CD8-GZMB/PRF1 cytotoxic cell state and the CD8-CXCL13 potentially dysfunctional state. The relationship of the CD8-MX1 IFN-I response cell state to the others appeared unclear from these trajectories, as they did not occur on one of the branches of the trajectory, suggesting they may be a separate state. Additionally, the contributions of CD8⁺ TIL from tumor and uninvolved samples overlapped, with the exception of CD8-CXCR6/IL7R (Fig. 2C). The inference of similarity and proximal state transition from the pseudotime analysis may suggest that a TIL population found outside the tumor (CD8-CXCR6/IL7R) could be contributing to a state inside the tumor (CD8-ZNF683).

Similar trajectory analysis was performed on the CD4⁺ TIL (Fig. 2B). These data showed that a CD4-CXCR4 cell state was adjacent to the T_CM cluster, with the suggested transition having a much higher confidence (higher density of cells) for CD4⁺ TIL in contrast to what the data suggested for CD8⁺ TIL. The trajectory analysis suggested that a small subset of effector cells may be related to the CD4-CXCL13 dysfunctional cluster. However, this portion of the trajectory had a relatively low number of cells whose transcriptional data supported the transition (as displayed by the low cell density surrounding the purported transition on the graph) and therefore represents an unlikely cell state transition. Likewise, the CD4⁺ T-cell trajectory also tenuously projected a transition from an effector state to a final Treg phenotype. However, this transition is again weakly supported based on the distance between the CD4-FOXP3 cell state and the other CD4⁺ cell states as well as the extremely sparse cellular density along the trajectory. This may indicate these were natural Tregs that infiltrated the tumor as opposed to induced Tregs that transitioned from another state in situ. Similar to the CD8-MX1 IFN-I cluster, the relationship between the CD4-MX1 IFN-I response population and the other clusters was unclear, likely due to its small sample size. Unlike the CD8⁺ TIL, there was homogeneity of the CD4⁺ TIL populations between tumor and the uninvolved normal samples (Fig. 2D).

To determine the degree to which T-cell states co-occur in the same PDAC sample, we performed Spearman correlations on cell-state frequencies of all tumor samples in a patient-by-patient manner, similar to an analysis performed by Bassez and colleagues (42). A significantly positive correlation between two cell states indicates that if a patient has a higher fraction of cluster “A,” they also tended to have a higher fraction of cluster “B” (Fig. 2E). This analysis showed that the frequencies of the highly related CD8-MX1 and CD4-MX1 IFN-I response cell states (r = 0.70, P = 8.17 × 10⁻⁸) were significantly correlated with each other. Although this cell state for both CD4 and CD8 shared many genes expressed by other cell states (Supplementary Fig. S2A), only CD8-MX1 significantly, positively correlated with CD4-FOXP3 (r = 0.45, P = 0.013), validating that these MX1 states stand on their own, but also that a TME rich with interferon, which influences these cell states, might be responsible for their significant correlation (43, 44). Finally, the CD8-CXCL13 and CD4-CXCL13 cell states (r = 0.53, P = 6.71 × 10⁻⁴) were significantly correlated with each other, plus a significant correlation between the CD8-CXCL13 cell state and the cycling population of cells (r = 0.40, P = 0.028) was also observed, indicating that the CXCL13 cell state might be a proliferative population of TIL despite potentially being dysfunctional.

Inverse relationships were also observed for several T-cell states (Fig. 2E). A significant negative correlation was found for the CD8/CD4-CXCL13 populations and the CD4-CXCR4 effector population (r = −0.37, P = 0.041 and r = −0.39, P = 0.031, respectively), suggesting that a TME with dysfunctional T cells has a lesser content of some but not all effector populations. On the other hand, when there was a richer content of the effector population CD8-GZMB/PRF1, it was associated with a scarcity of the CD8 T_CM cell state (CD8-CCR7/IL7R; r = −0.36, P = 0.04). Another effector cell state, previously shown as the most prevalent effector population (Fig.1D), CD8-GZMK, had a significant negative correlation with the T_RM-associated CD8-ZNF683 cluster (r = −0.37, P = 0.04; Fig. 2E). Generally, dysfunctional (CXCL13), interferon-response (MX-1), and suppressive (FOXP3) cell states clustered away from effector (ZNF683, GZMK, and GZMB/PRF1) or less-differentiated cell states (CCR7/IL7R; Fig. 2E).

To gain further insight into the potential function of T-cell states as well as the relationship between them, single-cell TCR sequencing and transcriptomic profiling on the same cells were performed on a subset of seven samples (the MDA1 cohort). As shown in Supplementary Fig. S1A, the MDA1 cohort is comparable to the combined data set given the degree of overlap between data sets. The similarity was further verified by quantifying the overlap of genes between the MDA1 data set and the larger data set (Fig. 3A). Each cell state had at least 70% overlap of the top genes (P ≤ 0.05 and fold change ≥0.2), except for the CD4-CXCL13 cell state, which was composed of only a few cells in the MDA1 cohort. Furthermore, unsupervised clustering of the MDA1 TIL showed that the cell state markers in this subset were similar to the larger data set (Supplementary Fig. S4A and S4B). For example, the top markers of the CD8⁺ cell states included GZMK, IL7R, ZNF683, PRF1, CXCL13, and MX1, whereas top markers in the CD4⁺ TIL cell states included FOXP3, CD40LG, SELL, CXCL13, and MX1. UMAP embedding of the smaller MDA1 cohort suggested the robustness of the clustering, with the original cell state labels retaining colocalization in the UMAP for both CD8⁺ and CD4⁺ TIL (Fig. 3B and C). To assess the possibility of sampling bias in the MDA1 sample set, the contribution of each sample to the different cell states was analyzed (Fig. 3D and E). Apart from the CD8-CXCL13 cell state deriving predominantly from one sample (MDA1_T05), the cell states were not dominated by a single patient (Fig. 3D). Together, these data suggest that conclusions drawn from the smaller MDA1 cohort could potentially be extended to the larger data set of patients with PDAC.

Next, TCR sequences and transcriptomic information from the same single cells were analyzed to determine the relationship between cell states, as clonotypes found within multiple cell states (referred to here as cell state sharing) may indicate a transition between those states. Visualization of the combination of the transcriptional cell states with their clonal frequency for CD8⁺ TIL and CD4⁺ TIL showed that most clones were detected at low frequencies, suggesting limited proliferation at the tumor site (Fig. 4A and B). In Fig. 4A and B, each bar of the Circos plots represents one T-cell clone where the height of the gray bar on the outer ring of the Circos plots is proportional to the number of cells bearing this TCR. The middle ring indicates the transcriptomic state(s) in which each TIL clonotype can be found, and the inner ring is colored by the patient to which this TCR belongs. When a clonotype is found in multiple transcriptomic clusters, the bar in the middle ring is comprised of several smaller colored segments. Focusing on the CD8 transcriptomic states in Fig. 4A, we found that the expanded TIL clonotypes (i.e., found in multiple copies in the tumor tissue) in the CD8-GZMK cell state (teal) had a high degree of cell state sharing (Fig. 4A, red trapezoidal segment), which is highlighted in the zoomed-in sections of the circular clonotype plots showing the contrast between expanded clones (Fig. 4C, left) and unexpanded clones (Fig. 4C, right). These observations were further emphasized by the Upset plot showing the amount of cell state sharing and the number of times a specific intersection was detected (Fig. 4D; Supplementary Fig. S5). For the CD8-GZMK cell state, 18/40 TIL clonotypes (with ≥2 cells detected) were also assigned to at least one other cell state, such as the six clonotypes shared between CD8-GZMK and CD8-ZNF683 (Fig. 4D). Furthermore, the CD8-ZNF683 cell state also exhibited many expanded clonotypes (Fig. 4A, purple bars) and a high number of cell state sharing (Fig. 4D), with 23 of its 41 clonotypes detected in at least one other cell state, particularly the CD8-CXCR6/IL7R state (blue bars) with 7 instances of overlap. Both the CD8-GZMK and CD8-ZNF683 clusters also exhibited sharing with the CD8-GZMB/PRF1 cell state (5 times and 3 times, respectively). The TCR overlap between these cell states shows a relationship and further strengthens the cell state transitions depicted by the pseudotime trajectory analysis in Fig. 2A. In contrast, the clones in the CD8-CXCL13 cell state (pinkish-red) were highly expanded (Fig. 4A) but were detected only once in other cell states, perhaps indicating a termination point on the trajectory (Fig. 4D). In the CD4⁺ TIL, only very limited cluster sharing was identified, pointing to a lack of proliferation in the tumor (Fig. 4B and E; Supplementary Fig. S6). This is particularly notable for the CD4-FOXP3 cell state (red), which was noted in the pseudotime plot (Fig. 2B) to have a very weak connection to the rest of the trajectory. This lack of TCR sharing between this cell state and other CD4 states further supports the idea that these are innate rather than induced Tregs.

Further information can be gleaned from the TCR-sequencing data by using clonal expansion and diversity of TIL clones at the tumor site as a surrogate for functional involvement in the antitumor response. In Fig. 4A, the effector cell states of CD8-GZMK, CD8-ZNF683, CD8-GZMB/PRF1, and CD8-CXCL13 showed the highest levels of clonal expansion, with the potentially dysfunctional CD8-CXCL13 cell state having the highest levels of clonal expansion. Although the CD8-CXCL13 cell state was found in only one patient (MDA1_T05), the observed high clonal expansion matches what has been previously reported by several other studies (33, 40, 45, 46). In contrast, very little CD4⁺ TIL expansion was observed, suggesting their involvement in the response to the tumor may have been suppressed (Fig. 4B). The lack of CD4⁺ TIL expansion reported here is also consistent with previous observations in the literature for other solid tumor types (29, 33, 34, 40, 47).

To quantify the degree of clonal expansion within each cluster, we measured clonal frequency evenness by computing the Gini index for CD4⁺ and CD8⁺ TIL across patients (Fig. 5A, left). The Gini coefficient is a measure of inequality, ranging from 0 to 1, where lower values indicate clonal frequency evenness and higher values signify more inequality in clonal expansion (i.e., expansion of only one or very few clones). With a median Gini coefficient value closer to 1, the CD8⁺ TIL population showed a higher degree of oligoclonality than the CD4⁺ TIL across patients (P = 0.0313, paired t test). This observation was validated using the related clonal evenness metrics DE50 (P = 0.0313) and the Pielou index (P = 0.0313), which both showed CD8⁺ TIL having a lower median value than CD4⁺ TIL (Supplementary Fig. S7A). When we evaluated the diversity of the repertoire using the diversity metric log(Chao1) index, we did not observe a significant difference between the CD4⁺ and CD8⁺ TIL, indicating the repertoires are composed of a similar number of clonotypes (P = 0.156; Fig. 5A, center). However, when we utilized the Shannon diversity index, another diversity metric that takes into account both clonal evenness and clonotype abundance, we observed that overall the CD8⁺ TIL have a lower diversity than the CD4⁺ TIL (P = 0.0313; Fig. 5A, right). This observation was again validated using the related Simpson index (P = 0.0313; Supplementary Fig. S7A, right). Based on previous analyses of T-cell repertoire metrics presented by Chiffelle and colleagues, we concluded that this indicated a greater number of clones comprise the expanded fraction of the CD8⁺ TIL repertoire than is the case for the CD4⁺ TIL repertoire (48).

When breaking down the T-cell repertoire by cell states, we observed a similar pattern with all CD4⁺ TIL states exhibiting a lower clonality than the majority of CD8⁺ TIL states (Fig. 5B, left; Supplementary Fig. S7B, left and center). Contributing the most to the higher clonality for the CD8⁺ TIL were the CD8-GZMK, CD8-GZMB/PRF1, and CD8-CXCL13 cell states. Notably, the median number of clones was lower for these oligoclonal CD8⁺ cell states (Fig. 5B, center). However, this was not necessarily the case with the CD4⁺ TIL as the most oligoclonally expanded CD4⁺ cell state (CD4-FOXP3) had nearly the highest value (i.e., number of clones in its repertoire; Fig. 5B, left and center). Finally, the oligoclonal CD8-GZMB/PRF1 and CD8-CXCL13 cell states had not only lower clonotype richness but the lowest median diversity as indicated by the Shannon and Simpson index (Fig. 5B, right; Supplementary Fig. S7B, right). Taken altogether, this indicates an expansion of select TIL clones in these cell states.

This degree of clonal expansion and the difference between the CD4⁺ and CD8⁺ repertoires in each cell state is further illustrated in Fig. 5C, now depicted by patient. The CD8-GZMK, CD8-ZNF683, CD8-GZMB/PRF1, and CD8-CXCL13 had the greatest proportion of clonally expanded cells within the TME, with at least half of the expanded (≥2) CD8⁺ TIL found ≥5 times, indicating a robust expansion. This contrasts with the CD8-CXCR6/IL7R cell states, where three fourths of the clones were detected only once (Fig. 5C). In concordance with our earlier findings regarding the clonality and diversity of the TCR repertoire, CD4⁺ cell states presented a minority of clonally expanded TIL as evidenced by the majority of the clones found in these cell states being detected only once (Fig. 5C). Furthermore, even within the expanded fraction, most CD4 clones were detected between 2 and 5 times. Taken together, these metrics indicate that although the CD8⁺ TIL population has a similar level of clonotype richness to the CD4⁺ TIL, the CD8⁺ effector populations are more clonal, suggesting the participation of select TIL clones in the immune response to the tumor.

Concurrently, as the tissue samples from the MDA1 cohort were processed for scRNA-seq, a piece from each tissue was used for ex vivo TIL growth. The cultured TIL were then used for transcriptomic and TCR scRNA-seq (n = 40,065; Fig. 6A). The propagated TIL expression data revealed seven major cell states comprising five CD8⁺ states (G1-CD8-CD27, G2-CD8-GNLY, G3-CD8-ZNF683, G4-CD8-CCL3, G4-CD8-MKI67), one CD4⁺ state (G6-CD4-CD40L), and one γδ⁺ T-cell state (G7-γδΤ-TRDC; Fig. 6B). As expected, a large fraction of the cultured TIL were cycling cells. Cells initially classified as cycling were reclassified into the noncycling clusters using a random forest classifier. CD8⁺ TIL states predominated the grown TIL products, which was expected because the culture method was optimized to enhance CD8⁺ TIL proliferation (49). The expression of signature genes suggested several activated CD8⁺ cell states, a proliferating CD8⁺ state, an effector CD4⁺ state, and a γδ⁺ T-cell state (Fig. 6C; Supplementary Fig. S8). Cell states were labeled after key defining genes that were highly expressed (Fig. 6C).

A major goal of the cellular-based cancer immunotherapy field has been to identify and select the best T-cell subtypes to provide the most efficacious ACT. To this end, we attempted to discern the fate of the different fresh TIL cell states identified in Fig. 1 to determine whether the TIL culture process expanded subpopulations of interest. The individual TCRs were used to track TIL from the fresh samples to their corresponding grown TIL product. We found that regardless of their initial transcriptomic cell state in the tissue, the TIL clonotypes were identified in multiple grown cell states with a majority in the G1-CD8-CD27 subtype in the cultured product (Fig. 6D). As opposed to the trajectory inferred between the different cell states in the tissue in our pseudotime analyses in Fig. 2 (where this cell state was not described), there was no apparent directionality from one specific fresh tissue T-cell state to a grown T-cell state. This difference can be attributed in part to the controlled activation environment provided by the in vitro expansion using the TIL 3.0 method, which is in contrast to the activation the TIL received in situ. Analysis of the TCR frequency of the grown TIL paired with their transcriptomic clusters showed a high degree of cluster sharing, where the same TCR sequence was found in multiple transcriptomic clusters (Supplementary Figs. S9 and S10A and S10B).

Despite the lack of connection between any final grown cell state and specific initial PDAC TIL populations, there was an indication that certain TIL populations from the tissue were preferentially grown over others (Fig. 6E and F). The T-cell clones originating from CD8-GZMK, CD8-ZNF683, and CD8-GZMB/PRF1 clusters were preferentially expanded ex vivo, whereas those from CD8-CXCL13 populations decreased in frequency relative to their initial frequency in the tissue (Fig. 6E). The CD4⁺ TIL showed a similar pattern, in which mainly CD4⁺ effector cells expanded while Tregs did not (Fig. 6F). In both CD8⁺ and CD4⁺ TIL, the activated cell states expanded, whereas cell states such as the potentially dysfunctional CD8/CD4-CXCL13 populations had limited expansion. Finally, the frequency of TIL clones in the tumor can be an indication of their involvement in an immune response, where high-frequency clones are assumed to have expanded at the tumor site in response to the disease (27). In considering the use of TIL as a therapeutic modality, the expansion of high-frequency clones from the tumor is highly desired. Thus, we investigated if the culture methodology used expanded high-frequency clones from the tissue. Comparing all TIL clones detected in fresh tissue and grown samples for six pairs, our data showed that the high-frequency clones detected in the fresh tissue were expanded, and they remained at a relatively high frequency in the grown product (Supplementary Fig. S11).

Here, we used scRNA-seq to investigate CD8⁺ and CD4⁺ T cells and defined 13 cell states in PDAC and uninvolved/normal pancreas tissue samples plus seven cell states resulting from ex vivo culture. This work was executed using an initial cohort of seven patients treated at The University of Texas MD Anderson Cancer Center, and then further strengthened by the addition of another cohort of 26 patients from The University of Texas MD Anderson Cancer Center as well as a cohort of 24 patients from a previous study at Peking Union Medical College Hospital (39). Based on the transcriptional data from all three cohorts, and concurrent single-cell TCR sequencing in the initial seven-patient cohort, we present a model for the organizational relationship and differentiation trajectory of PDAC TIL as well as their unrestricted expansion into each component population of the ex vivo culture (Fig. 7). The TIL cell states detected in this study are consistent with those reported in the literature for other cancer types and validate the lack of institutional and sample-size bias in the three data sets used in this study. For example, CD8⁺ TIL populations similar to the CD8-GZMK, CD8-ZNF683, and CD8-CXCL13 found here have been reported in other scRNA-seq studies in melanoma, triple-negative breast cancer, non–small cell lung cancer, colorectal cancer, and liver cancer (24, 29, 31, 33, 34).

In addition to GZMK, the CD8-GZMK cell state in our analysis showed expression of EOMES, which has notably been tied to a T-cell exhaustion phenotype when paired with the expression of immune-checkpoint molecules (50, 51). Moreover, a recent study in pancreatic cancer reported an exhausted CD8⁺ T-cell population expressing GZMK, EOMES, and TIGIT (41). However, the CD8-GZMK population described in our study lacked the expression of checkpoint molecules or other putative exhaustion markers. Its position on the pseudotime trajectory along with its sharing TCRs with multiple other cell states indicates it could be an intermediate cell state. Indeed, the CD8-GZMK population is similar to those reported in several other studies that classified these cells as a transitional “predysfunctional” population that gives rise to dysfunctional T cells (29, 31, 33, 34, 40, 52). The CD8-GZMK cell state found here and in these studies is characterized by intermediate clonality, low/intermediate expression of checkpoint markers, and shared TCRs with their respective cytotoxic (i.e., CD8-GZMB/PRF1) and/or dysfunctional (i.e., CD8-CXCL13) populations.

CD8-ZNF683 may be a T_RM-like cell state based on the expression of the T_RM transcription factor ZNF683 (53). Although this cell state did not express ITGAE (CD103) or have high expression of checkpoint markers, as is commonly reported (54), this cell state was similar to a state reported by Guo and colleagues in non–small cell lung cancer (29). The T_RM population is of particular interest due to its correlation with improved survival outcomes and increased cytotoxic activity in several epithelial cancers (25, 31, 55) and its involvement in response to checkpoint blockade therapy (26, 56). However, it is unclear whether this population plays the same role in PDAC due to the lack of checkpoint marker expression observed in this study.

Another population of interest is the CD8-CXCL13 cell state, which has been reported by other groups to be an exhausted/dysfunctional population based on its high ex­pression of checkpoint molecules. This population is also described as being highly clonally expanded and able to retain proliferative potential at the early stages of dysfunction (33, 40, 45, 46). Although we similarly report in PDAC that the CD8-CXCL13 cell state has high checkpoint expression and is clonally expanded, the population did not appear to retain great proliferative potential as it did not expand well in our TIL 3.0 culture method (simultaneous engagement of the TCR and 4-1BB combined with high-dose IL2). This suggests they could be at a later stage of dysfunction and/or require stimulation beyond that provided by our culture method. Nonetheless, studies in melanoma, ovarian cancer, and lung cancer respectively show that this population contains tumor-reactive TIL (40, 46, 57). Moreover, the dysfunctional appearance of this population could be a consequence of its ability to react against the tumor repeatedly, and therefore there is value in leveraging this TIL population for therapeutic application.

One of the big questions that remain unanswered is what factor(s) of the pancreatic cancer TME might be causing TIL to move to certain cell states or might be responsible for the overall composition of TIL. One factor that could have a role is malignant progression from intraductal papillary mucinous neoplasm (IPMN) precursor lesions to PDAC. It was previously reported that as lesions progress from IPMN to PDAC, there is a reduction in cytotoxic CD8⁺ TIL and CD69⁺ CD4⁺ TIL (58). Likewise, the same trend was observed for the conventional dendritic cell (cDC) subset cDC2, which is a critical mediator of tumor antigen cross-presentation, an important process in the activation of CD8⁺ and CD4⁺ TIL (58). Because of the use of scRNA-seq on CD3⁺ TIL, we were still able to detect both a cytotoxic CD8⁺ TIL cell state (CD8-GZMB/PRF1) and an activated CD4⁺ TIL cell state (CD4-CXCR4) in our study, although with the limitation of not being able to follow their presence through the tumor progression.

Transforming growth factor-beta (TGFβ) is a well-known factor with pleiotropic effects that is found at high levels in PDAC and not only is responsible for PDAC progression at late stage but can affect immune cell phenotype, which in turn can affect the progression (59). TGFβ is responsible for suppressing cytotoxic CD8⁺ and activated CD4⁺ T cells, either directly or indirectly through the suppression of cDC2. TGFβ is also partly responsible for the accumulation of Tregs and myeloid-derived suppressor cells in PDAC, the latter of which are responsible for cancer progression (60). Although the focus of our study was not to characterize the myeloid compartment, we did observe a strong signature for activated Tregs that appeared to be natural Tregs based on their lack of TCR sharing with other cell states and position on the pseudotime trajectory.

TGFβ is also a known factor for the development of T_RM, a population that closely matches with our CD8-ZNF683 cell state (54). In our study, this cell state and the CD8-CXCR6/IL7R cell state appear to be linked based on the pseudotime trajectory, their transcriptome, and TCR overlap. In fact, it has been reported in the literature that a CXCR6⁺ precursor cell can transition to a T_RM phenotype after being activated (61). The occurrence of this transition is supported here in this study by the observation that the CD8-ZNF683 cell state is more clonally expanded than the CD8-CXCR6/IL7R one. Perhaps some of the TIL that comprise the CD8-CXCR6/IL7R cell state were activated, expanded, and transitioned to the CD8-ZNF683 cell state. In contrast, the connection between the CD8-ZNF683 and the CD8-GZMK cell states that is suggested by the pseudotime and the shared TCRs is less clear. The relationship raises the question of whether the CD8-ZNF683 cells could transition to the CD8-GZMK state, which is postulated to be a predysfunctional state, after repeated antigen challenge (62).

Checkpoint blockade immunotherapy is a method that has been successful in altering the TME in some malignancies, but the benefit of this immunotherapy has not yet been demonstrated in patients with PDAC. As a potential alternative approach, previous work by our group focused on developing ACT using PDAC TIL showed it is feasible to generate large amounts of tumor-reactive CD8⁺ TIL from primary and metastatic tumors (38). The TCR scRNA-seq data presented here show that the culture process expands high-frequency clones from potentially important CD8⁺ TIL cell states such as CD8-GZMK, CD8-ZNF683, CD8-CXCL13, and CD8-GZMB/PRF1 while avoiding the detrimental CD4-FOXP3 Tregs. Clonotypes from the CD8-GZMK and CD8-ZNF683 cell states expanded favorably (frequency in situ < frequency ex vivo) during culture, and exploration of these cell states following adoptive transfer may be of interest due to their memory status and proliferative potential. Recently, scRNA-seq analysis by Lu and colleagues found that ZNF683 was associated with improved TIL persistence after adoptive transfer in a patient with metastatic colorectal cancer (63). Furthermore, if CD8-ZNF683 is truly a T_RM-like population, then it may have improved tissue-homing properties as previously reported for this type of T cell (64). In contrast, clonotypes from the CD8-CXCL13 cell state did not expand favorably (frequency in situ > frequency ex vivo). Given the fact that this population is reported to contain tumor-antigen reactive TIL, and therefore has therapeutic value, expansion techniques will need to be modified to take this into account as even the engagement of a powerful costimulatory molecule like 4-1BB used in this study did not optimally expand it.

Ultimately, regardless of the in situ cell states, the grown product consisted of several activated cell states with the largest cell state being the G1–CD8–CD27 cluster. This overall state of activation could be attributed in part to the controlled activation environment provided by the in vitro expansion using the TIL 3.0 method. Finally, our TIL culture process also expanded a population of γδ⁺ T cells that were present at a low frequency in the fresh tissue. Although the culture process was designed to favor CD8⁺ αβ T-cell expansion by using an agonistic 4-1BB antibody, γδ⁺ T cells can also express 4-1BB and respond to this stimulation. Even so, γδ⁺ T cells have antitumor potential, and thus their presence may add to the effectiveness of the TIL product (65).

In conclusion, this study defines the TIL cell states in patients with PDAC, their inferred differentiation trajectory in situ, and their fate during ex vivo TIL culture for application in TIL ACT. Although the observations presented here are consistent with other cancer types, it is important to note that the conclusions presented here are based on a snapshot of a dynamic process. Future studies in the context of immunotherapy-treated patients with PDAC with longitudinal combined transcriptomic and T-cell repertoire analysis will be needed to better understand the TIL cell states presented here and their clinical significance for this difficult cancer type. Finally, we expect that the granular knowledge gained from these data will help inform future T-cell targeted therapeutic strategies as well as manufacturing processes for TIL-based strategies.

After providing written informed consent, 33 patients with primary PDAC underwent surgical resection (Supplementary Table S2). Patients are referred to by their deidentified number. Tissue from surgical resections was used under protocols (PA 15-0176 and PA17-0793) approved by the Institutional Review Board of The University of Texas MD Anderson Cancer Center. This study was conducted in accordance with the Declaration of Helsinki. For the MDA1 and MDA2 data sets, tissue from normal pancreas is referred to as “uninvolved” because it came from patients with PDAC. For the PUMCH data set, tissue from control pancreases is referred to as “normal” because it came from patients without malignant pancreatic tumors.

Fresh tumor samples were cut into 1 to 3 mm³ fragments and disaggregated using a Medimachine (BD Biosciences) to create a single-cell suspension in 1× PBS according to the manufacturer's instructions. After disaggregation, CD3⁺ T cells were isolated for sequencing via magnetic bead separation using the EasySep Release Human CD3 Positive Selection Kit (STEMCELL Technologies) according to the manufacturer's instructions. The enriched cell suspension was collected in 1× PBS containing 2% FBS and 1 mmol/L EDTA, enumerated using an automatic cell counter, and sent for oil immersion on the 10X platform. For cultured TIL, previously cryopreserved samples were thawed in warm TIL-CM, washed with 1 × PBS, and resuspended at 1 × 10⁶ cells/mL in 1% BSA in 1× PBS before being sent for oil immersion on the 10X platform.

PDAC tumors were washed with sterile PBS and minced with sterile scalpels in human wash media (DMEM/F12, 0.01M HEPES, 1× GlutaMAX supplement, 10% FBS) until tumor fragments were <1 mm³. Tissue suspensions were centrifuged, and the supernatant was removed and resuspended in dissociation buffer [collagenase II (5 mg/mL), dispase (1 mg/mL) in human wash media]. The samples were placed on a rotation plate in a 37°C incubator for 45 minutes. The single-cell solution was then centrifuged and filtered. Single cells were collected as a pellet following centrifugation of the flow-through and resuspended in warm resuspension buffer (1% BSA in 1× PBS). Cells were enumerated using a hemocytometer and sent for oil immersion on the 10X platform.

Single-cell capture and library construction were performed with the 10X Genomics Chromium Single-Cell 5′ kits v.1.0 (product codes 1000014, 1000020, and 1000151) with TCR enrichment (1000005) for cohort MDA1 and 3′ V.2 kits (product codes120237, 120236, and 120262) for MDA2, according to the manufacturer's instructions. Briefly, cells were loaded into the Chromium Single-Cell Chip A for a recovery target of 10,000 cells. Reverse transcription was performed on a Bio-Rad T100 thermal cycler, and the barcoded cDNA was purified with Dynabeads (Thermo Fisher Scientific, 37002D) prior to 14 cycles of cDNA amplification. Of this transcriptome cDNA, 2 μL was used for TCR enrichment and subsequent TCR library construction in the MDA1 cohort. Per the manufacturer's protocol, up to 50 ng (or 20 μL) of transcriptome cDNA was used for single-cell library construction. Transcriptome libraries were sequenced on a HiSeq 400 (Read 1, 26 cycles; Index 1, 8 cycles; Read 2, 91 cycles), and TCR libraries were pooled and sequenced 150 cycles paired-end on a MiSeq (Illumina). For transcriptome libraries, the median sequencing depth for each sample was targeted for a median of 30,000 reads/cell; for TCR libraries, 1,000 reads/cell. Read counts can be found in Supplementary Table S3.

The 10X Genomics Cell Ranger pipeline (v. 3.1.0) was used to demultiplex and generate unique molecular identifier (UMI) matrices for all samples in the MDA1 cohort. UMI matrices for the MDA2 and PUMCH cohorts were also generated from the provided FASTQ files using this version of the Cell Ranger pipeline. The UMI data were then processed in R (v.3.5.3) using the Seurat package (V.3). Cells with more than 6,000 genes (doublet removal), fewer than 200 genes, greater than 20% UMI in mitochondrial genes, or a UMI:gene ratio greater than 10 or less than 1.3 were filtered. To filter out non–T cells from the data sets, all samples were combined, and the top 5,000 variable genes were reduced to 50 top principal components and clustered using Seurat's SNN clustering function. Clusters with significant average expression (Wilcoxon test) of CD3 genes (i.e., T cells) were then extracted for further analysis. This gene selection, principal component analysis, and clustering step were performed once more to further remove clusters not expressing CD3 genes.

To account for library chemistry, and other possible unknown differences between cohorts, T cells were integrated using the methods described in https://satijalab.org/seurat/v3.0/integration.html. Briefly, the top 5,000 variable features were selected with TCR variable genes (“^TR[ABGD][VDJ]”) excluded. These genes were used to compute 20 “anchors” with the FindIntegrationAnchors function, and subsequently “integrated” by cohort with IntegrateData to produce a batch-corrected data set. The first 30 principal components of the integrated expression data were used for subsequent UMAP embedding.

T cells with detectable expression of either CD4 or CD8 genes were assigned to their respective clusters. The first 30 principal components (excluding these genes) from a randomly selected training set of 75% of these cells were used to train a random forest classifier (RandomForest package in R) for CD4⁺ versus CD8⁺ T cells. The remaining 25% were used as a validation set. This classifier was then used to classify T cells with no detectible CD4 or CD8 expression (likely due to gene dropout) as either CD4⁺ or CD8⁺ T cells and then further clustered separately. After CD4/CD8 classification and separation, the top 3,000 variable features were reselected for each group, and principal components were computed. The top 30 principal components were used to perform clustering with Seurat's FindNeighbors and FindClusters functions. Top marker genes for each cluster were identified with Seurat's FindAllMarkers function, and P values were determined by the Wilcoxon test (Supplementary Tables S4 and S5).

We performed RNAscope on 5 μmol/L formalin-fixed, paraffin-embedded sections of PDAC and adjacent, uninvolved pancreas tissue using the ACD RNAscope Multiplex Fluorescent kit v2 and 4-Plex ancillary kit (cat nos. 323100 and 323120). Slides were prepared as described in the kit's user guide (Doc. No. 323100-USM), baked briefly at 60°C for 1 hour in an ACD HybEZ II oven (PN 321710/321720), and then deparaffinized in xylenes and 100% ethanol. Exogenous peroxides were blocked for 10 minutes, using the provided hydrogen peroxide, and target retrieval was performed for 10 minutes in the provided retrieval buffer in a 5.5qt Hamilton Beach digital steamer (MODEL: 37530Z). A hydrophobic barrier was applied around the tissue sections, and they were treated with RNAscope protease plus for 30 minutes and washed in water. Probes sets were then applied, and the slides were probed for 2 hours before washing (5 minutes ×2) in the kit's wash buffer and stored overnight in 5× SSC buffer. Hybridization Amplication, four rounds of HRP development, Opal Dye conjugation, and blocking were performed per the kit's instructions, washing twice, for 5 minutes with agitation in the provided wash buffer after each step. Slides were then stained with DAPI, washed in PBS, mounted in Invitrogen Prolong Diamond Antifade (cat. P36965), and dried overnight in the dark. Imaging was performed on a Nikon Eclipse Ti2E microscope (PN MEA54000) using the NIS Elements AR 5.30.05 software as recommended by the manufacturer.

For samples with T cells detected on the initial screen, further sections were then stained for T-cell subpopulations with the following 8 panels: panel 1 (CD4, CXCR4, CCR7, and IFIT1); panel 2 (CD4, CXCR4, FOXP3, and CXCL13); panel 3 (CD8A, CD4, GZMK, and ZNF683); panel 4 (CD8A, CXCR6, ZNF683, and IL7R); panel 5 (CD8A, GZMK, CXCL13, and PRF1); panel 6 (CD8A, CXCR6, CCR7, and IL7R); panel 7 (CD8A, GZMK, IFIT1, and CXCL13); tumor panel (KRT19/MUC1, CD3E/CD3D/CD3G, CD8A, and CD4). The panels were composed of the following probes from ACDBio: Hs-CD3-pool-C2 (426621-C2), Hs-CD4 (605601), Hs-CD4-C4 (605601-C4), Hs-CXCR4-C3 (310511-C3), Hs-CCR7-C2 (410721-C2), Hs-IFIT1 (415551), Hs-CXCL13-C4 (311321-C4), Hs-PRF1 (407381), Hs-GZMK (475901), Hs-GZMK-C2 (475901-C2), Hs-CXCR6 (468461), Hs-IL7R-C4 (408841-C4), Hs-ZNF683-C2 (586921-C2), Hs-MUC1 (310391), and Hs-KRT19 (310221). The probes were visualized using the following OPAL dye fluorophores from Akoya Biosciences: Opal 520 (FP1487001KT), filter: FITC, color: green, dilution: 1:500; Opal 570 (FP1488001KT), filter: Cy3, color: cyan, dilution: 1:750; Opal 620 (cat: FP1495001KT), filter: Texas Red, color: red, dilution: 1:750; Opal 690 (FP1497001KT), filter: Cy5, color: white, dilution: 1:750.

Publicly available scRNA-seq data from PBMCs from healthy donors and patients with PDAC were obtained from 10X Genomics [10K PBMCs from healthy donors (v3, 28 × 91), Single-Cell Immune Profiling data set by Cell Ranger 3.0.0, 10X Genomics (2018 November 19)] and Steele and colleagues (41). The data were filtered and processed as described previously in this article for the other single-cell data set. T cells were then integrated by cohort (10× PBMCs, Steele-Healthy PBMCs, and Steele-PDAC PBMCs) and clustered separately by CD4 and CD8 T cells as described for the other data sets in this article. The top 20 cluster marking genes by P value (Wilcox test) were used to compare with the similarly selected top marker genes from tissue T cells.

Pseudotime analysis was performed on cells in the CD4 and CD8 populations separately using Monocle 3 (version 0.2.2.0; https://github.com/cole-trapnell-lab/monocle3). Expression data were UMAP embedded using the Monocle function “reduce_dimension” with default parameters. The trajectory graph was inferred with the function “learn_graph” with the minimal branch length set to 15 and close_loop = FALSE. CD4 and CD8 trajectories were rooted in the clusters CD4-CCR7 and CD8-CCR7/IL7R, respectively.

For all tumor samples with at least 50 T cells, cluster fractions were calculated by dividing the number of cells in a cluster by the total number of T cells in the sample. For each cluster, the cluster fractions across samples were then pairwise Spearman correlated to all other cluster fractions using the cor() function in R. The Spearman Rho values obtained were then hierarchically clustered with Ward.D2 linkage for the graph and tested for being significantly different from 0 using a t test. P values were corrected by the Benjamini–Hochberg method for multiple testing.

TCR calls on the single-cell TCR-sequencing data were performed with the “cellranger vdj” function of the Cell Ranger software suite (v.3.1.0). Cells with only a TCR alpha or beta chain detectable whose TCR nucleotide sequence exactly matched that of another clonotype with both the alpha and beta chain detectable were reassigned to that clonotype. Clonotype data were then attached to single-cell metadata by exact matches in cell barcodes within a sample.

TCR Circos plots were created with the code adapted from: https://github.com/aislyn/TCR_circos. For each clonotype (defined by a unique CDR3 nucleotide sequence), the median expression was calculated across cells, and then hierarchically clustered using Euclidean distance between the medoids and ward.d2 linkage (using hclust function in R). The hierarchal graph is shown at the center of the plot with each leaf/spoke representing one clonotype. For each clonotype, the fraction of cells in that clonotype coming from each patient and each cluster was calculated and plotted as the inner and outer solid rings, respectively. Finally, the number of cells detected in each clonotype was plotted on the outermost ring of the plot.

For the Upset plots, the number of TCRs existing in single or multiple states was quantified by filtering for clonotypes with at least 2 cells detected (such that overlap was possible), and then the clonotype was classified by overlap state (e.g., CD8-GZMK only, CD8-ZNF683/CD8-CXCR6/IL7R overlap, etc.). The number of clonotypes observed in each overlapping state was then counted and plotted using the upsetR package in R. In order to quantify whether an overlap state existed more or less frequently than expected by chance (because high-frequency clonotypes have more chances to overlap between states by chance than lower-frequency ones), we performed a permutation test for our overlap counts by performing the same computation on permutated clonotypes and expression states within patients 1,000 times.

For each patient and cluster with at least 5 clonotypes present in the cluster, the number of cells in the clonotype was tallied, and the following indexes were calculated as follows: Chao1 species richness (computed with the Chao1 function from the “fossil” package in R), Gini coefficient (computed with the Gini function from the “DescTools” package in R), Diversity Evenness 50 (DE50) was the (frequency-ordered) number of clonotypes necessary to account for 50% of cells divided by the total number of clonotypes. The Simpson index was calculated as the sum of the clonotype fractions squared, the Shannon diversity index was calculated as the negative sum of the clonotype fractions multiplied by their natural log, and the Pielou evenness was calculated as the Shannon index divided by the natural log of the number of clonotypes.

Fresh tumor samples were cut into 1 to 3 mm³ fragments, and five fragments were placed in G-Rex10 flasks (Wilson Wolf) containing 20 mL of TIL culture media [TIL-CM; RPMI-1640 with GlutaMAX (Gibco/Invitrogen), 1 × penicillin and streptomycin (Gibco/Invitrogen), 50 μmol/L 2-mercaptoethanol (Gibco/Invitrogen), 20 μg/mL gentamicin (Gibco/Invitrogen), and 1 mmol/L sodium pyruvate (Gibco/Invitrogen)] with 6,000 IU/mL IL2, 10 μg/mL 4-1BB mAb, and 30 ng/mL anti-CD3 (OKT3) as previously described (49). Four to five days after culture initiation, 20 mL of additional TIL-CM with 6,000 IU/mL IL2 was added for a total volume of 40 mL. Half-media changes were done every 3 to 4 days with fresh TIL-CM containing 6,000 IU/mL IL2 for up to 35 days or until the cells formed a thick layer completely covering the bottom of the flask. Cell suspensions were collected and cryopreserved in FBS plus 10% DMSO.

A purified, human IgG4 monoclonal antibody (mAb) against human CD137/4-1BB, urelumab (663513), was kindly provided by Bristol-Myers Squibb. Human recombinant IL2 (Proleukin) was generously provided by Prometheus Therapeutics and Diagnostics. GMP-grade soluble anti-CD3 antibody (OKT3 clone) was obtained from Miltenyi Biotec (Bergisch Gladbach).

Due to the (intentional) highly proliferating state of the cultured T-cell samples, a large fraction of the T cells formed a cluster of cycling cells. Noncycling cell clusters were used to train a random forest classifier (similar to the CD4/CD8 classification in the fresh T cells), which was used to reclassify cycling cells to the most similar noncycling cluster.

The data generated in this study are publicly available on the Sequence Read Archive under BioProject ID PRJNA806978

D. Sakellariou-Thompson reports a patent for methods for expansion of TILs and use thereof licensed. M.-A. Forget reports a patent for methods for expansion of TILs and use thereof pending and licensed, and is an advisory board member for Proteios Technology. A. Reuben reports personal fees and other support from Adaptive Biotechnologies outside the submitted work. C.D. Tzeng reports other support from PanTher and Ethicon outside the submitted work. S. Pant reports other support from Xencor, 4D Pharma, Zymeworks, Ipsen, Novartis, and Janssen outside the submitted work. D.R. Fogelman reports other support from Merck outside the submitted work. A. Maitra reports royalties for a pancreatic cancer biomarker test from Cosmos Wisdom Biotechnology, and this financial relationship is managed and monitored by the UTMDACC Conflict of Interest Committee. A. Maitra is also listed as an inventor on a patent that has been licensed by Johns Hopkins University to ThriveEarlier Detection, and serves as a consultant for Freenome and Tezcat Biotechnology. C.L. Haymaker reports personal fees from Nanobiotix, other support from the Society for Immunotherapy of Cancer and Briacell, and grants from Iovance, BTG, Dragonfly Therapeutics, and Sanofi outside the submitted work, as well as a patent for WO2021 167908 A1, PCT/US2021/018252 licensed. C. Bernatchez reports grants from generous philanthropic contributions to the MD Anderson Cancer Moon Shots Program during the conduct of the study; grants from Iovance Biotherapeutics and Obsidian Therapeutics, and personal fees from Myst Therapeutics outside the submitted work; and a patent for B and T lymphocyte attenuator marker for use in adoptive T-cell therapy issued, a patent for methods for expansion of TILs and use thereof pending and licensed, and a patent for human 4-1BB agonist antibodies and methods of use thereof pending. No disclosures were reported by the other authors.

A. Schalck: Conceptualization, data curation, software, formal analysis, investigation, visualization, methodology, writing–original draft, writing–review and editing. D. Sakellariou-Thompson: Concep­tualization, investigation, visualization, methodology, writing–original draft, writing–review and editing. M.-A. Forget: Conceptualization, validation, visualization, writing–original draft, writing–review and editing. E. Sei: Data curation, investigation, project administration. T.G. Hughes: Investigation. A. Reuben: Writing–review and editing. S. Bai: Data curation, investigation. M. Hu: Investigation. T. Kumar: Investigation. M.W. Hurd: Resources, data curation, clinical support. M.H. Katz: Resources, clinical support. C.D. Tzeng: Resources, clinical support. S. Pant: Resources, clinical support. M. Javle: Resources, clinical support. D.R. Fogelman: Resources, clinical support. A. Maitra: Resources, supervision, funding acquisition, writing–review and editing, clinical support. C.L. Haymaker: Conceptualization, supervision, validation, writing–review and editing. M.P. Kim: Conceptualization, resources, supervision, funding acquisition, methodology, project administration, writing–review and editing, clinical support. N.E. Navin: Conceptualization, resources, supervision, funding acquisition, methodology, writing–original draft, project administration, writing–review and editing. C. Bernatchez: Conceptualization, resources, supervision, funding acquisition, writing–original draft, project administration, writing–review and editing.

The study was supported by philanthropic contributions to The University of Texas MD Anderson Moon Shots Program and the Sheikh Khalifa bin Zayed Foundation. N.E. Navin is a Grady Saunders Distinguished Professor, an AAAS Fellow, AAAS Wachtel Scholar, Damon-Runyon Rachleff Innovator, Andrew Sabin Fellow, and a Jack & Beverly Randall Innovator. This work was supported by grants to N.E. Navin from the NIH NCI (RO1CA240526 and RO1CA236864), the CPRIT Single Cell Genomics Center (RP180684), and the MD Anderson Sequencing Core Facility Grant (CA016672). This work was also supported by grants to M.P Kim from the NIH NCI (K08CA218690), and the Ben and Rose Cole Foundation, and to T.G. Hughes from the NIH (T32 CA009599). The authors thank Bristol-Myers Squibb for their generous contribution of the agonistic anti-4-1BB antibody urelumab (BMS-663513). Human recombinant IL2 (Proleukin) was generously provided by Prometheus Therapeutics & Diagnostics. Editorial support was provided by Bryan Tutt, Scientific Editor, Research Medical Library, The University of Texas MD Anderson Cancer Center. This paper is dedicated to the memory of Dr. Gauri Varadhachary, who is dearly missed by all her colleagues in the MD Anderson pancreatic cancer team.

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.