NF‐κB and stat3 transcription factor signatures differentiate HPV‐positive and HPV‐negative head and neck squamous cell carcinoma

Using high‐throughput analyses and the TRANSFAC database, we characterized TF signatures of head and neck squamous cell carcinoma (HNSCC) subgroups by inferential analysis of target gene expression, correcting for the effects of DNA methylation and copy number. Using this discovery pipeline, we determined that human papillomavirus‐related (HPV+) and HPV− HNSCC differed significantly based on the activity levels of key TFs including AP1, STATs, NF‐κB and p53. Immunohistochemical analysis confirmed that HPV− HNSCC is characterized by co‐activated STAT3 and NF‐κB pathways and functional studies demonstrate that this phenotype can be effectively targeted with combined anti‐NF‐κB and anti‐STAT therapies. These discoveries correlate strongly with previous findings connecting STATs, NF‐κB and AP1 in HNSCC. We identified five top‐scoring pair biomarkers from STATs, NF‐κB and AP1 pathways that distinguish HPV+ from HPV− HNSCC based on TF activity and validated these biomarkers on TCGA and on independent validation cohorts. We conclude that a novel approach to TF pathway analysis can provide insight into therapeutic targeting of patient subgroup for heterogeneous disease such as HNSCC.

Using high-throughput analyses and the TRANSFAC database, we characterized TF signatures of head and neck squamous cell carcinoma (HNSCC) subgroups by inferential analysis of target gene expression, correcting for the effects of DNA methylation and copy number. Using this discovery pipeline, we determined that human papillomavirus-related (HPV1) and HPV2 HNSCC differed significantly based on the activity levels of key TFs including AP1, STATs, NF-jB and p53. Immunohistochemical analysis confirmed that HPV2 HNSCC is characterized by co-activated STAT3 and NF-jB pathways and functional studies demonstrate that this phenotype can be effectively targeted with combined anti-NF-jB and anti-STAT therapies. These discoveries correlate strongly with previous findings connecting STATs, NF-jB and AP1 in HNSCC. We identified five top-scoring pair biomarkers from STATs, NF-jB and AP1 pathways that distinguish HPV1 from HPV2 HNSCC based on TF activity and validated these biomarkers on TCGA and on independent validation cohorts. We conclude that a novel approach to TF pathway analysis can provide insight into therapeutic targeting of patient subgroup for heterogeneous disease such as HNSCC.
Head and neck squamous cell carcinoma (HNSCC) is the fifth most common cancer worldwide. 1 HNSCC has been tra-ditionally associated with tobacco and alcohol exposure. A subset of high-risk human papilloma virus (HPV)-related HNSCC often found in tobacco non-exposed individuals has been described as a distinct clinicopathologic entity. 2 HPV2 HNSCC is typically associated with poorer outcomes, while HPV1 HNSCC has a more favorable prognosis. 3 For patients with locally advanced HNSCC, multimodality treatments including surgery, radiotherapy and/or cytotoxic chemotherapy have significantly improved, but five-year overall HNSCC survival remains poor. 3 Both genetic and epigenetic aberrations have been shown to play a role in HNSCC development, treatment response and survival. 4 However, the only targeted agent approved by the FDA for the treatment of HNSCC is cetuximab, a monoclonal antibody directed against epidermal growth factor receptor (EGFR). 5 The development of targeted therapy is challenged by diverse genetic and epigenetic alterations heterogeneously distributed in HNSCC. The data suggest that genetic abnormalities affect a limited number of mitogenic pathways, such as PI3K and NOTCH. [6][7][8] However, genome-wide expression array data demonstrate that HNSCC is characterized by expression dysregulation for a majority of genes. 9,10 Transcription factors are a primary determinant in the regulation of gene expression. Unfortunately, direct comprehensive analysis of the dysregulation of key TFs is challenging with existing screening techniques for several reasons: it is difficult to analyze changes in TFs expression due to their transient and low level expression; activation of most TFs requires posttranslational modifications, protein cleavage, formation of active transcriptional complex or protein translocation from cytoplasm to nucleus. 11 Therapeutic targeting of TF is complicated due to the low level and cycle-dependent expression of TF, diversity of posttranscriptional modifications and redundancy of the downstream gene regulation. Despite that, many laboratories have demonstrated efficiency of several anti-TF therapies, including agents such as bortezomib [targeting NF-jB through proteasome inhibition], 12 BEZ-235 [a dual PI3K/mTOR inhibitor], 6 and decoy oligonucleotides [STAT3 and NF-jB inhibitors]. 13,14 Using high throughput array datasets and the TRANSFAC database, we have developed an inferential method describing TF activity by the gene expression of their targets. 15 Using this method we demonstrate that HPV1 and HPV2 HNSCC are differentiated by the activity of key TFs including STATs and NF-jB. Approximately 25% of primary HNSCC, almost exclusively HPV2, have these TFs and their targets coordinately upregulated. After validation of this TF alteration pattern in separate clinical cohorts, we further confirmed that HNSCC with key TF pathway alteration are sensitive to therapy directed at these key TF pathways, providing evidence that TF pathway analysis can provide insight into therapeutic targeting.

Material and Methods
The detailed methods can be found in "Supporting Information".

Tissue samples
We used three independent cohorts of HNSCC and normal specimens, overall: 195 HNSCC and 63 noncancer affected patients (Supporting Information Tables S1-S3). Every participant signed a written informed consent before participating in this study. This study was approved by Johns Hopkins Medicine Internal Review Board (JHM IRB) and performed under approved research protocol NA_00036235. TMA tis-sues were collected under ECOG/RTOG-approved protocols. We have also used publicly available data (http://cancergenome.nih.gov/cancersselected/headandneck) for the TCGA-HNSCC cohort that includes 279 HNSCC and 50 control tissues.
DNA and RNA preparation from microdissected tissues DNA was isolated by phenol-chloroform extraction after proteinase K digestion 16 and RNA was isolated with the mir-Vana miRNA Isolation Kit (Ambion) per manufacturer's recommendations.

HPV analyses
HPV status of oropharyngeal SCC tumors was tested by in situ hybridization for high-risk HPV, 17 by quantitative PCR. 16 P16 expression was checked by immunohistochemical staining. 17 Reverse transcription (RT) and quantitative real time PCR (qRT-PCR) One mg of RNA from the validation cohort was reverse transcribed using the High Capacity cDNA RT Kit (Applied Biosystems). QRT-PCR was performed on Taqman 7900HT (Applied Biosystems, Foster City, CA) with gene-specific expression assays (Applied Biosystems) using standard qRT-PCR conditions.

Cell culture
Cell lines and cell culture conditions. Human HNSCC cell lines Ho1N1, HSC2 and SKN3 were purchased from the What's new? Changes in transcription factor (TF) expression contribute to genetic and epigenetic abnormalities in cancer. In order to capitalize on those changes and advance diagnostic and therapeutic strategies for cancer, researchers must first find a way to detect and quantify changes in TF activity. Here, TF activity was estimated globally in primary head and neck squamous cell carcinoma (HNSCC) using a novel inferential approach that accounted for gene silencing and loss of heterozygosity and homozygosity. Top-scoring pair biomarkers were identified and linked to human papillomavirus (HPV) status, enabling HPV1 and HPV-HNSCC to be distinguished based on TF activity. Japanese Collection of Research Bioresource. Each cell line was authenticated using a short tandem repeat analysis kit, GenePrint 10 (Promega), as directed at the Johns Hopkins University Core Facility. Cell grew at 37 C in 5% CO 2 .
Transient transfection and cell proliferation assay. The expression of STAT1, STAT3 and RELA genes were knocked down by ON-TARGETplus siRNA SMARTpool (Thermo Scientific) using RNAiMAX Transfection reagent (Life Technologies). Cell proliferation was measured by CCK-8 kit (Dojindo).

Statistical analyses
Preparation of TF target gene sets. We applied robust multiarray average (RMA) analysis to an expression array dataset for the discovery cohort. We annotated each TF with a list of its experimentally validated targets described in the Transcription Factor Database [TRANSFAC Professional]. 18 TFs with <5 targets were removed. From each TF gene set, we removed genes that were expected to have significantly reduced expression either due to increased methylation (b > 0.15) or DNA copy loss (CNV < 1.2) on a per tumor sample basis, creating tumor-specific TF gene sets. 15 Methylation arrays were preprocessed using R scripts that provided the estimated b value of percentage methylation. SNP-chip arrays were preprocessed using corrected robust linear models with maximum likelihood-based distances to provide copy number estimates. 19 Genes were removed from the TF gene set on a patient-specific basis using the resulting methylation and copy number estimates.
TF activity ranking. The resulting 1,325 tumor-specific TF target gene sets, independently corrected in the discovery and TCGA cohort for DNA-methylation and CNV-dependent expression loss, were used to compare samples from HPV1 and HPV2 groups to the pooled values of noncancer samples for each individual target gene by pseudo-t test, where each individual tumor sample was given a score based on the probability that it came from a normal distribution with mean and standard deviation equal to the expression of that gene in the normal samples. Wilcoxon rank-sum tests were applied to these pseudo-t statistics to compare the expression of targets of each TF relative to the expression of targets not regulated by the TF in each sample relative to normal. Wilcoxon p-values of all TF gene sets were used to rank TFs.
Top scoring pair (TSP). TSP was applied to the combined list of 72 target genes of STAT1, STAT3, NF-jB and AP1 pathways (Supporting Information Table S4). TSP aimed to discover five independent pairs of genes, whose changes in relative level best separate HPV1 and HPV2 samples from the discovery cohort. These relative changes were used to score patients in TCGA and an independent validation cohort, where measurements were made by sequencing and RT-PCR methods.
Correlation of IHC staining with clinical data. Protein staining scored by Aperio software was correlated with clinical data of ECOG-TMA cohort. Wilcoxon rank sum tests were used to compare two groups and Kruskal-Wallis tests were used to compare more than 2 groups. Markers' costaining was assessed by linear regression.

Cohort assembly and array analysis
To determine gene expression, methylation and copy number in primary HNSCC, we employed Affymetrix HuEx1.0 Gene-Chips Array, Illumina Infinium HumanMethylation27 Bead-Chips and Affymetrix Genome-wide SNP 6.0 Array for 44 HNSCC (multiple primary sites, broad stage distribution) and 25 normal samples from a discovery cohort described earlier (see Supporting Information Table S1 and Refs. 8 and 9).

Array data annotation and preparation of TF target gene sets
The information regarding the target genes for each TF was obtained from the TRANSFAC 2010.4 Professional database. 18 We retained TFs with a minimum of five experimentally validated targets, creating sets of target genes for 1,325 human TFs of the total 2,600 human TF described in TRANSFAC. The decreased expression signals from methylation and copy number-lost target genes would lower the calculated activity of TF, without necessarily lowering the net activity of the TF in a particular tumor sample. To correct for estimated changes in TF activity due to the effects of DNA methylation and copy number changes on target gene expression, we removed genes from each individual sample analysis that had a possibility of reduced expression regardless of TF activity by potential promoter methylation-driven silencing (b > 0.15) or by genetic loss of homo-or heterozygosity [Copy Number < 1.2 in Copy Number Variation (CNV) estimation]. These genes were identified by Methylation and SNP arrays used for the same discovery cohort of samples. Methylation and CNV-corrected 1,325 TF target gene sets for each sample were used for further analysis (Fig. 1).
compared for each individual samples to the pool of normal samples. In order to elucidate TF driven differences between HPV1 and HPV2 samples, we compared gene expression differences between HPV2 and HPV1 samples for the corrected list of gene targets of each TF. We used a Wilcoxon rank sum gene set test to rank the TFs based on the expres-sion changes in their targets. While a single target gene could be either upregulated or downregulated in the particular sample or patient group, the determination of TF activity significance is agnostic to direction of dysregulation in TF gene targets, assuming the presence of both upregulated and downregulated targets in its list. We isolated the 50 TFs with Figure 1. Experimental set up. Discovery cohort, 44 HNSCC (13 HPV1 and 31 HPV2) and 25 UPPP. Exon Array data for discovery cohort samples was normalized by RMA; 2,600 human TF were depicted from TRANSFAC database and were annotated by highly relevant experimentally validated target genes. The list of TF was reduced by removal of TF with <5 targets. For the remaining 1,325 TFs, the target genes from the TF gene sets were removed if their expression was silenced by hypermethylation or DNA copy loss (using methylation array and SNP assay data) on individual sample basis. Resulting TF target gene sets were then used to compare individual samples to pool of normal by pseudo-t test to evaluate patient specific TF activity. Control-normalized TF gene sets were compared in HPV1 vs. HPV2 groups. RMAnormalized exon array data alone was used in TSP to define HPV status from TF gene sets for AP1, STAT1, STAT3 and NF-jB pathways. TCGA cohort, 279 HNSCC (35 HPV1 and 244 HPV2) and 50 matched normal. RNA-Seq data was used to evaluate individual gene expression. Target genes from the 1,325 TF gene sets were used and were corrected for genes silenced by hypermethylation or DNA copy loss (using methylation array and SNP assay data) and compared for HPV1 and HPV2 patients on an individual sample basis similar to the discovery cohort. RNA-Seq was used to validate TSP analysis. Validation cohort, 61 HNSCC (18 HPV1 and 43 HPV2) and 28 Table S5). Please note, that while we detected canonical NF-jB pathway members, we did not detect any members of NF-jB alternative pathway (NFKB2 or RELB). 20 Our data suggest that HPV1 and HPV2 HNSCC demonstrate gene dysregulation related specifically to distinct alterations in TF activity (Supporting Information Table S5). These data were recapitulated using TCGA-HNSCC cohort (Supporting Information Table S6). These results support prior reports of overall dysregulation of several common TF pathways in head and neck cancer. 21 To determine robustness to different CNV and methylation thresholds, we repeated this analysis while varying threshold values, such as b > 0.25, b > 0.35 and b > 0.45 for DNA methylation and CNV < 1, CNV < 0.8 and CNV 0 for DNA copy number. The majority of the top TFs with diverse thresholds were found in the top 50 TF list with originally applied biologically relevant b > 0.15 and CNV < 1.2 conditions, with most TFs involved in STATs, NF-jB, AP1 and retinoic acid pathways (data not shown). However, we noted that p53 was lost during threshold variations. The p53 pathway is known to be affected in HNSCC via two mechanisms, onevia direct HPV E6 oncoprotein-induced degradation of p53 protein, secondvia multiple TP53 mutations. 22,23 Point mutations of TP53 genes as well as chromosomal loss 23 and lack of adjustment for CNV correction on p53 itself may have altered the significance of p53 TF pathway dysregulation in this analysis.

Coactivation of NF-jB and STAT3 in HPV2 HNSCC
To investigate the coordinated dysregulation and the direction of this dysregulation (upregulation or downregulation) for key TFs, we performed immunohistochemical (IHC) staining experiments on a primary HNSCC tissue microarray (TMA) cohort including 100 HNSCC (consisting of 90 HPV2 and 10 HPV1) from a prior clinical trial (ECOG E4393/RTOG 9614) and 13 control samples (Supporting Information Table S2). 23 We performed an analysis of STAT3 and NF-jB protein expression for these samples (Fig.  2). The protein staining was scored using Aperio software only in the tumor cells in each sample, where staining was independently quantified for the whole cell (reflecting overall protein expression) or for the nuclei (reflecting protein acti-vation and translocation to nuclei; Supporting Information Table S7 and Fig. S2). We have evaluated the correlation of protein staining with clinical data, such as age, gender, race, cancer site, as well as with disease, smoking and HPV statuses. HPV status has the strongest correlation with protein staining. In particular HPV2 tumors have increased protein expression and significant protein activation as defined by greater than average cellular and nuclear staining of both STAT3 and NF-jB, as compared to HPV1 (Fig. 3), supporting the conclusion that these two groups are differentiated by the activity of TFs such as NF-jB and STATs (Supporting Information Tables S5 and S6). Correlation of nuclear protein staining with HPV was supported by univariate and multivariate analyses (Supporting Information Table S8). The costaining of NF-jB and STAT3 was evaluated and showed significant coactivation of NF-jB and STAT3 in primary HNSCC (Supporting Information Table S9, Spearman correlation coefficient 0.30, p 5 0.003). We found significant correlation of STAT3 and NF-jB activation (costaining in nuclei; Supporting Information Table S10). We have calculated the median of nuclear and cellular staining of STAT3 and NF-jB for the entire TMA cohort of 100 HNSCC samples (Supporting Information Table S10). We have noticed that 50 of 90 HPV2 patients (55.6%) had NF-jB nuclear staining above the median staining intensity (12.85) while none of HPV1 samples had NF-jB staining above the median (Supporting Information Table S10). Similarly, 49 HPV2 patients (54.4%) had STAT3 nuclear staining above the median (23.17) while only one HPV1 sample had STAT3 staining above the median (p 5 0.0012 and 0.0157 for NF-jB and STAT3, Supporting Information Table S10); 30 of 90 HPV2 patients (33.3%) had both NF-jB and STAT3 nuclear staining above the median, while no HPV1 patient was found to have this signature (p 5 0.0302). See Supporting Information Table S10 for details.

Target genes of transcription factors STATs, NF-jB and AP1 distinguish HPV1 and HPV2 HNSCC
To define expression biomarkers related to TF activity, we defined the TF target genes that reflect the dysregulation of these pathways and examined the ability of these genes to distinguish HPV1 and HPV2 HNSCC. We combined the lists of all targets for NF-jB, STAT1, STAT3 and AP1 factors (Supporting Information Table S4) and applied top scoring pair [TSP] 24,25 analysis to this combined list. TSP was applied to the expression array data of 72 combined target genes from the discovery cohort. TSP allowed the discovery of 5 gene pairs among 72 targets that discriminate HPV1 and HPV2 (Supporting Information Fig. S1), namely CCND1, CEBPD, ICAM1, IFG1R, IL6ST, IRF1, JAG1, JAK3, NOS3 and SOCS3. Of note, 5 out of 10 of these paired genes belong to more than one of the above dysregulated pathways (Supporting Information Table S4). The five TSPs distinguished HPV1 and HPV2 samples from the discovery cohort with a  (Table 1).
To validate our discovery of TF differential dysregulation, we adopted RNA-Seq data available for the TCGA HNSCC (The Cancer Genome Atlas) cohort. This cohort consists of 279 tumor samples, including 35 HPV1 and 244 HPV2 samples. Validation of these genes as discriminators between HPV1 and HPV2 within the TCGA cohort also demonstrated significant separation with a p-value of 6.7 3 10 28 OR 5 8.9 [3.6 26, 95% CI] ( Table 1).

Validation of transcription factor target genes by qRT-PCR
To validate our results and evaluate the discriminative ability of the top five gene pairs from the TSP analysis, we assembled an independent validation cohort of 61 HNSCC (including 43 HPV2 and 18 HPV1 samples) and 28 control uvulopalatopharyngoplasty (UPPP) samples (Supporting Information Table S3). The TSPs based on expression of the top five gene pairs were analyzed by qRT-PCR, which was able to separate HPV1 and HPV2 HNSCC from this validation cohort with a p-value of 0.0006, OR 5 9.6 [2.2 60] ( Table 1). Using qRT-PCR data for 5 top scoring gene pairs from TSP, this validation cohort (Supporting Information  Table S3) was separated into three groups by unsupervised hierarchical clustering (Fig. 4 heat-map was built for the purpose of visualization). We found that 25 out of 61 samples were clustered into a separate group that was enriched by HPV1 patients. That group contains 16 out of 18 HPV1 patients and is distinct from two other primarily HPV2 groups by downregulation of CCND1 and partial downregulation of IRF1, ICAM1, IGF1R and CCND1 (p < 3.5 3 10 25 ). These data suggest that STATs, NF-jB and AP1 pathways are upregulated in HPV2 HNSCC patients as compared to HPV1 patients, in concordance with our immunohistochemical data. Two groups of primarily HPV2 patients can be distinguished by relative expression of TSP genes, with one group showing higher expression of CCND1, CEBPD, ICAM1, IRF1, IL6ST, JAG1, JAK3, NOS3 and SOCS3. These data suggest that dysregulation of STATs, AP1 and NF-jB pathways is coordinated and that these pathways are co-activated for a subset of HNSCC patients (25%, 15 out of total 61 samples), predominantly HPV2 patients (p < 0.05 by Fisher's exact test applied to hierarchically clustered groups, Fig. 4). As HPV2 tumors have been associated with poor clinical outcome, the detected gene signatures can potentially aid clinical assessment of HNSCC prognostic subgroups.

NF-jB and STAT coactivated HPV2 HNSCC cells are growth-inhibited by combined targeted therapy
We chose Ho1N1 and HSC2 from total 34 cell lines described in Ref. 26 after expression analysis of key TFs and their significant targets (detected via TSP analysis); these cell lines were selected as representative of STAT3-NF-jB overexpressing or base-line expressing tumors based on expres-sion pattern. To evaluate the clinical relevance and potential therapeutic susceptibility of NF-jB and STATs co-activation in HPV2 HNSCC, we knocked down the expression of these TFs in vitro by RNAi techniques. Ho1N1 cells are characterized by over-expression of STAT1, STAT3, RELA and NFKB1 and their signature targets as defined by our TSP analysis (IRF1, CEBPD, CCND1, ICAM1, JAG1, JAK3 and NOS3). Since our data demonstrate that expression of signature target genes discovered as five TSPs reliably describes the activity of key TF pathways, we conclude that Ho1N1 has hyperactivated STAT3 and NF-jB pathways. Conversely, HSC2 showed low baseline expression of STATs, NF-jB genes and their significant targets. Of note, JAK3 expression was not detected in HSC2 and the expression of JAG1 and NOS3 was not significantly different in these two cell lines (Supporting Information Fig. S2).
The siRNA SMARTpools for STAT1, STAT3 and RELA allowed significant inhibition of target gene expression in each cell line (Supporting Information Fig. S3). Growth inhibition was associated with induction of the cell death by downregulation of NF-jB, shown for wt TP53 head and neck cell lines (HPV2). 27 Interestingly, significant growth inhibition was noticed only in the STAT and NF-jB activated Ho1N1 cell line, but not in the HSC2 cell line (Fig. 5). Combined downregulation of two genes had additive effect on cell growth, suggesting that these pathways were downregulated independently. The combined inhibition of STAT1 and STAT3 expression had the strongest effect on Ho1N1 cell line proliferation (62% cell growth inhibition, as compared to control scrambled siRNA pool); and second strongest growth inhibition was archived with double STAT3-RELA inhibition (56% cell growth inhibition). The Ho1N1 cell line was the most sensitive to chemical NF-jB and STAT3 inhibitors (Bortezomib, Bay 11-7085, Cucurbitacin and SH-4-54; Supporting Information Figs. S4 and S5 and Refs. 12 and 28-30), supporting our siRNA experiments. NF-jB and STAT3 inhibitors were not used in combination due to the high cellular toxicity of the individual agents. No effect of siRNA inhibition was noted in the low STAT-NF-jB expressing HSC2 cell line (Fig. 5b). HSC2 was the most resistant to chemical NF-jB and STAT3 inhibitors (Supporting Information Figs. S4 and S5), supporting our siRNA experiments.  Fig. S6), with maximal cell growth inhibition with use of combined STAT3 and RELA siRNA (42% cell growth inhibition, Supporting  Information Fig. S6a). SKN3 has an intermediate response to both NF-jB inhibitors (Supporting Information Fig. S4). As for STAT3 inhibitors, SKN3 was sensitive to Cucurbitacin and resistant to SH-4-54 (Supporting Information Fig. S5), suggesting an overall intermediate response to STAT3 inhibitors. The diverse response of SKN3 to STAT3 inhibitors may  Cancer Genetics Gaykalova et al. be explained by high toxicity of Cucurbitacin and low specificity of SH-4-54.

Discussion
Over the past decade, there has been recognition that HNSCC includes favorable prognosis HPV1 and poor prognosis HPV2 patients. 31 The discovery of novel therapeutic agents for poor prognosis HPV2 HNSCC has been challenging, despite use of high throughput sequencing efforts to define novel therapeutic targets. TF are significant drivers of gene expression variation and characterization of TF alterations facilitates insight into HNSCC biology; however, reliable changes in TF activity are hard to directly detect due to the complexity of TF signaling. Using three high throughput platforms and an innovative integrative strategy that compensates for genetic and epigenetic changes in TF targets, we were able to annotate highly specific target genes to each known TF and infer whole-genome TF activity by the expression of those target genes in each particular tumor sample. We have found that HPV2 patients have strong coordinated dysregulation of several pathways, including coordinated activation of NF-jB and STATs pathways. Our results are supported by similar results shown by other investigators using a cDNA microarray analysis of 10 HPV2 HNSCC cell lines. 32,33 We confirmed our results on multiple validation datasets, including an ECOG TMA cohort, a TCGA-HNSCC cohort and independent HNSCC validation cohort, by protein and expression assays, as well as by confirmation of expression alteration for direct targets of specific TFs. The discovered coactivation of NF-jB and STAT3 in primary HNSCC tissue is also in agreement with prior published IHC data. 34 We have also shown that coactivation of NF-jB and STAT3 provide an avenue for selective therapeutic targeting in cell line models. We have demonstrated that these alterations occur in a substantial subset of HNSCCs. For example, 15 out of total 61 HNSCC patients (25%) or 14 out of total 43 HPV2 HNSCC patients (33%) have higher expression of direct targets of STATs, NF-jB and AP1 TFs, such as CCND1, CEBPD, ICAM1, IRF1, JAG1, JAK3 and NOS3 (Fig.  4). The direct dependence of these targets on TFs activity strongly correlates with the previous data showing NF-jB target activation in HNSCC. 27 Prior reports have described changes in the activity of the several TFs, like STAT1, STAT3, NF-jB and AP1, in HPV2 HNSCC, as well as in HPV mixed HNSCC population and HPV1 cervix patients. [35][36][37][38][39] For example, HPV16 E7 was shown to inhibit NF-jB activity. 40 Additional TF pathway activity alterations in HNSCC have been shown for p53, 22 RB1 23 and retinoic X receptor 41 cascades. There are many possibilities for these pathways to interact in HNSCC. Many of these TFs have common upstream regulators: including dependence of NF-jB, STATs and AP1 activation on IKKa, IKKb and on cytokine signals. 20,42 They also share common downstream targets: for example, STAT1 and NF-jB regulate CD40 expression 43,44 ; AP1 and NF-jB regulate IL6, IL8 and VEGF expression. 34 Finally, these TFs can also affect their own expression and expression of each other: NF-jB plays a role in paracrine-dependent activation of STAT3 pathway. 34  suggest that coactivation of NF-jB and STAT3, especially in the context of dysfunctional p53 will enhance BCL-XL expression and lead to HNSCC cell survival. 45 While direct ChIP-Seq experiments are not feasible to perform on the primary HNSCC samples due to the limited tumor availability, we expect overlap in TF binding to their target genes, consistent with the literature cited above.
The overlapping gene targets of p53, p63 and p73 TFs and their ability to distinguish HPV1 and HPV2 patients is summarized in Supporting Information Table S11. Notably, the mutation status of TP53 gene has a mixed downstream effect on the activity of its target genes (Supporting Information Table S12). The TP53 mutation status was available for 37 out of 44 HNSCC samples in the discovery cohort from previous deep-sequencing study, 46 where 25 samples were found to have wild-type TP53 and 7 samples were found to have disruptive TP53 mutations (2 samples with nonsense mutations and 5 samples with missense hot spot mutations: R175H, Y220C, R248Q, R273C and R282W). The diverse effect of disruptive TP53 mutations can be explained by downregulating effect of nonsense mutations and activating effect of gain-offunction hot spot mutations. Such results correlate with previously published data, summarized in Ref. 47.
In general, TFs were considered "undruggable" for many years, but with improved understanding of transcription factor biology and new insight into pathway analysis, as well as novel technologies in the development of targeting agents and agent delivery, there is a broad venue for treatment based on TFs and their pathways. 48,49 We have demonstrated that therapeutic agents directed against STATs and NF-jB, such as Cucurbitacin, 29 Bortezomib, 12 Bay 11-7085 28 and SH-4-54 30 have high toxicity and limited specificity (See Supporting Information Figs. S4 and S5 for details). However, small peptides and peptidomimetics have shown lower cellular toxicity and better specificity to decrease transcription activity of STAT3, NF-jB, MYC and others. 50 Oligonucleotide-based therapy, including RNA interference agents and singlestranded decoys have demonstrated specific downregulation of targeted pathways and are actively used in preclinical and early clinical studies to facilitate AP1, RB1, ERs, NF-jB and STATs downregulation. 13,14,51 The correlation of our findings with previous published data suggests that the inferential technique we have developed to define TF activity is robust and a useful tool for discovery of genome-wide changes in TF activity. While we evaluated only the most prominent TF pathways, such as STATs, NF-jB and AP1, many more TFs were found to be significantly and differentially activated in two HNSCC groups. The use of an inferential strategy to deduce activity of specific TFs also has the potential for application in therapeutic targeting, as we were able to use a list of ten target genes that helped us to validate TF signatures experimentally. This strategy could potentially be adapted for clinical practice to identify patients with dysregulated (upregulated or downregulated) pathways. Characterization of a limited list of TF target genes allowed us to accurately characterize TF signatures experimentally and could potentially be applied to define patients with tumors that are susceptible to specific TF pathway targeting agents as we have done in cell line systems, albeit with potential limitations in translating cell line responses to primary tumor responses to therapy. In addition, we used a correction strategy to infer TF signatures despite heterogeneity of genetic and epigenetic alterations that could potentially mask TF pathway alteration. This strategy has potential to clarify pathway alterations in specific individual tumors within tumor types that exhibit significant heterogeneity in terms of genetic and epigenetic alterations.
In summary, the use of high throughput integrated expression, methylation and genome copy number analysis to infer TF activity in HNSCC provides insight into biologic variability of behavior and treatment response for HPV1 and HPV2 HNSCC patients and may potentially be used to direct targeted therapy based on TF pathway analysis.