The UM-SCC-104 (referred to as SCC104) cell line was obtained from Sigma-Aldrich, and the UPCI : SCC152 cell line (referred to as SCC152) was obtained from ATCC. Both SCC104 and SCC152 cell lines have been reported to be HPV-16 positive (24, 25). SCC25 and SCC47 cell lines were purchased from ATCC and Millipore Sigma, respectively. Cell lines UM-SCC-11B, UM-SCC-74A, UM-SCC-29, UM-SCC-23, and UM-SCC-103 were obtained from Dr. Thomas Carey (University of Michigan). HSC-3 and CAL-27 cell lines were generously provided by Dr. Manish Bais (Boston University). All cell lines were grown and maintained in high-glutamine DMEM or DMED/F12 as recommended, with the following supplements: 10% FBS, 1% nonessential amino acids, and antibiotics. Other cell lines used in this study have been described before in Gluck et al. 2019 (26). The identities of the cell lines utilized in this study were confirmed via short tandem repeat profiling through services offered by Genetica. All cell lines were tested by the eMycoPlus Mycoplasma PCR Detection Kit (BulldogBio) to ensure that they were bereft of any mycoplasma infection.

Lentivirus-mediated depletion of p63 in SCC104 and SCC152 cells was performed using the pGIPZ system. GIPZ lentiviral shRNAs (clone IDs V2LHS_24248 [sh1] and V2LHS_24250 [sh2]) targeting TP63 were obtained from and virus was generated with the help of Gene Modulation Services Shared Core at Roswell Park Comprehensive Cancer Center. Viral infection and selection with Puromycin was performed as described before (27).

Protein extracts were prepared according to previously published protocol (27). Briefly 5 μL of protein lysates were loaded onto SDS-polyacrylamide gels and transferred to Immun-Blot PVDF membranes (Bio-Rad Laboratories). After blocking in 5% milk, the membranes were incubated in primary antibodies against the following: p63 (4A4, 1:20,000), ΔNp63 (E6Q3O; Cell Signaling Technology, 1:5000), ITGB1 (Proteintech, 1:10,000), ITGB4 (Proteintech, 1:10,000), cMYC (Santa Cruz Biotechnology, 1:5000), AKT1 (Proteintech, 1:10,000), mTOR (Proteintech, 1:10,000), Raptor (Proteintech, 1:10,000), S6 (Cell Signaling Technology, 1:5000), and p-S6 (Cell Signaling Technology, 1:5000). The MAB374 antibody (EMD Millipore) was used to detect GAPDH as a loading control at 1:20,000 dilution. HRP-conjugated secondary antibodies corresponding to the primary antibody host were incubated with each blot. Unbound antibodies were washed off in 0.05% Tween-20 in Tris-buffered saline. The LumiGLO peroxidase chemiluminescent substrate kit (SeraCare) was used to detect antibody-labeled proteins, and membranes were imaged using the Bio-Rad ChemiDoc imaging system.

The iDeal ChIP-seq kit for transcription factors (C01010055; Diagenode) or for histones (C01010051; Diagenode) and the associated protocols were used to perform ChIP-seq. SCC104 and SCC152 cells were grown to ~90% confluency and cross-linked in the supplied fixation buffer supplemented with 0.5% formaldehyde for 10 min. Lysates from the fixed cells were subsequently sonicated with a Diagenode Bioruptor to obtain sheared chromatin with an approximate fragment length of 150–400 bp. The ChIPs for p63 were carried out using 2 μg of p63 4A4 antibody (Santa Cruz Biotechnology) and 2 μg of ΔNp63-1.1 antibody (28). After cross-link reversal, proteinase-K/RNase A treatment, and DNA purification, libraries were prepared using the ThruPLEX DNA-seq kit (Rubicon Genomics). ChIP DNA and input controls were then subjected to 50-bp single-end sequencing on an Illumina HiSeq 2500, which resulted in 15–25 million reads per sample.

The raw ChIP-seq reads from all experiments were mapped to the Homo sapiens genome (hg19 build) using Bowtie v1.1.1 with the parameter m=1 to remove all reads mapping to multiple genomic loci (29). Peak calling was then performed using MACS2 v2.1.0 with a minimum FDR cutoff of 0.05 and sequenced Input used as control for each experiment, and resultant peaks were matched to the nearest gene using GREAT analysis with default settings (30, 31). For visualization of ChIP peaks, the package deeptools v3.3.2 was used to preprocess bam files to generate bigwig files which were then uploaded to IGV (32). Adobe illustrator was used for final image processing. Peak summits determined by MACS2 v2.1.0 were used as input to HOMER’s findMotifsGenome.pl program with the parameter “-size 200” (33).

Total RNA from cell lines was extracted using a Direct-zol RNA miniprep kit (Zymo Research). The extracted RNA was snap-frozen on dry ice and stored at −80°CC until library preparation. For each RNA sample, cDNA libraries were prepared using the TrueSeq RNA sample preparation kit (Illumina) and were then 50-bp single-end sequenced or paired-end sequenced on an Illumina HiSeq 2500. Quality control metrics were performed on raw sequencing reads using the FASTQC v0.11.9 application.

Reads were mapped to the appropriate reference genome, GRCh38/hg19 build, with HISAT2 v2.1.0 (34). Reads aligning to the reference genome were quantified with featureCounts v1.5.3 to generate a matrix of raw counts, which was then processed in R, to generate normalized expression values in transcripts per million according to the method proposed by Wagner et al. (35). Differential gene expression analysis comparing control to p63 knockdown was carried out using DESeq2 v1.24.0 (26). DEGS with an FDR value of ≤ 0.1 were considered statistically significant.

Total RNA from SCC104 and SCC152 knockdown cell lines was extracted using a Direct-zol RNA miniprep kit (Zymo Research). RNA was reverse transcribed with the Bio-Rad iScript cDNA synthesis kit according to the manufacturer’s instructions. The resulting cDNA was used for qPCR with Bio-Rad iQ SYBR green Supermix. A list of the qRT-PCR primers can be found in Supplementary Table 7 .

RNA-seq data from patient samples were obtained from GEO (GSE122512, GSE112026, GSE74927, and GSE72536) (11, 36–38). HPV+ tumors were assigned based on data presented in the original paper of each dataset. Alignment and quantification of counts for each dataset were performed as indicated by the original study. TCGA RNA-seq expression and CNA datasets were downloaded from cBioPortal (39, 40). Briefly, RNA-seq counts were extracted and normalized using the median-ratio method (DESeq2 v.1.24.0 [75]) and subsequently transformed to transcript per million values (35). For GSE112026, RSEM values were utilized for transcript quantification. HPV+ tumors were segregated into high and low p63 expression groups based on the median p63 expression level calculated from the RNA-seq data.

H3K27Ac ChIP-seq data from SCC104 and SCC152 cells were aligned to the human genome as described above. Narrow peaks were called using MACS2 v2.1.0 using the following parameters: -p 0.01, -nomodel, -extsize 150. The resulting narrowPeaks files were converted to gff format and used as inputs for the ROSE (rank order of super-enhancers) algorithm, which was run using default parameters along with appropriate input controls to generate typical enhancer and SE lists (41, 42).

The deepTools package was utilized to generate a signal matrix of histone modifications H3K27Ac, H3K4Me1, and H3K4Me3. The fluff python package was then utilized to generate heat maps showing the resulting signal of the histone modifications around a 2-kb window centered at each p63 ChIP-seq peak summit. The resulting histone signal enrichment was subjected to k-means clustering (k=3) (43).

The CEAS tool was used to annotate p63 ChIP-seq peaks to the nearest genomic feature of the hg19 genome assembly (44). The promoter region was considered up to 1,000 bp from a transcriptional start site, and the proximal enhancer was considered from 1,000 to 3,000 bp away. Any binding within a gene was considered intragenic, whereas any binding site greater than 3,000 bp upstream or downstream was considered distal intergenic.

To determine the top enriched DNA binding motifs of TFs found within nucleosome-free regions of SCC104 and SCC152 SEs, nucleosome-free regions were first determined using the HOMER findPeaks tool with the -nfr flag (33). The AME tool was used to determine enriched motifs found within the HOCOMOCO Human (v11 CORE) database. Motifs were ranked according to p value.

The GREAT tool was used to annotate binding loci to the nearest gene (31). Identified genes were then subjected to KEGG pathway analysis utilizing the DAVID functional annotation tool (45–47). For RNA-seq data, DEGs were subjected to both KEGG analysis utilizing the DAVID functional annotation tool and canonical pathway analysis by gene set enrichment analysis (48).

Statistical analyses were performed using R, a free software environment for statistical computing and graphics. A Shapiro-Wilk test was performed to check the normality of data, and then either a student’s t test or Wilcoxon signed-rank test was performed according to whether the data were normally distributed. A p value lower than 0.05 was considered significant.

Although a broad oncogenic role of p63 in HNSCC has been reported (21, 49, 50), its specific role in the HPV+ subtypes has not been fully explored. Thus, we first examined three independent RNA-seq datasets of HPV+ HNSCC tumors (GEO datasets GSE112026, GSE74927, and GSE72536) and observed a gradient in the pattern of p63 mRNA expression ( Supplementary Figures 1A–C ). Segregation of HPV+ HNSCC tumors according to median p63 expression revealed distinct p63^high and p63^low subtypes. This distinction was in agreement with previous unsupervised gene expression clustering analyses performed on HPV+ HNSCC tumors that had identified subtypes with distinct gene expression patterns, including different p63 levels (11, 14, 51).

To investigate the functional relevance of p63 in HPV+ HNSCC, we next focused on data from 67 HPV+ tumors that are available in The Cancer Genome Atlas (TCGA) patient datasets. Exploration of p63 expression across the HPV+ HNSCC tumors in this dataset revealed a similar pattern of p63 expression as observed in the GEO dataset ( Figure 1A ). Using the RNA-seq data from the TCGA datasets to segregate the HPV+ tumors according to p63^high and p63^low expression, we identified 6,459 differentially expressed genes (DEGs) between these two populations ( Figures 1B, C and Supplementary Table 1 ). We next examined GO biological processes that were enriched with upregulated and downregulated DEGs to identify pathways that are likely influenced by p63. Downregulated DEGs were significantly enriched in pathways involved in viral transcription and inflammatory immune responses, such as NF-κB and tumor necrosis factor signaling ( Figure 1D ). Upregulated DEGs were significantly enriched in pathways associated with cell adhesion and keratinization—processes linked to p63 ( Figure 1E ). These findings were of particular interest in lieu of prior HPV+ subtype studies in which gene expression-driven clustering analysis showed differential enrichment of immune response and cell adhesion pathways. Notably, the clustering of tumors according to p63 expression recapitulated the molecularly defined distinct subtypes of HPV+ HNSCC.

To verify that the pattern of p63 expression in tumors matches that in HNSCC cell lines to serve as suitable models for follow-up studies, we examined RNA-seq data generated from 9 HPV+ and 55 HPV− HNSCC cell lines ( Supplementary Figure 2A ) (38). Several of the HPV+ HNSCC cell lines had high p63 expression, and we verified this expression at the protein level by Western blotting with two anti-p63 antibodies. Similar to previous reports (23), four of the five HPV+ HNSCC cell lines consistently showed high levels of p63 protein expression, specifically the ΔNp63 isoform ( Supplementary Figure 2B ). Of these, we chose the SCC104 and SCC152 cell lines for follow-up mechanistic studies. These two well-characterized cell lines have been confirmed for HPV-positivity and shown to express viral factors and oncogenes E6 and E7, making them suitable for studies of HPV+ HNSCC in vitro (25). The SCC104 cell line was the primary choice for most of our experiments because of its robust growth and detailed phenotypic characterization compared to that for SCC152 (24); data from the SCC152 cell line were used to corroborate and/or validate the findings.

To identify the global network of p63 target genes, we performed ChIP-seq experiments in SCC104 cells with two anti-p63 antibodies. ChIP-seq of p63 with the widely used 4A4 antibody that recognizes all p63 isoforms identified 18,085 genomic sites, whereas ChIP-seq with a ΔNp63-specific antibody, ΔNp63-1.1, identified 10,028 p63-bound sites ( Figure 2A ); 9,724 sites were identified by both antibodies, which were deemed high-confidence p63 targets and utilized for subsequent analysis ( Supplementary Table 2 ). As expected, analysis of the p63 ChIP-seq peaks using HOMER revealed the consensus p63 motif (p = 1e-7867) as the most highly enriched motif, followed by the p53 motif (p = 1e-5462), which was independently confirmed by using MEME-ChIP ( Figure 2B and Supplementary Figure 3B ). Other enriched motifs were for TFs belonging to the AP-1 family, which cooperates with p63 to regulate target gene expression ( Figure 2B ) (52). The distribution of p63 peaks relative to the transcriptional start sites revealed that, in both cell lines, p63 preferentially targets intragenic and distal regulatory regions, which are likely to act as enhancer sites ( Figure 2C and Supplementary Figure 3C ). A DAVID-based pathway analysis of the genes associated with the top 2,500 p63 ChIP-seq peaks revealed several important pathways, including those deemed important in HPV-associated cancers, such as focal adhesion, p53, and Rap1 signaling pathways ( Figure 2D ) (53, 54). Interestingly, focal adhesion pathways are enriched in the subtype of HPV+ HNSCC tumors with high p63 expression levels (11, 14, 51). In parallel, we performed ChIP-seq of p63 in SCC152 utilizing the 4A4 antibody and identified 26,255 genomic peaks, with the p63 motif as the most enriched ( Supplementary Figures 3A, B and Supplementary Table 2 ). p63 binding-associated genes in SCC152 were enriched in pathways involving MAPK signaling, cell adhesion molecules, and Rap1, which all play a role in HPV infection ( Supplementary Figure 3D ). Of the 9,724 high-confidence p63 binding sites (identified with both antibodies), 5,933 sites were shared between SCC104 and SCC152 cell lines, providing a strong list of bona-fide p63 targets in HPV+ HNSCC.

We next explored the epigenomic landscape of SCC104 and SCC152 by ChIP-seq using histone marks H3K27Ac, H3K4Me1, and H3K4Me3 ( Supplementary Figures 4A , 5A ) (55), which identify gene regulatory features such as active enhancers (H3K27Ac^high and H3K4Me1^high), active promoters (H3K27Ac^high and H3K4Me3^high), and poised enhancers (H3K27Ac^low and H3K4Me1^high) (56).We performed k-means clustering on the histone marks centering around each H3K27Ac peak as described in Gluck et al. (27), which identified three distinct clusters of regulatory elements. In SCC104 cells, clusters 1 and 3 represented active enhancers, and associated genes were enriched in pathways such as mRNA processing, differentiation, and focal adhesion ( Supplementary Figures 4B ). Cluster 2 represented active promoters, and associated genes were enriched in viral processes and cell motility ( Supplementary Figures 4A, B ). To identify which TFs may regulate enhancer expression in these clusters, we performed a motif analysis and found enrichment of ZNF, IRF, KLF, and Ets family motifs as well as E2F motifs ( Supplementary Figure 4C ).

Similar results were obtained in clustering analysis of SCC152 cells, where clusters 1 and 3 also represented active enhancers, and genes associated with these sites were enriched in focal adhesion, cell junction assembly, and Notch signaling pathways ( Supplementary Figure 5B ). Cluster 2 was similarly associated with active promoters, and associated genes were enriched in mRNA processing and cell-cell adhesion pathways ( Supplementary Figures 5A, B ). Motif analysis of these regions in SCC152 also showed enrichment of E2F motifs within all identified clusters ( Supplementary Figure 5C ). The enrichment of E2F motifs across various gene regulatory elements in both HPV+ HNSCC cell lines is interesting and likely to be relevant given the known interaction of E2F TFs and HPV E7 and its effects on downstream pathways in HPV+ disease.

Multi-cluster enhancers, often referred to as super-enhancers (SEs), are associated with H3K27Ac^high marks and are densely occupied by key TFs (42, 57). These SEs are often associated with cell identity and lineage-driving genes and oncogenes, and in the context of HPV, viral oncoproteins have been found to play a central role in their activation (58, 59). It is postulated that SEs at the site of HPV integration likely upregulate the expression of HPV E6/E7, leading to activation of other SEs that facilitate disease progression (60). By applying the ROSE algorithm to the H3K27Ac ChIP-seq data, we identified 528 SEs in SCC104 and 317 SEs in SCC152 ( Figure 3A and Supplementary Figure 6A , Supplementary Table 3 ) (42, 61). As expected, several SEs were associated with master TFs, including p63, but we also observed SEs associated with genes important to HPV infection and HPV-induced carcinogenesis, such as WNT7A, ITGA2, and NOTCH1. ( Figure 3B and Supplementary Figure 6B ; Supplementary Table 3 ) (62, 63). Motifs for TFs with known roles in HNSCC, including FOSL1, E2F1, and E2F7 (10, 64), were also enriched at SCC104 SEs ( Figure 3C and Supplementary Figure 6C ). We found significant enrichment of the p63 motif in SCC104 SEs, suggesting p63 regulates the transcription of many SE-associated genes in HPV+ HNSCC ( Figure 3C and Supplementary Figure 6C ). The enrichment of the p63 motif in SEs is functionally relevant, because most (369/540) were occupied by p63 according to the ChIP-seq results from SCC104 cells.

Analysis of SE-marked genes revealed notable enrichment of pathways associated with cancer, including those for focal adhesion and those involving proteoglycans ( Figure 3D ). Interestingly, we found that SDC1, which encodes protein syndecan-1 and is a direct p63 target, was associated with SEs ( Supplementary Table 3 ). This is notable because syndecan-1 is the most abundant heparan sulfate proteoglycan in keratinocytes and serves as the primary HPV attachment receptor (65). HPV infection is known to affect the expression of genes involved with cell adhesion and cell motility, and several of the implicated genes were found in our data, including LAMA5, ITGA3, ITGA6, and LAMC2 (11, 36, 54, 66). SE-marked genes in SCC152 were similarly enriched in pathways important for cancer, including the Hippo signaling pathway, implicated in HPV-associated oropharyngeal SCC ( Supplementary Figure 6D ) (67).

We next explored the epigenomic state (as defined by histone marks) of the gene regulatory regions bound by p63 that were identified by ChIP-seq. For this purpose, we performed k-means clustering of the three histone marks, which again identified three distinct clusters ( Supplementary Figures 7A , 8A ). In SCC104 cells, cluster 1 represented active promoters, and genes associated with these sites were enriched in pathways for apoptosis, cell adhesion, and motility ( Supplementary Figure 7B ). Clusters 2 and 3 represented active enhancers, and the corresponding genes were associated with Notch, protein kinase B, and Rho signaling ( Supplementary Figures 7A, B ). Unsurprisingly, we observed enrichment of AP-1 motifs in p63-bound enhancer regions regulating such processes as signaling and differentiation, which has been reported in keratinocytes and breast cancer (68). Interestingly, we observed enrichment of E2F motifs in cluster 1, suggesting the interaction between p63 and E2F at active promoter regions regulates cellular movement and apoptosis ( Supplementary Figure 7C ).

Similar analyses in SCC152 cells showed identical patterns of clustering of gene regulatory elements, with cluster 1 representing active promoters and clusters 2 and 3 representing active enhancers ( Supplementary Figure 8A ). Cluster 1 showed enrichment of pathways associated with apoptosis and cellular organization ( Supplementary Figure 8B ). Unlike that in SCC104 cells, we did not find enrichment of E2F motifs in these regions but instead saw enrichment of FOX motifs ( Supplementary Figure 8C ). Genes found in clusters 2 and 3 were associated with tissue development, Notch signaling, and cell adhesion pathways ( Supplementary Figure 8B ). Similar to that for SCC104 cells, there was enrichment of ZNF, IRF, and ETS family motifs; however, there was no enrichment of AP-1 motifs ( Supplementary Figure 8C ). These findings suggest that p63 actively regulates genes and pathways considered important in HPV-induced carcinogenesis and that p63 may interact with cellular E2Fs at gene regulatory regions.

To identify p63 targets, we performed RNA-seq to profile the global transcriptomic changes resulting from loss of p63 expression. For these experiments, we stably knocked down p63 in SCC104 and SCC152 cells with two independent lentiviral mediated shRNAs. Western blotting confirmed that both shRNAs significantly reduced p63 expression, with sh2 showing markedly greater knockdown ( Figure 4A ). Loss of p63 substantially altered the transcriptomic landscape, with 6,607 and 5,809 genes exhibiting statistically meaningful changes in expression in SCC104 and SCC152 cells, respectively ( Supplementary Table 4 ). Of these, 2,189 and 1,461 DEGs showed statistically meaningful (p _adj < 0.1) changes with both p63 shRNAs in SCC104 and SCC152 cells, respectively ( Figure 4B and Supplementary Figure 9A ; Supplementary Table 4 ).

To explore pathways affected by loss of p63, we focused on a select group of 615 DEGs (≥log2 fold change of 1) common to both shRNA knockdowns in SCC104 cells. KEGG analysis of these DEGs showed enrichment of many pathways important in HPV-associated carcinogenesis, such as WNT, MAPK, and PI3K-Akt signaling ( Figure 4C and Supplementary Figure 9B ) (21, 53, 62, 69, 70). We also performed gene set enrichment analysis of canonical pathways associated with these DEGs, which identified cell cycle and retinoblastoma gene in cancer categories ( Figure 4D ). Genes involved in the pRb signaling pathway, including E2F1 and CCNA2, were significantly upregulated by the loss of p63, suggesting that p63 suppresses E2F-induced transcription and cell cycle activation.

The enrichment of pathways associated with HPV after p63 knockdown prompted us to look closer into the known HPV infection pathway. Of the 324 genes in the KEGG “human papillomavirus infection” pathway, 62 were differentially expressed upon p63 knockdown ( Supplementary Figure 10A ). HPV infection affects cell cycle regulation, which was also the case with p63 knockdown. Several factors involved in the cell cycle, including E2F1, RBL1, and cyclin A2, were upregulated upon loss of p63 ( Supplementary Figure 10A ). We also observed high enrichment of genes involved in focal adhesion and WNT and PI3K signaling ( Supplementary Figure 10B ).

To delineate our high-stringency p63-driven gene signature, we combined the gene signatures we identified from TCGA tumor data and from the p63 expression-defined cell lines ( Supplementary Figure 11A and Supplementary Table 1 ). To identify which genes are most reliant on p63 expression, we utilized our sh2 data from both cell lines to generate our HPV+ cell line-based signatures ( Supplementary Figure 11A and Supplementary Table 5 ). These analyses identified 1,052 genes shared between the TCGA and SCC104 datasets and 827 genes shared between the TCGA and SCC152 datasets ( Supplementary Figure 11B ). Then, to identify genes which were directly regulated by p63, we incorporated our p63 ChIP-seq data, which revealed 498 and 574 genes that were directly bound by p63 in SCC104 and SCC152 cells, respectively ( Supplementary Figure 11B ). We filtered all genes that were common between analyses to generate our combined p63 signature of 420 genes ( Supplementary Table 5 and Supplementary Figure 11B ). Finally, to identify genes of potential importance, we merged this gene signature with our super-enhancer landscape, which revealed 55 super-enhancer-associated genes ( Supplementary Figure 11B and Supplementary Table 5 ). These analyses provided a p63-driven gene expression signature in HPV+ HNSCC that is relevant in both cancer and HPV contexts for follow-up studies.

To explore how p63 regulates subtype-specific gene expression in HPV+ HNSCC, we compared our combined p63-driven signature of 420 genes with the aforementioned subtype-specific signatures. Our previous analyses of transcriptomic changes upon p63 knockdown revealed enrichment of cell adhesion and keratinization pathways, like the HPV-KRT subtype defined by Keck et al. (14) and Zhang et al. (11). Keck et al. (14) found HPV-KRT tumors upregulate genes involved in hypoxia, cell adhesion, and HER signaling as well as epithelial-associated genes, whereas genes involved in immune response and mesenchymal-associated genes are downregulated. Our p63^high signature displayed a similar pattern of upregulated gene expression (compared to expression in p63^low samples), supporting the notion that high p63 expression is a defining aspect of the HPV-KRT subtype ( Figure 5A ). Interestingly, we also found that 36 of our 420 p63 signature genes, including MAOA, SLC2A1, COL17A1, and KRT16, were associated with the reported HPV-KRT subtype signature ( Supplementary Table 6 ). Conversely, the p63^low expression signature was enriched with genes associated with the immune response and mesenchymal tissues ( Figure 5A ).

The HPV-KRT signature defined by Zhang et al. (11) had patterns of expression and pathway enrichment similar to the signature defined by Keck et al. (14). Accordingly, the gene signatures within our p63^high subgroups also showed comparable upregulation of genes associated with keratinization, redox processes, and cell adhesion previously seen by Zhang et al. (11) (relative to expression in the p63^low samples) ( Figure 5B ). The HPV-KRT-like subtype signature defined by Zhang et al. (11) was apparent in our p63-based DEGs (49 DEGs, including TP63, MAOA, PPARD, and KRT16) ( Figure 5B and Supplementary Table 6 ), whereas the p63^low samples had upregulated genes associated with the immune response and mesenchymal differentiation ( Figure 5B ). To validate the differences in expression observed between the p63^high and p63^low groups and in our p63-driven signature, we performed qRT-PCR for several genes identified in each subtype signature, including KRT16, SLC2A1, and MAOA, in the HPV+ cell line with p63 knockdowns, confirming our RNA-seq findings ( Figure 5C ). Altogether, these findings support a role for p63 in subtype-specific gene expression within HPV+ HNSCC and suggest that p63 directs the specific gene expression profiles that were discovered by hierarchal clustering of global gene expression in HPV+ HNSCC.

HPV+ HNSCC subtypes also have differences in copy number alteration (CNA) patterns and mutation frequencies, specifically in PIK3CA. Furthermore, HPV E6 and E7 oncoproteins are implicated in regulating the PI3K/AKT/mTOR network in cancer cells under both normoxic and hypoxic conditions, specifically by regulating AKT, a main effector of both PI3K and mTORC1 signaling (71). We utilized cBioPortal to examine the genomic properties of our p63 subgroups within the TCGA HNSCC tumor dataset, and compared our findings with those by Keck et al. (14) and Zhang et al. (11). Zhang et al. (11) found that HPV-KRT tumors had more amplifications in chr3q than HPV-IMU, the subtype with a strong immune response. Notably, TP63 is located on chr3q, providing more evidence toward p63 as a driver of this subtype. In accordance, p63^high samples had significantly more gains in chr3q than p63^low samples ( Figure 6A ). We also explored the frequency of Copy Number Alterations (CNA), specifically amplification events, of genes found within chr3q and found that TP63, SOX2, and PIK3CA had significantly more amplification events in the p63^high group (p = 0.0483, 0.0141, and 2.242e-3, respectively) ( Figure 6B ).

PI3K signaling plays a role in tumorigenesis, and activating mutations in PIK3CA have been found in various cancer types (72). Given the relationship between p63 and activated PIK3CA, we sought to explore if a p63/PI3K signaling axis is active within HPV+ HNSCC (73). First, we examined the mutational status of PIK3CA in the p63^high/p63^low TCGA samples. Strikingly, we found that PIK3CA was one of the most highly mutated genes (p = 9.522e-3) in p63^high samples, with 40% of tumors harboring PIK3CA mutations ( Figure 6C ), similar to what has been reported for the HPV-KRT subtype. Combined with the increased copy number for PIK3CA in p63^high samples, these data strongly suggest that PI3K activity is upregulated in p63^high tumors. The clinical data associated with these p63^high and p63^low groups revealed a significant difference in the Winter hypoxia score, with the p63^high tumors having significantly higher hypoxia scores (p = 0.0121) ( Figure 6D ). Hypoxia stimulates AKT signaling and downregulates E6/E7 expression, inducing reversible growth arrest that is a potential pathway by which HPV+ cancers, such as HPV-KRT tumors, evade the immune response and become resistant to therapy (9, 74). The p63 signature included several PI3K pathway members that were affected by p63 knockdown in HPV+ HNSCC cell lines and were part of the p63 expression-based DEGs from the TCGA datasets ( Figure 6E ). These data point to a likely role of p63 in regulating the PI3K signaling pathway within HPV+ HNSCC.

To systematically follow up on the p63-PI3K link, we examined the p63 knockdown datasets and found that many of the DEGs in the KEGG PI3K signaling pathway are known p63 targets, and importantly, some of these DEGs were associated with SEs ( Figure 7A ). Although the expression of many of the genes associated with PI3K signaling was decreased by p63 knockdown, the expression of AKT1 was modestly increased. We suspect this is due to the loss of repressive effects of other p63 targets. In support of this, qRT-PCR analysis showed that expression of PTEN, a negative regulator of phosphorylation and activation of AKT1, is decreased upon p63 knockdown in SCC104 and SCC152 cells ( Figure 7B ). By contrast, the expression of PIK3CA, the gene encoding the catalytic subunit of the PI3K complex, was decreased upon knockdown of p63 ( Figure 7B ). We also examined the mTOR pathway, which is downstream of PI3K, not only because it is dysregulated in HPV+ HNSCC but also because mTOR inhibitors show promising anticancer effects in HPV+ HNSCC mouse models (75). Similar to what we found for the PI3K pathway, we saw a general loss of downstream expression of mTORC1 targets upon p63 knockdown, suggesting a downregulation of mTORC1 signaling.

Finally, we performed Western blotting for some of the proteins involved in PI3K and mTORC1 signaling to confirm the transcriptomic findings. We found that the protein levels of upstream regulators of PI3K (ITGB1 and ITGB4) were decreased in SCC104 and SCC152 cells with p63 knockdown ( Figure 7C ), in line with the RNA-seq results. A downstream target of PI3K signaling, MYC, was also downregulated by p63 knockdown, suggesting a broader dampening of PI3K signaling upon loss of p63. However, the levels of phosphorylated AKT1 (pAKT1) were elevated in SCC104 cells with p63 knockdown, suggesting the PI3K pathway is activated despite the downregulated mRNA expression of some of the key signaling components. This suggests that a complex regulatory network controls AKT activation and that other regulators (aside from p63) indirectly play a role. Western blotting also revealed mixed results for components in the mTORC1 pathway ( Figure 7D ). Whereas mTOR expression was decreased in SCC104 cells with p63 knockdown, protein levels of Raptor, a key mTOR-interacting partner, was increased ( Figure 7D ), suggesting a more complex regulatory network for mTOR activation. However, the phosphorylation of S6 ribosomal protein (pS6), a key downstream mediator of mTOR activation was reduced in SCC104, indicating that overall mTORC1 signaling was dampened by the loss of p63 ( Figure 7D ). Taken together, these results point to a complex signaling network by which p63 regulates the PI3K and mTORC1 signaling pathway in HPV+ HNSCC, which may have implications for potential therapeutics.