Article Text

Original research
Single-cell analysis reveals cellular and molecular factors counteracting HPV-positive oropharyngeal cancer immunotherapy outcomes
  1. Junha Cha1,
  2. Da Hee Kim2,
  3. Gamin Kim3,
  4. Jae-Won Cho4,
  5. Euijeong Sung1,
  6. Seungbyn Baek1,
  7. Min Hee Hong3,
  8. Chang Gon Kim3,
  9. Nam Suk Sim2,
  10. Hyun Jun Hong2,
  11. Jung Eun Lee3,
  12. Martin Hemberg4,
  13. Seyeon Park5,
  14. Sun Ock Yoon6,
  15. Sang-Jun Ha5,
  16. Yoon Woo Koh2,
  17. Hye Ryun Kim3 and
  18. Insuk Lee1,7
  1. 1Department of Biotechnology, College of Life Science and Biotechnology, Yonsei University, Seoul, Republic of Korea
  2. 2Department of Otorhinolaryngology, Yonsei University College of Medicine, Seoul, Republic of Korea
  3. 3Division of Medical Oncology, Department of Internal Medicine, Yonsei Cancer Center, Yonsei University College of Medicine, Seoul, Republic of Korea
  4. 4The Gene Lay Institute of Immunology and Inflammation, Brigham and Women's Hospital, Massachusetts General Hospital, and Harvard Medical School, Boston, MA, USA
  5. 5Department of Biochemistry, College of Life Science and Biotechnology, Yonsei University, Seoul, Republic of Korea
  6. 6Department of Pathology, Yonsei University College of Medicine, Seoul, Republic of Korea
  7. 7POSTECH Biotech Center, Pohang University of Science and Technology (POSTECH), Pohang, Republic of Korea
  1. Correspondence to Professor Insuk Lee; insuklee{at}yonsei.ac.kr; Professor Hye Ryun Kim; nobelg{at}yuhs.ac; Professor Yoon Woo Koh; YWKOHENT{at}yuhs.ac; Professor Sang-Jun Ha; sjha{at}yonsei.ac.kr

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.

  • Immunotherapy
  • Head and Neck Cancer
  • co-inhibitory molecule
  • Immune Checkpoint Inhibitor
  • Tumor microenvironment - TME

Data availability statement

The scRNA-seq and scTCR-seq data generated in this study have been deposited in the Gene Expression Omnibus database (https://www.nicbi.nlm.nih.gov/geo/) under Accession No. GSE226620. The raw and processed data were deposited in the Korean Nucleotide Archive (KoNA) under Accession No. PRKJA230574. The code for analyzing and creating the figures in the manuscript is available upon reasonable request.

http://creativecommons.org/licenses/by-nc/4.0/

This is an open access article distributed in accordance with the Creative Commons Attribution Non Commercial (CC BY-NC 4.0) license, which permits others to distribute, remix, adapt, build upon this work non-commercially, and license their derivative works on different terms, provided the original work is properly cited, appropriate credit is given, any changes made indicated, and the use is non-commercial. See http://creativecommons.org/licenses/by-nc/4.0/.

Statistics from Altmetric.com

Request Permissions

If you wish to reuse any or all of this article please use the link below which will take you to the Copyright Clearance Center’s RightsLink service. You will be able to get a quick price and instant permission to reuse the content in many different ways.

WHAT IS ALREADY KNOWN ON THIS TOPIC

  • Transcriptional diversity of tumor is a significant factor associated with clinical outcomes.

  • Human papillomavirus (HPV)-associated oropharyngeal squamous cell carcinoma (OPSCC) manifests improved survival compared with HPV-negative counterpart, possibly due to intratumoral tertiary lymphoid structures.

  • KLRB1 is a candidate inhibitory receptor in CD8+ T cells promoting inhibitory functions in glioblastoma, while its expression in CD4+ T cells is associated with improved survival in head and neck squamous cell carcinoma.

WHAT THIS STUDY ADDS

  • CD161+ resident memory T cells specifically observed in virus-induced OPSCCs indicate suppression of follicular helper cells (Tfh) and B-cell interactions, and its density inversely correlates to immunotherapy response.

  • Transcriptional plasticity of malignant cells is associated with immunotherapy response regardless of HPV infection in OPSCC, and genes that correlate with HPV viral load suggest clinical significance.

  • Differentiation of CD4+ T naïve cells to Tfh distinguish HPV-positive OPSCC tumor microenvironments, implying Tfh as a critical factor in promoting enhanced clinical benefits of virus-associated OPSCCs.

HOW THIS STUDY MIGHT AFFECT RESEARCH, PRACTICE OR POLICY

  • Targeted inhibition of CD161+ resident memory T cells could upset tertiary lymphoid structure suppression and enhance the efficacy of immunotherapy in HPV-positive OPSCCs.

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.

Materials and methods

Collection of patient samples

Patients who were diagnosed with HNSCC and underwent surgical resection at Severance Hospital (Seoul, Republic of Korea) between 2018 and 2020 were prospectively enrolled in the study (online supplemental table S1). Of the 20 patients, 10 were treated with immunotherapy as part of a neoadjuvant clinical trial registered under the ClinicalTrials.gov. Fresh tumor tissues were mechanically and enzymatically dissociated using a gentleMACS dissociator (Cat# 130-093-235; Miltenyi Biotec, Gladbach Bergisch, Germany) and a Human Tumor Dissociation Kit (Miltenyi Biotec, Cat# 130-095-929), following the manufacturer’s instructions. After incubation for 1 hour at 37°C, the resuspended samples were filtered through a 70 µm MACS SmartStrainer (Cat# 130-098-462; Miltenyi Biotec) into RPMI-1640 medium (Corning, Corning, New York, USA) supplemented with 10% fetal bovine serum (Biowest, Riverside, Missouri, USA), and centrifuged at 300×g for 10 min. Tumor-infiltrating lymphocytes (TILs) were isolated by Ficoll (GE HealthCare, Chicago, Illinois, USA) density gradient centrifugation.

Supplemental material

Annotation of HPV status

To ascertain HPV positivity, p16 immunohistochemistry was performed according to established criteria. Samples were considered positive for p16 when strong and diffuse nuclear and cytoplasmic staining was observed in more than 70% of all HNSCC tumor cells; otherwise, they were recorded as negative.

Single-cell analysis of malignant and immune cells

Quality controlled single-cell data of 20 patient samples were integrated and corrected for batch effect when applicable. Cell types were annotated based on marker gene expressions and correlation to reference data sets, and largely divided into immune/stromal and malignant data set. Immune/stromal cells were analyzed for their gene expression pseudotime,7 gene velocity,8 and cell–cell interaction9 using previously published methods. Description of all bioinformatics approaches performed in this study are available in online supplemental methods.

Survival analysis

TCGA (The Cancer Genome Atlas)-HNSC data were downloaded from the GDC portal using the TCGAbiolinks R package.10 HTSeq counts were preprocessed using TCGAanlayze_Preprocessing(), with 0.6 as the cor.cut parameter. Data were subsequently normalized using TCGAanalyze_Normalization(). The preprocessed count data were normalized to sample-specific size factors calculated using DESeq2.11 70 oropharyngeal samples were divided into subsets based on the “primary site” column given in the clinical manifest token. Survival data were downloaded from TCGA Clinical Data Resource.

VirTect12 was used to annotate the HPV status of oropharyngeal cancers from TCGA, with the genecode.V.29.annotation.gtf file serving as input for the -ucsc_gene parameter, GRCh38.p12.genome serving as the -index parameter, and the virus_757.fasta file serving as an input for -index_vir, the default viral genome provided in the VirTect package. Samples containing at least one HPV16 strain were considered HPV-positive.

Results

Transcriptional 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 heterogeneous clinical outcomes of immunotherapy in HPV-positive OPSCC, we investigated oropharyngeal cancer via scRNA-seq and scTCR-seq. Cancer and intratumor immune cells are affected by viral infection, which alters the antitumor action of the latter. HPV infections provide clinical benefits in the treatment of OPSCC, although their role in immunotherapy remains controversial. While general quiescence of TIL cytotoxicity along with transcriptional heterogeneity of malignant cells may contribute to reduced immunotherapy efficacy, we found that virus-induced CD161 expression on Trm prevented antitumor activity. Specifically, the single-cell analysis suggested reduced interactions between TLS B cells and T cells in patients with HPV-positive OPSCC with high KLRB1 expression compared with Trm with low KLRB1 levels. These interactions are important for TLS formation and antitumor activity. Consistently, tumors with KLRB1-high Trm showed low expression of CLEC2D which encodes a ligand of CD161 in TLS B cells, resulting in reduced TLS activity. Furthermore, we observed a higher probability of transition from stem-like Trm to effector-like Trm in KLRB1-low Trm than in KLRB1-high Trm in HPV-positive OPSCC. Together, these results suggest that CD161 on Trm counteracts antitumor activity. Immunofluorescence further confirmed the association between CD161+ Trm and poor immunotherapy outcomes. Therefore, we propose that CD161 in Trm acts as an inhibitory immune checkpoint molecule36 40 41 and could serve as a potential new target for immunotherapy in patients with HPV-infected oropharyngeal cancer. Blocking CD161 on Trm via monoclonal antibodies bispecific for CD161 and CD103 could have clinical relevance as it might boost the antitumor effect of the exhausted oropharyngeal TME. Evidence of CD161 as a co-inhibitory receptor molecule requires further validation in virus-induced cancer before it is accepted as a novel and effective immunotherapy target.

We examined HPV-infected HNSCC tumors occurring at anatomical sites other than the oropharynx19 and found no difference in the expression of KLRB1 and CLEC2D. This implies that the effect exerted by CD161 against HPV-associated clinical benefits may be restricted to oropharyngeal cancers. We could not find adequate public data for assessing the roles of CD161+Trm in major HPV-associated cervical cancers whereby a similar function of CD161 should be investigated in the future. Similarly, no analogous roles of CD161+Trm were detected in other types of virus-associated cancers, such as hepatocellular carcinoma induced by the hepatitis B virus.42

Using bulk RNA sequencing, Wei et al.43 recently reported that CD161 was associated with prolonged survival in HPV-positive OPSCC. In contrast, our single-cell analysis of cellular heterogeneity in tumor tissues from patients with HPV-positive OPSCC revealed a dampened antitumor effect of KRLB1-high Trm. Given that CD161 is a marker for NK and MAIT cells, the expression of KLRB1 in bulk tissue may not be a surrogate for the abundance of CD161+ Trm. Interestingly, co-expression of CLEC2D and KLRB1 was associated by Wei et al43 with improved overall survival. Our results suggest that this co-expression is based on CLEC2D expression in B cells and KLRB1 expression in Tfh of tumors with KLRB1-low Trm (figure 5J).

CD161 expression is significantly higher in both CD4+ and CD8+ T cells in HPV-positive OPSCC,24 28 and CD4+ CD161+ T cells are associated with improved clinical outcomes. Consistent with this, our single-cell analysis suggested that CD4+ T cells strongly expressing KLRB1 displayed enhanced development into effector Tfh and increased interactions between CD4+ T cells and B cells, thereby promoting antitumor activity in an HPV-infected context. These results suggest that CD161 on CD4+ T cells promotes antitumor activity, whereas CD161 on CD8+ T cells inhibits it. However, our study did not resolve the underlying mechanisms by which the same gene counteracted the two different types of lymphocytes in the context of virus-associated antitumor effects. In addition, we could not demonstrate the causality of CD161 in Trm as an inhibitor of CLEC2D in B cells, although the exclusive distribution of CD161+ Trm and CLEC2D+ B cells across patients strongly suggests an inhibitory relationship. Functional analysis using blockade antibodies specific for CD161 on Trm will provide further insights into the causality and mechanisms underlying the regulation of antitumor activity.

In summary, our research not only provides an insight into the immune landscape of HPV-driven OPSCC but also sheds light on a special subset of Trm with prognostic and therapeutic potential. We provide compelling evidence of CD161+ Trm playing a significant role in counteracting the efficacy of immunotherapy in HPV-positive OPSCCs. Through our research, we have identified that the antitumor effect is attenuated due to reduced TLS activity and impaired T-cell effector functions, which are associated with elevated expression of KLRB1, the gene encoding CD161, in Trm cells. Notably, the density of CD161+ Trm within the tumor exhibits an inverse correlation with the immunotherapy outcome in HPV-positive OPSCC. Moreover, single-cell analysis revealed that the transcriptional diversity of malignant cells also contributes significantly to the negation of immunotherapy response. These findings shed light on crucial factors that hinder the effectiveness of immunotherapy in HPV-positive OPSCCs and will facilitate the developing novel strategies to enhance immunotherapeutic approaches for this specific group of patients.

Data availability statement

The scRNA-seq and scTCR-seq data generated in this study have been deposited in the Gene Expression Omnibus database (https://www.nicbi.nlm.nih.gov/geo/) under Accession No. GSE226620. The raw and processed data were deposited in the Korean Nucleotide Archive (KoNA) under Accession No. PRKJA230574. The code for analyzing and creating the figures in the manuscript is available upon reasonable request.

Ethics statements

Patient consent for publication

Ethics approval

The studies were approved by the Institutional Review Board of Yonsei University Severance Hospital with IRB No 4-2018-1210 and IRB No 4-2020-0945. Written informed consent was obtained prior to enrollment and sample collection at Yonsei University Severance Hospital. The research conformed to the principles of the Helsinki Declaration. Participants gave informed consent to participate in the study before taking part.

References

Supplementary materials

  • Supplementary Data

    This web only file has been produced by the BMJ Publishing Group from an electronic file supplied by the author(s) and has not been edited for content.

Footnotes

  • JC, DHK and GK contributed equally.

  • Contributors JC, HRK and IL conceived and designed the study. JC performed the bioinformatics analysis of single-cell omics data and formulated the study hypotheses. DHK, YWK and HRK collected and annotated the clinical data. GK performed the experimental validation. J-WC, ES and SB assisted with the bioinformatics analysis. MHH, CGK, NSS, HJH and JEL contributed to clinical annotations. MH provided advice on single-cell data analysis. SP assisted with the experiments. SOY performed pathological annotations. YWK, HRK, S-JH and IL supervised the clinical, experimental, and bioinformatic analyses, respectively. All authors contributed to the writing of the manuscript.

    Guarantor author: IL.

  • Funding This study was supported by the National Research Foundation funded by the Ministry of Science and ICT (grants 2018R1A5A2025079 and 2022R1A2C1092062 to IL; HA22C001102 and 2021R1A2C2094629 to HRK), the Korean Health Technology R&D Project through the Korean Health Industry Development Institute (KHIDI) funded by the Ministry of Health & Welfare (HV23C0090), and the Technology Innovation Program (20022947) funded by the Ministry of Trade Industry & Energy. This work was supported in part by the Brain Korea 21 (BK21) FOUR Program, Korea Drug Development Fund (RS-2021-DD121415), and the Korean Foundation for Cancer Research (9 CB-2022-C-1). This work was partially supported by the Yonsei Fellow Program, funded by Lee Youn Jae and Severance Hospital Research Fund for Clinical Excellence (SHRC).

  • Competing interests None declared.

  • Provenance and peer review Not commissioned; externally peer reviewed.

  • Supplemental material This content has been supplied by the author(s). It has not been vetted by BMJ Publishing Group Limited (BMJ) and may not have been peer-reviewed. Any opinions or recommendations discussed are solely those of the author(s) and are not endorsed by BMJ. BMJ disclaims all liability and responsibility arising from any reliance placed on the content. Where the content includes any translated material, BMJ does not warrant the accuracy and reliability of the translations (including but not limited to local regulations, clinical guidelines, terminology, drug names and drug dosages), and is not responsible for any error and/or omissions arising from translation and adaptation or otherwise.