Integrative Single-Cell RNA-Seq, Single Cell ATAC-Seq, and Spatial RNA-Seq Analysis of Heterogeneity in HPV-Related Head and Neck Squamous Cell Carcinoma

Purpose/Objective(s) Human papillomavirus (HPV) integration into the human genome is an established driver event in head and neck squamous cell carcinoma (HNSCC). Functionally, HPV oncogenic proteins inhibit tumor suppressor proteins and activate proliferation of infected epithelial cells. Despite seeming homogeneous, individual epithelial cells among a tumor may have significant cellular and molecular heterogeneity which could affect treatment outcomes. This study aimed to better characterize the single cell molecular heterogeneity within HPV-positive HNSCC. Materials/Methods Using a well-characterized HPV-positive HNSCC cell line (UM-SCC-47) we performed integrated single cell multiome RNA (scRNA-seq) and ATAC sequencing using 10x and sequencing platforms. Integration site analysis was done using multiple sequence alignment. Multimodal analysis was then performed using the Seurat package. Molecular profiling from sections of four HPV-positive primary tumors was also performed using the Visium Spatial Platform from 10X Genomics. Cell type identification analysis was performed using SpaceXR. HPV-associated phenotype enrichment was resolved relative to the HPVon signature derived from Puram et al., Nature Genetics, 2023. Results Using single cell ATAC-seq data we performed multiple sequence alignment with the transitive consistency score to demonstrate the HPV16 integration at a single cell level. This showed that there are two integration sites for HPV (3:189879025, 3:189895061). Using scRNAseq we also demonstrated that TP63 expression is positively correlated with E7 expression (R2 = 0.676). Other HPV-associated genes include KTN1 (R=0.43), ALDH1A3 (0.42), PTMA (0.41), and VMP1 (0.40). Spots containing an enriched HPVon signature also have enrichment of multiple gene sets including, MYC signaling (R=0.77), E2F transcription (R=0.65), and G2-M DNA damage checkpoint (R=0.63). Within the cell type identification analysis, after normalizing by tumor proportion greater than 50%, additional gene set analysis showed enrichment of cytoskeletal regulation and cell adhesion pathways. Conclusion Typically, HPV-associated HNSCC has better prognosis than other forms of HNSCC. Recently, there has been a trend towards treatment de-escalation for patients with the goal to maintain excellent prognosis and minimize treatment side-effects. However, our study supports the emerging hypothesis that high levels of intra-tumoral heterogeneity may be a critical feature of HPV-associated HNSCC that will fail therapy. Our multi—omics single cell analyses show there is significant heterogeneity in HPV expression which should be considered in future biomarker studies. Human papillomavirus (HPV) integration into the human genome is an established driver event in head and neck squamous cell carcinoma (HNSCC). Functionally, HPV oncogenic proteins inhibit tumor suppressor proteins and activate proliferation of infected epithelial cells. Despite seeming homogeneous, individual epithelial cells among a tumor may have significant cellular and molecular heterogeneity which could affect treatment outcomes. This study aimed to better characterize the single cell molecular heterogeneity within HPV-positive HNSCC. Using a well-characterized HPV-positive HNSCC cell line (UM-SCC-47) we performed integrated single cell multiome RNA (scRNA-seq) and ATAC sequencing using 10x and sequencing platforms. Integration site analysis was done using multiple sequence alignment. Multimodal analysis was then performed using the Seurat package. Molecular profiling from sections of four HPV-positive primary tumors was also performed using the Visium Spatial Platform from 10X Genomics. Cell type identification analysis was performed using SpaceXR. HPV-associated phenotype enrichment was resolved relative to the HPVon signature derived from Puram et al., Nature Genetics, 2023. Using single cell ATAC-seq data we performed multiple sequence alignment with the transitive consistency score to demonstrate the HPV16 integration at a single cell level. This showed that there are two integration sites for HPV (3:189879025, 3:189895061). Using scRNAseq we also demonstrated that TP63 expression is positively correlated with E7 expression (R2 = 0.676). Other HPV-associated genes include KTN1 (R=0.43), ALDH1A3 (0.42), PTMA (0.41), and VMP1 (0.40). Spots containing an enriched HPVon signature also have enrichment of multiple gene sets including, MYC signaling (R=0.77), E2F transcription (R=0.65), and G2-M DNA damage checkpoint (R=0.63). Within the cell type identification analysis, after normalizing by tumor proportion greater than 50%, additional gene set analysis showed enrichment of cytoskeletal regulation and cell adhesion pathways. Typically, HPV-associated HNSCC has better prognosis than other forms of HNSCC. Recently, there has been a trend towards treatment de-escalation for patients with the goal to maintain excellent prognosis and minimize treatment side-effects. However, our study supports the emerging hypothesis that high levels of intra-tumoral heterogeneity may be a critical feature of HPV-associated HNSCC that will fail therapy. Our multi—omics single cell analyses show there is significant heterogeneity in HPV expression which should be considered in future biomarker studies.