p63 Directs Subtype-Specific Gene Expression in HPV+ Head and Neck Squamous Cell Carcinoma

The complex heterogeneity of head and neck squamous cell carcinoma (HNSCC) reflects a diverse underlying etiology. This heterogeneity is also apparent within Human Papillomavirus-positive (HPV+) HNSCC subtypes, which have distinct gene expression profiles and patient outcomes. One aggressive HPV+ HNSCC subtype is characterized by elevated expression of genes involved in keratinization, a process regulated by the oncogenic transcription factor ΔNp63. Furthermore, the human TP63 gene locus is a frequent HPV integration site and HPV oncoproteins drive ΔNp63 expression, suggesting an unexplored functional link between ΔNp63 and HPV+ HNSCC. Here we show that HPV+ HNSCCs can be molecularly stratified according to ΔNp63 expression levels and derive a ΔNp63-associated gene signature profile for such tumors. We leveraged RNA-seq data from p63 knockdown cells and ChIP-seq data for p63 and histone marks from two ΔNp63high HPV+ HNSCC cell lines to identify an epigenetically refined ΔNp63 cistrome. Our integrated analyses reveal crucial ΔNp63-bound super-enhancers likely to mediate HPV+ HNSCC subtype-specific gene expression that is anchored, in part, by the PI3K-mTOR pathway. These findings implicate ΔNp63 as a key regulator of essential oncogenic pathways in a subtype of HPV+ HNSCC that can be exploited as a biomarker for patient stratification and treatment choices.


INTRODUCTION
Head and neck squamous cell carcinoma (HNSCC) is the sixth most common cancer worldwide and has a 5 year-mortality rate of nearly 50%, making it a leading cause of cancer-related death (1). HPV infection has overtaken alcohol and tobacco consumption as the predominant risk factor in the majority of newly diagnosed HNSCC cases (2)(3)(4). Intriguingly, patients with HPV+ HNSCC have better overall survival and progression-free survival than those with HPV− HNSCC. However, current treatment options for both HPV+ and HPV− HNSCCs consist of standard care regimens of chemoradiotherapy concurrent with cisplatin (2,5). More attention has been given toward deescalation of therapy for HPV+ HNSCC, and a clearer understanding of the underlying biology may aid in identifying patients who would benefit from new treatment modalities (6).
HPV infection of epithelial cells, primarily in the oropharynx, can result in the integration of the viral genome into the host genome, leading to dysregulated expression of viral and cellular oncoproteins and carcinogenesis (7). HPV E6 and E7 oncoproteins are the primary drivers of the pathogenesis of HPV and function by degrading tumor suppressor p53 and retinoblastoma protein (pRb), respectively, leading to activation of the cell cycle-promoting E2F family of transcription factors (TFs) (8,9). Integration of the HPV genome alters the expression and DNA methylation profiles of a broad range of host genes (10). Although the HPV genome can integrate throughout the human genome, it occurs with a higher incidence in some regions, including the 3q region surrounding TP63 (11)(12)(13)(14). TP63 encodes p63, a member of the p53 family of transcription factors, which plays an essential role in the development and maintenance of the stratified squamous epithelium (15)(16)(17)(18). DNp63a is the most prevalent p63 isoform in tissues of epithelial origin and acts predominantly as an oncogene in several cancers, including HNSCC, while TAp63 has much more restricted expression and shows tumor-suppressor features (18)(19)(20)(21)(22). Ectopic expression of HPV oncoproteins in human keratinocytes leads to upregulation of p63 at both the mRNA and protein levels (23). Conversely, silencing of E6/E7 expression in HPV+ cell lines leads to a loss of p63 expression (23). Despite these known functional interactions between p63 and HPV, very few studies have examined the specific role of p63 in modulating gene expression in HPV+ HNSCC.
Hierarchical clustering analyses of HPV+ HNSCCs revealed two distinct subtypes based on gene expression profiles, copy number alterations, mutational profiles, and patient outcomes (11,14). One subtype characterized by the amplification of the 3q chromosomal region, including the TP63 locus, was shown by two independent studies to be enriched in pathways involved in keratinization and cell adhesion (11,14). Interestingly, patients with this subtype tend to have worse outcomes and respond more poorly to treatment compared to patients with tumors belonging to the other HPV+ subtype (11,14). These studies suggest that there are HPV-dependent mechanisms that affect p63 function in HPV+ HNSCC and the pathology of this disease.
To explore the oncogenic role of p63 in HPV+ HNSCC, we established a p63-driven gene regulatory network based on both preclinical cell culture models and tumor datasets. Our in-depth examination of p63 in the broader transcriptomic and genomic context revealed that p63 regulates critical sets of genes and pathways in the HPV infection pathway and HPV-associated malignancy, including PI3K signaling, WNT signaling, and cell cycle control which may inform the clinical differences between the HPV+ HNSCC subtypes. Importantly, we found that p63 expression correlates with the more aggressive HPV+ HNSCC subtype and that it directs the associated gene expression programs. Finally, we identified a potentially important role for p63 in regulating PI3K signaling and mTOR signaling in HPV+ HNSCC, which may have implications for future treatment choices. Our studies suggest that p63 is an important driver of the subtype-specific gene expression program in HPV+ HNSCC, and can serve as a biomarker to identify patients with more aggressive disease.

Cell Culture Studies
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.

Knockdown of p63
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).

ChIP of p63 and Histone Marks
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 crosslinked 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 mg of p63 4A4 antibody (Santa Cruz Biotechnology) and 2 mg of DNp63-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.

ChIP-Seq Analysis
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).

RNA Isolation and Library Preparation for RNA-Seq
Total RNA from cell lines was extracted using a Direct-zol RNA miniprep kit (Zymo Research). The extracted RNA was snapfrozen 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.

RNA-Seq Analysis
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.

qRT-PCR Analysis
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.

HNSCC Dataset Analysis
RNA-seq data from patient samples were obtained from GEO (GSE122512, GSE112026, GSE74927, and GSE72536) (11,(36)(37)(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.

Determination of Enhancers and SEs According to H3K27Ac Marks
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).

Histone Modification Enrichment at p63 Binding Sites
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).

Genomic Feature Assignment
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.

Motif Enrichment Analysis of Enhancers
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.

Gene Ontology/Pathway Enrichment Analysis
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)(46)(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).

Statistics
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.

Generation of a p63 Gene Signature From the HPV+ HNSCC TCGA Tumor Dataset
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-kB 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.

Mapping the Genomic Targets of p63 in Representative HPV+ HNSCC Cell Lines
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 DNp63 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 DNp63-specific antibody, DNp63-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 HPVassociated 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.
Characterizing the Enhancer Landscape of p63 high HPV+ Cells identify gene regulatory features such as active enhancers (H3K27Ac h i g h and H3K4Me1 h i g h ), 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.
p63 Is Super-Enhancer Marked and Regulates Expression of Super-Enhancer-Associated Genes in HPV+ HNSCC 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 HPVinduced 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 p63bound 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.

Loss of p63 Expression Dysregulates Signaling Pathways Involved in HPV-Associated Carcinogenesis
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 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).

Generation of an Overall p63-Driven Gene Signature in HPV+ HNSCC
To delineate our high-stringency p63-driven gene signature, we combined the gene signatures we identified from TCGA tumor   Figure 11A and Supplementary Table 5). These analyses identified 1,052 genes shared between the TCGA and SCC104 datasets and 827 genes s h a r e d b e t w e e n t h e T C G A a n d S C C 1 5 2 d a t a s e t s (Supplementary Figure 11B). Then, to identify genes which were directly regulated by p63, we incorporated our p63 ChIPseq 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-enhancerassociated genes (Supplementary Figure 11B and Supplementary Table 5). These analyses provided a p63driven 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 subtypespecific 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 epithelialassociated 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 p63 Exerts Broad Control of the PI3K-Signaling Pathway 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.

DISCUSSION
HNSCC associated with high-risk human HPV infection is a growing problem that is clinically and biologically distinct from HPV-HNSCC. Molecular studies of HPV+ HNSCC in the past focused primarily on tumor suppressor pathways that are targeted by viral oncoproteins such as E6, which inactivates the p53 tumor suppressor protein by instigating its degradation (9). However, the operation of oncogenic drivers in the underlying complex genomic and epigenomic milieu of HPV+ HNSCC remains unclear. Studies of p63, especially oncogenic DNp63, have established its role in directing broad transcriptional programs in SCCs in various anatomical sites, including the oral cavity (21,49,76). However, the specific role of p63 in modulating the transcriptomic landscape of HPV+ HNSCC has received less attention. To address this shortcoming, we leveraged genomic, transcriptomic, and epigenomic data from HPV+ HNSCC tumors and preclinical cell line models and identified p63 as a critical regulator that affects multiple facets of HPV+ HNSCC biology, including pathways essential in HPVmediated carcinogenesis, and HNSCC subtype-specific gene expression.
One notable finding is the wide range of p63 expression levels across the HPV+ tumors, which is clearly evident in several independent datasets. This reinforces results from previously defined subtypes of HPV+ HNSCC based on hierarchical clustering of gene expression (11,13,14,51). We provide evidence that the molecular and phenotypic attributes of the more aggressive p63 hig h HPV-KRT subtype, such as keratinization and cell adhesion, are likely controlled by a p63-driven direct transcriptional program. Indeed, this is well supported by our integrated analysis of several complementary models that include HPV+ patient tumors and by an epigenetically refined DNp63 cistrome assembled from RNAseq after p63 knockdown and ChIP-seq of p63 and histone marks in SCC104 and SCC152 cell lines. We suspect that the heterogeneity of HPV+ tumors is shaped not only by the expression and the transcriptional output of p63 in the tumor epithelium but also by the enhanced influence of p63 in the tumor microenvironment, as evidenced in triple-negative breast cancer (77). This notion is supported by our identification of enriched immune pathways and immune-related genes as p63 targets, some of which were associated with SEs in SCC104 and SCC152 cells.
The tumor microenvironment of HPV+ HNSCC is distinct from that of its HPV-counterpart, with greater immune infiltration, T-cell activation, and immunoregulation (51). Furthermore, two recent landmark single-cell-based studies indicated that HPV-specific infiltrating lymphocytes may mount an immune-based response to HPV+ HNSCC (78,79). In these studies, HPV oncoprotein E2 was a target of particular interest, as E2 expression is often maintained in HPV+ HNSCC unlike in HPV+ cervical cancer. Interestingly, E2 expression is lower in the p63 high HPV-KRT subtype, most likely because of higher levels of HPV genomic integration, which disrupts HPV early gene expression (11,51). We suspect that the relatively lower immune cell infiltration and activation in HPV-KRT tumors may be attributable to their p63 high and E2 low state. Taken together, these studies suggest that more attention should be paid to HPV+ subtypes, because their inherent differences, especially in their tumor microenvironments, likely affect their clinical outcomes.
Our integrated analysis revealed many interesting insights into p63 and its link to specific aspects of HPV biology. One novel potential regulatory mechanism of p63 involves E2F family members, notably E2F1 and E2F7, which regulate cell cycle gene expression in HPV infection and HPV-associated cancers (10,80). Indeed, E2F7 was consistently enriched throughout our analyses as a direct target of p63 and as a TF whose motif was enriched around many p63-bound genomic sites and in SEs. E2F7 is an atypical E2F factor and a transcriptional repressor of genes involved in genomic stability, cell proliferation, and migration (80,81). Knockdown of p63 in HPV+ HNSCC cell lines upregulated E2F7 expression, leading us to speculate that the increased p63 expression in HPV+ HNSCC negates the repressive function of E2F7. It is possible that p63 and E2F7 physically interact and coregulate genes that have joint p63/E2F motifs in their regulatory elements, an interesting notion that could be experimentally tested pending the availability of ChIP-grade E2F7 antibodies. Further investigation into the p63-E2F network is needed given the p63-dependent enrichment of pathways related to pRB, the cell cycle, and cell cycle checkpoints that are relevant for HPV+ HNSCC. Along the same line, it is likely that p63-p53 interactions are a key component in the cell cycle circuitry given p63's ability to regulate some p53 targets in the presence or absence of p53 (15,82,83).
The results from our studies reemphasize the need to further characterize the drivers and molecular attributes of HPV+ HNSCC subtypes given the potential differences in overall patient survival between the subtypes. This is particularly important in light of recent efforts to identify treatment deescalation strategies for HPV+ HNSCC to reduce adverse events while maintaining better oncologic outcomes. Unfortunately, two large phase III clinical trials have shown inferior overall survival and progression-free survival as well as increased rates of locoregional failure (6,84,85), prompting reevaluation of ongoing deintensification trials for HNSCC. Although current HNSCC treatment options take HPV status into account, we posit that a personalized genomics approach that considers HPV+ subtypes would better inform treatment options and prevent failure of treatment de-escalation. The activated PI3K signaling that we uncovered in the p63 high HPV-KRT subtype also suggests a potential avenue for therapeutic intervention. Compensatory activation of downstream signaling pathways, including PI3K, has been suggested as one of the major mechanisms of resistance to EGFR inhibitors, including cetuximab (86). The addition of a mTOR/PI3K inhibitor effectively controls cell growth in EGFR inhibitor-resistant HNSCC, suggesting that combination therapy may increase treatment efficacy (87,88). In addition, mTOR inhibitors show promising anticancer effects in HPV+ HNSCC xenograft mouse models (75). We suspect that patients with the p63 high tumor subtype would benefit from a combination of EGFR and PI3K inhibitor treatment and radiotherapy.
Altogether, the results from our studies suggest p63 and its key downstream effectors can be used as stratification markers for HPV+ HNSCC patients. However, this requires validation of our preclinical genomic and epigenomic data in tumors from HPV+ HNSCC patients.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi. nlm.nih.gov/geo/query/acc.cgi?acc=GSE182133.

AUTHOR CONTRIBUTIONS
AG designed and performed experiments, acquired and analyzed data, prepared the figures, and wrote the manuscript. AO and CG performed experiments and acquired data. JB analyzed data. SS supervised the project and experimental design and wrote the manuscript. All authors contributed to the article and approved the submitted version.
FUNDING AG was supported partly by NYSTEM contract no. C30290GG. this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Copyright © 2022 Glathar, Oyelakin, Gluck, Bard and Sinha. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.