4-1BB+ Tregs and inhibitory progenitor exhausted T cells confer resistance to anti-PD-L1 and anti-CTLA-4 combination therapy

Abstract

Predictors of immune checkpoint inhibitor response in cancer remain elusive. From a previous phase 2 neoadjuvant immunotherapy window-of-opportunity study, we present the single-cell RNA and T cell receptor (TCR) sequencing analysis of 57 pre- and post-treatment tumor biopsies from head and neck cancer patients treated with durvalumab (anti-PD-L1) alone or with tremelimumab (anti-CTLA-4), identifying key cellular and molecular predictors of immune checkpoint inhibitor (ICI) response. Malignant cells and neutrophil senescence promote ICI response. While CXCL13⁺ exhausted T (Tex) cells enhance response through 4-1BB signaling, anti-CTLA-4 induces 4-1BB⁺ regulatory T cells (Tregs) restricting ICI efficacy. These opposing roles of 4-1BB in different cellular contexts may explain the limited benefit of combinatorial immunotherapy observed in clinical trials. We identify two subsets of tumor-reactive progenitor Tex (Tpex): ICI-responsive Tpex1 and ICI-resistant Tpex2, a subset characterized by KLRB1 and IL17R. The balance of Tpex1 and Tpex2 associates with ICI response across multiple cancers, offering insights into sustaining response. This study was registered at ClinicalTrials.gov (NCT03737968).

Keywords: dual ICB; durvalumab; head and neck squamous cell carcinoma; immune checkpoint inhibitor; immunotherapy; progenitor exhausted T; single-cell RNA; single-cell TCR; tremelimumab.

Graphical abstract

Figure 1

Overall research design of HNSCC immunotherapy clinical trial Study design and sampling scheme of the window-of-opportunity study. Study timeline for enrolled samples and their description of treatment, pathological response, and the generated data type are illustrated.

Figure 2

Malignant senescence program associated with ICI sensitivity (A) Uniform manifold approximation and projection (UMAP) of identified malignant cells pre-treatment (left) and post-treatment (right) colored by patient origin. (B) Overlap of selected patient-specific NMF modules. Rows and column clustered with hierarchical clustering. Color indicates −log10 of p value calculated via one-sided hypergeometric test (overlap of gene sets). (C) Pie chart proportion of nine meta-programs in pre- and post-ICI responders and non-responder malignant cells. Each cell is labeled a meta-program with the highest gene set module score of the nine. (D) Scatterplot of arithmetic mean of epithelial senescence meta-program score for each pre-ICI patient sample and their tumor regression values post-ICI. Pearson correlation coefficient and its p value are denoted. (E) Variance of senescence module score for malignant cells in pre- or post-ICI samples. p value is calculated with Levene’s test. (F) Variable importance score for each predictor (calculated as the mean decrease of Gini impurity when a variable is chosen to split a node) of 8 random forest models. Data are represented as mean ± SEM. (G) Area under reciever operating characteristic (AUROC) curve for public signatures in NSCLC dataset, melanoma dataset, and GBM dataset. Dashed line indicates random expectation value of 0.5. Gene signature scores are calculated via gene set variation analysis (GSVA) for each sample. (H) ROC curve for ensemble random forest model without senescence module as a predictor (black, immune model) and with it (red). For both models, tumor purity, stromal score, and immune signature (Ayer et al.²¹) were used as predictors. FPR, false positive rate; TPR, true positive rate

Figure 3

ICI response predictors in CD8⁺ T states (A) UMAP of CD8⁺ T cells extracted, re-normalized, and re-clustered. Tex, exhausted T cell; Tmem, memory T cell; Teff, effector T cell; Tstr, stressed T cell; Tprolif, proliferating T cell. (B) Network representation of enriched KEGG pathway for pre-ICI CD8⁺ T cells upregulated in responders (R) compared to non-responders (NR), visualized via emapplot() of clusterProfiler package. Colors indicate the enriched subcluster within CD8⁺ T cells. Size represents the overlapping genes in each term. Color within the term node represents the contribution of that cluster to the term based on the number of overlapping genes from each cluster. (C) Boxplot of stress-associated heat shock protein (HSP) signature identified to be upregulated in Tstr, significantly less so in responders at baseline. p values were calculated via Wilcoxon rank-sum test. (D) The number of expanded T cell clonotypes in post- versus pre-treatment tumor biopsies, classified by expansion of frequency (above 2 red and above 5 green) and expansion of proportion and frequency (blue). Overlaid scatterplot indicates tumor regression value post-therapy. (E) Top 50 intracellular ligand-receptor pairs across CD8⁺ T subclusters that significantly increased interaction post-therapy for each response group as measured by MultiNicheNet. TNFSF9 and KLRB1 are highlighted in red. (F) Feature plot of TNFSF9 expression in CD8⁺ T cell UMAP dimension and violin plot of TNFSF9 expression for each subclusters (left) and feature plot of TNFRSF9 expression in CD8⁺ T cell UMAP dimension and violin plot of TNFRSF9 expression for each subclusters (right). p values were calculated via Wilcoxon rank-sum test for each cluster vs. all others, adjusted by Benjamini-Hochberg method. (G) Cell-type-specific ligand-receptor pseudo-bulk (by patient) product value from MultiNicheNet of top 8^th to 11^th that showed statistically significant increase in responder vs. non-responders post-therapy compared to baseline. Two-sided t test was additionally performed and denoted at the top of each group comparison (pre- vs. post-ICI). p value significance: ∗p < 0.05, ∗∗p < 0.01, ∗∗∗p < 0.001, and ∗∗∗∗p < 0.0001.

Figure 4

CXCL13⁺ Tex network topology as a predictor for ICI response (A) Conceptual depiction of the network analysis using scHumanNet performed for CXCL13⁺ Tex cells. (B) Similarity of patient-specific network nodes calculated with pairwise Euclidean distance from adjacency matrix of a union gene set. p values were calculated via two-sided t test. Additionally, difference of variance was tested with Fligner-Killeen (F-K) test, and its p values are denoted. (C) Cell-type-specific network of responder CXCL13⁺ Tex from pre-ICI treatment group. Gene nodes are colored according to each different subcommunity determined by Louvain clustering with top 3 nodes depicted for each community. The 30 genes ordered by centrality in a responder community that decreased in centrality in the corresponding non-responder network are summarized as a table below (green) and termed responder network signature. (D) Conceptual depiction of the in silico perturbation analysis using a fine-tuned (all immune cells) Geneformer foundation model. For each gene perturbed, the embedding shift of responder cells (CXCL13⁺ Tex only) toward non-responder cells was measured through cosine similarity. (E) Two-dimensional UMAP representation of fine-tuned Geneformer embeddings (512 dimensions) with immune cells as input, labeled for ICI response status. (F) Cosine similarity of shifted embedding for each gene set. 100 random genes were sampled from perturbed genes as control. The 21 genes of the immune hub signature show positive shift toward non-response when perturbed. p value significance: ∗p < 0.05, ∗∗p < 0.01, ∗∗∗p < 0.001, and ∗∗∗∗p < 0.0001; ns, non-significant.

Figure 5

Identification of ICI-resistant Tpex subpopulation (A) Conceptual depiction of TCR-based selection and identification of tumor-specific exhausted T cells (Tex) and progenitor exhausted T cells (Tpex). (B) Dot plot of selected genes significantly upregulated in specific subsets of tumor-specific CD8⁺ T cells. (C) Average expression genes in Tex and Tpex (Tpex1 + Tpex2). Red dots indicate significantly upregulated genes in each group calculated by two-sided t test. Key genes associated with memory/exhaustion/dysfunction are highlighted in green. (D) Scatterplot of ratio differences (post-ICI minus pre-ICI) for patients with TCR information. Pearson correlation coefficient and its p value is depicted. Ratio were calculated as Tpex1/(Tpex2 + Tex). (E) RNA velocity streamline plot for tumor-reactive CD8⁺ T cells in non-responder group (left) and responder group (right). (F) Cell transition probability density for Tpex1 to Tpex2 (left) and persistence of Tpex2 (right). (G) Reverse transition of Tex to Tpex1 (left) and Tpex2 (right), for non-responder and responder group. p values are calculated via Kolmogorov-Smirnov test. (H) Proportion of tumor-specific CD8⁺ T subclusters identified via Seurat label transfer algorithm. Datasets were divided by cancer type, ICI treatment, and response group as annotated by the original authors. p value significance: ∗p < 0.05, ∗∗p < 0.01, ∗∗∗p < 0.001, and ∗∗∗∗p < 0.0001; ns, non-significant.

Figure 6

Combinatorial ICI effect inhibition in the CD4⁺ T compartments (A) UMAP of CD4⁺ T cells extracted, re-normalized, and re-clustered. (B) CTLA-4 expression colored in the CD4⁺ T cell UMAP dimension. (C) Percentage of clones that are expanded (n > 1) and that are single (n = 1) for three CD4⁺ T cell types. p values are calculated via Wilcoxon rank-sum test. (D) Mosaic plot of post-ICI CD4⁺ T cells that express CTLA-4 (expression above 0). The yellow line indicates the expected ratio. Pearson residual p values are colored red or blue if significantly depleted or enriched, respectively. (E) Transition probability density of selected cell-to-cell transitions for monotherapy (D) and combination group (D + T). p values are calculated via Kolmogorov-Smirnov test. (F) RNA velocity streamline plot for CD4⁺ T cells in monotherapy group (left) and combination group (right). (G) Volcano plot of differentially expressed genes in 4-1BB⁺ Tregs with positive log fold values for combination group. Genes above 0.25 log2 fold change and under adjusted p value (Benjamini-Hochberg) of 0.05 are colored in red. Green text indicates genes associated with IL-2 signaling. (H) Split violin plot for module score of “IL-2-STAT5 signaling” term (collected from BioPlanet 2024) in CD4⁺ T cell subclusters divided by therapy group. (I) Top 10 enriched terms (from BioPlanet database) of 4-1BB⁺ upregulated DEGs in combination group, sorted by q value, adjusted via Benjamini-Hochberg method. The red dashed line indicates adjusted p value of 0.05 p value significance: ∗p < 0.05, ∗∗p < 0.01, ∗∗∗p < 0.001, and ∗∗∗∗p < 0.0001; ns, non-significant.