Background

Head and neck squamous cell carcinoma (HNSCC) is a common cancer1 that develops in the mucous membranes of the oral cavity, oropharynx, larynx or hypopharynx. Human papillomavirus (HPV) induces HNSCC, especially in the oropharynx, giving rise to oropharyngeal squamous cell carcinoma (HPV-positive OPSCC). Its clinical features differ from those of HPV-negative OPSCC, which is induced by non-viral factors (e.g., alcohol and smoking) and is characterized by a relatively high mutational burden.2

HPV status is a significant prognostic factor in OPSCC, with HPV-positive patients exhibiting more favorable overall survival outcomes.3 PD-1 (programmed cell death protein 1) inhibitors, such as pembrolizumab and nivolumab, prolong overall survival in patients with HNSCC who failed platinum-based therapies, leading to their approval as second-line or first-line (e.g., pembrolizumab) therapeutics.4–6 However, these PD-1 inhibitors have demonstrated clinical benefits in a limited number of patients4–6 (overall response rate ~10%–20%). Although HPV-positive OPSCC individuals are expected to respond favorably to PD-1 inhibitors by promoting de novo or pre-existing antitumor immune responses to viral antigens,3 subgroup analysis indirectly revealed that the harzard ratio in HPV-positive patients was comparable to that in HPV-negative patients, suggesting a similar efficacy in both groups.4 5 These findings point to the possible existence of unknown factors that counteract HPV-associated clinical benefits in oropharyngeal cancer immunotherapy.

In this study, we sought to identify the cellular and molecular components responsible for the inconsistent clinical advantages of HPV infection in oropharyngeal cancer immunotherapy by conducting single-cell RNA and T-cell receptor sequencing (scRNA-seq and scTCR-seq, respectively) on tumors from HPV-positive and HPV-negative OPSCCs. We hypothesized that distinct tumor microenvironment (TME) features would offset the favorable outcomes of HPV infection, possibly hinting at potential therapeutic targets. Single-cell transcriptome analysis suggested that the clinical benefit of HPV infection was voided by the heterogeneity of malignant cells. Specifically, intratumoral immune cells, including CD4+ follicular helper T cells (Tfh) and CD8+ resident memory T cells (Trm), were modulated by HPV infection. Moreover, Trm exhibiting elevated expression of KLRB1 (encoding the CD161 receptor) in HPV-positive patients was associated with TME features responsible for reduced antitumor activity. Additionally, immunofluorescence of tumors from HPV-positive patients treated with immunotherapy showed an inverse correlation between the density of CD161+ Trm and changes in tumor size. Collectively, these results suggest that CD161+ Trm counteracts the clinical benefits of HPV infection in OPSCC immunotherapy.

ResultsTranscriptional heterogeneity of malignant cells inversely correlates with immunotherapy outcome in OPSCC

To investigate cellular and molecular changes in the TME of OPSCC caused by HPV infection, we collected tonsil or base of tongue tumor biopsies from 14 patients with HPV-positive and 6 patients with HPV-negative OPSCC (figure 1A, online supplemental table S1) and determined the expression of the p16 gene by immunohistochemistry (Methods). Biopsies from 10 patients, who were subsequently treated with immunotherapy, were collected before the commencement of treatment. 10 samples were analyzed by 3’ end sequencing; whereas the other 10 were analyzed by 5’ end sequencing, which also enabled scTCR-seq. Initially, we employed 3’ end sequencing for our analysis but faced challenges with the cell sorting settings, which led to the inability to retrieve malignant cells from dissociated tumor samples. Subsequently, with the advent of the more recent 5’ end sequencing, we adjusted our cell sorting approach, enabling the successful isolation of malignant cells from the samples analyzed. This adjustment was pivotal in obtaining the desired cell profiles, highlighting the importance of proper set-up in cell sorting methodologies. Eight HPV-positive and two HPV-negative tumor samples collected from patients with OPSCC treated with immunotherapy as part of a neoadjuvant clinical trial were analyzed using both scRNA-seq and scTCR-seq. The immunotherapy response status was determined based on the pathological response. A responder was defined as those who showed tumor regression of at least 30%. The analyzed cells were categorized as either malignant or immune/stromal based on reference annotation transfer and marker gene expression (online supplemental methods). Malignant cells were confirmed to have abnormal copy number variation, and those from HPV-negative tumors displayed significantly more amplification and aberration compared with HPV-positive tumors, consistent with a previous report (online supplemental figure S1).13 Moreover, we observed a strong amplification of the 3q region in malignant cells from HPV-positive tumors. The genes encoded in this region include PIC3KA14 and SOX2,15 which have been implicated in HPV-positive oropharyngeal and vulvar cancers, respectively. These two genes are located in the 3q26 OncoCassette,16 suggesting their involvement in driving tumorigenesis during HPV infection.

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.

Following quality control, we obtained the transcriptional profiles of 21,271 malignant cells, not applying batch correction due to their patient-specific nature. These cells, obtained solely via 5’ end sequencing, displayed transcriptome phenotypes influenced more by biological than technical variations. Consequently, their transcriptomes were clustered according to patients (figure 1B) and were partitioned by HPV infection status in the Uniform Manifold Approximation and Projection (UMAP) space (figure 1C). This finding indicated that viral infection dictated the transcriptional phenotype of malignant cells. The partitioning of malignant cells by copy number variation profiles largely explained their dichotomous transcriptional groups according to HPV infection status. Transcriptome diversity of malignant cells has a prognostic value, with high “tumor biodiversity” being associated with worse treatment outcome.17 18 Therefore, we hypothesized that malignant cells with high transcriptional diversity would exacerbate the clinical outcomes of OPSCC. Based on the number of transcriptional clusters in malignant cells per patient (online supplemental methods) and the nearest-neighbor graph approach, we defined 25 distinct clusters of malignant cells (figure 1D). The 10 patients treated by immunotherapy displayed varying cancer transcriptome plasticity, with transcriptionally more diverse malignant cells in the non-responder (NR) group than in the responder (R) group (figure 1E). There were significantly more clusters of malignant cells in the NR group than in the R group, regardless of HPV infection status (figure 1F). As expected, we found that the HPV16 reads detected in malignant cells significantly correlated with each patient’s average expression of CDKN2A (coding p16, figure 1G, online supplemental methods) and that these reads were mostly found in malignant rather than immune cells (figure 1H). These results suggest that the heterogeneity of malignant cells may counteract the clinical advantages of immunotherapy for HPV-positive OPSCC.

Differentiation of CD4+ Tfh with antitumor activity is facilitated by HPV infection in OPSCC

We clustered 89,339 immune/stromal cells from tumor samples based on single-cell transcriptome profiles and annotated them (online supplemental methods). The data revealed a severe batch effect between sequencing chemistry version by 10x Genomics platforms (online supplemental figure S2A), which was mitigated via sample-wise normalization and regressing out cell cycle effects, followed by data integration (online supplemental methods). We identified 20 distinct immune cell clusters and annotated them based on the expression of canonical markers (figure 2A, online supplemental figure S3A, online supplemental table S4). Following batch correction, cell clusters showed a balanced distribution between HPV-positive versus HPV-negative groups and between 3’ end versus 5’ end scRNA-seq, with no severe patient bias (online supplemental figure S2B,C).

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.

By assessing the composition of immune/stromal cells, we found that the proportion of Tfh cells was significantly higher in HPV-positive than in HPV-negative tumors (p<0.05, two-sided Mann-Whitney U test). The same trend was observed for CD20+ B cells and tertiary lymphoid structure (TLS) B cells, although the difference was not statistically significant (figure 2B). An increase in Tfh, activated B cells, and TLS was previously reported in HNSCC HPV-positive compared with HPV-negative tumors, which correlated positively with prognosis.19 20 Thus, our results suggest that immune modulation by HPV infection may be similar in oropharyngeal cancer and other HNSCC subtypes. The TLS, which is composed of Trm, Tfh, and B cells, is an ectopic lymphoid organ that develops at sites of chronic inflammation including tumors. The extent of the TLS has been associated with a favorable clinical outcome in various cancers,21 suggesting that a higher proportion of intratumor Tfh and B cells contributes to favorable clinical outcomes in HPV-positive OPSCC.

We hypothesized that the higher proportion of Tfh cells in HPV-positive tumors might be due to the enhanced differentiation of Tfh from CD4+ T naïve cells following viral infection. Thus, we subclustered CD4+ T cells into four subsets: CD4+ naïve T cells (CD4 Tn), CD4+ memory T cells (CD4 Tm),22 CD4+ effector Tfh cells (Tfh eff), and CD4+ exhausted Tfh cells (Tfh ex)23 (figure 2C, online supplemental figure S3B, online supplemental table S5). Pseudotime analysis revealed two lineages from the naïve population that differentiate toward follicular helper and memory states, confirming our annotated CD4+subsets (figure 2D). KLRB1 in CD4+ T effector cells is associated with favorable clinical outcomes in HPV-positive HNSCC.24 Consistent with previous reports, we observed the upregulation of KLRB1 along the Tfh trajectory, particularly in Tfh eff. Differentially expressed gene (DEG) analysis of Tfh eff (figure 2E, online supplemental table S7) revealed the upregulation of genes involved in both effector and exhaustion functions. In humans, KLRB1 inhibits the cytotoxicity of natural killer (NK) cells.25 Here, we observed high expression of KLRB1 in NK and Tfh cells (figure 2F,G). CLEC2D (LLT1), which encodes a ligand of CD161, was expressed in diverse immune cell types. In CD20+ B and TLS B cells, the expression of CLEC2D was restricted to HPV-positive tumors (figure 2H,I). In addition, KLRB1 was upregulated in the Trm of HPV-positive tumors (figure 2G, online supplemental table S8). Together, these results suggest that CD161 receptors on T cells play roles in Trm-B-Tfh cooperation26 to promote TLS activity in HPV-positive OPSCC, which translates into a stronger antitumor response and favorable clinical outcome.20

To determine differences in cellular differentiation or homeostasis of the identified cell types, we applied RNA velocity analysis. HPV-positive tumors presented higher RNA velocity from CD4 Tn toward Tfh (figure 2J). By comparing the transition probability between cell states (online supplemental methods), we found that CD4 Tn had a lower probability of persistence, but a higher probability of transitioning to Tfh eff in HPV-positive than in HPV-negative tumors (figure 2K). Given that Tfh are associated with antitumor activity, a higher probability of transition into Tfh eff may contribute to HPV-associated clinical benefits in OPSCC.

Differentiation of CD8+ Trm for antitumor activity is facilitated by HPV infection in OPSCC

To investigate the modulation of tumor infiltrating CD8+ T cells by HPV infection and their impact on clinical outcomes, we selected three CD8+ T-cell subsets from the 20 immune cell subsets (figure 3A) and subclustered them into 4 new subsets representing different cellular states: stem-like CD8+ T cells (CD8 Tstem), effector-like CD8+ T cells (CD8 Teff), interferon (IFN)-α-stimulated CD8+ T cells (CD8 Tstim), and exhausted CD8+ T cells (CD8 Texh) (figure 3B, online supplemental figure S3C, online supplemental table S9). ENTPD1 (encoding CD39), ITGAE (encoding CD103), GZMB, and PRF1 were among the DEGs in CD8 Texh, indicating that these cells were cytotoxic. The initially annotated Trm cells of the 20 immune cell subsets were mostly composed of CD8 Texh, with some CD8 Tstem and CD8 Teff. Given the highly immunosuppressive TME, the majority of cytotoxic Trm were exhausted, upregulating the expression of genes encoding co-inhibitory immune checkpoint molecules, such as HAVCR2, TIGIT, and LAG3. The expression of signature genes for progenitor and terminally differentiated CD8+ T cells derived from a previous study27 was assessed for individual cells and visualized in the UMAP space. As expected, we observed a localized pattern of the progenitor signature in CD8 Tstem and a terminally differentiated signature in CD8 Texh (online supplemental figure S4A). Moreover, we observed a declining progenitor score and an increasing score for terminal differentiation along the subsets, ordered from stem-like to effector-like to exhausted (figure 3C). In the HPV-positive group, we observed a significant decrease in progenitor scores for CD8 Teff and CD8 Texh, but a significant increase in the terminal differentiation scores for CD8 Tstem, CD8 Tstim, and CD8 Texh. These results imply a faster differentiation of CD8+ T cells in the HPV-positive than in the HPV-negative group. To quantify the tendency for the transition from CD8 Tstem to CD8 Texh, we performed RNA velocity analysis (online supplemental methods). Marked differences in cell state transitions were observed between HPV-positive and HPV-negative groups (figure 3D). Specifically, CD8+ T cells showed a significantly higher probability of transitioning from CD8 Tstem to CD8 Teff and then to CD8 Texh in the HPV-positive compared with the HPV-negative group (figure 3E). Given that CD8 Texh are mostly composed of terminally differentiated cytotoxic Trm cells, this result suggests that HPV infection induces favorable clinical outcomes by promoting differentiation into cytotoxic Trm in OPSCC.

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.

Next, we compared the clonal expansion of T cells between HPV-positive and HPV-negative groups. As expected, we observed a large proportion of clonally expanded CD8+ Teff and CD8+ Texh cells (figure 3F). Moreover, CD8+ T cells in all four subsets showed greater clonal expansion in the HPV-positive group (figure 3G, online supplemental methods). T-cell receptors specific for virus antigens may partly explain the larger population of clonally expanded T cells. Given the enhanced clinical benefits provided by HPV-specific CD8+ T cells in HPV-positive OPSCCs,28 it is reasonable to suggest that T cells, once infiltrated and are activated by HPV antigens on cancer cells, enhance antitumor activities in HPV-positive OPSCC compared with the negative counterpart. Consistent with this hypothesis, clonally expanded Trm showed significantly higher signature scores for PD-L1 (programmed death-ligand 1) signaling, IFN-γ, major histocompatibility complex class 2, and T cell-mediated immunity in the HPV-positive group (figure 3H), which may explain the stronger antitumor activity of these cells and, ultimately, better outcomes.

Better clinical outcomes due to HPV infection may also be due to changes in cell–cell interactions, which are important for the antitumor activity of immune cells. For example, B-Tfh and B-Trm interactions are critical for TLS formation and activity. Assuming that the number of detected ligand-receptor pairs translated to the probability of interaction between cells, we estimated the number of significant ligand-receptor pairs (online supplemental methods). We observed a dramatic increase in the interactions between B-cell subsets (TLS B and CD20+ B cells) and T-cell subsets (Tfh, CD8 Tstem, CD8 Teff, and CD8 Texh) (figure 3I, online supplemental table S10). These results suggest that the favorable outcomes of HPV-positive OPSCC can be partly explained by stronger interactions of B cells with Tfh and Trm, which are important for TLS formation and activity.

CD161 of Trm can be induced and inhibit antitumor activity in HPV-positive OPSCC

To identify the key molecules in CD8+ T cells responsible for modulating HPV-associated clinical outcomes, we first performed DEG analysis between HPV-positive and HPV-negative groups for each cell subset. KLRB1, a well-known marker of NK and mucosal-associated invariant T (MAIT) cells, was upregulated in three major CD8+ T-cell subsets: CD8 Tstem, CD8 Teff, and CD8 Texh (figure 4A,B, online supplemental figure S4A,B, online supplemental table S11). Induction of CD161 (also known as NKR-P1) on T cells has been reported in virus-infected mice, where it might play a role in antiviral immunity.29–31 However, the function of CD161+ T cells in tumors remains controversial, as there is disagreement about whether it has a co-inhibitory32 33 or co-stimulatory34 35 antitumor effect. A recent study on brain tumors suggested that CD161 expressed on CD8+ T cells was a co-inhibitory receptor, with CD161 blockade enhancing T cell-mediated killing of glioma cells in vitro and antitumor functions in vivo.36 In our cohort KLRB1 was found to be upregulated specifically in the HPV-positive group in all CD8 subsets, most prominent in the exhausted state.

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.

Immune checkpoint molecules generally show changes in expression across cell subsets. To investigate the dynamics of gene expression during CD8+ T-cell differentiation, we performed trajectory analysis using Slingshot7 and examined the proportional changes in each CD8+ T-cell subset along pseudotime (figure 4C). Then, we also examined changes in the expression of genes involved in the regulation of cytotoxic CD8+ T cells, such as co-inhibitory (PDCD1, LAG3, HAVCR2, TIGIT) and co-stimulatory (CD27, CD82, TNFRSF4, TNFRSF9, CXCR3) immune checkpoint molecules, effector molecules (GZMA, GZMB, PRF1, IFNG, FASLG), regulators of CD8+T cell stemness (TCF7, IL7R, STAT4, FOXP1, GPR183) and exhaustion (TOX, ENTPD1, BATF, PRDM1), and cell signaling factors (CCL3, CCL5, CXCL13) along pseudotime (figure 4D). As expected, co-inhibitory immune checkpoint molecules and factors promoting T-cell exhaustion showed increased expression during the transition from CD8 Tstem to CD8 Teff to CD8 Texh. Interestingly, the potential immune checkpoint molecule KLRB1 also showed increased expression during pseudotime, but only in the HPV-positive group (online supplemental figure S4C). Together with the upregulation of KLRB1 in Trm and expression dynamics of CD8+T subsets along pseudotime, our findings suggest that CD161 expression of Trm is induced by HPV infection.

Notably, KLRB1 expression in Trm varied among HPV-positive OPSCCs. KLRB1-high tumors (n=7) were identified by selecting the upper half of HPV-positive tumors sorted by the median expression of KLRB1 in Trm. Along the pseudotime of CD8+ T cells, KLRB1-low Trm showed a higher expression of genes involved in antiviral activity (IRF1, IRF7, OASL, IFI27, IFI6) and other possible IFN-related genes (IFIT1 and ISG20) that have not yet been reported in the context of HPV infection (figure 4E). These results suggest that the CD161 level in Trm dictates the immunogenicity of the TME, partly explaining the favorable prognostic effect in HPV-positive OPSCC.

To investigate whether KLRB1 expression induced by HPV infection altered the differentiation into cytotoxic Trm in tumors, we compared RNA velocity between KLRB1-high and KLRB1-low groups of CD8+ T cells from HPV-positive patients (figure 4F). Among HPV-positive patients, KLRB1-low group showed a significantly higher probability of transition from CD8 Tstem to CD8 Teff and CD8 Texh than the KLRB1-high group (figure 4G). Given that CD8 Texh include mostly terminally differentiated cytotoxic Trm, this result suggests that in HPV-positive OPSCC, the CD161 receptor inhibits the differentiation of cytotoxic Trm, thereby reducing antitumor activity.

Next, we observed a marked reduction in the interactions of B-cell subsets with T cells in the KLRB1-high group compared with KLRB1-low group of HPV-positive OPSCC (figure 4H). Additionally, we examined the differences in immunity between KLRB1-high and KLRB1-low groups of HPV-positive OPSCC through clonal expansion using scTCR-seq data. In HPV-positive patients, a larger number of expanded clonotypes were observed in the KLRB1-low group (figure 4I), which suggests that CD161 on Trm may reduce the antitumor activity of HPV-positive OPSCC.

To obtain a mechanistic clue explaining the lower antitumor activity in CD161+ Trm, we examined the expression of KLRB1 and CLEC2D, a ligand of CD161, in three major cellular components of TLS in patients with OPSCC: Tfh, Trm, and B cells (CD20+ B and TLS B). Notably, while KLRB1 was highly expressed in the Tfh of both KLRB1-low Trm and KLRB1-high Trm tumors of HPV-positive OPSCC, expression of CLEC2D in B cells was restricted to tumors with KLRB1-low Trm (figure 4J). This suggests that HPV infection induces CD161 expression on some Trm and that these cells potentially inhibit interactions with B cells by suppressing CLEC2D expression on B cells, eventually resulting in reduced TLS activity. While KLRB1 was mainly expressed in Tfh and Trm, CLEC2D was expressed across various lymphoid cell types (figure 2H,I). Therefore, CD161 may play an important role in the regulation of antitumor activity in HPV-positive OPSCC.

CD161+ Trm density inversely correlates with immunotherapy efficacy in HPV-positive OPSCC

To validate HPV-induced KLRB1 expression in Tfh and Trm (see figure 2G) in an independent cohort of patients with oropharyngeal cancer, we selected 70 HNSCC tumor samples derived from the oropharyngeal regions in TCGA. By assessing HPV infection status using viral genome alignment as proposed by Khan et al12 (online supplemental methods), we obtained 50 HPV-positive and 20 HPV-negative oropharyngeal cancer samples from the TCGA cohort. We observed significantly higher expression correlations between Tfh markers (CD4 and CXCR537) and KLRB1 in HPV-positive tumor samples than in HPV-negative ones (figure 5A), tested using the cocor package.38 Similarly, we observed a stronger correlation between Trm markers (CD8 and ITGAE) and KLRB1 in HPV-positive tumors than in HPV-negative ones (figure 5B). These results support the idea that HPV induces KLRB1 expression in intratumoral Tfh and Trm in oropharyngeal cancer. In addition, immunofluorescence of archival tissue samples from 13 patients with HPV-positive OPSCC, of which 8 had been treated with neoadjuvant immunotherapy (online supplemental methods, online supplemental table S12), revealed colocalization of CD161+, CD103+, and CD8+cells (figure 5C), further confirming the expression of CD161 on Trm in HPV-positive OPSCC.

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/mm2 (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.

A strong correlation was found also between Tfh markers and the B-cell marker CD79A in HPV-positive tumors compared with HPV-negative ones (figure 5D). A similar yet not as significant correlation was observed between Trm and B-cell markers (figure 5E). These results suggest that CD161+ Tfh and CD161+ Trm in HPV-positive tumors have a higher probability of TLS formation by recruiting B cells via CD161-CLEC2D interactions, resulting in enhanced antitumor activity.

Gene expression analysis showed that HPV-positive tumors with KLRB1-high Trm included TLS B cells with low-level expression of CLEC2D, suggesting reduced TLS activity. Given the critical role of TLS in the response to immune checkpoint blockade treatment,39 we hypothesized that CD161+ Trm negatively affected immunotherapy efficacy. To validate the effect of CD161+ Trm on immunotherapy, we tested the relationship between the intratumoral density of CD161+Trm (marked by CD161+CD103+ CD8+ cells) and changes in tumor size in eight patients with HPV-positive OPSCC treated with neoadjuvant immunotherapy. Supporting our hypothesis, we observed a significant inverse correlation between the density of intratumoral CD161+Trm and changes in tumor size (figure 5F,G). This result suggests that the density of CD161+Trm in the TME can be potentially used to predict the efficacy of immunotherapy in HPV-positive OPSCC.

Discussion

To identify the cellular and molecular factors that explain the