A Single-Cell Transcriptome Atlas of Epithelial Subpopulations in HPV-Positive and HPV-Negative Head and Neck Cancers

Abstract

Persistent infection with HPV causes nearly 5% of all cancers worldwide, including cervical and oropharyngeal cancers. Compared to HPV-negative (HPV−) head and neck squamous cell carcinomas (HNSCCs), HPV-positive (HPV+) HNSCCs exhibit a significantly improved treatment response; however, established treatment regimens were largely developed for HPV− disease. Effectively de-escalating therapy and optimizing treatment protocols to minimize toxicity for both HPV+ and HPV− tumors has been variably successful, in part due to the heterogeneity of cellular subpopulations. Single-cell RNA sequencing (scRNAseq) has primarily been used to define immune cell populations rather than the cell type of origin, epithelial cells. To address this, we analyzed published scRNAseq data of HPV+ and HPV− HNSCCs to distinguish epithelial tumor cell populations as a function of HPV status. We identified the transcriptome signatures, ontologies, and candidate biomarkers of newly identified epithelial subpopulations with attention to those that are shared or enriched in HPV+ or HPV− HNSCCs. We hypothesize that distinct epithelial cell populations and reprogramming in HPV− versus HPV+ HNSCC represent important components of the pro-tumor environment. These are described here as a foundation for the identification of new epithelial-cell-specific biomarkers, effectors, and candidate targets for optimizing the treatment of HNSCC.

1. Introduction

Although a causal relationship between HPV infection and head and neck squamous cell carcinoma (HNSCC) development was first reported over two decades ago [1] and a preventive vaccine is available, the incidence of HPV+ HNSCCs has continued to increase and is predicted to represent the majority of HNSCCs in the future [2,3]. This trend reflects the continued prevalence of oral HPV infection (i.e., ~7% of the US population, higher amongst men) [4] that is associated with a 50-fold increased risk for developing HPV+ HNSCC [5]. It is well established that HPV+ and HPV− HNSCCs differ in anatomical distribution, molecular drivers, patient demographics, treatment response, and clinical course, and as such, can be considered distinct cancer types that warrant distinct treatment algorithms [6]. For example, HPV-positive and HPV-negative HNSCCs tend to arise in different anatomical locations: HPV-negative HNSCCs primarily occur in the oral cavity, larynx, and hypopharynx, whereas HPV-positive HNSCCs commonly occur in the oropharynx (e.g., tonsils, base of tongue, soft palate). This distinction aligns with epidemiological data showing that approximately 90% of oral cavity HNSCCs are HPV-negative [7], while around 70% of oropharyngeal HNSCCs are HPV-positive [8]. Altogether, these differences support ongoing efforts to tailor HNSCC treatments based on patient and/or tumor characteristics. For example, window-of-opportunity studies in HPV− HNSCCs utilizing PD1 or PDL1-targeting immune therapy have shown improved outcomes in intermediate-risk HNSCC patients [9] and in locally advanced resectable HNSCC (KEYNOTE-689). However, additional studies are needed to understand why some patients will respond favorably to de-escalation and how to optimize treatment paradigms based on tumor characteristics, for example, by HPV status.

With regards to HPV+ HNSCCs, various deintensification strategies have been tested in clinical trials to maximize response while minimizing toxicity. These approaches have included (1) replacement [10,11], reduction [12,13,14], or omission of cytotoxic chemotherapy [15], (2) de-escalated adjuvant (chemo)radiotherapy following ablative surgery [16,17], (3) reduction in dose of radiation following induction chemotherapy [18,19,20,21,22], and (4) reduction in dose of radiation during definitive (chemo)radiotherapy [23]. Unfortunately, recent phase III clinical trials aimed at de-escalation approaches in HPV+ tumors (De-ESCALate HPV [10], RTOG 1016 [11]) resulted in inferior outcomes and have demonstrated the need for caution in deviation from standard of care [24]. In the wake of these studies, there is a renewed interest in predictive and prognostic biomarkers that may help guide patient selection for treatment tailoring [24,25]. For example, ongoing research on the utility of serum HPV circulating tumor DNA assays [26,27,28,29] has led to promising initial results. Other studies have correlated genomic biomarkers representative of tumor biology subtypes with patient outcomes. These have included p53 status [30], PIK3CA mutation [31], MATH score/ER-alpha expression [32], and TRAF/CYLD loss [33]. Thus, there is a clear need to better understand the biology of HPV-driven carcinogenesis and the profiles and/or reprogramming of key cell subpopulations specific to HPV+ disease.

Given that HPV exclusively infects the mucosal epithelium, and that squamous epithelial cells are the cells of origin for HNSCCs, defining the effects of infection on epithelial subpopulations is essential to advance new insights into HPV-reprogrammed host biology, candidate biomarkers, and targeted therapeutic strategies. Recent advances in single-cell technologies such as single-cell RNA sequencing (scRNAseq) have enabled the bioinformatic identification of distinct subpopulations within heterogenous tissues. With scRNAseq, the transcriptomic signature of distinct epithelial subpopulations can be identified and analyzed to identify their respective gene and pathway signatures. A previous publication [34] described scRNAseq data on HPV+ and HPV− HNSCCs to define distinct cell types, and explored immune-cell-related processes and interactions that shape the tumor immune microenvironment. Herein, we systematically analyzed the corresponding epithelial compartment to identify cellular subpopulations and ontologies enriched in HPV+ vs. HPV− HNSCCs as a foundation for the discovery of novel biomarkers and subtype-specific therapeutic targets.

2. Materials and Methods

2.1. scRNAseq Data Processing

Publicly available preprocessed single-cell RNA sequencing (scRNAseq) data were downloaded from the Gene Expression Omnibus database: accession ID GSE164690 [34] (Table S1). The R v4.1.1 library Seurat v4.4 was used for cell-type clustering and marker gene identification [35]. Each sample was normalized by SCTransform, using the glmGamPoi method and the number of RNA molecules per cell was regressed out. Samples were integrated with common anchor genes using the rPCA method to minimize sample-to-sample variation. Cell clusters were determined by the Louvain algorithm. UMAP dimension reduction was performed using the first 30 principal components. Marker genes for each cell type were calculated using the Wilcoxon rank-sum test returning only genes that are present in a minimum of 25% of the analyzed cluster. Epithelial cells were subset by isolating the CD45- cells, reintegrating the data and curating the epithelial cells of interest by using gene expression markers. Original annotations of HPV status, defined based on clinical p16 IHC testing and confirmed by detection of viral genes, were retained for all samples. Global expression of all epithelial cells originating from HPV− versus HPV+ tumors was compared to define globally differentially expressed genes (DEGs). Each scRNAseq cluster was considered a distinct epithelial subpopulation.

2.2. Ontology Analysis and Cytoscape Visualization

The transcriptomes and DEGs of select clusters were analyzed using the gProfiler web server [36] for Gene Ontology (GO) terms, which includes Gene Ontology (GO) terms, curated by the Ensembl database [37], and pathways from KEGG (https://www.genome.jp/kegg/, accessed on 1 January 2025), WikiPathways (https://www.wikipathways.org/index.php/WikiPathways, accessed on 3 January 2025), and Reactome (https://reactome.org/, accessed on 4 January 2025). Complex enriched gene sets were visualized using the Enrichment Map App [38] within Cytoscape 3 using parameters previously defined [39]. In brief, enrichment files and gene sets outputs from g:Profiler were processed using Enrichment Map and terms with a false discovery rate less than 0.01 were considered significantly enriched. Nodes represent pathways, and they are connected by edges (lines) between nodes with genes in common. Node size is dependent on the enrichment score associated with that pathway, while edges thickness is relative to the number of genes shared between nodes.

3. Results

3.1. Identification of Epithelial Cells by Joint scRNAseq Analysis of HPV+ and HPV− HNSCCs

We sought to define the epithelial subpopulations present in HNSCC tumors and their gene signatures, which have not been previously reported. To this end, published scRNAseq data for six HPV+ and nine HPV− HNSCC tumors were processed with batch correction (Figure 1A). UMAPs of all cells split by the HPV status of the parent tumor confirmed a high level of overlap (Figure 1B). Unbiased clustering revealed 20 distinct populations whose ontologies (Figure 1C) were assigned based on the corresponding distinct transcriptomic signatures visualized by heatmap (Figure 1D). The identified cellular populations included immune, epithelial, and stromal cells, demonstrating the heterogeneous cell types present in the tumors.

3.2. Identification of Global Transcriptomes and Pathways Distinguished by HPV Status

Keratin-expressing clusters were jointly re-clustered to isolate epithelial cells from other cell types (UMAP Figure 2A). The global signature of epithelial cells in HPV+ vs. HPV− HNSCCs was first analyzed (heatmap, Figure S1). The cytoscape visualization of GO biological processes upregulated in HPV− vs. HPV+ HNSCCs revealed pathways related to immune activation, peptidase activity, and processes such as adhesion, angiogenesis, and apoptosis (Figure 2B). Conversely, a similar analysis for processes upregulated in HPV+ vs. HPV− HNSCCs included pathways related to metabolism, migration, and differentiation (Figure 2C). Additionally, KEGG pathway (Figure 2D) and reactome pathway (Figure 2E) enrichment was assessed. Relative to HPV− HNSCCs, there was an upregulation in HPV+ HNSCCs in processes related to metabolism, tissue development, migration, and signal transduction. In HPV− HNSCCs, there was an upregulation of immune response, cell motility, and peptidase activity, which notably indicated a relative downregulation in immune response activation, cytokine signaling, and T-cell proliferation.

3.3. Distinct Epithelial Cell Subpopulations Are Present in HPV− vs. HPV+ HNSCC Tumors

Epithelial subpopulations were identified by unbiased clustering (UMAP Figure 3A, heatmap Figure 3B, dot plot Figure 3C. Cells from both HPV+ and HPV− HNSCCs were represented in all subpopulations. Enrichment (defined as ≥30% difference by cell number) of subpopulations in HPV+ versus HPV− HNSCCs was determined by analyzing the distribution of tumor cells across clusters (Figure 3D, arrows indicating enrichment). C1, C4, C5, C7 were found to be enriched in HPV+ HNSCCs, C2 and C3 in HPV− HNSCCs, and C0 and C6 were roughly equally represented in HPV+ and HPV− HNSCCs. C8 had too few cells to be confidently assessed. Each cluster/subpopulation was defined by a distinct transcriptomic profile, although in some subpopulations, there were differences in the signature between HPV+ and HPV− specimens (heatmap, Figure S2).

3.4. HPV+ and HPV− Epithelial Subpopulations Are Defined by Differentiation Status and Cell-Cycle Phase

We determined the differentiation status and cell-cycle phase of epithelial subpopulations as previously described [39]. To this end, we probed for markers of embryonic/basal (COL17A+/KRT5+) (Figure 4A) versus differentiated (KRT13+/DMKN+) (Figure 4B) character, and markers of S (PCNA) (Figure 4C) and G2/M (MKI67) (Figure 4D) phase based on a well-defined set of markers for each [40,41]. This allowed the predominant differentiation state and cell-cycle phase(s) of each epithelial subpopulation to be defined and provide context into major biological processes that may appear in pathway analyses. We find that cells in S phase were predominately in C5, cells in G2/M phase in C3, and cells in C6 were in a mix of S and G2/M phases (Figure 4E). When broken down by HPV status, an increase in cells in S phase was observed for HPV+ vs. HPV− HNSCCs (Figure 4F).

3.5. Genes Ontologies of Epithelial Subpopulations Enriched in HPV+ HNSCCs

The C1, C4, C5, and C7 subpopulations were proportionately increased in HPV+ HNSCCs. We have previously discussed C7, which most closely matches the previously reported HIDDEN cell signature, as being upregulated in HPV+ HNSCCs [39]. For each of these subpopulations, we sought to determine the defining processes and gene ontologies to elucidate the identity of these epithelial cell subtypes. We analyzed each of these subpopulations using Cytoscape visualization (Figure 5A, Figure 6A, Figure 7A and Figure 8A) and defined the top hits for GO biological processes (Figure 5B, Figure 6B, Figure 7B and Figure 8B), KEGG pathways (Figure 5C, Figure 6C, Figure 7C and Figure 8D), and reactome pathways (Figure 5D, Figure 6D, Figure 7D and Figure 8C). To identify candidate drivers of each subpopulation, transcription factor binding site enrichment analysis was performed through Transfac analysis (Figure 5E, Figure 6E, Figure 7E and Figure 8E). Altogether, we found that C1 was predominated by mitochondrial-related processes, e.g., cellular respiration and oxidative phosphorylation, C4 by cell migration, p53 signaling, and cancer cell metabolism, C5 by processes related to the S phase of the cell cycle such as DNA replication, and C7 by migration and tissue differentiation.

3.6. Ontologies of Epithelial Subpopulations Enriched in HPV− HNSCCs

Similarly, we examined the subpopulations proportionately increased in HPV− HNSCCs, namely, C2 and C3, using Cytoscape (Figure 9A and Figure 10A) and by defining top hits for GO biological processes (Figure 9B and Figure 10B), KEGG pathways (Figure 9C and Figure 10C), and reactome pathways (Figure 9D and Figure 10D). We also determined transcription factor binding site enrichment by Transfac analysis (Figure 9E and Figure 10E). C2 was dominated by mitochondrial-related pathways such as biosynthesis and cellular respiration, signatures shared with HPV+ enriched C1 cells. Lastly, C3 was defined by genes involved in the G2/M cell-cycle regulation, including processes of organelle and chromosome organization, nuclear division, cytokinesis, and nucleic acid metabolism.

3.7. Ontologies of Epithelial Subpopulations Equally Represented in HPV− and HPV+ HNSCCs

Two epithelial subpopulations (C0 and C6) were found to be present in HNSCC tumors regardless of HPV status. The same analysis for these subpopulations was performed to define ontologies: Cytoscape (Figure 11A and Figure 12A), GO biological processes (Figure 11B and Figure 12B), KEGG pathways (Figure 11C and Figure 12C), reactome pathways (Figure 11D and Figure 12D), and transcription factor binding site enrichment by Transfac (Figure 11E and Figure 12E). C0 processes predominately related to immune responses, biosynthesis, and cell differentiation, while C6 featured a range of pathways related to the S and G2/M phases of the cell cycle, including ATP synthesis, translation, biosynthesis, and p53 class mediator signaling.

3.8. Summary of Newly Identified Transcriptome Signatures

Signatures discovered by Cytoscape visualization, GO biological processes, KEGG pathways, and reactome pathways for each epithelial subpopulation and global expression comparison are summarized in Table 1. Potential avenues for targeting cellular subpopulations based on defining signatures were added.

4. Discussion

Ongoing deintensification efforts and the study of new therapies in the treatment of HPV+ and HPV− HNSCCs will benefit from the definition of differences within specific cellular populations and signatures. These offer opportunities for targeted treatment and/or prognostic biomarkers towards improved treatment stratification, especially with the increasing interest in neoadjuvant immunotherapy. Single-cell technologies enable such studies by discovering cellular subpopulations and their respective transcriptomes for greater granularity in our understanding of tumor- and microenvironment-specific HNSCC heterogeneity. Recent studies using scRNAseq have largely focused on differences in immune cell and stromal subpopulations and transcriptomic profiles, which may impact immunotherapy development [34,42]. For example, transcriptome profiles of tumor-infiltrating leukocytes isolated from immunotherapy naïve HPV− vs. HPV+ HNSCC tumors revealed that certain cell types, such as B cells, myeloid cells, and conventional CD4⁺ T cells, have divergent signatures by HPV status [43]. Thus, immunotherapeutic strategies targeting CD8⁺ T cells or CD4⁺ regulatory cells may be equally efficacious between HPV+ and HPV− HNSCCs since these cell types have similar signatures independent of HPV status; however, strategies targeting B cells, myeloid cells, and conventional CD4⁺ T cells could require tailored designs to match the divergent signatures between HPV+ vs. HPV− HNSCCs. These findings highlight the clinical importance of understanding divergent cell subpopulations and signatures of HNSCCs by HPV status.

Because HNSCCs originate in epithelial tissues, we, therefore, reasoned that elucidating the signatures of epithelial subpopulations enriched in HPV+ vs. HPV− tumors may identify biomarkers, as well as drivers of malignancy that might be therapeutic targets. First, the global signature of epithelial cells was defined in HPV+ and HPV− tumors. HPV+ HNSCCs harbored relatively decreased immune response activation signatures, including decreased antigen presentation and interferon signaling. This signature is consistent with numerous reported immunosuppressive activities of the high-risk HPV E5, E6, and E7 oncoproteins [44,45], including the modulation of antigen presentation [46], MHC surface expression [47], and interferon signaling [48]. Additionally, this finding in epithelial cells complements recent scRNAseq data in tumor-infiltrating immune cells, wherein conventional CD4⁺ T cells were found to differ in their differentiation trajectory by HPV status, having downregulated interferon responses and effector memory phenotypes [43]. There was also an upregulation in processes related to epithelial differentiation, consistent with the published disruption of the epithelial differentiation program by HPV towards the formation of an HPV-induced differentiation-dissonant epithelial nonconventional (HIDDEN) surface compartment [39]. Distinct differential signatures related to motility (HPV− HNSCCs) and migration (HPV+ HNSCCs) may reflect the HPV-dependent epithelial-to-mesenchymal transition (EMT) phenotypes. For example, EMT signature genes predicted worse five-year overall survival in HPV−, but not HPV+, HNSCCs [49,50], highlighting underlying mechanisms dependent on HPV. These differences may relate to distinct mechanisms whereby HPV reprograms epithelial cells and drives carcinogenesis, as the E6 and E7 oncogenes have been shown to downregulate E-Cadherin [51] and induce EMT-like processes [52,53]. Lastly, peptidase activity implicated in cancer progression was globally upregulated in epithelial HPV− but not HPV+ cells, indicating therapeutic potential for peptidase inhibitors tailored to HPV− tumors.

The particular advantage of scRNAseq technologies was next leveraged to delineate transcriptomically distinct epithelial subpopulations. Eight subpopulations were subsequently subjected to in-depth analysis. Two of these (C0 and C6) were roughly equally represented; two were enriched in HPV− HNSCCs (C2, C3), and four in HPV+ HNSCCs (C1, C4, C5, C7). First, we examined shared, and, therefore, likely HPV-independent, subpopulations. The C0 subpopulation was defined by metabolic biosynthetic processes while C6 contained cells expressing markers of the G2/M and S phases of the cell cycle, indicative of the HPV-independent presence of cells defined by metabolic and proliferative signatures. Studies of cancer cell metabolism have long highlighted the Warburg effect, unique metabolic rewiring that allows persistent aerobic glycolysis and fuels rapid biosynthesis, growth, and proliferation [54,55]. This signature was indeed observed in the shared C0 subpopulations, and also in the HPV+ enriched C4 subpopulation. In the C6 subpopulation harboring proliferative cells, signaling by p53 class mediators was observed, perhaps representing signaling downstream of the p53 pathway that permits cellular survival and uncontrolled growth. Given that the p53 tumor suppressor is the most commonly mutated gene in HPV− HNSCCs [56] and that high-risk HPV E6 causes the degradation of p53 [57,58,59], this ontology may reflect a shared functional loss of wild-type p53 followed by unrestricted cellular growth, despite the difference in p53 mutational status between HPV− and HPV+ tumors. This ontology was also shared with that of the HPV+ enriched C4 subpopulation. Altogether, standard therapeutic techniques such as chemo/radiation targeting shared proliferation and cancer cell metabolism signatures would be expected to demonstrate efficacy against these HNSCC subpopulations regardless of HPV status.

Next, we examined epithelial subpopulations relatively enriched in HPV− HNSCCs (C2 and C3). The C2 subpopulation was dominated by mitochondrial processes, such as cellular respiration and electron transport chain functions, while C3 consisted of cells in the G2/M phase of the cell cycle, featuring processes involved in organelle and chromosomal organization, nuclear division, and cytokinesis. Future functional analysis of C3 might yield therapeutic approaches targeting mitotic machineries in HNSCC specifically. This is particularly relevant in light of the recent literature elucidating mechanisms and agents to target cancer cell mitochondria as a novel therapeutic [60,61] and/or to improve the response to existing radiotherapy protocols [62]. Interestingly, a potentially analogous subpopulation for C2 was identified in the HPV+ enriched C1 population, which was also defined by mitochondrial and respiration processes. Several mechanisms of mitochondrial reprogramming by HPV have been reported [63], such as the induction of proapoptotic proteins, reactive oxygen species (ROS), and oxidative stress [64]. A comparative analysis of these subpopulations may further define mechanisms whereby HPV reprograms and dysregulates mitochondrial biology and/or ATP production. Differences and similarities between C3 and C1 may then guide similar and/or distinct mitochondrial targeting strategies for both HPV− and HPV+ HNSCCs.

Lastly, the remaining subpopulations enriched in HPV+ HNSCCs were considered. The defining ontologies of C1 (mitochondrial processes), C4 (cancer cell metabolism), and C7 (aberrant differentiation) were previously discussed, leaving the C5 subpopulation consisting of cells in the S phase of the cell cycle. This enrichment reflected a global shift in the distribution of cells defined by cell-cycle markers between HPV+ and HPV− tumors, wherein we observed a notable increase in cells in the S phase (C5) in HPV+ tumors. This finding was expected given that reprogramming of the cell cycle in high-risk HPV-infected SE by the E6 and E7 oncogenes is classically defined by dysplasia, hyperplasia, and hyperproliferation, and infected cells remain dependent on E6/E7 expression upon transformation into SCC. Thus, the enrichment of cells in the S phase of the cell cycle is anticipated for HPV+ tumors and underscores the validity of novel cell populations and their possible functions identified here based on gene expression patterns.

Altogether, this work represents the first step towards a detailed understanding of distinct epithelial subpopulations in HNSCCs, and their defining transcriptomes and associated biological processes as a function of HPV status. In-depth analysis and functional testing of differentially expressed genes, e.g., those regulating mitochondrial processes, may advance the translational goal of developing novel biomarkers and approaches for the tailored or joint suppression of HPV-driven and -independent effectors. Future analysis may permit the discovery of novel biomarkers and therapeutic targets to eliminate epithelial cells that are crucial for sustaining and promoting the development of each tumor subtype, and new treatment approaches that are shared or uniquely dependent on HPV status. A larger sample size of HPV+ tumors will be needed in the future to address the high clinical variability amongst HPV+ HNSCCs and to identify sub-classes of HPV+ tumors that can be used for stratification in de-escalation attempts.

5. Conclusions

We identify differentially expressed genes and pathways that define and distinguish HPV+ and HPV– epithelial subpopulations in HNSCC. These can now be used as a foundation for the discovery of biomarkers and, through functional laboratory testing, of effectors of HNSCC progression and therapy resistance. New treatment strategies shared by or tailored to HPV tumor status, both HPV+ and HPV−, may emerge from these data.

Author Contributions

Conceptualization, M.C.B., R.F., M.A. and S.I.W.; Data curation, M.C.B., C.M.R.-C., R.F. and M.A.; Formal analysis, M.C.B., C.M.R.-C., T.C., A.V., M.T.W. and M.A.; Funding acquisition, M.C.B., M.T.W. and S.I.W.; Investigation, M.C.B., M.A. and S.I.W.; Methodology, M.C.B., T.C. and M.A.; Resources, S.I.W.; Software, M.A.; Supervision, S.I.W.; Writing—original draft, M.C.B.; Writing—review and editing, M.C.B., A.L.T., C.Z., M.T.W., R.F., T.M.W.-D. and S.I.W. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Data are contained within the article and Supplementary Materials.

Acknowledgments

We thank all members of the Wells laboratory for reagents and insightful discussions.

Conflicts of Interest

The authors declare no conflicts of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Table

Figure 1 Joint analysis of HPV+ and HPV− HNSCC tumors with batch correction. (A) Overview of the bioinformatic pipeline to identify transcriptomically distinct epithelial subpopulations in scRNAseq data of HPV+ (n = 6) and HPV− (n = 9) HNSCCs. (B) Resulting UMAP of all cells separated by patient HPV status. (C) Cell clustering demonstrates a mix of immune, epithelial, and stromal cell types. (D) Heatmap of transcriptomically distinct clusters with select genes labeled.

Figure 2 Signatures upregulated in epithelial cells from HPV+ and HPV− HNSCCs. (A) UMAP of keratinocytes isolated from the published scRNAseq data and re-clustered, separated by HPV status to show extensive overlap. (B–D) Global differential expression analysis by HPV status. (B) Cytoscape visualization of GO biological processes upregulated in HPV− vs. HPV+ HNSCCs reveals pathways related to immune activation, peptidase activity, and processes such as adhesion, angiogenesis, and apoptosis. (C) Cytoscape visualization of GO biological processes upregulated in HPV+ HNSCCs includes pathways related to metabolism, migration, and differentiation. Top 10 KEGG pathways (D) and reactome pathways (E) upregulated in HPV− and HPV+ HNSCCs.

Figure 3 Enrichment of distinct epithelial subpopulations in HPV+ and HPV− HNSCC. (A) UMAP showing transcriptomically distinct clusters representing distinct epithelial subpopulations, separated by HPV status of HNSCC tumor. (B) Heatmap corresponding to clusters with select genes labeled. (C) Dot plot of top genes expressed per cluster. (D) Distribution of cells from HPV+ and HPV− HNSCCs across scRNAseq clusters, with those enriched ≥30% by cell counts highlighted by arrows.

Figure 4 Analysis of differentiation status and cell-cycle phase in epithelial subpopulations. Feature plots demonstrating distribution of select basal (A), differentiated (B), S phase (C), and G2/M phase (D) markers. (E) Distribution of G1, G2/M, and S-phase cells. C5 harbors predominately cells in S, C3 cells in G2/M, and C6 a mix of S and G2/M cells. (F) Distribution of epithelial cells across cell-cycle phases by HPV status.

Figure 5 Gene ontologies of the HPV+ enriched C1 HNSCC epithelial subpopulation. (A) Cytoscape visualization of GO biological processes highlight genes related to the electron transport chain and aerobic respiration. Summary of top GO biological processes (B), KEGG pathways (C), and reactome pathways (D) demonstrates that the transcriptome of C1 centers on mitochondrial processes and ATP generation. (E) Transfac analysis for transcription factor drivers of the C1 subpopulation.

Figure 6 Gene ontologies of the HPV+ enriched C4 HNSCC epithelial subpopulation. (A) Cytoscape visualization of GO biological processes highlights pathways related to migration, cytoskeleton organization, and metabolic processes. Summary of top GO biological processes (B), KEGG pathways (C), and reactome pathways (D). (E) Transfac analysis for transcription factor drivers of the C4 subpopulation.

Figure 7 Gene ontologies of the HPV+ enriched C5 HNSCC epithelial subpopulation. (A) Cytoscape visualization of GO biological processes highlights pathways related to proliferation. Summary of top GO biological processes (B), KEGG pathways (C), and reactome pathways (D) centers on cell-cycle processes, particularly S phase. (E) Transfac analysis for candidate transcription factor drivers of the C5 subpopulation.

Figure 8 Gene ontologies of the HPV+ enriched C7 HNSCC epithelial subpopulation. (A) Cytoscape visualization of GO biological processes highlights pathways related to migration, development, and metabolism. Summary of top GO biological processes (B), reactome pathways (C), and KEGG pathways (D) include keratinocyte differentiation related processes. (E) Transfac analysis discover candidate transcription factor drivers of the C7 subpopulation.

Figure 9 Gene ontologies of the HPV− enriched C2 HNSCC epithelial subpopulation. (A) Cytoscape visualization of GO biological processes highlights pathways related to biosynthesis and cellular respiration. Summary of top GO biological processes (B), KEGG pathways (C), and reactome pathways (D) centers on mitochondrial and ATP metabolic processes. (E) Transfac analysis for transcription factor drivers of the C2 subpopulation.

Figure 10 Gene ontologies of the HPV− enriched C3 HNSCC epithelial subpopulation. (A) Cytoscape visualization of GO biological processes highlights pathways related to cell cycle and proliferation. Summary of top GO biological processes (B), KEGG pathways (C), and reactome pathways (D) centers on the G2/M phase of the cell cycle. (E) Transfac analysis for transcription factor drivers of the C3 subpopulation.

Figure 11 Gene ontologies of the shared C0 HNSCC epithelial subpopulation. (A) Cytoscape visualization of GO biological processes highlights pathways related to immune response and proliferation. Summary of top GO biological processes (B), KEGG pathways (C), and reactome pathways (D) centers on antigen processing and translation. (E) Transfac analysis for transcription factor drivers of the C0 subpopulation.

Figure 12 Gene ontologies of the shared C6 HNSCC epithelial subpopulation. (A) Cytoscape visualization of GO biological processes highlights pathways related to biosynthesis and p53. Summary of top GO biological processes (B), KEGG pathways (C), and reactome pathways (D) centers on translation and biosynthesis. (E) Transfac analysis for transcription factor drivers of the C6 subpopulation.

Table 1

Summary of select signatures and potentially targetable phenotypes of cellular subpopulations identified by scRNAseq analysis.

	Enrichment	Select Signatures	Candidate Targetable Phenotypes
Global	HPV−	- Immune response - Cell motility - Peptidase activity	- Peptidases
Global	HPV+	- Biosynthesis - Epithelial differentiation - Migration - Signal transduction	- Metabolic pathways - Immune reactivation - Plasticity
C0	Shared	- Biosynthesis - Antigen processing and presentation - Cell differentiation	- - Metabolic pathways - - Immune reactivation
C6	Shared	- S and G2/M phase of cell cycle - Biosynthesis - Oxidative phosphorylation - Signaling by p53 class mediators	- Cell-cycle checkpoints - Mitochondrial processes - p53 signaling
C2	HPV−	- Cellular respiration - Electron transport chain	- Mitochondrial processes
C3	HPV−	- G2/M phase of cell cycle - Organelle and chromosome organization - Nuclear division and cytokinesis - RNA/DNA metabolism	- Cell division - Mitotic machinery - Nucleotide metabolism
C1	HPV+	- Oxidative phosphorylation - Electron transport chain	- Mitochondrial processes
C4	HPV+	- Cell migration - Protein metabolism - p53 signaling - Glycolysis/gluconeogenesis	- Invasion - p53 signaling - Cancer cell metabolism
C5	HPV+	- S phase of cell cycle - DNA replication - Chromosome and cytoskeletal organization - Translation	- DNA replication machinery - Translation machinery
C7	HPV+	- Cell migration - Epithelial differentiation	- Aberrant HPV-driven differentiation (HIDDEN cells)

Supplementary Materials

The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/v17040461/s1, Figure S1: Heatmap of genes differentially expressed in epithelial cells in HPV− vs HPV+ HNSCC tumors; Figure S2: Heatmap of genes differentially expressed in epithelial cells in HPV− vs HPV+ HNSCC tumors organized by clusters, with select clusters with distinct differences in gene profiles highlighted with arrows; Table S1: Data—scRNAseq gene lists.

References

1. Gillison, M.L.; Koch, W.M.; Capone, R.B.; Spafford, M.; Westra, W.H.; Wu, L.; Zahurak, M.L.; Daniel, R.W.; Viglione, M.; Symer, D.E. . Evidence for a Causal Association Between Human Papillomavirus and a Subset of Head and Neck Cancers. JNCI J. Natl. Cancer Inst.; 2000; 92, pp. 709-720. [DOI: https://dx.doi.org/10.1093/jnci/92.9.709] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/10793107]

2. Chaturvedi, A.K.; Engels, E.A.; Pfeiffer, R.M.; Hernandez, B.Y.; Xiao, W.; Kim, E.; Jiang, B.; Goodman, M.T.; Sibug-Saber, M.; Cozen, W. . Human papillomavirus and rising oropharyngeal cancer incidence in the United States. J. Clin. Oncol.; 2011; 29, pp. 4294-4301. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/21969503][DOI: https://dx.doi.org/10.1200/JCO.2011.36.4596]

3. Sung, H.; Ferlay, J.; Siegel, R.L.; Laversanne, M.; Soerjomataram, I.; Jemal, A.; Bray, F. Global Cancer Statistics 2020: GLOBOCAN Estimates of Incidence and Mortality Worldwide for 36 Cancers in 185 Countries. CA A Cancer J. Clin.; 2021; 71, pp. 209-249. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/33538338][DOI: https://dx.doi.org/10.3322/caac.21660]

4. Gillison, M.L.; Broutian, T.; Pickard, R.K.L.; Tong, Z.-Y.; Xiao, W.; Kahle, L.; Graubard, B.I.; Chaturvedi, A.K. Prevalence of Oral HPV Infection in the United States, 2009–2010. JAMA; 2012; 307, pp. 693-703. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/22282321][DOI: https://dx.doi.org/10.1001/jama.2012.101]

5. Gillison, M.L.; D’Souza, G.; Westra, W.; Sugar, E.; Xiao, W.; Begum, S.; Viscidi, R. Distinct Risk Factor Profiles for Human Papillomavirus Type 16–Positive and Human Papillomavirus Type 16–Negative Head and Neck Cancers. JNCI J. Natl. Cancer Inst.; 2008; 100, pp. 407-420. [DOI: https://dx.doi.org/10.1093/jnci/djn025] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/18334711]

6. Sabatini, M.E.; Chiocca, S. Human papillomavirus as a driver of head and neck cancers. Br. J. Cancer; 2020; 122, pp. 306-314. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/31708575][DOI: https://dx.doi.org/10.1038/s41416-019-0602-7]

7. Katirachi, S.K.; Grønlund, M.P.; Jakobsen, K.K.; Grønhøj, C.; von Buchwald, C. The Prevalence of HPV in Oral Cavity Squamous Cell Carcinoma. Viruses; 2023; 15, 451. [DOI: https://dx.doi.org/10.3390/v15020451] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/36851665]

8. Gooi, Z.; Chan, J.Y.K.; Fakhry, C. The epidemiology of the human papillomavirus related to oropharyngeal head and neck cancer. Laryngoscope; 2016; 126, pp. 894-900. [DOI: https://dx.doi.org/10.1002/lary.25767]

9. Wise-Draper, T.M.; Gulati, S.; Palackdharry, S.; Hinrichs, B.H.; Worden, F.P.; Old, M.O.; Dunlap, N.E.; Kaczmar, J.M.; Patil, Y.; Riaz, M.K. . Phase II Clinical Trial of Neoadjuvant and Adjuvant Pembrolizumab in Resectable Local-Regionally Advanced Head and Neck Squamous Cell Carcinoma. Clin. Cancer Res.; 2022; 28, pp. 1345-1352. [DOI: https://dx.doi.org/10.1158/1078-0432.CCR-21-3351]

10. Mehanna, H.; Robinson, M.; Hartley, A.; Kong, A.; Foran, B.; Fulton-Lieuw, T.; Dalby, M.; Mistry, P.; Sen, M.; O’Toole, L. . Radiotherapy plus cisplatin or cetuximab in low-risk human papillomavirus-positive oropharyngeal cancer (De-ESCALaTE HPV): An open-label randomised controlled phase 3 trial. Lancet; 2019; 393, pp. 51-60. [DOI: https://dx.doi.org/10.1016/S0140-6736(18)32752-1]

11. Gillison, M.L.; Trotti, A.M.; Harris, J.; Eisbruch, A.; Harari, P.M.; Adelstein, D.J.; Jordan, R.C.K.; Zhao, W.; Sturgis, E.M.; Burtness, B. . Radiotherapy plus cetuximab or cisplatin in human papillomavirus-positive oropharyngeal cancer (NRG Oncology RTOG 1016): A randomised, multicentre, non-inferiority trial. Lancet; 2019; 393, pp. 40-50. [DOI: https://dx.doi.org/10.1016/S0140-6736(18)32779-X] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/30449625]

12. Nguyen-Tan, P.F.; Zhang, Q.; Ang, K.K.; Weber, R.S.; Rosenthal, D.I.; Soulieres, D.; Kim, H.; Silverman, C.; Raben, A.; Galloway, T.J. . Randomized Phase III Trial to Test Accelerated Versus Standard Fractionation in Combination With Concurrent Cisplatin for Head and Neck Carcinomas in the Radiation Therapy Oncology Group 0129 Trial: Long-Term Report of Efficacy and Toxicity. J. Clin. Oncol.; 2014; 32, pp. 3858-3867. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/25366680][DOI: https://dx.doi.org/10.1200/JCO.2014.55.3925]

13. Chera, B.S.; Amdur, R.J.; Green, R.; Shen, C.; Gupta, G.; Tan, X.; Knowles, M.; Fried, D.; Hayes, N.; Weiss, J. . Phase II Trial of De-Intensified Chemoradiotherapy for Human Papillomavirus–Associated Oropharyngeal Squamous Cell Carcinoma. J. Clin. Oncol.; 2019; 37, pp. 2661-2669. [DOI: https://dx.doi.org/10.1200/JCO.19.01007] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/31411949]

14. Chera, B.S.; Amdur, R.J.; Tepper, J.E.; Tan, X.; Weiss, J.; Grilley-Olson, J.E.; Hayes, D.N.; Zanation, A.; Hackman, T.G.; Patel, S. . Mature results of a prospective study of deintensified chemoradiotherapy for low-risk human papillomavirus-associated oropharyngeal squamous cell carcinoma. Cancer; 2018; 124, pp. 2347-2354. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/29579339]

15. Yom, S.S.; Torres-Saavedra, P.; Caudell, J.J.; Waldron, J.N.; Gillison, M.L.; Truong, M.T.; Jordan, R.; Subramaniam, R.; Yao, M.; Chung, C. . NRG-HN002: A Randomized Phase II Trial for Patients With p16-Positive, Non-Smoking-Associated, Locoregionally Advanced Oropharyngeal Cancer. Int. J. Radiat. Oncol. Biol. Phys.; 2019; 105, pp. 684-685. [DOI: https://dx.doi.org/10.1016/j.ijrobp.2019.08.038]

16. Ferris, R.L.; Flamand, Y.; Weinstein, G.S.; Li, S.; Quon, H.; Mehra, R.; Garcia, J.J.; Chung, C.H.; Gillison, M.L.; Duvvuri, U. . Phase II Randomized Trial of Transoral Surgery and Low-Dose Intensity Modulated Radiation Therapy in Resectable p16+ Locally Advanced Oropharynx Cancer: An ECOG-ACRIN Cancer Research Group Trial (E3311). J. Clin. Oncol.; 2021; 40, pp. 138-149. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/34699271]

17. Ma, D.J.; Price, K.A.; Moore, E.J.; Patel, S.H.; Hinni, M.L.; Garcia, J.J.; Graner, D.E.; Foster, N.R.; Ginos, B.; Neben-Wittich, M. . Phase II Evaluation of Aggressive Dose De-Escalation for Adjuvant Chemoradiotherapy in Human Papillomavirus–Associated Oropharynx Squamous Cell Carcinoma. J. Clin. Oncol.; 2019; 37, pp. 1909-1918. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/31163012]

18. Fakhry, C.; Westra, W.H.; Li, S.; Cmelak, A.; Ridge, J.A.; Pinto, H.; Forastiere, A.; Gillison, M.L. Improved Survival of Patients With Human Papillomavirus–Positive Head and Neck Squamous Cell Carcinoma in a Prospective Clinical Trial. JNCI J. Natl. Cancer Inst.; 2008; 100, pp. 261-269.

19. Marur, S.; Li, S.; Cmelak, A.J.; Gillison, M.L.; Zhao, W.J.; Ferris, R.L.; Westra, W.H.; Gilbert, J.; Bauman, J.E.; Wagner, L.I. . E1308: Phase II Trial of Induction Chemotherapy Followed by Reduced-Dose Radiation and Weekly Cetuximab in Patients With HPV-Associated Resectable Squamous Cell Carcinoma of the Oropharynx—ECOG-ACRIN Cancer Research Group. J. Clin. Oncol.; 2016; 35, pp. 490-497. [DOI: https://dx.doi.org/10.1200/JCO.2016.68.3300]

20. Villaflor, V.M.; Melotek, J.M.; Karrison, T.G.; Brisson, R.J.; Blair, E.A.; Portugal, L.; De Souza, J.A.; Ginat, D.T.; Stenson, K.M.; Langerman, A. . Response-adapted volume de-escalation (RAVD) in locally advanced head and neck cancer. Ann. Oncol.; 2016; 27, pp. 908-913. [DOI: https://dx.doi.org/10.1093/annonc/mdw051]

21. Seiwert, T.Y.; Foster, C.C.; Blair, E.A.; Karrison, T.G.; Agrawal, N.; Melotek, J.M.; Portugal, L.; Brisson, R.J.; Dekker, A.; Kochanny, S. . OPTIMA: A phase II dose and volume de-escalation trial for human papillomavirus-positive oropharyngeal cancer. Ann. Oncol.; 2019; 30, pp. 297-302. [DOI: https://dx.doi.org/10.1093/annonc/mdy522] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/30481287]

22. Chen, A.M.; Felix, C.; Wang, P.-C.; Hsu, S.; Basehart, V.; Garst, J.; Beron, P.; Wong, D.; Rosove, M.H.; Rao, S. . Reduced-dose radiotherapy for human papillomavirus-associated squamous-cell carcinoma of the oropharynx: A single-arm, phase 2 study. Lancet Oncol.; 2017; 18, pp. 803-811. [DOI: https://dx.doi.org/10.1016/S1470-2045(17)30246-2] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/28434660]

23. Yom, S.S.; Harris, J.; Caudell, J.J.; Geiger, J.L.; Waldron, J.; Gillison, M.; Subramaniam, R.M.; Yao, M.; Xiao, C.; Kovalchuk, N. . nterim Futility Results of NRG-HN005, A Randomized, Phase II/III Non-Inferiority Trial for Non-Smoking p16+ Oropharyngeal Cancer Patients. Int. J. Radiat. Oncol. Biol. Phys.; 2024; 120, pp. S2-S3. [DOI: https://dx.doi.org/10.1016/j.ijrobp.2024.08.014]

24. Rosenberg, A.J.; Vokes, E.E. Optimizing Treatment De-Escalation in Head and Neck Cancer: Current and Future Perspectives. Oncologist; 2021; 26, pp. 40-48. [DOI: https://dx.doi.org/10.1634/theoncologist.2020-0303]

25. Farah, C.S. Molecular landscape of head and neck cancer and implications for therapy. Ann. Transl. Med.; 2021; 9, 915. [DOI: https://dx.doi.org/10.21037/atm-20-6264] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/34164549]

26. Routman, D.M.; Chera, B.S.; Jethwa, K.R.; Van Abel, K.; Kumar, S.; DeWees, T.A.; Garcia, J.J.; Price, D.L.; Kasperbauer, J.L.; Laack, N.N., II . Detectable HPV ctDNA in Post-Operative Oropharyngeal Squamous Cell Carcinoma Patients is Associated With Progression. Int. J. Radiat. Oncol. Biol. Phys.; 2019; 105, pp. 682-683. [DOI: https://dx.doi.org/10.1016/j.ijrobp.2019.08.033]

27. Chera, B.S.; Kumar, S.; Shen, C.; Amdur, R.; Dagan, R.; Green, R.; Goldman, E.; Weiss, J.; Grilley-Olson, J.; Patel, S. . Plasma Circulating Tumor HPV DNA for the Surveillance of Cancer Recurrence in HPV-Associated Oropharyngeal Cancer. J. Clin. Oncol.; 2020; 38, pp. 1050-1058. [DOI: https://dx.doi.org/10.1200/JCO.19.02444] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32017652]

28. Naegele, S.; Efthymiou, V.; Das, D.; Sadow, P.M.; Richmon, J.D.; Iafrate, A.J.; Faden, D.L. Detection and Monitoring of Circulating Tumor HPV DNA in HPV-Associated Sinonasal and Nasopharyngeal Cancers. JAMA Otolaryngol.–Head Neck Surg.; 2023; 149, pp. 179-181. [DOI: https://dx.doi.org/10.1001/jamaoto.2022.4107] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/36520425]

29. Yang, R.; Li, T.; Zhang, S.; Shui, C.; Ma, H.; Li, C. The effect of circulating tumor DNA on the prognosis of patients with head and neck squamous cell carcinoma: A systematic review and meta-analysis. BMC Cancer; 2024; 24, 1434. [DOI: https://dx.doi.org/10.1186/s12885-024-13116-6]

30. Carlos de Vicente, J.; Junquera Gutiérrez, L.M.; Zapatero, A.H.; Fresno Forcelledo, M.F.; Hernández-Vallejo, G.; López Arranz, J.S. Prognostic significance of p53 expression in oral squamous cell carcinoma without neck node metastases. Head Neck; 2004; 26, pp. 22-30. [DOI: https://dx.doi.org/10.1002/hed.10339]

31. Beaty, B.T.; Moon, D.H.; Shen, C.J.; Amdur, R.J.; Weiss, J.; Grilley-Olson, J.; Patel, S.; Zanation, A.; Hackman, T.G.; Thorp, B. . PIK3CA Mutation in HPV-Associated OPSCC Patients Receiving Deintensified Chemoradiation. JNCI J. Natl. Cancer Inst.; 2020; 112, pp. 855-858. [DOI: https://dx.doi.org/10.1093/jnci/djz224] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/31747025]

32. Patel, K.B.; Mroz, E.A.; Faquin, W.C.; Rocco, J.W. A combination of intra-tumor genetic heterogeneity, estrogen receptor alpha and human papillomavirus status predicts outcomes in head and neck squamous cell carcinoma following chemoradiotherapy. Oral Oncol.; 2021; 120, 105421. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/34198234]

33. Hajek, M.; Sewell, A.; Kaech, S.; Burtness, B.; Yarbrough, W.G.; Issaeva, N. TRAF3/CYLD mutations identify a distinct subset of human papillomavirus-associated head and neck squamous cell carcinoma. Cancer; 2017; 123, pp. 1778-1790. [DOI: https://dx.doi.org/10.1002/cncr.30570] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/28295222]

34. Kürten, C.H.L.; Kulkarni, A.; Cillo, A.R.; Santos, P.M.; Roble, A.K.; Onkar, S.; Reeder, C.; Lang, S.; Chen, X.; Duvvuri, U. . Investigating immune and non-immune cell interactions in head and neck tumors by single-cell RNA sequencing. Nat. Commun.; 2021; 12, 7338. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/34921143]

35. Hao, Y.; Hao, S.; Andersen-Nissen, E.; Mauck, W.M., 3rd; Zheng, S.; Butler, A.; Lee, M.J.; Wilk, A.J.; Darby, C.; Zager, M. . Integrated analysis of multimodal single-cell data. Cell; 2021; 184, pp. 3573-3587.e29. [DOI: https://dx.doi.org/10.1016/j.cell.2021.04.048] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/34062119]

36. Raudvere, U.; Kolberg, L.; Kuzmin, I.; Arak, T.; Adler, P.; Peterson, H.; Vilo, J. g:Profiler: A web server for functional enrichment analysis and conversions of gene lists (2019 update). Nucleic Acids Res.; 2019; 47, pp. W191-W198. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/31066453]

37. Cunningham, F.; Achuthan, P.; Akanni, W.; Allen, J.; Amode, M.R.; Armean, I.M.; Bennett, R.; Bhai, J.; Billis, K.; Boddu, S. . Ensembl 2019. Nucleic Acids Res.; 2019; 47, pp. D745-D751. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/30407521]

38. Merico, D.; Isserlin, R.; Stueker, O.; Emili, A.; Bader, G.D. Enrichment Map: A Network-Based Method for Gene-Set Enrichment Visualization and Interpretation. PLoS ONE; 2010; 5, e13984.

39. Bedard, M.C.; Chihanga, T.; Carlile, A.; Jackson, R.; Brusadelli, M.G.; Lee, D.; VonHandorf, A.; Rochman, M.; Dexheimer, P.J.; Chalmers, J. . Single cell transcriptomic analysis of HPV16-infected epithelium identifies a keratinocyte subpopulation implicated in cancer. Nat. Commun.; 2023; 14, 1975.

40. Kowalczyk, M.S.; Tirosh, I.; Heckl, D.; Nageswara Rao, T.; Dixit, A.; Haas, B.J.; Schneider, R.; Wagers, A.J.; Ebert, B.L.; Regev, A. Single cell RNA-seq reveals changes in cell cycle and differentiation programs upon aging of hematopoietic stem cells. Genome Res.; 2015; 25, pp. 1860-1872.

41. Tirosh, I.; Izar, B.; Prakadan, S.M.; Wadsworth, M.H.; Treacy, D.; Trombetta, J.J.; Rotem, A.; Rodman, C.; Lian, C.; Murphy, G. . Dissecting the multicellular ecosystem of metastatic melanoma by single-cell RNA-seq. Science; 2016; 352, pp. 189-196. [DOI: https://dx.doi.org/10.1126/science.aad0501] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/27124452]

42. Li, H.; Zandberg, D.P.; Kulkarni, A.; Chiosea, S.I.; Santos, P.M.; Isett, B.R.; Joy, M.; Sica, G.L.; Contrera, K.J.; Tatsuoka, C.M. . Ferris, Distinct CD8+ T cell dynamics associate with response to neoadjuvant cancer immunotherapies. Cancer Cell; 2025; ahead of print [DOI: https://dx.doi.org/10.1016/j.ccell.2025.02.026]

43. Cillo, A.R.; Kürten, C.H.L.; Tabib, T.; Qi, Z.; Onkar, S.; Wang, T.; Liu, A.; Duvvuri, U.; Kim, S.; Soose, R.J. . Immune Landscape of Viral- and Carcinogen-Driven Head and Neck Cancer. Immunity; 2020; 52, pp. 183-199.e9. [DOI: https://dx.doi.org/10.1016/j.immuni.2019.11.014] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/31924475]

44. Kondoh, N.; Mizuno-Kamiya, M. The Role of Immune Modulatory Cytokines in the Tumor Microenvironments of Head and Neck Squamous Cell Carcinomas. Cancers; 2022; 14, 2884. [DOI: https://dx.doi.org/10.3390/cancers14122884]

45. Conarty, J.P.; Wieland, A. The Tumor-Specific Immune Landscape in HPV+ Head and Neck Cancer. Viruses; 2023; 15, 1296. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/37376596][DOI: https://dx.doi.org/10.3390/v15061296]

46. de Freitas, A.C.; de Oliveira, T.H.A.; Barros, M.R.; Venuti, A. hrHPV E5 oncoprotein: Immune evasion and related immunotherapies. J. Exp. Clin. Cancer Res.; 2017; 36, 71. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/28545552][DOI: https://dx.doi.org/10.1186/s13046-017-0541-1]

47. Westrich, J.A.; Warren, C.J.; Pyeon, D. Evasion of host immune defenses by human papillomavirus. Virus Res.; 2017; 231, pp. 21-33. [DOI: https://dx.doi.org/10.1016/j.virusres.2016.11.023] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/27890631]

48. Miyauchi, S.; Kim, S.S.; Jones, R.N.; Zhang, L.; Guram, K.; Sharma, S.; Schoenberger, S.P.; Cohen, E.E.W.; Califano, J.A.; Sharabi, A.B. Human papillomavirus E5 suppresses immunity via inhibition of the immunoproteasome and STING pathway. Cell Rep.; 2023; 42, 112508. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/37171962][DOI: https://dx.doi.org/10.1016/j.celrep.2023.112508]

49. Jung, A.R.; Jung, C.-H.; Noh, J.K.; Lee, Y.C.; Eun, Y.-G. Epithelial-mesenchymal transition gene signature is associated with prognosis and tumor microenvironment in head and neck squamous cell carcinoma. Sci. Rep.; 2020; 10, 3652. [DOI: https://dx.doi.org/10.1038/s41598-020-60707-x]

50. van der Heijden, M.; Essers, P.B.M.; Verhagen, C.V.M.; Willems, S.M.; Sanders, J.; de Roest, R.H.; Vossen, D.M.; Leemans, C.R.; Verheij, M.; Brakenhoff, R.H. . Epithelial-to-mesenchymal transition is a prognostic marker for patient outcome in advanced stage HNSCC patients treated with chemoradiotherapy. Radiother. Oncol.; 2020; 147, pp. 186-194. [DOI: https://dx.doi.org/10.1016/j.radonc.2020.05.013]

51. Matthews, K.; Leong Cheng, M.; Baxter, L.; Inglis, E.; Yun, K.; Bäckström, B.T.; Doorbar, J.; Hibma, M. Depletion of Langerhans Cells in Human Papillomavirus Type 16-Infected Skin Is Associated with E6-Mediated Down Regulation of E-Cadherin. J. Virol.; 2003; 77, pp. 8378-8385. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/12857907]

52. Jung, Y.-S.; Kato, I.; Kim, H.-R.C. A novel function of HPV16-E6/E7 in epithelial–mesenchymal transition. Biochem. Biophys. Res. Commun.; 2013; 435, pp. 339-344. [DOI: https://dx.doi.org/10.1016/j.bbrc.2013.04.060] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/23628416]

53. Rezaei, M.; Mostafaei, S.; Aghaei, A.; Hosseini, N.; Darabi, H.; Nouri, M.; Etemadi, A.; Neill, A.O.; Nahand, J.S.; Mirzaei, H. . The association between HPV gene expression, inflammatory agents and cellular genes involved in EMT in lung cancer tissue. BMC Cancer; 2020; 20, 916. [DOI: https://dx.doi.org/10.1186/s12885-020-07428-6] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32972386]

54. Hsu, P.P.; Sabatini, D.M. Cancer Cell Metabolism: Warburg and Beyond. Cell; 2008; 134, pp. 703-707.

55. Liberti, M.V.; Locasale, J.W. The Warburg Effect: How Does it Benefit Cancer Cells?. Trends Biochem Sci; 2016; 41, pp. 211-218. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/26778478][DOI: https://dx.doi.org/10.1016/j.tibs.2015.12.001]

56. Nathan, C.-A.; Khandelwal, A.R.; Wolf, G.T.; Rodrigo, J.P.; Mäkitie, A.A.; Saba, N.F.; Forastiere, A.A.; Bradford, C.R.; Ferlito, A. TP53 mutations in head and neck cancer. Mol. Carcinog.; 2022; 61, pp. 385-391.

57. Thomas, M.; Pim, D.; Banks, L. The role of the E6-p53 interaction in the molecular pathogenesis of HPV. Oncogene; 1999; 18, pp. 7690-7700. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/10618709]

58. Werness, B.A.; Levine, A.J.; Howley, P.M. Association of human papillomavirus types 16 and 18 E6 proteins with p53. Science; 1990; 248, pp. 76-79.

59. Scheffner, M.; Werness, B.A.; Huibregtse, J.M.; Levine, A.J.; Howley, P.M. The E6 oncoprotein encoded by human papillomavirus types 16 and 18 promotes the degradation of p53. Cell; 1990; 63, pp. 1129-1136. [DOI: https://dx.doi.org/10.1016/0092-8674(90)90409-8]

60. Wen, S.; Zhu, D.; Huang, P. Targeting Cancer Cell Mitochondria as a Therapeutic Approach. Future Med. Chem.; 2013; 5, pp. 53-67. [DOI: https://dx.doi.org/10.4155/fmc.12.190]

61. Guo, X.; Yang, N.; Ji, W.; Zhang, H.; Dong, X.; Zhou, Z.; Li, L.; Shen, H.-M.; Yao, S.Q.; Huang, W. Mito-Bomb: Targeting Mitochondria for Cancer Therapy. Adv. Mater.; 2021; 33, 2007778. [DOI: https://dx.doi.org/10.1002/adma.202007778] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/34510563]

62. McCann, E.; O’Sullivan, J.; Marcone, S. Targeting cancer-cell mitochondria and metabolism to improve radiotherapy response. Transl. Oncol.; 2021; 14, 100905. [DOI: https://dx.doi.org/10.1016/j.tranon.2020.100905] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/33069104]

63. Cruz-Gregorio, A.; Aranda-Rivera, A.K.; Roviello, G.N.; Pedraza-Chaverri, J. Targeting Mitochondrial Therapy in the Regulation of HPV Infection and HPV-Related Cancers. Pathogens; 2023; 12, 402. [DOI: https://dx.doi.org/10.3390/pathogens12030402] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/36986324]

64. Cruz-Gregorio, A.; Aranda-Rivera, A.K.; Pedraza-Chaverri, J. Human Papillomavirus-related Cancers and Mitochondria. Virus Res.; 2020; 286, 198016. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32445871][DOI: https://dx.doi.org/10.1016/j.virusres.2020.198016]