Single-cell transcriptomic analysis of HPV-related multiphenotypic sinonasal carcinoma uncovers MYB-HPV association

(192 of 200 words) Human papillomavirus (HPV)-related multiphenotypic sinonasal carcinoma (HMSC) is a rare tumor that morphologically resembles high grade adenoid cystic carcinoma (ACC), yet exhibits indolent clinical behavior. Both demonstrate MYB proto-oncogene upregulation, but HMSC lacks the MYB translocation typically seen in ACC. Transcriptional changes in HMSC tumors remain uncharacterized. We performed single-cell RNA sequencing (scRNA-seq) on a human HMSC tumor and compared expression profiles with published ACC and oropharyngeal squamous cell carcinoma (OPSCC) scRNA-seq datasets. Primary malignant cells from HMSC (n=134) and ACC (n=980) clustered separately, and HMSC lacked bicellular differentiation into luminal and myoepithelial cells, distinguishing it from ACC. A greater proportion of HMSC cells expressing HPV-related genes (HPVon) expressed MYB (83% vs. 62%, p=0.022) and MYB targets (p=6.4×10-6), suggesting an HPV-MYB association. This finding was validated in HPV-positive OPSCC, with 7/10 tumors showing MYB upregulation in HPVon versus HPVoff cells (p<0.05). A 264- gene signature from HPVon HMSC cells was also associated with worse prognosis in HPV+ OPSCC (p<0.003), suggesting an alternate role for HPV that has not been well characterized. Further validation of the HPV-MYB association and prognostically relevant HPV gene signature may improve patient stratification and therapeutic strategies in HPV-related malignancies. .CC-BY-NC 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 25, 2025. ; https://doi.org/10.1101/2025.04.21.649795doi: bioRxiv preprint 4

Human papillomavirus (HPV)-related multiphenotypic sinonasal carcinoma (HMSC) represents a subset of HPV-related carcinomas that arise within the sinonasal tract and histologically resemble salivary adenoid cystic carcinoma (ACC). First described 10 years ago,(1) HMSC remains poorly characterized, as it is an extremely rare tumor entity. HMSC exhibits a broad morphologic spectrum, reflecting features of both salivary-type and squamous cell carcinomas.(1–3) Despite its typical high grade histology, HMSC has low potential for metastasis and is rarely fatal.(1) Like ACC, myeloblastosis (MYB) genes are typically upregulated in HMSC.(4) However, ACC harbors MYB or MYBL1 chromosomal translocations in over 70% of cases,(5–7) which juxtapose super- enhancers next to the MYB locus to drive expression.(5) These translocations are notably absent in HMSC,(1) thus suggesting an alternate mechanism of MYB upregulation. Moreover, MYB may facilitate paracrine interactions between myoepithelial and luminal cells in ACC, driving oncogenic NOTCH signaling in luminal cells.(5,8) Conversely, HMSC lacks evidence of bicellular differentiation, suggesting MYB may act differently in these tumors.(2) HMSC is also unique from ACC in harboring high-risk HPV33.(1–3) HPV is known to drive oropharyngeal carcinogenesis through viral protein E6/E7-induced cell cycle entry and dysregulation of p53 and Rb. This leads to unchecked DNA damage, excessive DNA replication and genomic instability,(9) but the predominant HPV subtypes driving these oropharyngeal cancers are HPV 16 and 18.(10) Recent studies have also suggested the importance of HPV expression heterogeneity across malignant cells within HPV-positive (HPV+) oropharyngeal squamous cell carcinoma (OPSCC).(11) Specifically, a subset of OPSCC tumor cells may lose HPV expression and enter a senescent phenotype, which drives therapy resistance, invasion, and worsened prognosis.(11) In HMSC, though, the role of HPV and the impact of its potential expression heterogeneity remain unclear. Despite their differences, distinguishing HMSC from ACC remains a diagnostic challenge,(3) and some question whether this rare tumor is, in fact, a separate molecular entity. We employed single cell RNA-sequencing (scRNA-seq)(12,13) to characterize transcriptional heterogeneity in HMSC with an eye .CC-BY-NC 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 25, 2025. ; https://doi.org/10.1101/2025.04.21.649795doi: bioRxiv preprint 5 towards understanding MYB expression and the role of HPV in these tumors, while providing a comparison with salivary ACC and HPV-positive oropharyngeal SCC. This work represents the first detailed transcriptional characterization of intra-tumoral heterogeneity in this rare disease entity.

Human tumor sample. This study received approval from the institutional review board. Prior to surgery, a patient at Massachusetts Eye and Ear with HMSC provided informed consent to participate in this project. A fresh tumor biopsy was obtained during the procedure, and the diagnosis of HMSC was confirmed by pathology. Tumor dissociation, cell sorting, and SMART-Seq2. Tumor dissociation, cell sorting, library preparation, and sequencing were performed as previously described.(8) Briefly, the fresh tumor sample was mechanically and enzymatically dissociated to a single cell suspension using a Human Tumor Dissociation Kit (Miltenyi Biotec, Bergisch Gladbach, Germany). The cell suspension was filtered and resuspended, and cell viability was confirmed to be > 90% using trypan blue. Before fluorescence- activated cell sorting (FACS), cells were stained with 1 μM calcein AM (ThermoFisher Scientific, Waltham, MA, USA) and 0.22 μM TO-PRO-3 iodide (ThermoFisher Scientific) for viability and CD45- vioblue (Miltenyi Biotec) for immune cell depletion. Viable, singlet, CD45- cells were FACS sorted using 488 nm (calcein AM, 530/30 filter), 640 nm (TO-PRO-3, 670/14 filter), and 405 nm (Vioblue, 450/50 filter) lasers, standard forward scatter height versus area criteria, and calceinhigh/TO-PROlow gates. Cells were sorted into 96-well plates containing TCL-buffer (QIAGEN, Hilden, Germany) with 1% β- mercaptoethanol and snap frozen prior to library prep. Full length cDNA libraries were generated using a modified SMART-seq2 protocol.(12–14) RNA purification was performed with Agencourt RNAClean XP beads (Beckman Coulter, Brea, CA, USA), reverse transcription with Superscript II or Maxima reverse transcriptase (ThermoFisher Scientific), and .CC-BY-NC 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 25, 2025. ; https://doi.org/10.1101/2025.04.21.649795doi: bioRxiv preprint 6 whole transcriptome amplification using KAPA HiFi HotStart ReadyMix (KAPA Biosystems, Wilmington, MA, USA). Tagmentation of cDNA libraries was performed with the Nextera XT Library Prep Kit (Illumina, San Diego, CA, USA) and sequencing was conducted as paired-end 38-base reads on a NextSeq 500 (Illumina). Processing of HMSC and ACC data and removal of non-malignant cells. HMSC data were aligned to the GRCh38.d1.vd1 reference with STAR version 2.5.2(15) and counted with featureCounts.(16) scRNA- seq data of ACC cells were taken from Parikh et al.(8) Cells with less than 2 000 transcripts with at least 1 read, or more than 40% mitochondrial reads were removed from the analysis. Genes that were expressed in less than 3 cells across the combined HMSC-ACC cohort were also removed. Log2(TPM+1) values were calculated and used for downstream analysis. We then used the Seurat 4.1.0(17) pipeline to perform dimensional reduction and clustering. We first scaled the data using the ScaleData function with the 15 000 most variable genes across the cohort, regressing out the number of detected reads and the percent of mitochondrial reads. Principal component analysis (PCA) was calculated with the same genes, and uniform manifold approximation and projection UMAP was done using the first 30 principal components (PC). Cells were grouped using the Louvain algorithm in Seurat, considering the first 30 principal components and setting the resolution to 1. Each cluster was classified as malignant or non-malignant based on ACC score, defined as mean expression of known ACC markers,(18) MYB expression, and expression of non-malignant markers - COL3A1 and COL1A2 as CAFs markers; CD79A and PTPRC as white blood cells markers; and VWF and PECAM1 as endothelial markers. All cells in the non-malignant clusters were removed from further analysis. Quality control and batch effect correction of HMSC data. To remove potential artifacts, we focused only on HMSC malignant cells with at least 2 500 transcripts. Potential batch effects between different plates were corrected with Seurat integration: 1 000 integration features were selected based on the 1 000 .CC-BY-NC 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 25, 2025. ; https://doi.org/10.1101/2025.04.21.649795doi: bioRxiv preprint 7 most variable features using the SelectIntegrationFeatures function. Integration anchors were found using the FindIntegrationAnchors function with k.filter = 50, and data was integrated using the IntegrateData function considering the 50 nearest neighbors using k.weight = 50. Integrated data from HMSC malignant cells were then rescaled using the ScaleSeurat function with default arguments, and dimensionality reduction and clustering were performed using the Seurat pipeline, considering the first 10 principal components and a resolution of 0.5. Scoring for canonical pathways in HMSC and ACC. Mean expression of the genes JAG1, JAG2, and DLL1 were used for the NOTCH ligands score. For Myoepithelial score, Luminal score, and Notch targets score, we used the mean expression of genes as used in Parikh et al.(8) We used MYB targets as identified in Drier et al.(5) A cell was classified as MYB+ if MYB expression was > 1 TPM. An HMSC cell was classified as HPV+ if the number of reads mapped to the HPV33 genome was greater than 10. Cell cycle score was calculated using the Z-score mean of genes in GO mitotic cell cycle (GO:0000278). A cell was considered cycling if it had a cell cycle score of more than 0.3. Analysis of HMSC heterogeneity. cNMF (Consensus Non-negative Matrix Factorization on single-cell RNA-Seq data)(19) version 1.5.0 python package with python 3.8.15 was used to determine expression programs from HMSC log2(TPM+1) values. We used the preprocess_for_cnmf function for batch correction specifying the plate variable under the “harmony_vars” parameter. We also set “max_scaled_thresh = 50” and “theta = 2”. Then we ran cNMF for 3 to 10 programs, and 20 iterations per run. For differential analysis between HPVon and HPVoff HMSC cells, we used log2(TPM+1) values of the 3 000 most highly expressed genes in HMSC (of the 15 000 most variable genes in HMSC cells) and corrected for plate-based batch effects using linear regression with the argument test.use = "LR". .CC-BY-NC 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 25, 2025. ; https://doi.org/10.1101/2025.04.21.649795doi: bioRxiv preprint 8 Gene set enrichment analysis (GSEA) was calculated with the fgsea 1.20.0 R package,(20) and visualized using hypeR version 2.0.1.(21) Gene symbols with average log2 fold change, sorted in decreasing order, were used as inputs for GSEA. MSigDB Hallmarks V7.0(22) and GO mitotic cell cycle (GO:0000278) were used in the enrichment analysis. Comparison between HMSC and ACC. The corrected HMSC expression data was then combined with the primary ACC cells with log2(TPM+1) values above. The combined HMSC-ACC data was rescaled with Seurat ScaleData regressing out the number of detected reads and percent of mitochondrial reads, using the 15 000 most variable genes of the combined dataset. PCA was calculated using the same number of variable genes and UMAP was calculated using the first 20 PC. To identify distinct cell populations within the dataset, we performed a clustering analysis using the Louvain algorithm with a resolution of 0.01. For differential expression analysis, we initially used the log2(TPM+1) values of the 15 000 most variable genes across the combined dataset as an input to FindMarkers. Subsequently, we identified cycling cells as described above, removed them and then re-calculated the top 15 000 most variable genes. Log2(TPM+1) values of these genes across combined dataset of non-cycling cells were used as the input for FindMarkers, with the same parameters as before filtering. Copy number variation estimation of HMSC and ACC. Infercnv version 1.10.1(23) was used to estimate copy number variation from integrated normalized expression values of the combined HMSC–ACC (primary and non-primary) cohort. Since the data were already normalized, step 3 (normalization by sequencing depth) and step 4 (log transformation of data) were skipped. Denoising was done with default Infercnv settings. Copy number aneuploidy score was defined as the square root of the average of deviation of diploidy squared. .CC-BY-NC 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 25, 2025. ; https://doi.org/10.1101/2025.04.21.649795doi: bioRxiv preprint 9 OPSCC single cell analysis. Raw data was downloaded from Gene Expression Omnibus with GEO accession GSE182227. Cells with less than 1 000 expressed genes and genes that were expressed in less than 3 cells were filtered out. Data was normalized using the Seurat NormalizeData function with log normalization and a scale factor of 10 000. For downstream analysis, we only used HPV+ tumors. A cell was identified as HPVon if it had a count of 3 or more for any HPV33 gene. For differential analysis between HPVon and HPVoff OPSCC cells, we used the Seurat pipeline with the 5 000 most highly expressed genes (of the 15 000 most variable genes) and log2(TPM+1) values were used along with linear regression for patient identity to account for inter-patient variability. Statistical analysis and data visualization. Data analysis was performed in R (version 4.4.1)(24) and ggplot2 version 3.5.1(25) was used to generate figures. To generate heatmaps, we used ComplexHeatmap version 2.10.0,(26) and we generated stacked barplots with the ggbarstats function from ggstatsplot package, version 0.12.5.(27) We calculated Fisher tests using fisher_test and T-tests with t_test functions from the rstatix 0.7.0 R package(28) with default arguments. False discovery rate (FDR) was calculated using the p.adjust function in stats package. TCGA data analysis. Transcriptomic and clinical data were downloaded using TCGAbiolinks(29) with Data Release 42.0, while patients’ follow-up data were downloaded from the GDC portal. RNA-seq data was log2 normalized and only samples originating from one of the following were included in the analysis: "Base of tongue, NOS", "Oropharynx, NOS", "Posterior wall of oropharynx", "Tonsil, NOS". HPV status was taken from cBioPortal(30) and the TCGA pan cancer atlas 2018.(31) For survival analyses, we split patients by median expression, and the signature was calculated with the mean of log2(TPM+1). Analysis and visualization were done using the TCGAanalyze_survival function from TCGAbiolinks. Creating the signature without cell cycle genes was done by excluding genes that were annotated as GO mitotic cell cycle (GO:0000278). .CC-BY-NC 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 25, 2025. ; https://doi.org/10.1101/2025.04.21.649795doi: bioRxiv preprint 10

HMSC is genomically and transcriptionally distinct from ACC To profile intratumoral heterogeneity in HMSC, we harvested a fresh tumor biopsy from a 56-year-old woman who underwent a left partial maxillectomy for recurrent poorly differentiated basaloid carcinoma. Pathology revealed solid sheets of cells with focal ductal differentiation (Figure 1A) and multifocal positivity for both P16 (Figure 1B) and MYB (Figure 1C). HPV PCR was positive for HPV type 33, supporting the diagnosis of HMSC and distinguishing the lesion from high grade, solid-type ACC. The freshly resected specimen was dissociated into a single-cell suspension (see Methods),(8) flow-sorted to deplete nonviable and CD45-positive cells, and profiled by the SMART-seq2 protocol,(12) which provided full-length transcripts for sequencing and analysis. Transcriptomes from a total of 155 HMSC single cells were retained after initial quality control (see Methods). We first sought to compare HMSC with ACC and combined our dataset with our previously published single cell transcriptomic analysis, which resulted in 980 primary malignant cells from seven head and neck ACC tumors after quality control.(8) Expression analysis of known cell markers (see Methods) revealed that 143/155 cells in the HMSC dataset were malignant cells, with 134/143 cells that passed additional quality control and batch correction (see Methods). Moreover, analysis of inferred copy number variations (CNV; see Methods) across the two datasets confirmed the distinction between malignant cells and non-malignant cells by demonstrating more CNVs for malignant cells. In addition, this analysis also revealed that malignant HMSC cells harbored a significantly higher degree of CNVs compared to malignant ACC cells (Figure 2A, S1A). Specifically, HMSC cells exhibited losses within chromosomes 4, 6, 9, 12-15, and 17, and gains within chromosomes 3, 18-20, and 22 (Figure 2A). HMSC cells clustered separately from ACC cells based on global gene expression (Figure 2B, S1B-C), supporting the hypothesis that these represent distinct cellular entities. We then scored cells based on MYB expression and found that malignant HMSC cells and ACC cells demonstrated high MYB expression (Figure 2C), validating previous studies that demonstrated .CC-BY-NC 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 25, 2025. ; https://doi.org/10.1101/2025.04.21.649795doi: bioRxiv preprint Figure 1 B C 100 μm 400 μm 100 μm 400 μm 100 μm 400 μmH&E H&E p16 p16 MYB MYB A .CC-BY-NC 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 25, 2025. ; https://doi.org/10.1101/2025.04.21.649795doi: bioRxiv preprint 25 Figure 1. HMSC shows strong p16 and MYB staining. Histopathologic images of HMSC taken at 100X (top row, scale bar represents 400 µm) and 400X (bottom row, scale bar represents 100 µm) show (A) H&E staining, (B) strong and diffuse p16 immunohistochemical (IHC) staining, and (C) strong and diffuse MYB IHC staining. .CC-BY-NC 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 25, 2025. ; https://doi.org/10.1101/2025.04.21.649795doi: bioRxiv preprint ACC vs. HMSC Malignant Cells UMAP_1 UMAP_2 4 0 -4 5-5 150 10 B C D E HMSC ACC CAF Endothelial + WBC Red= amplification Blue= deletions CNVs in ACC and HMSC Single CellsA ACC vs. HMSC Malignant Cells UMAP_1 UMAP_2 4 0 -4 5-5 150 10 Kaye ACC Score UMAP_1 UMAP_2 0.25 0.00 -0.25 4 0 -4 5-5 150 10 Up in HMSC Up in ACC FDR FDR Pathway Pathway 10-510-0 10-1010-2. 5 10-7. 5 10-12 10-210-0 10-1 10-3 count Modified Expression Distribution of Expression 0.7 0.9 1 1.1 1.30 2e+06 4e+06 6e+06 chr1 chr2 chr3 chr4 chr5 chr6 chr7 chr8 chr9 chr10 chr11 chr12 chr13 chr14 chr15 chr16 chr17 chr18 chr19 chr20 chr21 chr22 Figure 2 MYB Expression 0.0 2.5 5.0 7.5 .CC-BY-NC 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 25, 2025. ; https://doi.org/10.1101/2025.04.21.649795doi: bioRxiv preprint 26 Figure 2. Single-cell RNA sequencing reveals that HMSC is genomically and transcriptionally distinct from ACC. A) Heatmaps show CNVs in HMSC and ACC non-malignant (top) and malignant (bottom) single cells. Columns represent chromosomal regions; red represents amplification, while blue represents deletion. B) Uniform manifold approximation and projection (UMAP) shows that primary malignant cells from HMSC cluster separately from those of ACC tumors. All cells are colored by patient. C) UMAP shows MYB expression in primary malignant cells from HMSC and ACC tumors. D) UMAP shows Kaye ACC score in primary malignant cells from HMSC and ACC tumors. E) Bar plots show gene set enrichment analysis (GSEA) for pathways upregulated in HMSC relative to ACC (left panel, blue bars, dashed line represents false discovery rate (FDR) <0.05) and pathways upregulated in ACC relative to HMSC (right panel, red bars, dashed line represents FDR<0.05). .CC-BY-NC 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 25, 2025. ; https://doi.org/10.1101/2025.04.21.649795doi: bioRxiv preprint 11 moderate to strong MYB staining in a majority of HMSC cases.(4) Given the histologic similarities of these two entities, we also scored malignant cells for a previously defined ACC expression signature consisting of 1 160 mRNAs and 22 miRNAs,(18) to examine the degree to which HMSC resembled ACC at the transcriptional level (Figure 2D). Both ACC and HMSC cells scored highly for this signature. Thus, despite differences in global clustering and CNV profile, malignant HMSC cells shared common transcriptional features with ACC cells. Differential gene expression analysis comparing HMSC and ACC malignant cells revealed 94 genes that were more highly expressed in HMSC and 49 genes more highly expressed in ACC (FDR 1.25 fold-change, Figure S1D). Gene set enrichment analysis revealed that genes with higher expression in HMSC were enriched for E2F targets (Figure S1E), consistent with well described changes that occur with HPV infection.(32) Additionally, genes upregulated in HMSC demonstrated an enrichment of mitotic cell cycle genes and G2M checkpoint genes (Figure S1E), while genes with lower expression in HMSC were enriched for epithelial-to-mesenchymal transition (EMT) genes (Figure S1E). To estimate if this was due to a higher rate of cycling cells in HMSC, we filtered out cycling cells from both HMSC and ACC and repeated the analysis (Figure 2E and S1F-S1G). Even after the removal of cycling cells, genes higher in HMSC were still enriched in E2F targets and cell cycle genes (Figure 2E). For genes higher in ACC, filtering out cycling cells also revealed enrichment of KRAS signaling and EMT genes (Figure 2E). Intratumoral expression heterogeneity in HMSC We next explored malignant cell expression heterogeneity in HMSC. We utilized non-negative matrix factorization (NMF, see Methods) to uncover three distinct expression programs within malignant cells (Figure 3A-B and S2A). Program one was enriched with cell cycle genes, program two with TNF-a signaling via NF-kB and apoptosis genes, and program three with estrogen response, coagulation, and complement related genes (Figure S2B). Most cells showed either a strong program one or program two signature, with only a few cells strongly expressing program three (Figure 3A-B). .CC-BY-NC 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 25, 2025. ; https://doi.org/10.1101/2025.04.21.649795doi: bioRxiv preprint 0.000 0.025 0.050 0.075 0.100 −10 −5 0 5 10 Luminal−Myoepithelial spectrum Density r =−0.51 0.00 0.05 0.10 0.15 0.20 −10 −5 0 5 10 Luminal−Myoepithelial spectrum Density r =−0.08 B C D E Program 1 Program 2 Metagene.1 topMetagene.2 topMetagene.3 top RBM5 SYCP2 PDZD2 BCL11A SMARCA2 MAT2A RP5−857K21.8 MALAT1 HNRNPH1 RP1−222E13.10 NPIPA5 CALCOCO1 HMGB1 CCNL2 ANKRD36B MIR205HG RP11−429G19.2 KIAA0907 CDCA7 SMC4 SERPINB5 RCAN1 NFKBIA PMAIP1 CDKN1A ZC3H12A IER2 DUSP1 ZFP36 JUNB NR4A1 NR4A2 MAFF NCOA7 EGR1 ATF3 BTG2 BHLHE40 CLDN4 MYC WFDC2 C4B RP11−139E24.1 MMP20 PLIN5 RP11−528N21.3 CFI C5orf63 MLPH UNC5C RP11−393I2.2 C2orf74 XIST KRT19 FAM129A SHC2 FMO2 HSD17B2 PIGR LTF metagene.1 metagene.2 metagene.3 Z−score expression −4 −2 0 2 4 metagene.1 0 0.5 1 metagene.2 0 0.5 1 metagene.3 0 0.5 1 Metagene.1 topMetagene.2 topMetagene.3 top RBM5 SYCP2 PDZD2 BCL11A SMARCA2 MAT2A RP5−857K21.8 MALAT1 HNRNPH1 RP1−222E13.10 NPIPA5 CALCOCO1 HMGB1 CCNL2 ANKRD36B MIR205HG RP11−429G19.2 KIAA0907 CDCA7 SMC4 SERPINB5 RCAN1 NFKBIA PMAIP1 CDKN1A ZC3H12A IER2 DUSP1 ZFP36 JUNB NR4A1 NR4A2 MAFF NCOA7 EGR1 ATF3 BTG2 BHLHE40 CLDN4 MYC WFDC2 C4B RP11−139E24.1 MMP20 PLIN5 RP11−528N21.3 CFI C5orf63 MLPH UNC5C RP11−393I2.2 C2orf74 XIST KRT19 FAM129A SHC2 FMO2 HSD17B2 PIGR LTF metagene.1 metagene.2 metagene.3 Z−score expression −4 −2 0 2 4 metagene.1 0 0.5 1 metagene.2 0 0.5 1 metagene.3 0 0.5 1 Program 1Program 3 Program 2 Metagene.1 topMetagene.2 topMetagene.3 top RBM5 SYCP2 PDZD2 BCL11A SMARCA2 MAT2A RP5−857K21.8 MALAT1 HNRNPH1 RP1−222E13.10 NPIPA5 CALCOCO1 HMGB1 CCNL2 ANKRD36B MIR205HG RP11−429G19.2 KIAA0907 CDCA7 SMC4 SERPINB5 RCAN1 NFKBIA PMAIP1 CDKN1A ZC3H12A IER2 DUSP1 ZFP36 JUNB NR4A1 NR4A2 MAFF NCOA7 EGR1 ATF3 BTG2 BHLHE40 CLDN4 MYC WFDC2 C4B RP11−139E24.1 MMP20 PLIN5 RP11−528N21.3 CFI C5orf63 MLPH UNC5C RP11−393I2.2 C2orf74 XIST KRT19 FAM129A SHC2 FMO2 HSD17B2 PIGR LTF metagene.1 metagene.2 metagene.3 Z−score expression −4 −2 0 2 4 metagene.1 0 0.5 1 metagene.2 0 0.5 1 metagene.3 0 0.5 1 Program 3 Z Score Expression A UMAP_1 Program 1 UMAP_1 Program 2 UMAP_1 Program 3 2 0 -2 -4 -2 0 2 4 2 0 -2 -4 -2 0 2 4 2 0 -2 -4 -2 0 2 4 0 2 4 6 8 0.00 0.25 0.50 0.75 pearson count data acc hmsc Myo genes correlations Pearson Score Count Myoepithelial Gene Correlation 0 2 4 6 0.0 0.2 0.4 pearson count data acc hmsc Luminal genes correlations ACC HMSC 4 0 2 0.500.00 0.25 0.75 6 8 0 2 4 6 0.0 0.2 0.4 pearson count data acc hmsc Luminal genes correlations ACC HMSC Luminal Gene Correlation 0 2 4 6 0.0 0.2 0.4 pearson count data acc hmsc Luminal genes correlations Pearson Score Count 4 0 2 0.400.00 0.20 6 Myoepithelial-Luminal Spectrum −4 0 4 0 10 UMAP_1 UMAP_2 −5 0 5 luminal_over_myo −4 0 4 0 10 UMAP_1 UMAP_2 −5 0 5 luminal_over_myo Luminal Myoepithelial UMAP_1 UMAP_2 4 0 -4 0 10 r = -0.51 r = -0.08 ACC HMSC Myoepithelial cells Luminal cells Myoepithelial-Luminal Spectrum Density Density 0.1 00 0.0 75 0.0 50 0.0 25 0.0 00 -10 -5 0 5 10 -10 -5 0 5 10 0.2 00 0.1 50 0.1 00 0.0 50 0.0 00 Figure 3 .CC-BY-NC 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 25, 2025. ; https://doi.org/10.1101/2025.04.21.649795doi: bioRxiv preprint 27 Figure 3. HMSC lacks bicellular differentiation. A) Heatmap shows top genes within three programs in HMSC. Cells are arranged with hierarchical clustering. B) UMAPs show expression of the three programs in malignant HMSC cells. C) Bar plots show counts of binned Pearson correlations among myoepithelial (left) and luminal (right) genes in ACC (red) and HMSC (blue). Correlations are generally weaker in HMSC than ACC. D) UMAP shows HMSC and ACC malignant cells on a spectrum from myoepithelial (blue) to luminal (red). HMSC cells demonstrate an intermediate phenotype. E) Line plots show smoothed densities of cells at points along the myoepithelial-luminal spectrum in ACC (top) and HMSC (bottom). ACC cells exhibit a bimodal distribution while HMSC cells exhibit an intermediate phenotype. .CC-BY-NC 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 25, 2025. ; https://doi.org/10.1101/2025.04.21.649795doi: bioRxiv preprint 12 Notably, distinct myoepithelial and luminal expression programs were not uncovered in an unbiased fashion, as they were in ACC malignant cells.(8) As a result, we next assessed the presence of distinct myoepithelial and luminal genes in HMSC. We annotated HMSC malignant cell clusters based on the expression of known myoepithelial (Figure S2C) and luminal (Figure S2D) markers.(8) We then assessed the correlations within these myoepithelial and luminal gene sets for ACC and HMSC (Figure 3C). There were substantially higher correlations across luminal and myoepithelial genes in ACC than in HMSC, supporting the lack of distinct cell types in the latter. We then scored the HMSC and ACC malignant cells on a myoepithelial-luminal spectrum (Figure 3D-E, see Methods) and found that while ACC cells ranged from highly myoepithelial to highly luminal, HMSC cells largely scored toward the middle of the spectrum, suggesting that cells did not fall clearly into one or the other category. Thus, as previously described,(8) ACC showed a bimodal, bicellular differentiation distribution; by contrast, HMSC cells displayed a single intermediate phenotype along this axis. Given the prominent role of Notch signaling in ACC pathology,(33–35) we also scored all malignant cells based on Notch signaling using known Notch target genes (Figure S3A) and found that Notch was only activated in a handful of HMSC cells. We investigated the correlation between myoepithelial score and both Notch ligands and Notch targets in malignant HMSC cells (Figure S3B). Unlike ACC, there was no significant correlation between the myoepithelial score and Notch ligand expression (r=-0.088, p=0.31), likely because myoepithelial score does not capture essential biological features in the context of HMSC. There was a weak but statistically significant positive correlation between myoepithelial score and Notch target expression (r=0.21, p=0.015), but this association was dependent on a single cell and was no longer significant once it was removed (r= -0.028, p=0.75). Taken together, these results suggest that Notch signaling is unlikely to play an essential role in HMSC oncogenesis and heterogeneity as it does in ACC. .CC-BY-NC 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 25, 2025. ; https://doi.org/10.1101/2025.04.21.649795doi: bioRxiv preprint 13 HPV33 is associated with MYB expression in HMSC Given that HPV infection is a distinguishing trait of HMSC, we next characterized HPV expression in HMSC malignant cells. Plotting HMSC cells by number of reads per cell that map to the HPV33 genome revealed a bimodal distribution of HPV33 expression, with most cells expressing either 0 reads (henceforth HPVoff cells) or at least 24 reads (henceforth HPVon cells, see Figure 4A). Most malignant cells (105/134, 78.4%) were HPVon (Figure 4A-B). HPVon cells were distributed similarly across different clusters and on the uniform manifold approximation and projection (UMAP) (Figure 4B), suggesting HPV status may not drive a primary axis of transcriptional variability across cells. 264 genes were differentially expressed between HPVon and HPVoff cells, all of which were upregulated in HPVon cells (Figure 4C, Supplemental Table 1). Moreover, genes higher in HPVon cells were enriched for E2F targets, suggesting HPV activates E2F factors in HMSC (Figure 4D). One notable gene that was upregulated in HPVon cells was MYB (Figure 4C, fold-change=1.37, p<0.012, FDR<0.1). A significantly higher fraction of HPVon cells had overexpression of MYB (MYB+ cells) (Figure 4E, 83% vs. 62%, Fisher’s exact test p=0.022), and MYB targets were significantly upregulated in HPVon cells (Figure 4F, p=6.4 x 10-6). Taken together, these findings suggest HPV33 infection may drive MYB expression in HMSC tumors, instead of the canonical MYB translocations responsible for MYB upregulation in most ACC(5). Given that MYB may drive alternate cell fates in ACC,(5) including JAG1 and Notch1 expression in myoepithelial and luminal cells, respectively,(8) we next explored expression of Notch ligands and targets in HPVon and HPVoff cells in HMSC (Figure S3C-F). While there was higher expression of JAG1 in HPVon cells (p=0.01), no other significant differences were observed in expression of Notch ligands or targets between cell types, consistent with the notion that Notch signaling may not play a significant role in HMSC oncogenesis. HPV-MYB association is also observed in HPV+ oropharyngeal SCC. To validate the observed HPV33-MYB association in HMSC, we next assessed whether this association was retained in another HPV-related tumor of the head and neck, OPSCC. We reanalyzed scRNA-seq .CC-BY-NC 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 25, 2025. ; https://doi.org/10.1101/2025.04.21.649795doi: bioRxiv preprint B FE A C 0 reads 1 read 2 reads 3-23 reads 24+ reads Count HPV33 Reads Per Cell 100 75 50 25 0 28 0 1 0 105 HPV33 Expression UMAP_1 UMAP_2 -2 0 2 -2 -1 0 1 2 D Figure 4 HPVon HPVoff Average log2 fold-change ANLN TUBB6 AXL IFT122 GFM1SEPT2 PBRM1 INSIG1CAPG POP7 MYC MYB JAG1 10-0 10-2 10-4 10-6 −0.5 0.0 0.5 1.0 1.5 Average log2 fold−change p−value Significant DEG (FDR1.25) Genes up in HPV33 positive Same Differential Gene Expression by HPV StatusHMSC DEGs by HPV Status p value 10-6 10-4 10-2 10-0 0.0 0.5 1.0 1.5 Genes up in HPVon Significant DEG (FDR 1.25) -0.5 (N=29) HPVoff (N=105) HPVon 90% 80% 60% 50% 40% 30% 10% 0% MYB Expression within HMSC 100% 70% 20% MYB+ MYB- p-value = 0.0224 HPVonHPVoff MYB target score Downstream MYB targets within HMSC 6.4e-06 4 3 2 HPVoff HPVon p = 6.4e-06 HALLMARK_PI3K_AKT_MTOR_SIGNALING HALLMARK_WNT_BETA_CATENIN_SIGNALING GO_MITOTIC_CELL_CYCLE HALLMARK_MYC_TARGETS_V2 HALLMARK_INTERFERON_ALPHA_RESPONSE HALLMARK_CHOLESTEROL_HOMEOSTASIS HALLMARK_E2F_TARGETS 0.150 0.087 0.047 0.026 0.014 0.008 FDR FDR 0.05 0.10 NES 1.7 1.6 1.5 1.4 up 0.087 HALLMARK_PI3K_AKT_MTOR_SIGNALING HALLMARK_WNT_BETA_CATENIN_SIGNALING GO_MITOTIC_CELL_CYCLE HALLMARK_MYC_TARGETS_V2 HALLMARK_INTERFERON_ALPHA_RESPONSE HALLMARK_CHOLESTEROL_HOMEOSTASIS HALLMARK_E2F_TARGETS 0.150 0.087 0.047 0.026 0.014 0.008 FDR FDR 0.05 0.10 NES 1.7 1.6 1.5 1.4 up 0.05 GSEA for HPVon HMSC Cells FDR 0.150 FDR 1.7 1.4 0.10 Up 1.5 NES 0.047 0.026 0.014 0.008 1.6 HALLMARK_PI3K_AKT_MTOR_SIGNALING HALLMARK_WNT_BETA_CATENIN_SIGNALING GO_MITOTIC_CELL_CYCLE HALLMARK_MYC_TARGETS_V2 HALLMARK_INTERFERON_ALPHA_RESPONSE HALLMARK_CHOLESTEROL_HOMEOSTASIS HALLMARK_E2F_TARGETS 0.1500.0870.0470.0260.0140.008 FDR FDR 0.05 0.10 NES 1.7 1.6 1.5 1.4 up HALLMARK_E2F_TARGETS HALLMARK_CHOLESETEROL HOMOEOSTASIS HALLMARK_INTERFERON ALPHA_RESPONSE HALLMARK_MYC TARGETS_V2 GO_MITOTIC CELL_CYCLE HALLMARK_WNT_BETA CATENIN_SIGNALING HALLMARK_PI3K_AKT MTOR_ SIGNALING .CC-BY-NC 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 25, 2025. ; https://doi.org/10.1101/2025.04.21.649795doi: bioRxiv preprint 28 Figure 4. HPV is associated with MYB in HMSC. A) Bar plot shows counts of HMSC malignant cells with the specified numbers of HPV33 reads detected. Most cells had either 0 reads (28/134) or >24 reads (105/134). B) UMAP shows the distribution of HPVon cells across all HMSC malignant cells. HPVon cells are shown in blue. C) Volcano plot shows differentially expressed genes (DEGs) in HPVon HMSC malignant cells with upregulated genes shown in red (FDR1.25 fold-change). MYB is upregulated in HPVon cells. D) Dot plot shows GSEA for pathways that are upregulated in HPVon relative to HPVoff HMSC cells (dotted line marks FDR<0.05, color indicates significance, size indicates Normalized Enrichment Score (NES)). E) Stacked bar plots show proportion of MYB+ cells by HPV status. HPVon cells show greater MYB positivity than HPVoff cells (83% vs 62%, p=0.0224). F) Violin plot shows upregulation of downstream MYB targets in HPVon cells relative to HPVoff cells in HMSC (t-test, p=6.4x10-6). .CC-BY-NC 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 25, 2025. ; https://doi.org/10.1101/2025.04.21.649795doi: bioRxiv preprint 14 data from a previously published cohort of 10 patients with HPV+ OPSCC (11) to identify HPVon and HPVoff cells using similar approach. A significantly higher fraction of HPVon cells were MYB+ in 7/10 tumors (Figure 5A, Fisher’s exact test, p<0.05). Moreover, MYB targets were upregulated in HPVon cells (relative to HPVoff cells) in all 10 tumors assessed (Figure 5B), supporting the HPV-MYB association seen in HMSC. Finally, we evaluated MYB expression across all OPSCC tumors in The Cancer Genome Atlas (TCGA), finding significant upregulation of MYB in HPV+ tumors relative to HPV- tumors (Figure 5C). We next sought to assess the degree to which HPV-related gene expression in these two tumor types was similar. In the cohort of HPV+ OPSCC, 16 genes were significantly differentially upregulated in HPVon cells (FDR1.25) (Figure 5D, Supplemental Table 2). These genes did not overlap with genes significantly upregulated in HMSC HPVon cells. However, as in HMSC, genes upregulated in HPVon cells were enriched for E2F targets in OPSCC as well (Figure 5E). We next devised a gene signature based on differentially expressed genes in HMSC and examined its significance in OPSCC tumors. All 10 tumors exhibited a statistically significant upregulation of this gene signature within their HPVon cells relative to the HPVoff cells (Figure 6A), suggesting this signature may carry importance in OPSCC as well. To understand its prognostic significance, we examined survival associations of this gene signature in TCGA patients with OPSCC. As previously mentioned, some of the genes in the signature were associated with the cell cycle and E2F targets, so we specifically excluded these genes to demonstrate the association with survival was not merely a reflection of increased proliferation rate. While there was no survival association in patients with HPV- OPSCC (Figure S4A- S4B), we found a negative association between the signature and overall survival for all OPSCC patients as well as specifically amongst HPV+ OPSCC patients (Figure 6B, Figure S4C-S4D). In contrast with the notion that HPV expression in OPSCC is typically associated with a favorable prognosis,(11,36) this signature may reflect an alternate, negative role of HPV in HPV+ OPSCC that is not well described and may contribute to variability in patient outcomes.(36,37) .CC-BY-NC 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 25, 2025. ; https://doi.org/10.1101/2025.04.21.649795doi: bioRxiv preprint 1.61e−31 3.68e−06 4.44e−11 1.21e−99 9.27e−07 6.2e−116 5.71e−16 3.8e−08 4.01e−13 0.0135 OP17 OP20 OP33 OP34 OP35 OP4 OP5 OP6 OP9 OP14 HPV− HPV+ HPV− HPV+ HPV− HPV+ HPV− HPV+ HPV− HPV+ 0.1 0.2 0.3 0.4 0.5 0.1 0.2 0.3 0.4 0.5Myb Target Score HPV Status HPV− HPV+ MYB targets by HPV Status D B E A 100% 90% 80% 70% 60% 50% 40% 30% 20% 10% 0% 100% 90% 80% 70% 60% 50% 40% MYB + p-value=1.61e -05 OP14 OP17 OP20 OP33 OP34 MYB status p-value=0.00514 p-value=1.26e -07 p-value=0.000475p-value=0.00796 MYB - p-value=0.336 OP35 OP4 OP5 OP6 OP9 p-value=0.0456 p-value=0.778 p-value=0.367p-value=5.93e -06 HPV- HPV+ (n=1243) (n=646) HPV- HPV+ (n=169) (n=132) HPV- HPV+ (n=169) (n=132) HPV- HPV+ (n=1675) (n=250) HPV- HPV+ (n=113) (n=44) HPV- HPV+ (n=74) (n=78) (n=977) (n=1,462) HPV- HPV+ (n=1,103) (n=2,060) HPV- HPV+ (n=1,420) (n=159) HPV- HPV+ (n=131) (n=19) HPV- HPV+ 30% 20% 10% 0% 95% 5% 89% 11% 0.1% 99% 5% 95% 3% 4% 97% 96% 61% 58% 42%39% 17% 83% 52% 48% 5%2% 98% 95%89%96% 4% 11%12%6% 94% 88%95% 5%3% 97% 17% 83%87% 3% AA MYB Expression within OPSCC Figure 5 A OPSCC DEGs by HPV Status Genes up in HPVoff Significant DEG (FDR 1.25) Genes up in HPVon GSEA for HPVon OPSCC Cells 2.2 2.0 1.8 1.6 NES FDR 0.050 0.075 0.100 0.125 8.0e-04 FDR 2.4e-01 G2M_CHECKPOINT E2F_TARGETS GO_MITOTIC_CELL_CYCLE OXIDATIVE_PHOSPHORYLATION FATTY_ACID_METABOLISM MYC_TAR GETS_V1 ESTRO GEN_RESPONSE_LATE SPERMATOGENESIS C 4.04e−11 0 2 4 HPV− HPV+ hpv_status MYB expression: log2(TPM) HPV Status HPV− HPV+ MYB expression by HPV statusMYB expression by HPV status HPVoff HPVon 4 2 0 MYB expression: log ₂(TPM) p = 4.04e -11 HPV Status HPVoff HPVon HALLMARK_ESTROGEN_RESPONSE_LATE HALLMARK_OXIDATIVE_PHOSPHORYLATION HALLMARK_SPERMATOGENESIS HALLMARK_FATTY_ACID_METABOLISM HALLMARK_MYC_TARGETS_V1 GO_MITOTIC_CELL_CYCLE HALLMARK_G2M_CHECKPOINT HALLMARK_E2F_TARGETS 2.4e−01 1.4e−02 8.0e−04 4.5e−05 2.6e−06 1.5e−07 FDR NES 2.2 2.0 1.8 1.6 FDR 0.025 0.050 0.075 0.100 0.125 up 2.6e-04 HALLMARK_ESTROGEN_RESPONSE_LATE HALLMARK_OXIDATIVE_PHOSPHORYLATION HALLMARK_SPERMATOGENESIS HALLMARK_FATTY_ACID_METABOLISM HALLMARK_MYC_TARGETS_V1 GO_MITOTIC_CELL_CYCLE HALLMARK_G2M_CHECKPOINT HALLMARK_E2F_TARGETS 2.4e−01 1.4e−02 8.0e−04 4.5e−05 2.6e−06 1.5e−07 FDR NES 2.2 2.0 1.8 1.6 FDR 0.025 0.050 0.075 0.100 0.125 up 0.025 HPVoff HPVon HPVoff HPVon HPVoff HPVon HPVoff HPVon HPVoff HPVon HPVoff HPVon HPVoff HPVon HPVoff HPVon HPVoff HPVon HPVoff HPVon MYB Targets by HPV Status OP4 OP5 OP6 OP9 OP14 OP17 OP20 OP33 OP34 OP35 HPV Status HPVoff HPVon Myb Target Score Myb Target Score 0.5 0.1 0.4 0.3 0.2 HPVonHPVoff HPVonHPVoff HPVonHPVoff HPVonHPVoff 5.71e-16 HPVonHPVoff 0.5 0.1 0.4 0.3 0.2 HPVonHPVoff HPVonHPVoff HPVonHPVoff HPVonHPVoff HPVonHPVoff HALLMARK_ESTROGEN_RESPONSE_LATE HALLMARK_OXIDATIVE_PHOSPHORYLATION HALLMARK_SPERMATOGENESIS HALLMARK_FATTY_ACID_METABOLISM HALLMARK_MYC_TARGETS_V1 GO_MITOTIC_CELL_CYCLE HALLMARK_G2M_CHECKPOINT HALLMARK_E2F_TARGETS 2.4e−011.4e−028.0e−044.5e−052.6e−061.5e−07 FDR NES 2.2 2.0 1.8 1.6 FDR 0.025 0.050 0.075 0.100 0.125 up .CC-BY-NC 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 25, 2025. ; https://doi.org/10.1101/2025.04.21.649795doi: bioRxiv preprint 29 Figure 5. HPV is associated with MYB in OPSCC. A) Stacked bar plots show MYB+ cells by HPV status in OPSCC tumors. Higher MYB positivity (Fisher’s exact test, p<0.05) was seen in HPVon cells for 7/10 individual OPSCC tumors. B) Violin plots demonstrate upregulation of MYB targets in HPVon relative to HPVoff cells (t-test, p<0.05) in all 10 individual OPSCC tumors. C) Violin plot shows higher MYB expression in HPVon (right) versus HPVoff (left) OPSCC tumors from TCGA (t-test, p=4.04x10-11). D) Volcano plot shows DEGs in HPVon OPSCC malignant cells. Red represents genes upregulated in HPVon cells (FDR1.25 fold-change), and green represents genes upregulated in HPVoff cells (FDR1.25 fold-change). E) Dot plot shows GSEA for pathways upregulated in HPVon relative to HPVoff OPSCC cells (dotted line marks FDR<0.05, color indicates significance, size indicates NES). .CC-BY-NC 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 25, 2025. ; https://doi.org/10.1101/2025.04.21.649795doi: bioRxiv preprint 3.4e−31 1.3e−07 8e−06 7.9e−42 2e−09 4.1e−110 6.9e−26 2.1e−10 1.2e−16 0.0083 OP17 OP20 OP33 OP34 OP35 OP4 OP5 OP6 OP9 OP14 HPV− HPV+ HPV− HPV+ HPV− HPV+ HPV− HPV+ HPV− HPV+ 0.2 0.4 0.6 0.8 0.2 0.4 0.6 0.8 hpv Signature Score hpv HPV− HPV+ HMSC HPV signature B AAA HMSC HPV Signature in OPSCC by HPV Status OP4 Signature Score OP5 OP6 OP9 OP14 Signature Score HPV Status HPVoff HPVon Log-rank p = 0.0028 Log-rank p = 0.0079 Figure 6 0.8 0.6 0.4 HPVoff HPVon 0.2 HPVoff HPVon HPVoff HPVon HPVoff HPVon HPVoff HPVon HPVoff HPVonHPVoff HPVonHPVoff HPVonHPVoff HPVon 0.8 0.6 0.4 HPVoff HPVon 0.2 OP17 OP20 OP33 OP34 OP35 Survival Plots for HMSC HPV signature survival probability All OPSCC HPV+ OPSCC + +++ +++++++++++++++ ++++ + ++ +++ +++ ++++++++ ++ ++++ +++++ + ++ p = 0.0024Log−rank 0.00 0.25 0.50 0.75 1.00 0 12 24 36 48 60 Time since diagnosis (months) Survival probability Legend + +High HMSC HPV signature no CC (n = 39) Low HMSC HPV signature no CC (n = 38) OPSCC samples 39 31 8 3 1 0 38 33 22 16 10 0Low HMSC HPV signature no CC (n = 38) High HMSC HPV signature no CC (n = 39) 0 12 24 36 48 60 Time since diagnosis (months) Legend Number at risk + ++ ++++++++++ +++ + + ++ +++ + ++ +++ + ++ +++ + ++ p = 0.007Log−rank 0.00 0.25 0.50 0.75 1.00 0 12 24 36 48 60 Time since diagnosis (months) Survival probability Legend + +High HMSC HPV signature no CC (n = 25) Low HMSC HPV signature no CC (n = 24) OPSCC HPV+ samples 25 21 6 3 1 0 24 19 16 12 10 0Low HMSC HPV signature no CC (n = 24) High HMSC HPV signature no CC (n = 25) 0 12 24 36 48 60 Time since diagnosis (months) Legend Number at risk .CC-BY-NC 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 25, 2025. ; https://doi.org/10.1101/2025.04.21.649795doi: bioRxiv preprint 30 Figure 6. HMSC HPV signature is upregulated in HPV+ OPSCC and is associated with poor prognosis. A) Violin plots show upregulation of HMSC HPV signature in HPVon relative to HPVoff cells (t- test, p<0.05) for all ten individual OPSCC tumors (each p<0.01). B) Kaplan-Meier curves of TCGA all OPSCC (left) and HPV+ OPSCC (right) patients stratified by expression of HMSC HPV signature. Red survival curves represent patients with high expression of HMSC HPV genes and demonstrate poorer survival in both analyses (p<0.01). .CC-BY-NC 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 25, 2025. ; https://doi.org/10.1101/2025.04.21.649795doi: bioRxiv preprint 15

HMSC represents an extremely rare entity, with less than 100 total cases described in the English literature. In this study, we present the first scRNA-seq analysis of this unique tumor type and provide a comparative analysis with scRNA-seq data from ACC and OPSCC. Our findings highlighted several important transcriptional differences between HMSC and ACC, including the distinct clustering of HMSC cells from ACC cells, the lack of bicellular differentiation in HMSC without defined myoepithelial or luminal programs, and the lack of oncogenic Notch signaling. These findings are consistent with previous histologic studies of HMSC(1–3) and strongly support the notion that this is a distinct molecular entity, not to be characterized as a subset of ACC. Notable expression differences between HMSC and ACC included lower expression of EMT genes and KRAS signaling in HMSC compared to ACC. EMT is hypothesized to be associated with invasion and metastasis in several epithelial tumor types,(38–40) and in ACC we previously found upregulation of mesenchymal genes in myoepithelial cells.(8) Decreased EMT in HMSC may contribute to the more indolent clinical behavior of these tumors, especially in comparison to their typical high grade histology at presentation.(3) Downregulated KRAS signaling may have a similar impact, as activating mutations in the KRAS gene have been associated with poorer survival in head and neck squamous cell carcinoma.(41) The most prominent non-cycling population of cells in this tumor was enriched for TNFa/NFkB signaling. Among its many effects, HPV infection has been shown to induce TNFa/NFkB signaling,(42) which may have tumor suppressive activity.(43,44) Upregulation in this pathway is associated with improved survival for patients with HPV+ tumors,(45) again consistent with the slow- growing nature of HMSC tumors. These observations may begin to account for some of the distinct clinical behavior observed in HMSC compared to ACC and other similar histologic subtypes. Despite these differences between HMSC and ACC, HMSC did display concordant expression of a previously described ACC transcriptomic signature,(18) suggesting some expression similarities between these entities. These similarities may be related to upregulation of the MYB proto-oncogene in .CC-BY-NC 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 25, 2025. ; https://doi.org/10.1101/2025.04.21.649795doi: bioRxiv preprint 16 both tumor types,(4,7) though by distinct mechanisms, as HMSC lacks the canonical MYB translocations described in ACC.(5–7) Here, we uncovered an association between HPV33 infection and MYB expression that suggests an alternate, HPV-driven, mechanism of MYB upregulation in HMSC. Furthermore, this HPV-MYB association was validated in single cell and bulk transcriptional data in OPSCC, another HPV-related tumor type in the head and neck, which suggests it may be a pan-cancer association. Interestingly, the downstream impacts of both MYB and HPV in HMSC may be different than those previously described in other tumor types. Specifically, regarding MYB, the lack of bicellular differentiation and Notch signaling or Notch ligand expression, as observed in ACC,(8) suggests that MYB may have a distinct downstream role in HMSC. Further analysis of MYB binding sites in HMSC may help to elucidate these alternate mechanisms. In addition, regarding HPV, we uncovered 264 genes differentially upregulated in HPVon cells in HMSC. Although these genes were upregulated in HPVon OPSCC cells in all 10 tumors assessed, they were separate from genes differentially expressed by OPSCC HPVon cells in an unbiased analysis. Moreover, this gene signature was associated with poor survival in the TCGA OPSCC cohort, both when assessing all tumors and when filtering specifically to HPV+ tumors. In conjunction with data suggesting HPV positivity typically portends improved survival and outcomes in OPSCC,(11,36,37) the poor survival association of this HPV-related signature suggests an alternate role of HPV that is less well characterized. Developing a better understanding of this signature may help to elucidate reasons for treatment failure in HPV+ OPSCC, a tumor that typically responds well to radiation therapy.(46) Further investigation of these genes and this program may be instrumental in identifying future opportunities for new therapeutics for HPV+ OPSCC. Refinement of this signature was limited by the small number of cells and the inclusion of only a single patient sample. However, HMSC is an exceedingly rare tumor entity, and both the HPV-MYB association and the significance of the HMSC HPV gene signature were validated in OPSCC. Prospectively capturing more HMSC samples is additionally complicated by the fact that these tumors are often not identified until a definitive diagnosis is made on the final pathology report. As a rare disease, .CC-BY-NC 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 25, 2025. ; https://doi.org/10.1101/2025.04.21.649795doi: bioRxiv preprint 17 HMSC also lacks cell line or PDX models to validate findings, underscoring the need for further model development in this rare tumor type. Finally, the identified HMSC HPV gene signature requires further validation in larger patient cohorts and across diverse HPV+ tumor types to better understand its biological relevance, prognostic implications, and value as a therapeutic target. In conclusion, this study reports the first scRNA-seq analysis of a rare HMSC tumor and highlights the transcriptional differences between HMSC and ACC. Whereas ACC cells typically exhibit either myoepithelial or luminal differentiation, HMSC cells demonstrated a single cellular phenotype. Our findings also highlighted a potential HPV-related mechanism of MYB upregulation within HMSC, distinguishing it from the MYB-NFIB fusion-driven overexpression in ACC. Further, we identified HPV- related genes in HMSC that were distinct from an HPV-related signature in OPSCC, suggesting alternate effects of HPV in these two tumor types. These HMSC HPV-related genes were associated with poor overall survival in TCGA OPSCC tumors, suggesting this alternate HPV effect may be relevant across HPV-related tumors. Further work is warranted to explore the biologic and prognostic utility of this program across HPV-related malignancies. .CC-BY-NC 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 25, 2025. ; https://doi.org/10.1101/2025.04.21.649795doi: bioRxiv preprint 18

Funding: Supported by V Foundation (S.V.P.), Cancer Research Foundation (S.V.P.), Barnes Jewish Hospital Foundation (S.V.P.), NCI 1K08CA237732 (S.V.P), the European Research Council Horizon 2020 grant 949029 (Y.D.), and NIDCR 1K08DE033093 (A.S.P.). The funding sources had no involvement in the design, conduct, and reporting of the research. Author Contributions: AW, MDAS, ASP, and YD were responsible for designing and executing the study as well as writing/editing the manuscript. DTL and SVP provided patient samples and edited the manuscript. MM, WB, OO, VY, IC, SC, SHT, WCF, and IT contributed to data interpretation, data analysis, and editing the manuscript. ASP, YD, SVP, and IT were also responsible for supervision. .CC-BY-NC 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 25, 2025. ; https://doi.org/10.1101/2025.04.21.649795doi: bioRxiv preprint 19 Competing Interests I.T. is a co-founder and adviser of Cellyrix Therapeutics and an adviser of Immunitas Therapeutics. All other authors declare no potential conflicts of interest.

1. Bishop JA, Andreasen S, Hang JF, Bullock MJ, Chen TY, Franchi A, et al. HPV-related Multiphenotypic Sinonasal Carcinoma: An Expanded Series of 49 Cases of the Tumor Formerly Known as HPV-related Carcinoma With Adenoid Cystic Carcinoma-like Features. American Journal of Surgical Pathology. 2017 Dec;41(12):1690–701. 2. Andreasen S, Bishop JA, Hansen TVO, Westra WH, Bilde A, Von Buchwald C, et al. Human papillomavirus‐related carcinoma with adenoid cystic‐like features of the sinonasal tract: clinical and morphological characterization of six new cases. Histopathology. 2017 May;70(6):880–8. 3. Bishop JA, Westra WH. Human papillomavirus-related multiphenotypic sinonasal carcinoma: An emerging tumor type with a unique microscopic appearance and a paradoxical clinical behaviour. Oral Oncology. 2018 Dec;87:17–20. 4. Shah AA, Oliai BR, Bishop JA. Consistent LEF-1 and MYB Immunohistochemical Expression in Human Papillomavirus-Related Multiphenotypic Sinonasal Carcinoma: A Potential Diagnostic Pitfall. Head and Neck Pathol. 2019 Jun;13(2):220–4. 5. Drier Y, Cotton MJ, Williamson KE, Gillespie SM, Ryan RJH, Kluk MJ, et al. An oncogenic MYB feedback loop drives alternate cell fates in adenoid cystic carcinoma. Nat Genet. 2016 Mar;48(3):265–72. .CC-BY-NC 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 25, 2025. ; https://doi.org/10.1101/2025.04.21.649795doi: bioRxiv preprint 20 6. Mitani Y, Liu B, Rao PH, Borra VJ, Zafereo M, Weber RS, et al. Novel MYBL1 Gene Rearrangements with Recurrent MYBL1–NFIB Fusions in Salivary Adenoid Cystic Carcinomas Lacking t(6;9) Translocations. Clinical Cancer Research. 2016 Feb 1;22(3):725–33. 7. Brayer KJ, Frerich CA, Kang H, Ness SA. Recurrent Fusions in MYB and MYBL1 Define a Common, Transcription Factor–Driven Oncogenic Pathway in Salivary Gland Adenoid Cystic Carcinoma. Cancer Discovery. 2016 Feb 1;6(2):176–87. 8. Parikh AS, Wizel A, Davis D, Lefranc-Torres A, Rodarte-Rascon AI, Miller LE, et al. Single-cell RNA sequencing identifies a paracrine interaction that may drive oncogenic notch signaling in human adenoid cystic carcinoma. Cell Reports. 2022 Nov;41(9):111743. 9. Lechner M, Liu J, Masterson L, Fenton TR. HPV-associated oropharyngeal cancer: epidemiology, molecular biology and clinical management. Nat Rev Clin Oncol. 2022 May;19(5):306–27. 10. Sabatini ME, Chiocca S. Human papillomavirus as a driver of head and neck cancers. Br J Cancer. 2020 Feb 4;122(3):306–14. 11. Puram SV, Mints M, Pal A, Qi Z, Reeb A, Gelev K, et al. Cellular states are coupled to genomic and viral heterogeneity in HPV-related oropharyngeal carcinoma. Nat Genet. 2023 Apr;55(4):640–50. 12. Picelli S, Faridani OR, Björklund ÅK, Winberg G, Sagasser S, Sandberg R. Full-length RNA-seq from single cells using Smart-seq2. Nat Protoc. 2014 Jan;9(1):171–81. 13. Picelli S, Björklund ÅK, Faridani OR, Sagasser S, Winberg G, Sandberg R. Smart-seq2 for sensitive full-length transcriptome profiling in single cells. Nat Methods. 2013 Nov;10(11):1096–8. 14. Puram SV, Tirosh I, Parikh AS, Patel AP, Yizhak K, Gillespie S, et al. Single-Cell Transcriptomic Analysis of Primary and Metastatic Tumor Ecosystems in Head and Neck Cancer. Cell. 2017 Dec;171(7):1611-1624.e24. .CC-BY-NC 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 25, 2025. ; https://doi.org/10.1101/2025.04.21.649795doi: bioRxiv preprint 21 15. Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics. 2013 Jan 1;29(1):15–21. 16. Liao Y, Smyth GK, Shi W. featureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics. 2014 Apr 1;30(7):923–30. 17. Stuart T, Butler A, Hoffman P, Hafemeister C, Papalexi E, Mauck WM, et al. Comprehensive Integration of Single-Cell Data. Cell. 2019 Jun;177(7):1888-1902.e21. 18. Gao R, Cao C, Zhang M, Lopez MC, Yan Y, Chen Z, et al. A unifying gene signature for adenoid cystic cancer identifies parallel MYB-dependent and MYB-independent therapeutic targets. Oncotarget. 2014 Dec 30;5(24):12528–42. 19. Kotliar D, Veres A, Nagy MA, Tabrizi S, Hodis E, Melton DA, et al. Identifying gene expression programs of cell-type identity and cellular activity with single-cell RNA-Seq. eLife. 2019 Jul 8;8:e43803. 20. Korotkevich G, Sukhov V, Budin N, Shpak B, Artyomov MN, Sergushichev A. Fast gene set enrichment analysis [Internet]. Bioinformatics; 2016 [cited 2025 Mar 26]. Available from: http://biorxiv.org/lookup/doi/10.1101/060012 21. Federico A, Monti S. hypeR: an R package for geneset enrichment workflows. Wren J, editor. Bioinformatics. 2020 Feb 15;36(4):1307–8. 22. Liberzon A, Birger C, Thorvaldsdóttir H, Ghandi M, Mesirov JP, Tamayo P. The Molecular Signatures Database Hallmark Gene Set Collection. Cell Systems. 2015 Dec;1(6):417–25. 23. Tickle T, Tirosh I, Georgescu C, Brown M, Haas B. inferCNV: Visualizing large-scale copy number variation in single cell RNA-seq expression data. [Internet]. 2019. Available from: https://www.bioconductor.org/packages/devel/bioc/vignettes/infercnv/inst/doc/inferCNV.html .CC-BY-NC 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 25, 2025. ; https://doi.org/10.1101/2025.04.21.649795doi: bioRxiv preprint 22 24. R Core Team. R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing; 2024. 25. Wickham H. ggplot2: elegant graphics for data analysis. Second edition. Cham: Springer international publishing; 2016. 1 p. (Use R!). 26. Gu Z, Eils R, Schlesner M. Complex heatmaps reveal patterns and correlations in multidimensional genomic data. Bioinformatics. 2016 Sep 15;32(18):2847–9. 27. Patil I. Visualizations with statistical details: The “ggstatsplot” approach. JOSS. 2021 May 25;6(61):3167. 28. Kassambara A. rstatix: Pipe-Friendly Framework for Basic Statistical Tests [Internet]. 2019 [cited 2025 Mar 26]. p. 0.7.2. Available from: https://CRAN.R-project.org/package=rstatix 29. Colaprico A, Silva TC, Olsen C, Garofano L, Cava C, Garolini D, et al. TCGAbiolinks: an R/Bioconductor package for integrative analysis of TCGA data. Nucleic Acids Research. 2016 May 5;44(8):e71–e71. 30. Cerami E, Gao J, Dogrusoz U, Gross BE, Sumer SO, Aksoy BA, et al. The cBio Cancer Genomics Portal: An Open Platform for Exploring Multidimensional Cancer Genomics Data. Cancer Discovery. 2012 May 1;2(5):401–4. 31. Liu J, Lichtenberg T, Hoadley KA, Poisson LM, Lazar AJ, Cherniack AD, et al. An Integrated TCGA Pan-Cancer Clinical Data Resource to Drive High-Quality Survival Outcome Analytics. Cell. 2018 Apr;173(2):400-416.e11. 32. Galloway DA, Laimins LA. Human papillomaviruses: shared and distinct pathways for pathogenesis. Current Opinion in Virology. 2015 Oct;14:87–92. .CC-BY-NC 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 25, 2025. ; https://doi.org/10.1101/2025.04.21.649795doi: bioRxiv preprint 23 33. Ho AS, Ochoa A, Jayakumaran G, Zehir A, Valero Mayor C, Tepe J, et al. Genetic hallmarks of recurrent/metastatic adenoid cystic carcinoma. Journal of Clinical Investigation. 2019 Sep 4;129(10):4276–89. 34. Hanna GJ, Bae JE, Lorch JH, Schoenfeld JD, Margalit DN, Tishler RB, et al. Long-term outcomes and clinicogenomic correlates in recurrent, metastatic adenoid cystic carcinoma. Oral Oncology. 2020 Jul;106:104690. 35. Wang S, Yu Y, Fang Y, Huang H, Wu D, Fang H, et al. Whole-exome sequencing reveals genetic underpinnings of salivary adenoid cystic carcinoma in the Chinese population. Journal of Genetics and Genomics. 2020 Jul;47(7):397–401. 36. Sinha P, Karadaghy OA, Doering MM, Tuuli MG, Jackson RS, Haughey BH. Survival for HPV- positive oropharyngeal squamous cell carcinoma with surgical versus non-surgical treatment approach: A systematic review and meta-analysis. Oral Oncology. 2018 Nov;86:121–31. 37. Fischer CA, Kampmann M, Zlobec I, Green E, Tornillo L, Lugli A, et al. p16 expression in oropharyngeal cancer: its impact on staging and prognosis compared with the conventional clinical staging parameters. Annals of Oncology. 2010 Oct;21(10):1961–6. 38. Ting DT, Wittner BS, Ligorio M, Vincent Jordan N, Shah AM, Miyamoto DT, et al. Single-Cell RNA Sequencing Identifies Extracellular Matrix Gene Expression by Pancreatic Circulating Tumor Cells. Cell Reports. 2014 Sep;8(6):1905–18. 39. Yu M, Bardia A, Wittner BS, Stott SL, Smas ME, Ting DT, et al. Circulating Breast Tumor Cells Exhibit Dynamic Changes in Epithelial and Mesenchymal Composition. Science. 2013 Feb;339(6119):580–4. .CC-BY-NC 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 25, 2025. ; https://doi.org/10.1101/2025.04.21.649795doi: bioRxiv preprint 24 40. Ye X, Weinberg RA. Epithelial–Mesenchymal Plasticity: A Central Regulator of Cancer Progression. Trends in Cell Biology. 2015 Nov;25(11):675–86. 41. Tinhofer I, Budach V, Saki M, Konschak R, Niehr F, Jöhrens K, et al. Targeted next-generation sequencing of locally advanced squamous cell carcinomas of the head and neck reveals druggable targets for improving adjuvant chemoradiation. European Journal of Cancer. 2016 Apr;57:78–86. 42. James MA, Lee JH, Klingelhutz AJ. Human Papillomavirus Type 16 E6 Activates NF-κB, Induces cIAP-2 Expression, and Protects against Apoptosis in a PDZ Binding Motif-Dependent Manner. J Virol. 2006 Jun;80(11):5301–7. 43. Tilborghs S, Corthouts J, Verhoeven Y, Arias D, Rolfo C, Trinh XB, et al. The role of Nuclear Factor-kappa B signaling in human cervical cancer. Critical Reviews in Oncology/Hematology. 2017 Dec;120:141–50. 44. Basile JR, Zacny V, Münger K. The Cytokines Tumor Necrosis Factor-α (TNF-α) and TNF-related Apoptosis-inducing Ligand Differentially Modulate Proliferation and Apoptotic Pathways in Human Keratinocytes Expressing the Human Papillomavirus-16 E7 Oncoprotein. Journal of Biological Chemistry. 2001 Jun;276(25):22522–8. 45. Schrank TP, Prince AC, Sathe T, Wang X, Liu X, Alzhanov DT, et al. NF-κB over-activation portends improved outcomes in HPV-associated head and neck cancer. Oncotarget. 2022 May 24;13(1):707–22. 46. Wang H, Wang B, Wei J, Meng L, Zhang Q, Qu C, et al. Molecular mechanisms underlying increased radiosensitivity in human papillomavirus-associated oropharyngeal squamous cell carcinoma. Int J Biol Sci. 2020;16(6):1035–43. .CC-BY-NC 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 25, 2025. ; https://doi.org/10.1101/2025.04.21.649795doi: bioRxiv preprint