Tumor-infiltrated double-negative regulatory T cells predict outcome of T cell-based immunotherapy in nasopharyngeal carcinoma

Abstract

Adoptive cell therapy (ACT) using tumor-infiltrating lymphocytes (TILs) has demonstrated clinical success in solid tumors. We analyze 47 TIL infusion products and 62 pretreatment tumor microenvironments (TMEs) from a randomized phase 2 clinical study of concurrent chemoradiotherapy plus TIL-ACT (NCT02421640). Using single-cell and bulk RNA sequencing along with flow cytometry, we identify 14 CD3⁺ T cell clusters within 26 TIL infusion products: 11 CD3⁺CD8⁺ TILs, 2 CD3⁺CD4⁺ TILs, and 1 CD3⁺CD8^-CD4^- double-negative (DN) TIL. (DN) TILs, significantly associated with poor TIL-ACT outcomes, exhibit an activated regulatory T cell-like phenotype and include two CD56⁺ and four CD56^- subsets. Among them, CD56^-KZF2⁺ (DN) TILs are predominantly suppressive. (DN) TILs inhibit CD8⁺ TIL expansion via Fas-FasL, transforming growth factor β (TGF-β), and interleukin (IL)-10 signaling. Distinct CD8⁺ T subsets differentially impact on TIL-ACT outcomes, while 9 baseline TME gene signatures and 14 intracellular T cell genes hold prognostic value. Our findings identify predictive TIL subsets and biomarkers for TIL-ACT outcomes.

Keywords: TIL infusion products; adoptive cell therapy; double-negative T cells; nasopharyngeal carcinoma; single-cell RNA sequencing.

Graphical abstract

Figure 1

Study design and single-cell transcriptomic landscape of TIL infusion products (A) Graphical overview of the experimental design and bioinformatics workflow. In our clinical trial cohort, bulk RNA sequencing was performed on tumor tissue from NPC patients before CCRT+TIL (n = 22) or CCRT alone (n = 40), and TIL infusion products were isolated from pretreatment tumor tissue from these patients, expanded ex vivo, and subsequently prepared for single-cell (sc)RNA sequencing (n = 26) and flow cytometry (n = 47). (B) Uniform manifold approximation projection (UMAP) plot of 65,978 quality-controlled T cells colored by 14 T cell subsets with the distinct G1 phase, S phase, and G2M phase. The cell-cycle phase of 14 T cell subsets was labeled in the UMAP plot and determined by the cell-cycle phase scores calculated by the Seurat package. (C) Heatmap showing the expression of the top 50 differentially expressed genes among cell subsets obtained by cell type annotation utilizing the expression of canonical marker genes in 14 cell subgroups. Information on the 14 cell subsets is displayed on the right. (D) UMAP plots showing the distribution of indicative T cell subsets between the non-progression group and progression group. (E) Dot plot showing the composition of CD3⁺CD4⁻CD8⁻, CD3⁺CD4⁺, and CD3⁺CD8⁺ T cells in NPC TIL infusion products according to flow cytometry (n = 47). Scatter dot plots show individual values and mean ± SEM. (F) Bar plot comparing the frequency of 14 cell subgroups in the TIL infusion products from NPC patients in the non-progression (n = 15) and progression groups (n = 11). The data are presented as mean ± SEM. The Wilcoxon rank-sum test was used to determine the significance. ∗p < 0.05; ns, not significant; adjusted for multiple comparisons using the Benjamini-Hochberg procedure. (G) Boxplot comparing the frequency of (DN) TILs in the non-progression (n = 28) and progression groups (n = 19) according to flow cytometry. Median and interquartile were shown in the boxplots. Wilcoxon rank-sum test was used to determine the significance. (H) Kaplan-Meier survival curves for progression-free survival (PFS) and overall survival (OS) in NPC patients stratified according to the abundance (high abundance vs. low abundance) of the C8_DN_ T_cells cluster inferred from (sc)RNA sequencing data (n = 26). p values were determined based on the two-sided log rank test. Optimal cutoff values were defined using the “survminer” R package. See also Figure S1 and Table S1.

Figure 2

Transcriptional characterization and developmental trajectory of CD3⁺CD8⁻CD4⁻ (DN) TILs (A) Heatmap showing the expression levels of genes encoding NK cell markers, Treg-related molecules, immune checkpoint proteins, cytokines, and (DN) T cell-related markers among CD8⁺, CD4⁺, and (DN) TIL clusters. (B) Gene set enrichment analysis (GSEA) showing significant upregulation of TGF-β (left) and IL-10 (right) signaling pathways in (DN) TILs compared with other TIL subsets. (C) Violin plot showing the expression of a regulatory T cell signature across (DN) TIL, CD4⁺ TIL, and CD8⁺ TIL clusters. A dashed line indicates the median of the signature score of the (DN) TIL cluster. Median and interquartile were shown in the boxplots. Non-parametric Kruskal-Wallis test was used to determine the significance. (D) Representative overlap histogram for the expression of IL-10, TGF-β, IKZF2, FOXP3, and CTLA4 among (DN) TILs, CD4⁺ TILs, and CD8⁺ TILs determined by flow cytometry. (E) Developmental trajectory of 16,978 (DN) T cells from tumor tissues and peripheral blood of naive treated NPC patients (n = 10) inferred by Monocle 2 and CytoTRACE algorithm. Solid and dotted lines denote distinct cell fates based on expression profiles, with colors indicating the origin of (DN) T cells, pseudo-time, and CytoTRACE score (GEO: GSE162025). (F) Dot plot showing the correlation between (DN) TIL signature scores of TMEs in this study and each state of (DN) T cells (GEO: GSE162025; PB, Cell Fate_1, and Cell Fate_2). The sizes and colors of the circles represent the strength of the relationship, assessed using Spearman’s correlation test. PB, peripheral blood. (G) Representative overlap histogram for the expression of IL-10, TGF-β, IKZF2, FOXP3, and CTLA4 in the (DN) T cells from TIL infusion products, peripheral blood of NPC patients, or healthy donors determined by flow cytometry. See also Figures S2 and S3.

Figure 3

Transcriptomic characterization and association of CD56⁻ (DN) TILs with clinical outcome (A) UMAP plot depicting 2,609 quality-controlled (DN) TILs colored by six cell subsets. (B) UMAP plot of (DN) TILs, with cells colored based on the relative normalized expression of NCAM1 (encoding CD56), FCGR3A (encoding CD16), and IKZF2. (C) Volcano plot showcasing differentially expressed genes (DEGs) between CD56⁻ (DN) TILs and CD56⁺ (DN) TILs. Red dots and blue dots indicate up-regulated genes (n = 770) and down-regulated genes (n = 316) in CD56⁻ (DN) TILs, respectively. Selected genes were labeled. (D) Heatmap showing the expression of genes encoding TCR proteins, NK cell markers, immune checkpoint proteins, cytokines, and (DN) T cell-related markers among CD56⁻ and CD56⁺ (DN) TIL clusters. (E) GSEA plots showing that Treg signaling pathways were significantly up-regulated in CD56⁻ (DN) TILs compared with CD56⁺ (DN) TILs. (F) Violin plot showing the expression of an activated DN Treg signature across four (DN) TIL clusters (CD56⁻IKZF2⁺, CD56⁻IKZF2^-, CD56⁺IKZF2^-, and CD56⁺IKZF2⁺ (DN) TILs). A dashed line indicates the median of the signature score of the CD56⁻IKZF2⁺ (DN) TIL. Median and interquartile were shown in the boxplots. Kruskal-Wallis test was used to determine the significance. (G) Kaplan-Meier survival curves for PFS (top) and OS (bottom) in NPC patients stratified according to the abundance (high abundance vs. low abundance) of CD56⁻IKZF2⁺ (DN) TILs inferred from (sc)RNA sequencing data (n = 26). p values were determined based on the two-sided log rank test. See also Figure S4.

Figure 4

Associations of CD8⁺ TIL subsets with (DN) TILs and clinical outcome (A) Dot color represents the communication probability of the specific ligand-receptor pairs, including TGF-β, IL-10, and Fas-FasL signaling between the indicated sender cluster and receiver clusters. (B) Representative histograms and statistical graph showing carboxyfluorescein succinimidyl ester (CFSE) dilution (left) and proliferation inhibition rates (right) of (DN) TILs with CD4⁺ and CD8⁺ naive T cells at a ratio of 4:1, 2:1, 1:1, and 1:2. The data are presented as mean ± SD, and three biological replicates were included. (C) Representative histograms and statistical graph showing CFSE dilution (left) and proliferation inhibition rates (right) of (DN) TILs with CD4⁺ and CD8⁺ naive T cells at a 1:1 ratio with neutralizing antibodies against TGF-β and IL-10, as well as the presence or absence of a Fas-FasL signaling antagonist. The data are presented as mean ± SD, and three biological replicates were included. p values were evaluated by t test (two-sided). ∗∗p < 0.01; ∗∗∗p < 0.001. (D) Bubble plot showing the correlation between the frequency of (DN) TILs and other TIL subsets in TIL infusion products (n = 26). p values were determined by Spearman correlation analysis. ∗p < 0.05. (E) Dot plot showing the correlation between the frequency of (DN) TILs and CD8⁺ TILs in TIL infusion products (n = 26). p values were determined by Spearman correlation analysis. (F) Representative histograms and statistical graph showing CFSE dilution (left) and proliferation inhibition rates (right) of (DN) TILs with autologous expanded CD8⁺ TILs at ratios of 4:1, 2:1, and 1:1. Conventional CD4⁺ (i)Tregs were included as a control. The data are presented as mean ± SD, and three biological replicates were included. (G–K) Kaplan-Meier survival curves of PFS and OS in NPC patients stratified according to the abundance (high vs. low) of total CD8⁺ TILs (G), C2_proliferation_CD8_T_cells-PCNA (H), C4_MHC_II_CD8_T_cells (I), C6_proliferation_CD8_T_cells-TOP2A (J), and C12_tissue_resident_memory_CD8_T_cells (K) inferred from (sc)RNA sequencing data (n = 26). p values were determined based on the two-sided log rank test. See also Figures S5 and S6.

Figure 5

Relationship between transcriptome features in the baseline tumor microenvironment and clinical outcome and TIL subset composition (A) Heatmap showing the expression of the indicative gene signatures of the baseline TME in 62 NPC patients. (B) Bubble plots showing the Cox regression analysis of gene expression signatures in patients treated with CCRT+TIL-ACT or CCRT alone. The p values and hazard ratios were estimated from a stratified Cox model with the low-gene-signature-score group as the reference group. (C) Kaplan-Meier survival curves for PFS in NPC patients treated with CCRT+TIL-ACT (n = 22) stratified according to the indicative gene signature expression (high score vs. low score) of NK cell, MHC class II, CD8⁺ cytotoxic T cell, CD39⁻CD69⁻ stem-like T cell, NeoTCR8, BCR, LAMP3⁺ DC, interferon receptor, and mismatch repair signatures. (D) Heatmap showing the correlation between gene signature scores and the frequency of 14 TIL subsets in TIL infusion products, with significance assessed by the Spearman correlation. ∗p < 0.05, ∗∗p < 0.01. See also Figure S7 and Table S2.

Figure 6

Prognostic value analysis of 14 differentially expressed genes in TIL infusion products from non-progression and progression groups (A) Euler diagrams represent the overlapping genes enriched in the non-progression group and associated with PFS and OS. Thirteen overlapped genes were associated with non-progression, PFS, and OS. The right diagram represents 13 overlapped genes significantly associated with favorable PFS and OS. The colors of dots represent the influence of the gene on prognosis: red indicates a protective factor and green indicates a risk factor. (B) Pathway enrichment analysis of the 13 overlapped genes. The top 8 significant pathways are shown. (C) Dot plot showing the expression of 13 overlapped genes for each T cell subset in NPC TIL infusion products. For each gene, dot colors represent the average expression in each cell type (scaled and log-normalized), whereas size reflects the percentage of cells with detectable expression in each TIL subset. (D) Dot plot showing the expression of TNFRSF18 for each T cell subset in NPC TIL infusion products. (E) Kaplan-Meier survival curves for PFS (top) and OS (bottom) in NPC patients stratified according to the expression of TNFRSF18 based on a median split inferred from (sc)RNA sequencing data (n = 26). p values were determined based on the two-sided log rank test. See also Figure S8.