HPV⁺ (UM-SCC-47, UM-SCC-104) and HPV^- HNSCC (UM-SCC-6, UM-SCC-17A, UM-SCC-1) cell lines were purchased from Sigma. UPCI : SCC090 were from ATCC. Cells were incubated in DMEM medium supplemented with 1% L−glutamine, 10% FBS and 100 U/mL penicillin/streptomycin for 24h in an atmosphere of 5% CO₂:95% air. The glioma cells in the co-culture system were collected and RNA was extracted by using Trizol. Transcriptome sequencing was performed by Novogene and differential gene expression analyses were performed by using the R software package, DESeq2, according to the manufacturer’s instructions.

The bulk RNA sequencing and the clinical data were obtained from the HNSCC dataset of TCGA database, which consisted of 487 samples (https://portal.gdc.carcinoma.gov/). Each probe ID received an annotation with respect to the gene from the corresponding platform annotation profile of the GDC website and the raw matrix data received the quantile normalization and log2 conversion. As smoking is the independent risk factor for HNSCC, we excluded the patients who smoked. HPV groupings for HNSCC patients were based on the HPV condition which provided by the Pan-Cancer Microbiome (cBioportal.org). 72 of the 249 remaining patients were HPV⁺ and 177 were HPV^-. The raw gene expression datasets were processed. Samples with missing data were excluded from this study.

The information for the presence of HPV originated from the pan cancer dataset. These data were registered to the clinical data from TCGA accordingly (7). Correlations between the presence of HPV and the survival of HNSCC patients were studied as described in the following methods.

The scRNA-seq data of patients with HNSCC were obtained from the GEO database (GSE164690) (8). The data were integrated using the SCTransform method of R software Seurat package and analysed using mutual principal component analysis (PCA). Clustering was performed in a resolution of the FindClusters that set to 0.1 and visualized with the uniform manifold approximation and projection (UMAP). Unique molecular modifier (UMI) counts of less than 500, doublets and cells with >5% of their mitochondrial genes were filtered out.

Viral mapping was performed as previously described (7). Raw scRNA fastq data were processed with UMI-tools (https://github.com/CGATOxford/UMI-tools).Virus infected cells were defined with the Viral-Track approach (6). In brief, the sequencing data containing the single cell index were mapped to the virus genome reference database and the status of the single cell was added to the expression matrix so as to correlate with the presence of HPV infection and the corresponding transcriptome.

For the cell line, we calculated mRNA expression-based stemness indices (mRNAsi) using one-class logistic regression machine-learning algorithm (OCLR). For the scRNA data, we employed the trajectory analysis using Monocle 2 package.

GO ontology and KEGG pathway analyses were performed using the WebGestalt webserver. Cell- communication analyses was performed using the CellphoneDB R package (9). The Monocle2 R package was used for the trajectory analysis (10).

Primers were used to amplify the HPV 16 E6 region as follows: E6 forward primer, 5’-TCAGGACCCACAGGAGCG-3′, reverse primer, 5’-CCTCTCACGTCGCAGTAACTGTTG-3′. β-actin, which was used as housekeeping gene, was amplified by using forward primer 5’-TCACCCACACTGTGCCCATCTACGA-3′, reverse primer, 5’-CAGCGGAACCGCTCATTGCCAATGG-3′. HPV⁺ refers to the mean fluorescence intensity values above 1.5-fold negative control tissues.

Immunofluorescence staining was performed as described previously (11). In brief, after embedding with a frozen section compound (Leica, #3801480), biopsies were sectioned into 4-μm on a microtome (Leica CM1950). For immunofluorescent staining, the sections were fixed in pre-cooled methanol (-20°C) for 5 minutes, after washing twice with PBS. The sections were then blocked with PBS/5%BSA/Fcγ blocker at 4°C for 1 hour. Primary antibodies were incubated with sections at 4°C overnight. After washing twice with PBS, the sections were incubated with fluorescent-coupled secondary antibodies for 1 hour at room temperature. After two further washes, the sections were mounted and imaged on an immunofluorescence microscope (Leica DMI3000B).

Antibodies used for staining included: Mouse-anti-human CD3 monoclonal antibody (Invitrogen, # MA1-21454), Rabbit-anti-human CTLA4 polyclonal antibody (Invitrogen, # PA5-115060), Rabbit-anti-human PD1 polyclonal antibody (Invitrogen, # PA5-20350), Rabbit-anti-human CCL4 polyclonal antibody (Invitrogen, #PA5-114961), Rabbit-anti-human FGFR2 polyclonal antibody (Invitrogen, #PA5-14651), A594 goat-anti-rabbit IgG (Affinity, #S0006), A488 donkey-anti-mouse IgG (Invitrogen, #A-21202) and DAPI (Invitrogen, #D21490).

Statistical studies were performed with the R software. Kaplan-Meier (KM) and receiver operating characteristic (ROC) tests were carried out as previously described (12–14), which was done with the “survivor” and “survROC” packages (15). The ‘survminer’ package (16) was used to calculate the optimal cut-off data points. We assessed the prognostic correlates of potentially interesting genes with univariate and multivariate Cox regression correlations. For prognostic correlates, we used the hazard ratios and 95% confidence intervals. Differences between groups were analysed using GraphPad Prism 8.0 software. Comparisons between the two groups were made using the Student’s t-test and P<0.05 was considered to be statistically significant.

HPV groupings for HNSCC patients were based on the HPV condition which provided by the Pan-Cancer Microbiome (cBioportal.org). Kaplan-Meier curves showed shorter survival times for the HPV^- patients when compared to those who were HPV⁺ ( Figure 1A ). As smoking is the independent risk factor for HNSCC, the patients who smoked were excluded so that 72 HPV⁺ and 177 HPV^- cases remained. The survival time was shorter in the HPV^- group compared to HPV⁺ patients after the smokers were excluded ( Figure 1B ). In order to be able to study the effect of HPV on HNSCC tumour cells, the scRNA-seq data were used to study HPV infection in vivo. During data mining, the HNSCC cells were divided into nine clusters, of which clusters 0 and 2 were tumour cells ( Figure 1C , sFigures 1A, B ). Approximately half of the cells in cluster 2 tumour cells were found to be infected with HPV, while more than 75% of cluster 0 tumour cells were infected ( Figure 1D ).

Differentially expressed genes (DEGs) in HPV⁺ and HPV^- tumour cells and GSEA functional analysis of these genes clearly indicated that HPV^- tumour cells upregulated the genes associated with the IL17 signalling pathway ( Figure 1E , sFigure 1B ). This correlates with increasing evidence that supports the pro-tumourigenic mechanisms of IL17 signalling in the immunosuppressive endogenous microenvironment within the growing tumour (17). To further investigate the intrinsic effects of HPV on tumour cells, three different types of HPV⁺ or HPV^- tumour cells were cultured. Transcriptomic analyses were performed on HPV⁺ and HPV^- tumour cells obtained from these experiments. PCA analysis revealed significant differences between HPV⁺ and HPV^- tumour cells ( Figure 1F ). Gene expression analysis identified 1341 DEGs between HPV⁺ and HPV^- tumour cells and 90 of these overlapped with the DEGs found between HPV⁺ and HPV^- cells in vivo ( Figure 1G ). A heatmap using these genes with a log fold change of more than 1 was constructed ( Figure 1H ). Although the presence of HPV in cells can lead to significant differences in gene expression, the data obtained from the cell lines and in vivo also show significant variations, suggesting the importance of the tumour microenvironment.

The stemness of tumour cells is closely related to the malignancy and recurrence of tumours (18). In the transcriptomic data from cell lines, we found a high level of stemness in HPV^- tumour cells, despite the lack of significance in the data due to the small sample size ( Figure 2A ). We therefore studied the level of stemness of tumour cells using scRNA data, in which we can identify genes that are in a transitional state between two different cell fates. Trajectory analysis was performed on HPV⁺ and HPV^- tumour cells, and we identified nine trajectory states in HNSCC tumour cells, with states 1 and 9 pointing to early stages of pseudo-time and therefore correlating with stemness of the tumour ( Figures 2B, C ). HPV^- tumour cells showed a high proportion of cancer stem cells (state1) compared to HPV⁺ cells ( Figure 2D ). While the expression seen in cluster 2 of both HPV⁺ and HPV^- tumour cells was upregulated in the later states, a higher proportion of HPV⁺ cells in the cluster 0 displayed a developed manner compared to the HPV^- cells ( Figure 2E ). This indicated an enhanced level of stemness in HPV^- tumour cells. DEGs of the early state (state1) showed a reduction in cytokine and chemokine signalling pathways ( Figures 2F, G ), while the HPV⁺ state 1 tumour cells showed a reduced level of IL17 signalling. This was similar to the DEGs of HPV⁺ cluster 0 cells ( Figure 1E , sFigures 2A, B ).

To further investigate the genes associated with stemness in the tumour cell lines, we compared the stemness-related genes that were differentially expressed in the HPV⁺ tumour cell lines cultured in vitro and found that 134 of them overlapped with the state 1 stemness genes ( Figure 2H ). A heatmap was constructed using the genes with a fold change of more than 1, in which we found that many DEGs in either HPV⁺ or HPV^- tumour cells remained in the same category ( Figure 2I ). Both in vitro and in vivo data implied a tendency of high stemness in the HPV^- tumour cells.

As mentioned, we studied the level of stemness of HNSCC tumour cells using both cell lines and scRNA data. We filtered then the DEGs of the low stemness cells using the TCGA bulk-seq data. 34 risk negatively-related stemness genes were filtered for a univariate Cox regression study and these were found to significantly correlate with a reduced risk of HNSCC (p < 0.05; Figure 3A ). By using the random survival forest algorithm, which is a powerful non-parametric method that can construct predictive models with time-to-event outcome, the top significant 6 genes were screened. These were ATP1A1, KMT2E, RAD23A, TRR, WASF2 and DNAJC9 ( Figure 3B ). The HNSCC patients in the TCGA dataset were then divided into groups with high- or low-risk, based on the expression levels of the gene signature. Kaplan–Meier curves showed that patients in the high-risk group had longer survival times than those in the low-risk group ( Figure 3C ). ROC curves obtained from the cases of HNSCC were plotted by estimating the predictive power of the genetic characteristics, and AUCs of 0.913 and 0.984 were obtained for 3- and 5-year survival rates, respectively ( Figure 3D ). The expression of the 6 signature genes were at a low level in state 1, but at a high level in the HPV⁺ group ( Figures 3E, F ). This was observed for both the scRNA and TCGA data ( Figures 3F, G ). Up to this point, we defined the stemness-based low and high risk TCGA samples using signature genes and compared the level of these signature genes in the HPV⁺ tumour cells and TCGA HNSCC cohort. Both data showed that the HPV situation is associated with the level of stemness-based signature genes.

We also performed a similar analysis on the positively-related stemness genes and found that the majority were correlated with an increased risk of HNSCC ( sFigures 3A, C ). 3 of these (ALDOA, CTNNA1 and TMBIM6) were significantly different, although the AUCs for the 3- and 5-year survival rates were only 0.705 and 0.76, respectively ( sFigures 3B, D ). These genes were highly expressed in state 1 but were similar in both the HPV⁺ and HPV – groups ( sFigures 3E, F ).

In agreement with the previous studies, CIBERSORT estimated that the infiltration of immune cell subsets showed a greater proportion of CD8⁺ T cells in the HPV⁺ HNSCC patients when compared to those who were HPV^- ( Figure 4A ). Also, the immune checkpoint genes, including CTLA4, LAG3, PDCD1, HAVCR2, PDCD1LG2 and TIGIT, were found to be higher in the HPV⁺ HNSCC patients ( Figure 4B ). Similar results were also observed with the scRNA data ( sFigures 4A, B ). To confirm the checkpoint protein expression, we performed immunostaining on the surgical sections from the HPV⁺ and HPV^- HNSCC patients ( Figure 4C ). We observed significantly higher levels of CTLA4 and PD1 in the HPV⁺ HNSCC samples and a higher proportion of dysfunctional CD8⁺ T cells were found in HPV⁺ HNSCC patients ( Figure 4C ). Risk scores of HNSCC data and subtypes were evaluated according to the 6-gene signature that was established to be stemness-related genes. CIBERSORT estimated infiltration showed that the low-risk group of HNSCC patients had a higher proportion of CD8⁺ T cells ( Figure 4D ). In addition, the gene profiles showed that the immune checkpoint genes, CTLA4, LAG3, PDCD1 and TIGIT, were significantly increased in these patients ( Figure 4E ), suggesting a correlation between the stemness-related genes and the immune microenvironment.

To investigate the effect of HPV⁺ or HPV^- tumour cells on immune cells, we constructed an intercellular communication network and analysed the interactions between tumour cells and immune cells ( Figure 5A ). HPV⁺ tumour cells displayed 25 extra cell-cell interactions when compared to HPV^- cells, and of these BMPR1A showed the most promise as a favorable gene for the prognosis of HPV⁺ HNSCC patients. However, this was not the case for the HPV^- HNSCC patients ( Figures 5B, C ). When the B cells, T cells and macrophages interacted with tumour cells, it was found that interactions with the FGFR2 gene were upregulated, which was not detected in the endothelial cell-tumour interactions ( Figures 5D, E , sFigures 5A-C ). The FGFR2 gene could also be used as a possible candidate for the prognosis of HPV⁺ HNSCC ( Figure 5F ). To confirm the FGFR2 expression, we performed immunostaining on the surgical sections from the HPV⁺ and HPV^- HNSCC patients and observed a significantly higher level of FGFR2 in the samples of those who were positive ( Figure 5G ).

More importantly, in this study, HPV⁺ tumour cells appeared to have a significantly higher level of cellular communication with CD8⁺ T cells, which correlated with the high proportion of tumour infiltrated CD8⁺ T cells in HPV⁺ HNSCC patients ( sFigure 5C ). The chemokine, CCL4, and the T cells activator, CD83, were also found to be upregulated in the HPV⁺ tumour-T cell interactions ( Figure 6A ). These were also genes that were upregulated in the HPV⁺ group and were significantly associated with a favorable prognosis in the HPV⁺ HNSCC patients ( Figures 6B, C ). Gene correlation assays showed that both CCL4 and CD83 not only correlated with the amount of recruited T cells, but also had a substantial correlation with the immune checkpoint genes such as CTLA4, LAG3, PDCD1, HAVCR2 and TIGIT ( Figure 6D ). Furthermore, the expression of CCL4 was validated and quantified using immunostaining, and this showed a significantly higher level of CCL4 in the HPV⁺ HNSCC samples ( Figure 6E, F ). The data indicated an important role of the interactions between HPV⁺ tumour cells and infiltrated T cells, which might account for the different prognoses observed in HPV⁺ and HPV^- HNSCC patients.