Single-Cell Analyses Reveal a Functionally Heterogeneous Exhausted CD8+ T-cell Subpopulation That Is Correlated with Response to Checkpoint Therapy in Melanoma

Abstract

PD-1 pathway inhibitors have revolutionized cancer therapy. However, most patients do not durably benefit, highlighting the need for biomarkers to stratify patients as responders or nonresponders. Although CD8+ tumor-infiltrating lymphocytes (TIL) have been associated with immune checkpoint therapy response, there is no consensus on which CD8+ TIL subpopulations have the most prognostic value. Preclinical studies have focused on progenitor-like exhausted CD8+ T cells (TPEX) because TPEX proliferate more in response to PD-1 inhibitors than other exhausted T-cell (TEX) subpopulations. However, immune checkpoint inhibitor treatment drives TPEX differentiation into other TEX populations that can mediate antitumor immunity. These data complicate the ability to identify prognostically important T-cell populations in patients that predict immune checkpoint inhibitor treatment response. In this study, we found that patients with advanced melanoma with ≥20% of CD8+ TILs coexpressing PD-1 and CTLA4 (termed CPHi TIL) had better objective response rates and survival following PD-1 monotherapy than those below this threshold. Characterization of the CPHi TIL subset using bulk and single-cell RNA sequencing showed that although TPEX-like cells were present within the CPHi subset, they were in the minority of these cells. Rather, the CPHi population was numerically dominated by other subsets, including cycling, terminally exhausted-like, cytotoxic-like, and/or resident memory-like TEX populations, and a subset enriched for glycolytic genes. Collectively, these data show that CPHi TILs correlate with response in melanoma, but this TIL subset is a heterogeneous mix of different subpopulations that may differentially contribute to antitumor immunity following checkpoint blockade. Significance: The PD-1+ CTLA4+ CD8+ tumor-infiltrating lymphocyte population correlating with immunotherapy response is a heterogeneous mix of subpopulations, which has important implications for optimizing checkpoint-based immunotherapy.

Figure 1:. The frequency of CP^Hi CD8⁺ TILs correlates with response to PD-1 checkpoint blockade and survival.

Frequency of (A) CP^Hi CD8⁺ tumor-infiltrating lymphocytes (TILs co-expressing PD-1 and CTLA-4) and (B) CD4⁺ Foxp3⁺ regulatory T cells (Tregs) from pre-treatment melanoma biopsies from anti-PD-1 therapy naïve patients shown in Table 1 grouped by patient age, sex, clinical stage, and serum LDH level. (C) Frequency of CP^Hi TILs grouped by tissue biopsy site (skin metastasis versus involved lymph node). (D) Frequency of CP^Hi TILs grouped based upon treatment response. Shown are all patients, p=0.0040 (left), PD-1 monotherapy patients, p=0.0477 (middle), and combination therapy patients, p=0.1679 (right). (E) Frequency of Foxp3⁺ Tregs grouped based upon treatment response (all patients). For (A)-(E), CP^Hi TIL and Treg frequencies were assessed by flow cytometry and statistical significance was determined by the Mann-Whitney Test (**, p < 0.01; ns, not significant). (F) Progression-free survival (PFS) and (G) overall survival (OS) of patients who received anti-PD-1 monotherapy grouped by pre-treatment biopsy CP^Hi CD8⁺ TIL frequency (<20% (red) or ≥20% (blue)). Vertical dashes indicate censored data. For (F)-(G), survival statistics were assessed using the Gehan-Breslow-Wilcoxson Test and the Mantel-Haenszel method. ***, p < 0.001.

Figure 2:. Transcriptional features of CP^Hi CD8⁺ TILs in advanced melanoma.

(A) Schematic for study design incorporating bulk and single cell RNA seq data to study CP^Hi CD8⁺ TILs in advanced melanoma patients. (B)-(I) Clustering and UMAP visualization of CD8⁺ TILs sorted from advanced melanoma patients (n=14,855 cells), integrated from 8 samples from 7 patients. Details of each sample provided in Table S1. (B) Shown is the distribution of PDCD1 transcript expression (encoding PD-1) (left) and CTLA4 transcript expression (encoding CLTA-4) (right). (C) UMAP showing distribution of CD8⁺ T cells classified as “CP^Hi” or “CP^Lo” based on PDCD1 and CTLA4 transcript expression (see Methods). (D) Clustering and functional annotations of the CD8⁺ TILs from advanced melanoma patients shown in Fig. 2B. Colors denote transcriptional clusters, labeled with functional annotations (see Methods). Full list of up regulated genes per cluster can be found in Table S2. (E) Stacked bar plot showing the distribution of CP^Hi vs. CP^Lo cells across the different clusters. (F)-(H) UMAP visualization (top) and violin plots quantifying enrichment in CP^Hi vs. CP^Lo cells (bottom) of signatures derived from the literature associated with (F) terminal T cell exhaustion (left, CD8-B signature; p = 4.92×10⁻¹⁶⁶) or the progenitor T_EX (T_PEX) subset (right, CD8-G signature; p = 5.86×10⁻⁵⁷) [27], (G) human T_RM [17] (p = 1.81×10⁻⁴¹), and (H) cell cycle [56] (p = 2.73×10⁻²²). (I) UMAP visualization (top) and box plots quantifying clonal expansion levels in CP^Hi vs. CP^Lo cells (bottom) for CD8⁺ TILs that also paired TCR seq data (p = 1.2×10⁻¹⁵). P values for (F)-(I) calculated using the Wilcoxon Rank-Sum Test.

Figure 3:. Modeling identified transcriptional modules associated with glycolysis and T_RM/dysfunction that mark unique subpopulations within CP^Hi CD8⁺ TILs.

Topics from LDA modeling shown that are associated with exhaustion/effector features (p = 8.93×10⁻¹⁷⁴) (A) vs. glycolytic metabolism features (p = 2.96×10⁻⁶⁵) (B) Wilcoxon Rank Sum Test. Displayed are the top 15 genes by weight contributing to each Topic. Inset into the bar plot is a UMAP visualization of the Topic enrichment. Below each bar plot are plots of Empirical Cumulative Distribution Functions (ECDF) showing the distribution of Topic weights. A full list of Topics and the top 50 genes contributing to each Topic by weight can be found in Table S4. (C) UMAP visualization (left) and violin plots quantifying enrichment in CP^Hi vs. CP^Lo cells (right) of signatures associated with a glycolysis signatures (obtained from HALLMARKs database, http://www.gsea-msigdb.org/gsea/msigdb/cards/HALLMARK_GLYCOLYSIS.html). P values for (A)-(C) calculated using Wilcoxon Rank-Sum Test. (D) GSEA analysis quantifying enrichment of the glycolysis signature in the bulk RNA seq data set comparing CP^Hi TILs vs. CP^Lo TILs from advanced melanoma patients. NES = 1.6559837, q = 0.03973635. (E) (Top) Average module score for the T_RM signature[17] vs. the glycolysis signature (HALLMARKs) for each individual transcriptional cluster shown in Fig. 2D. (Middle) Average module score for the CD8-B signature[27] vs. the glycolysis signature (HALLMARKs) for each individual transcriptional cluster shown in Fig. 2D. (Bottom) Average module score for the Cell Cycle signature[56] vs. the glycolysis signature (HALLMARKs) for each individual transcriptional cluster shown in Fig. 2D. (F) Heat map showing Z scores for individual metabolism-related transcripts across the different transcriptional clusters shown in Fig. 2D.

Figure 4:. Distribution of CD8⁺ TILs co-expressing different combinations of immunotherapy targets based on transcript expression.

(A) UMAP showing the distribution of transcripts for coinhibitory (HAVCR2 encoding TIM-3, LAG3 encoding LAG3, TIGIT encoding TIGIT) and costimulatory (TNFRSF4 encoding OX40R, TNFRS9 encoding 4–1BB, TNFRSF18 encoding GITR, ICOS encoding ICOS, CD40LG encoding CD40L) receptors in the CD8⁺ TIL single cell dataset. (B) Bhattacharya coefficient for indicated genes. (C) Heat map showing Z scores for individual genes across the different transcriptional clusters shown in Fig. 2D. (D) Stacked bar plot showing the distribution of cells classified as “high” for the indicated combination of immunotherapy targets based on a set threshold (see Methods) across the different Seurat clusters shown in Fig. 2D.