Single-cell analysis reveals cellular and molecular factors counteracting HPV-positive oropharyngeal cancer immunotherapy outcomes

Abstract

Background: Oropharyngeal squamous cell carcinoma (OPSCC) induced by human papillomavirus (HPV-positive) is associated with better clinical outcomes than HPV-negative OPSCC. However, the clinical benefits of immunotherapy in patients with HPV-positive OPSCC remain unclear.

Methods: To identify the cellular and molecular factors that limited the benefits associated with HPV in OPSCC immunotherapy, we performed single-cell RNA (n=20) and T-cell receptor sequencing (n=10) analyses of tonsil or base of tongue tumor biopsies prior to immunotherapy. Primary findings from our single-cell analysis were confirmed through immunofluorescence experiments, and secondary validation analysis were performed via publicly available transcriptomics data sets.

Results: We found significantly higher transcriptional diversity of malignant cells among non-responders to immunotherapy, regardless of HPV infection status. We also observed a significantly larger proportion of CD4⁺ follicular helper T cells (Tfh) in HPV-positive tumors, potentially due to enhanced Tfh differentiation. Most importantly, CD8⁺ resident memory T cells (Trm) with elevated KLRB1 (encoding CD161) expression showed an association with dampened antitumor activity in patients with HPV-positive OPSCC, which may explain their heterogeneous clinical outcomes. Notably, all HPV-positive patients, whose Trm presented elevated KLRB1 levels, showed low expression of CLEC2D (encoding the CD161 ligand) in B cells, which may reduce tertiary lymphoid structure activity. Immunofluorescence of HPV-positive tumors treated with immune checkpoint blockade showed an inverse correlation between the density of CD161⁺ Trm and changes in tumor size.

Conclusions: We found that CD161⁺ Trm counteracts clinical benefits associated with HPV in OPSCC immunotherapy. This suggests that targeted inhibition of CD161 in Trm could enhance the efficacy of immunotherapy in HPV-positive oropharyngeal cancers.

Trial registration number: NCT03737968.

Keywords: Head and Neck Cancer; Immune Checkpoint Inhibitor; Immunotherapy; Tumor microenvironment - TME; co-inhibitory molecule.

Figure 1

Heterogeneity of OPSCC malignant cells and its influence on immunotherapy outcome. (A) Summary of the procedure used in this study. Both 3’ (GEX) and 5’ (GEX and TCR) 10x Genomics technology were used to analyze 20 patients with OPSCC: 6 HPV-negative (1 R and 1 NR) and 14 HPV-positive (3 R and 5 NR). (B, C) Heterogeneity of malignant cells from the 10 patients analyzed by 5’ GEX. Samples were labeled per patient (with an indication of HPV infection status by either + or –) (B) or per HPV infection status (C) and visualized in UMAP dimensions. (D) Clustering of malignant cells into 25 cell subsets according to the Louvain algorithm. (E) Proportion of clusters of malignant cells for the 10 patients divided by immunotherapy outcome. (F) Bar plot of malignant cluster count values for R and NR (p=0.04, two-sided Wilcoxon rank-sum test). (G) Scatterplot showing the average expression of p16 (CDKN2A) in relation to average HPV16 read counts of patients. The PCC and its p value are indicated. (H) Bar plot showing the count of malignant and immune cells with at least one HPV16 read. GEX, gene expression; HPV, human papillomavirus; NR, non-responder; OPSCC, oropharyngeal squamous cell carcinoma; PCC, Pearson correlation coefficient; R, responder; scRNA-seq, single-cell RNA sequencing; scTCR-seq, single-cell T-cell receptor sequencing; TCR, T-cell receptor; UMAP, Uniform Manifold Approximation and Projection.

Figure 2

Modulation of intratumor CD4⁺ T cells in OPSCC by HPV infection. (A) UMAP visualization of immune cells after normalization and batch correction. 20 immune and stromal cell subsets were clustered and annotated. (B) Composition of immune cell types for individual patients represented as grouped bar plots. The significance of compositional differences between groups was assessed using a two-sided Wilcoxon rank-sum test (*p<0.05). (C) UMAP visualization of subclusters of CD4⁺ T cells. (D) Two different lineages of CD4⁺ T cells were identified by trajectory analysis using Slingshot and alignment along pseudotime. Tfh lineage (left, Tn → Tfh eff → Tfh ex) and CD4⁺ T memory lineage (right, Tn → Tm). (E) Volcano plot showing significant DEGs (q<0.05, |log2FC|>0.25) in red. (F) UMAP of immune cells with shoring expression of KLRB1. (G) Violin plot showing KLRB1 expression across immune cell types in HPV-positive and HPV-negative groups. (H) UMAP of immune cells with shoring expression of CLEC2D. (I) Violin plot showing CLEC2D expression across immune cell types in HPV-positive and HPV-negative groups. (J) RNA velocity of differentiating CD4⁺ T cells in HPV-positive (left) and HPV-negative (right) groups. (K) Probability density plot of persistent CD4⁺ Tn (left, Tn → Tn) and transition from Tn to Tfh eff (right) in HPV-positive (left) and HPV-negative (right) groups. The median transition probability is indicated by a dotted vertical line for each group. Significant differences between probability distributions were assessed using the two-sided Kolmogorov-Smirnov test. DEGs, differentially expressed genes; HPV, human papillomavirus; OPSCC, oropharyngeal squamous cell carcinoma; Teff, effector-like T cells; Tfh, follicular helper T cells; Tm, memory T cells; Tn, naïve T cells; UMAP, Uniform Manifold Approximation and Projection.

Figure 3

Modulation of intratumor CD8⁺ T cells in OPSCC by HPV infection. (A) UMAP of CD8⁺ T cells with initial cell-type annotations. (B) UMAP of subclustered CD8⁺ T cells annotated by distinct cell states. (C) Split violin plot of GSVA scores for progenitor (left) or terminal differentiation (right) signature genes between HPV-positive and HPV-negative groups for each cell state. (D) RNA velocity of CD8⁺ T lineage cells visualized in the UMAP dimension for HPV-positive and HPV-negative groups. (E) Probability density plot of transition from stemness to effector T cells (left) and from effector to exhausted T cells (right) in HPV-positive and HPV-negative groups. The median probability is indicated by a dotted vertical line for each group. Significant differences between probability distributions were assessed using the two-sided Kolmogorov-Smirnov test. (F) Clonal expansion of CD8⁺ T cells visualized in the UMAP dimension. Cells with no expansion (n=1), expanded (2≤n <5), or hyper-expanded (5≤n) TCR clones were labeled using color codes. NA: cells with GEX barcodes but no TCR information. (G) Cumulative cell count with different clonal expansion status for each cell state normalized to the sum of all HPV-negative or HPV-positive cells per group. (H) Box plot of the GSVA score for selected MSigDB hallmark gene sets in hyperexpanded Trm. Significant differences between HPV-positive and HPV-negative groups were assessed using a two-sided Wilcoxon rank-sum test. (I) Chord diagram for selected immune cells (CD8⁺ Tstem, CD8⁺ Teff, CD8⁺ Texh, CD20⁺ B, TLS B, and CD4⁺ Tfh) and significant ligand-receptor interactions. Circumference size denotes the number of statistically significant interactions (p<0.05). ****p<0.0001, two-sided Wilcoxon rank-sum test. GSVA, gene set variation analysis; GEX, gene expression; HPV, human papillomavirus; IFN, interferon; MHC, major histocompatibility complex; OPSCC, oropharyngeal squamous cell carcinoma; PD-L1, programed death-ligand 1; TCR, T-cell receptor; Teff, effector-like T cells; Texh, exhausted T cells; Tfh, follicular helper T cells; TLS, tertiary lymphoid structure; Tstem, stem-like T cells; Trm, resident memory T cells; UMAP, Uniform Manifold Approximation and Projection.

Figure 4

Two subsets of HPV-positive Trm in OPSCC present distinct KLRB1 expression. (A, B) Volcano plot showing significant DEGs (q<0.05, |log2FC|>0.4) depicted in red, following HPV infection in CD8⁺ Texh (A) and CD8⁺ Teff (B). (C) Trajectory of Trm, from Tstem to Teff and then to Texh, identified by Slingshot and aligned along pseudotime. (D) Heatmap showing the average expression of selected genes associated with CD8⁺ T-cell immunity for the three Trm states in HPV-positive and HPV-negative groups. (E) Average expression of selected genes along the pseudotime depicted in (C) for the HPV-positive group, whose Trm presented high or low KLRB1 expression. (F) RNA velocity for HPV-positive Trm with low or high KLRB1 expression. (G) Probability density plots of transition from stem cells to effector cells (left, Tstem → Teff) and from effector cells to exhausted cells (right, Teff → Texh) in HPV-positive patients with high or low KLRB1 expression. The median probability is indicated by a dotted vertical line for each group. Significant differences between probability distributions were assessed using the two-sided Kolmogorov-Smirnov test. (H) Chord diagram for the detected ligand-receptor interaction via CellPhoneDB analysis in selected immune cell subsets (CD8⁺ Tstem, CD8⁺ Teff, CD8⁺ Texh, CD20⁺ B, TLS B, and CD4⁺ Tfh) for HPV-positive patients, whose Trm present low or high KLRB1 expression. Circumference size denotes the number of statistically significant interactions (p<0.05). (I) Normalized number of clonally expanded T cells for HPV-positive patients with high or low KLRB1 expression. (J) Violin plots showing KLRB1 expression in Tfh or Trm cells, as well as CLEC2D expression in TLS B or CD20⁺ B cells of OPSCC patients. DEGs, differentially expressed genes; HPV, human papillomavirus; OPSCC, oropharyngeal squamous cell carcinoma; Teff, effector-like T cells; Texh, exhausted T cells; Tfh, follicular helper T cells; TLS, tertiary lymphoid structure; Tstem, stem-like T cells; Trm, resident memory T cells.

Figure 5

Intratumor density of CD161⁺ Trm is associated with immunotherapy outcome for HPV-positive OPSCC. (A, B) Pearson correlation between average expression of CD4⁺ Tfh markers (CD4 and CXCR5) (A) or CD8⁺ Trm markers (CD8 and ITGAE) (B) and expression of KLRB1 in patients with HPV-negative and HPV-positive oropharynx cancer compiled from TCGA database. (C) Multiplex immunofluorescence staining for CD8 (orange), CD161 (red), CD103 (green), PD-1 (yellow), cytokeratin (white), and DAPI (blue); scale bar: 20 µm. White squares enclose CD161⁺ CD103⁺ CD8⁺ cells in the HPV-positive HPV01 sample. (D, E) Pearson correlation between average expression of CD4⁺ Tfh markers (CD4 and CXCR5) (D) or CD8⁺ Trm markers (CD8 and ITGAE) (E) and expression of the B-cell marker CD79A in patients with HPV-negative and HPV-positive oropharynx cancer compiled from TCGA database. (F) Changes in tumor size (blue) and density of CD161⁺ Trm for six patients with HPV-positive OPSCC treated by immunotherapy and with more than five CD161⁺ CD103⁺ CD8⁺ cells/mm² (pink) are depicted as a connected scatterplot. (G) Scatterplot linking changes in tumor size and the density of CD161⁺ Trm for the same six patients with OPSCC. The significance of the Spearman correlation is denoted. DAPI, 4′,6-diamidino-2-phenylindole; HPV, human papillomavirus; OPSCC, oropharyngeal squamous cell carcinoma; PD-1, programed cell death protein 1; TCGA, The cancer genome atlas; and PD-1Tfh, follicular helper T cells; Trm, resident memory T cells.