Noncanonical HPV carcinogenesis drives radiosensitization of head and neck tumors

Significance Human papillomavirus–associated (HPV+) head and neck squamous cell carcinoma (HNSCC) is now the most common HPV-associated cancer with increasing incidence. HPV-mediated oncogenesis is generally thought to rely on the integration of the viral DNA into the host genome, loss of HPV early gene 2 (HPV E2) expression, activation of Phosphatidylinositol-4,5-Bisphosphate 3-Kinase Catalytic Subunit Alpha (PIK3CA), and apolipoprotein B mRNA editing catalytic polypeptide (APOBEC)-mediated mutagenesis. We report the identification of a subclass of HPV+ carcinomas comprising ~45% of HPV+ HNSCC that is not associated with any of these classic features. Patients in this subgroup have robustly improved clinical outcomes, and cell models with genomic and transcriptomic features of this class have increased sensitivity to radiation. The recognition of biologically distinct tumor subclasses of HPV+ HNSCC with differential responses to radiotherapy may fundamentally alter how patients are treated.

We analyzed transcriptional data from 104 HPV+ (Human papillomavirus) HNSCC (head and neck squamous cell carcinoma) tumors together with two publicly available sources to identify highly robust transcriptional programs (modules) which could be detected consistently despite heterogeneous sequencing and quantification methodologies. Among 22 modules identified, we found a single module that naturally subclassifies HPV+ HNSCC tumors based on a bimodal pattern of gene expression, clusters all atypical features of HPV+ HNSCC biology into a single subclass, and predicts patient outcome in four independent cohorts. The subclass-defining gene set was strongly correlated with Nuclear factor kappa B (NF-κB) target expression. Tumors with high expression of this NF-κB module were rarely associated with activating PIK3CA alterations or viral integration, and also expressed higher levels of HPHPV E2 and had decreased APOBEC mutagenesis. Alternatively, they harbored inactivating alterations of key regulators of NF-κB, TNF receptor associated factor 3 (TRAF3), and cylindromatosis (CYLD), as well as retinoblastoma protein (RB1). HPV+ HNSCC cells in culture with experimental depletion of TRAF3 or CYLD displayed increased expression of the subclass-defining genes, as well as robust radio-sensitization, thus recapitulating both the tumor transcriptional state and improved treatment response observed in patient data. Across all gene sets investigated, methylation to expression correlations were the strongest for the subclass-defining, NF-κBrelated genes. Increased tumor-infiltrating CD4+ T cells and increased Estrogen receptors alpha (ERα) expression were identified in NF-κB active tumors. Based on the relatively high rates of cure in HPV+ HNSCC, deintensification of therapy to reduce treatment-related morbidity is being studied at many institutions. Tumor subclassification based on oncogenic subtypes may help guide the selection of therapeutic intensity or modality for patients with HPV+ HNSCC.

HPV | head and neck cancer | radiosensitization | tumor microenvironment | patient outcome
The incidence of Human papillomavirus-associated (HPV+) head and neck squamous cell carcinoma (HNSCC) has dramatically increased over the last few decades and continues to rise (1,2). Due to the long latency from infection to tumor formation, HPV vaccination is not expected to impact this epidemic until 2060 (3). Despite the higher cure rate associated with HPV positivity, 20 to 30% of patients with HPV+ HNSCC suffer recurrence and have limited curative options (4)(5)(6)(7).
Aggressive multimodality therapy developed for HPV-negative HNSCC is used for HPV+ patients and results in better survival compared to HPV-negative patients. Treatment-related morbidity coupled with relatively good outcomes in HPV+ HNSCC has sparked interest in therapeutic de-escalation for these patients (1,8). Presently, deintensification of therapy to reduce morbidity is being frequently studied in patients with locally advanced disease, and preliminary results indicate that in nonsurgical patients, the radiation dose or field, but not chemotherapy, can be reduced (4,(9)(10)(11). However, these deintensification studies exclude many patients because of disease extent or tobacco history and likewise may include patients with aggressive tumors despite early stage and minimal smoking history. Since standard intensive therapy fails to cure ~30% of patients (4,6), appropriate patient selection has been a persistent barrier for the expansion of de-escalation strategies. As such, there is a growing demand for molecular tools to guide personalized treatment.
HPV+ HNSCC research has focused largely on HPV oncogenes E6 and E7 and how viral integration promotes carcinogenesis, tacitly assuming that HNSCC biology would mirror the historically more common HPV-induced cervical cancer. However, several key differences between cervical cancer and HPV+ HNSCC contradict the validity of this assumption, including near-complete absence of HPV type 18 in HNSCC (12,13), and OPEN ACCESS distinct gene and protein expression (14,15), as well as mutation (16)(17)(18)(19), profiles. Moreover, HNSCCs more commonly maintain expression of all HPV16 early genes and lack viral integrations (20)(21)(22). These differences suggest a broader understanding of oncogenesis in this relatively understudied disease is warranted (19).
Therefore, we hypothesized that tumors lacking classical features of HPV-driven carcinogenesis (loss of E2, viral integration, PIK3CA alteration, and APOBEC mutagenesis) might represent a group of tumors with alternate requirements for tumor development. To address this question in an unbiased fashion, we turned to unguided analysis of transcriptomic data. Transcriptome-wide profiling has transformed the classification of some cancer types impacting prognostication and early detection (23). Gene expression signatures have been widely used to identify important characteristics of genotypes or phenotypes of interest (24)(25)(26)(27)(28)(29)(30)(31). The development of statistical tools to address gene set abundance in the context of the larger transcriptome has solidified the utility of gene expression signatures in the study of biology (32)(33)(34).
Given the high dimensionality of transcriptomic data, it is not surprising that gene sets often can be identified with strong correlation to phenotypes using single datasets (35), but many of these signatures become less cohesive and lose the power to predict phenotype when applied to external datasets (36). Increasing amounts of publicly available transcriptomic data will hopefully accelerate the identification of robust gene expression signatures that are transferrable across labs, experimental techniques, and cohorts of patients for a disease of interest (36). Therefore, we analyzed our transcriptional data from 104 HPV+ HNSCCs in tandem with two publicly available sources to identify gene sets that were detected consistently, despite heterogeneous sequencing and quantification methodologies.
Surprisingly, we identified a single transcriptional program that naturally subclassifies HPV+ HNSCC tumors based on a bimodal pattern of gene expression, clusters all atypical features of HPV+ HNSCC biology into a single subclass, and predicts patient outcome. The gene expression signature and related radiosensitivity were recapitulated in HPV-positive head and neck cancer cells in culture, validating our translational findings.

Results
Defining Robust Transcriptional Programs in HPV+ HNSCC. To identify robust gene expression modules indifferent to sequencing strategy and technical variability, we analyzed three independent datasets (16,37). After normalization and filtering of very low expression genes, 11,843 genes were available in all three datasets. Internally correlated gene sets were identified from the three cohorts using weighted gene correlation network analysis (WGCNA). Dissimilarity matrices are displayed in Fig. 1A. In each dataset, genes were binned into groups (modules) with similar expression patterns using the WGCNA dynamic tree cutting algorithm (38). Consensus modules were then created by selecting genes that grouped together in WGCNA analyses from all three independent cohorts. To eliminate spurious assignment of consensus in the case of large modules, consensus modules were also required to comprise at least 5% of the largest single dataset module. Consensus module size ranged from 5 to 464 genes (Fig. 1B). As expected, intermodule correlation testing revealed relatively few consistent relationships between consensus modules (Fig. 1C), whereas many were apparent when examining single datasets (SI Appendix, Fig. S1 WGCNA modules are defined independent of polarity, that is anticorrelated genes may be grouped together; however, the polarities are important for biological interpretation. The gene numbers and polarities of comprising genes for each module are summarized in Fig. 2A. Since module compactness (intramodule correlation) is a key quality metric for gene expression signatures (36), the compactness of modules was compared to that of randomly selected genes showing that consensus module genes were correlated (intramodule correlation) far above background (Fig. 2B). Interestingly, modules with higher levels of negatively correlated genes (chiffon, dark goldenrod, chocolate, cornflower blue, and goldenrod) had weaker intramodule correlation despite the use of the absolute value of Spearman's Rho. In a few cases (chocolate, cornflower blue, and goldenrod), intramodule correlations were no better than randomly selected genes, suggesting that these modules were of lower quality. The compactness of gene expression modules can also be assessed with principle component analysis (36), specifically examining how much variance of the module genes can be explained by one principal component. Although polarity independent, the bidirectional modules were also less well described by one principal component (SI Appendix, Table S1).
We also assessed modules for evidence of multimodality distribution of gene expression across tumor samples using excess mass-based testing (39). Interestingly, only the sky blue module displayed multimodality in all cohorts with tumors grouping to high and low expression ( Fig. 2 C-E and SI Appendix, Table S1). Quantification metric independence is also a key feature of high-quality gene signatures (36). Quantification of the sky blue module was also highly metric independent, as seen by the correlation between the PC1 score and median of expression ( Fig. 2 C-E). Similar metrics for each consensus module are presented in SI Appendix, Table S1.
Gene Set Enrichment Testing. To begin determining the biological context of identified modules in HPV+ HNSCC, hypergeometric enrichment analysis was performed for MSigDB Hallmark gene sets, (Fig. 3A). Estrogen receptor alpha expression, defects in NF-κB regulators, and immune infiltration have been correlated with outcomes in HPV+ HNSCC. Interestingly, the sky blue module associated with NF-κB signaling and early estrogen response genes, and because of this association, we dubbed the sky blue module the NF-κB related module. The wheat tan module associated with various signatures related to inflammation, and therefore may correlate with tumor immune infiltration, and the thistle grey module associated with early and late estrogen response genes. A total of 30 Hallmark gene sets had a significant association with the consensus WGCNA modules; of these, 19 were associated with only one module (Fig. 3A). Only the Hallmark Allograft Rejection signature (1/30) was associated with more than 2 WGCNA modules, suggesting that these consensus modules are relatively specific for associated biological functions.
We identified three groups of consensus modules containing strong intermodule correlations across all cohorts ( Fig. 1C No. 1 to 3) and interrogated whether they might also have biological relevance using a similar enrichment analysis. While these module groups (meta-modules) retained associations identified by analysis of the individual modules, no additional biological associations arose when modules were analyzed as a group.
Tumor Microenvironment (TME). Since the modules were developed from bulk tumor sequencing, both cancer cells and components of the TME were sequenced. Using EcoTyper, a non-negative matrix factorization-based approach, we estimated the proportion of various cell types in the TME (40). The average correlation between the consensus module expression and TME components across the three cohorts is displayed in Fig. 3B. Because of the bimodal distribution of tumors in the NF-κB module, we were particularly interested in the correlation between the sky blue module and CD4+ T cells. To begin determining whether CD4+ T cell abundance in tumors correlated with identified subclasses, we queried CD4 signatures in high and low expressers (Fig. 2 C-E). Across all three cohorts, CD4+ T cells in tumors were elevated in the NF-κB (sky blue) high samples compared to NF-κB low samples (Fig. 3C). A graphical summary of the biological profiling associated with gene expression modules, including both gene set enrichment and TME analyses, is presented in Fig. 3D.
Using a subset of tumors with available FFPE from the UNC cohort (n = 50), we validated the RNA-based findings related to CD4+ cell infiltration (Fig. 3 E-J). Interestingly, the CD4+ cells correlated with the NF-κB module exclusively in the epithelial component defined by cytokeratin staining (Fig. 3G), but not in mesenchymal areas of the tumors. The nine tumors with the most tumor-invading CD4+ cells were all in the high sky blue (NF-κB) expressing group, Fig. 3H. The relative proportion of CD4+ cells in the epithelial component versus the mesenchymal component was also increased in the NF-κB high group, consistent with the correlation being specific to the regions of epithelial (cancer) cells (Fig. 3I). The ratio of CD4+ to CD8+ cells in the tumor mesenchymal tumor component was also distinct, with the high NF-κB expressers having a higher mesenchymal CD4+/CD8+ ratio (Fig. 3J). Based on these initial findings, we also performed multiplex immunofluorescence microscopy on the same samples staining for FoxP3 and CD4. We found that NF-κB high tumors had a higher proportion of FoxP3+/ CD4+ T cells, suggesting that they may be enriched in regulatory T cells ( Fig. 3K and SI Appendix, Fig. S2). Considering the correlation between the sky blue (NF-κB) module and early estrogen response genes and the reported prognostic value of ERα expression (41), we examined estrogen receptor alpha (ERα) expression and found that the ERα level in tumor cells was highly associated with expression of the NF-κB module, (Fig. 3L).
Somatic Variant Analyses. Next, we examined which genes were most frequently altered in the tumors with high versus low expression in the NF-κB module. Publicly available TCGA mutation and copy number calls were downloaded and used without  modification (42). Mutations in large genes (Titan, MUC16), known to be primarily associated with tumor mutational burden, were excluded from analysis. Interestingly, alterations in a known head and neck cancer driver, PIK3CA, were strongly enriched in the low NF-κB group. Conversely, NF-κB regulatory genes (CYLD and TRAF3) were nearly exclusively altered in the high NF-κB group (Fig. 4A). Notably, RB1 mutations and deletions were also enhanced in the NF-κB high group. An association between RB1 and CYLD alterations has been previously reported (43).
Considering that single copy (GISTIC score = −1) losses in TRAF3 and CYLD were quite common in the TCGA cohort, we asked whether low level but simultaneous loss of these NF-κB regulators associated with expression of the sky blue module. Indeed, tumors with single copy loss of TRAF3 and CYLD had higher expression of the sky blue module genes, (Fig. 4B). To validate this finding, we performed copy number profiling based on the UNC cohort using the CNVkit RNA pipeline, which identifies large-scale copy number aberrations based on regionally correlated patterns of RNA expression. In the UNC cohort, tumors with simultaneous chromosomal arm level losses (log 2 Ratio < −0. 35

) in both TRAF3
and CYLD also demonstrated increased expression of the NF-κB module (Fig. 4B).
To determine whether other NF-κB pathway components were altered in the high NF-κB group, we performed pathway-oriented mutational analysis, examining a set of 9 NF-κB genes reported to be altered in other NF-κBdriven cancers (Fig. 4C). Interestingly, this analysis revealed that many tumors in the NF-κB active group harbored multiple pathway alterations and that ~80% of NF-κB active tumors had some NF-κB pathway change. Excluding TRAF3 and CYLD, other NF-κB pathway alterations were also detectably enhanced in tumors with high expression of the sky blue module, Fig. 4C. NF-κB pathway mutations were predictive of recurrence-free survival in TCGA data, Fig. 4D.
To determine whether differences in mutagenic processes underlie or parallel differences in gene defects, we analyzed single nucleotide polymorphisms (SNPs) and their trinucleotide context, Fig. 4 E-H. NF-κB low tumors had significantly more SNPs, (Fig. 4E). C > G mutations were the most common in both groups, but NF-κB low tumors demonstrated more trinucleotide context specificity (SI Appendix, Fig. S3). Related to this, the fraction of APOBEC-related variants (as      High Low  quantified by Maftools) was marginally higher in the NF-κB low group (Fig. 4F). Evaluation of COSMIC signatures in the subgroups revealed that APOBEC mutagenesis was associated with more than 50% of mutations in the NF-κB low group, compared to only 24% in the NF-κB high group (Fig. 4 G and H). Taken together, these results reveal distinct patterns of somatic alteration and mutagenic processes in high vs. low gene expression subgroups defined by the sky blue (NF-κB) module.

Viral Gene Expression and Integration Analyses.
We previously reported that TRAF3/CYLD alterations in HPV+ HNSCC are associated with increased NF-κB activity and lack of viral integration (44). Integration of the HPV genome correlates with viral gene expression patterns (45), and increased expression of E6 and E7 is a hallmark of viral genomic integration in HPV+ HNSCC (45). To determine whether subsets of HPV+ HNSCC identified by high or low expression of the NF-κB module had differences in viral integration or viral gene expression, we investigated the TCGA and UNC cohorts. Viral integration was determined by read pairs that displayed discordant mapping to the human and HPV16 genomes (SI Materials and Methods). In our analyses, the ratio of E6 and E7 expression to E2 and E5 was most related to integration, and thresholds for high E6/E7 expression relative to E2/E5 were selected to best agree with integration calls.
Unguided clustering based on viral gene expression demonstrated gross relationships to sky blue module, as well as expected relationships to viral integration status ( Fig. 5A and SI Appendix, Fig. S4A). Tumors with high E6E7/E2E5 gene expression ratios were associated with lower expression of sky blue module genes ( Fig. 5B and SI Appendix,  Fig. S4B). Tumors with discordant read pairs supporting viral integration were associated with lower expression of sky blue (NF-κB) module genes, (Fig. 5C and SI Appendix, Fig. S4C) Taken together, these results demonstrate that tumor subclasses defined by the sky blue (NF-κB) module are associated with a lack of viral genomic integration and a distinct pattern of viral gene expression.

Genomic Methylation Analysis.
Considering that HPV+ HNSCC tumors have hypermethylated genomes compared to HPV-HNSCC tumors, we compared global methylation differences across the tumor sample groups delineated by the sky blue (NF-κB) module and correlated expression of module genes located close to CpG islands. Using TCGA data, we performed consensus clustering across the samples using the 2,000 methylation probes with the highest variance among HPV+ HNSCC. This approach resulted in two groups of HPV+ HNSCC tumors, with relatively high and low levels of methylation across the included probes, see Fig. 6A. NF-κB (sky blue) module expression was higher in tumors in the more methylated group (Fig. 6B) with a corresponding enrichment for NF-κB high tumors in the high methylation group (Fig. 6C). Although the average beta-value across all probes was not significantly different between the NF-κB active and inactive groups, there was a trend toward hypermethylation in the NF-κB active tumors (P = 0.097, Wilcoxon test). This trend is consistent with the clustering on high all variance probes, as well as mutational signature analysis which revealed increased 5-methylcytidine deamination in NF-κB active tumors (Fig. 4F). Consensus clustering on the CpG probes identified a smaller (~130) group of probes (Figs. 6A and 7D, dark purple) with a tighter association to NF-κB (sky blue) module expression upon repeat clustering (Fig. 6 D-F). Interestingly, the probes identified in this unguided manner were largely anticorrelated with the global differences in methylation levels, with these probes showing lower methylation in the NF-κB active group (Fig. 6 D-F).
Given the relationship between the NF-κB module and gross patterns of genomic methylation, we analyzed the other consensus modules for methylation sites that could influence module gene expression. Methylation probe-gene pairs for each module were identified if probes were proximal in the genomic space and had a high correlation between expression and methylation. We then used robust rank aggregation to control for the number of probes associated with each gene and to produce a gene-level P-value (SI Appendix). We found that ~80% of the genes in the sky blue (NF-κB) module had significant expression-methylation association (Fig. 6G). To assess the relative degree to which the expression-methylation correlations in identified modules were stronger than expected at random, we calculated Gene set enrichement analysis scores for each module compared to the distribution of gene-level methylation-expression correlations across all genes. We found that the sky blue (NF-κB) module was by far the most enriched compared to all other modules (Fig. 6H). In summary, the methylation analysis demonstrated that global epigenetic programming is different between the NF-κB high and low tumors and highlighted a strong relationship between methylation of sky blue (NF-κB) module GpC probes and RNA expression of this gene module.   Fig. S4 for similar analysis of the TCGA cohort. *P value < 5*10^−2. **P value < 5*10^−3.

Survival Analysis.
Considering the general interest in developing predictive biomarkers for guiding therapy in HPV+ HNSCC, we looked for correlations with outcome. For modules that did not intrinsically separate tumors into 2 classes (all but sky blue), the median module principal component analysis score was used as the cut point to divide tumors into 2 groups (high and low). For the bimodal sky blue module, tumors were stratified by an empiric threshold based on the nadir between the two expression peaks (Fig. 2 C-E). Interestingly, the sky blue (NF-κB) module was the only gene set that segregated patients into groups significantly associated with differences in recurrence-free survival (RFS) in all three cohorts (Fig. 7). Finding that the NF-κBrelated module intrinsically segregated patients into two groups with better or worse prognosis confirms the significance of the biological differences between these tumor groups.
To provide external validation of these clinical findings, we performed RNA sequencing of 57 patients with available Formalin-Fixed Paraffin-Embedded (FFPE) tumor tissue from the ECOG-ACRIN E1308 clinical trial. Tumor groups were defined based on the expression of sky blue (NF-κB) module genes based on an empirical threshold. We found that in these patients treated with induction chemotherapy followed by chemoradiation, both progression-free survival and overall survival were improved in the group expressing higher levels of sky blue (NF-κB) module genes(- Fig. 7 D and E and SI Appendix, Fig. S5). To understand the relationship between disease stage, tobacco smoke exposure, and NF-κB status, we pooled all patients treated with chemoradiation, for whom quantitative smoking data were available (TCGA, UNC, E1308). Using this combined cohort of 166 patients, multivariate analysis was performed and demonstrated that NF-κB active tumors were associated with less tobacco smoke exposure and less tumor stage 4 disease (SI Appendix, Tables S2-S4). Multivariate event-free survival analysis of this pooled cohort demonstrated that NF-κB status was not only prognostic while controlling for stage and tobacco smoke exposure (P = 0.005) but was also the most prognostic factor in all models we examined (SI Appendix, Table S5). Interestingly, although tobacco smoke exposure is an accepted adverse prognostic factor in HPV+ HNSCC, our data suggest that this association is driven exclusively by NF-κB inactive tumors (Fig. 7F). To determine whether NF-κB status was relevant to outcomes in patients with low or no tobacco smoke exposure, who might be considered candidates for deintensified therapy, we compared all patients treated with chemoradiation with no or less than 10 pack-years of smoking history to all NF-κB active patients treated with chemoradiation. We found that the NF-κB active group had significantly improved (HR = 0.4, P = 0.02) event-free survival (RFS or progression-free survival, based on availability) (Fig. 8F).
Oncogenic NF-κB Activation Promotes Radiation Sensitivity in HPV+ HNSCC. Locoregionally advanced HPV+ HNSCCs are usually treated with radiation, either in the primary or postoperative setting. Our data show that patients whose tumors display high expression of the sky blue module have significant survival benefits that are most likely explained by increased sensitivity of tumors with constitutively active NF-κB to radiation. Therefore, we hypothesized that NF-κB activity may contribute to intrinsic tumor cell characteristics that increase radiation sensitivity in HPV+ HNSCC. Indeed, we found that experimental TRAF3 and CYLD inactivation was sufficient to strikingly sensitize HPV+ head and neck cancer cells to radiation, as well as promote increased expression of well-known NF-κB targets and sky blue (NF-κB) module genes ( Fig. 8 A-E). Thus, our data support the premise that NF-κB signaling, which drives an alternative mechanism of HPV carcinogenesis, sensitizes cells to radiation through cancer cell-intrinsic effects.

Discussion
Using three independent HPV+ HNSCC cohorts, an unbiased approach identified robust gene expression profiles existing in these tumors (Figs. 1 and 2) and identified the major molecular pathways and cell types associated with each of these expression profiles (Fig. 3). Our group previously identified an etiologically distinct subset of HPV+ HNSCC characterized by the inactivation of TRAF3 or CYLD leading to constitutively active NF-κB (44). This prior work correlated the inactivation of TRAF3 and CYLD with improved patient survival and a lack of HPV integration, suggesting that NF-κB activation may enable cells to maintain HPV episomes as a distinct mechanism of carcinogenesis (44). In the current study, 22 expression modules were identified, but only one had ideal statistical characteristics for tumor subclassification, namely a bimodal distribution in all three independent cohorts (Fig. 2); this sky blue module was strongly associated with activated NF-κB signaling and contained all tumors with TRAF3 and CYLD inactivating defects (Figs. 3 and 4). Herein, we also report that tumors with simultaneous shallow (single copy) losses at both, CYLD and TRAF3, loci were robustly associated with high expression of the sky blue (NF-κB) module (Fig. 4), suggesting that mechanisms in addition to deep deletion of TRAF3 or CYLD modulate NF-κB activation in HPV+ HNSCC. Applying pathway-based analysis, we find that while individually rare, when viewed as a group, mutations in other NF-κB pathway genes were also more common in tumors with NF-κB activation (Fig. 4). We also found that RB1 deletions and mutations were significantly enriched in the NF-κB active group (Fig. 4). Interestingly, constitutively active NF-κB has been shown to reduce p53 (46), as well as Rb activity, while phosphorylated Rb can inhibit NF-κB (47), suggesting that the loss of Rb function may serve an important oncogenic role in tumors dependent on highly active NF-κB. An association between RB1 mutations and CYLD mutations in HPV+ HNSCC has recently been reported by other groups (43).
APOBEC3 cytidine deamination hypermutates viral genomes, including HPV (48)(49)(50)(51), and deep sequencing of human tumors revealed that APOBECs are major sources of mutations in several cancers, including HNSCC, where they create driver mutations (e.g., PIK3CA) (49). We recently reported that HPV+ tumors containing APOBEC signatures had markedly higher mutational burden compared to tumors lacking APOBEC mutations (49). The data presented here further reveal that APOBEC mutation signatures, as well as APOBEC-dependent PIK3CA mutations, are significantly more frequent in NF-κB inactive tumors (Fig. 4), which is consistent with the observed relative importance of APOBEC mutagenesis across the entire genome in the NF-κB inactive group (Fig. 4). We also find that missense mutations and focal amplifications of PIK3CA are present predominantly in the NF-κB inactive group and are relatively uncommon in the NF-κB active group (Fig. 4). Consistent with our results, PIK3CA mutations have been associated with poor prognosis in HPV+ HNSCC (52). Considering that activation of PIK3CA may limit the utility of related targeted therapies (cetuximab), it is intriguing to consider that the NF-κB module may represent a biomarker for a subset of patients that are more responsive to this therapy.
In line with our prior findings, here we show that NF-κB active tumors were associated with different patterns of viral gene expression and were less associated with HPV viral genomic integration (Fig. 5). Mechanisms explaining how a nonintegrated form of HPV may induce head and neck carcinogenesis are not well described. Considering that NF-κB is activated in various tumor types, and its activity promotes carcinogenesis through cellular proliferation, transformation, and angiogenesis, while protecting tumor cells from apoptosis (55), we hypothesize that activation of NF-κB due to genetic alterations in NF-κB regulators (e.g., TRAF3 or CYLD, Fig. 4) in the context of permissive epigenetic landscape promoting expression of NF-κB target genes (Fig. 6H) represents a distinct pathway for HPV16 mediated oncogenesis in the oropharynx (SI Appendix, Fig. S6).
We further explored differences between subtypes based on NF-κB signaling and found that patients with tumors harboring constitutively active NF-κB had improved survival in all three independent cohorts (Fig. 7). This finding was validated by performing RNAseq of samples from the ECOG-ACRIN1308 clinical trial (Fig. 7). Furthermore, in multivariate analysis, NF-κB status was more prognostic than established clinical factors and also remained a significant predictor of outcome while controlling for these factors (SI Appendix, Table S5). HNSCCs are often treated with radiation; rates of response and cure are higher for patients with HPV+ compared to HPV-negative HNSCC; however, our data suggest that the majority of survival benefit is attributable to the subtype of HNSCC with highly active NF-κB (Fig. 4). Our data suggest that constitutively active NF-κB, determined by the presence of inactivating mutations in NF-κB regulators or by a specific RNA signature (56), may serve as prognostic biomarkers to help clinicians with therapeutic decisions. Other groups have suggested ERα expression as a prognostic marker for this purpose, and we herein show that groups defined by high ERα expression and NF-κB activity are largely identical (Fig. 3L) (41). Since NF-κB activity is most commonly associated with resistance to therapy, we searched for drivers of therapeutic sensitivity and found an increased number of tumor-infiltrating CD4+ cells (Fig. 3). Better prognosis and elevated sensitivity to radiotherapy of HPV+ HNSCC, as compared to HPV-negative tumors, has been recently correlated with increased infiltrating immune cells (CD8+ T and B cells) (57), and NF-κB activity has been shown to increase T cell infiltration (58,59). CD4+ T cells have been implicated as both enhancers and inhibitors of antitumor immune response based on CD4+ T regulator or T helper status (60). A higher proportion of T regulatory cells in NF-κB active tumors ( Fig. 3K and SI Appendix, Fig. S2) indicates a potential for immunotherapy. Importantly, activation of NF-κB via depletion of TRAF3 or CYLD resulted in elevated radiosensitivity of HPV+ head and neck cancer cells (Fig. 8). Thus, increased NF-κB activity in HPV+ HNSCC contributes to intrinsic tumor cell (Fig. 8) and TME characteristics (Fig. 3) that promote increased radiation sensitivity.
In summary, the data presented here reveal the existence of two intrinsically different subtypes of HPV-positive head and neck cancer. The NF-κB active subset of HPV+ HNSCC displays not only a distinct pattern of NF-κB related gene expression but also markedly different mutational and methylation profiles. The striking lack of PIK3CA variants, lack of HPV integration, and decreased APOBEC mutagenesis in the subtype of tumors with active NF-κB suggest that they are driven by an alternative mechanism of HPV-related carcinogenesis (SI Appendix, Fig. S6).
Our data also have important clinical implications. Clinicians treating HPV+ HNSCC are actively searching for actionable biomarkers to guide therapeutic intensity (61). These data, as well as our prior work (44,56), demonstrate that two prognostic subclasses with different levels of NF-κB activity can be identified by many approaches related to their striking differences in many aspects of the cancer genomes. These groups are robustly prognostic and integrate prior observations related to the prognostic value of viral integration (62), PIK3CA mutations (52), and ERα expression (41). We believe that our data suggest that these two tumor groups harbor fundamentally different requirements for tumorigenesis and tumor maintenance. As such, we hope that the strong biological differences will also present opportunities to explore targeted therapies, which might be effective when applied with precision to these highly distinct tumor classes.

Materials and Methods (For Detailed Materials and Methods, Please See SI Appendix)
Three different HPV+ HNSCC cohorts were used in our study to determine and characterize transcriptional modules and perform survival analysis: TCGA (has both, DNAseq and RNAseq), UNC (RNAseq and immunofluorescent staining), and Vanderbilt (RNAseq). To validate clinical benefits of a sky blue module, we performed RNAseq of an independent, prospectively collected and uniformly treated cohort from the E1308 trial. To investigate whether overactive NF-κB contributes to intrinsic tumor cell characteristics that increase radiation sensitivity in HPV+ HNSCC, we used HPV+ HNSCC cells UMSCC47 and CRISPR/Cas9 system or siRNAs to delete TRAF3/CYLD. Data Acquisition. Deidentified publicly available clinical and genomic data from the Cancer Genome Atlas (TCGA) were utilized for this study. In this work, we consider a Gistic score of −2 synonymous with deep deletion. Gistic uses a dynamic segmentation algorithm to define chromosomal arm level (−1) and deeper focal deletions (−2) based on per tumor thresholds (63) Clinical data for the TCGA HNSCC cohort were acquired from Liu et al. (64) Variant calls were downloaded using the R TCGAbiolinks package; (65) calls performed with VarScan (66) were used for all analyses. TCGA RNA sequencing BAMs were downloaded from dbGaP, with NIH request #99293-1 for project #27853. Data, Materials, and Software Availability. UNC and E1308 RNAseq data were deposited to dbGaP (dbGaP Study Accession: phs002935.v1.p1, "Transcriptomic Profiling of Oropharyngeal Squamous Cell Carcinoma" (67); and dbGaP Study Accession: phs003320.v1.p1, "RNA sequencing of ECOG-E1308" (68), respectively). All other data are included in the article and/or supporting information.