HPV E6 regulates therapy responses in oropharyngeal cancer by repressing the PGC-1α/ERRα axis

HPV⁺ cancers comprise over 5% of malignancies worldwide, and rising incidence of HPV⁺ oropharyngeal squamous cell carcinoma (OPSCC) has made this the most common HPV-related cancer in the United States (1). HPV⁺ OPSCCs are typically sensitive to the combination of radiation plus cisplatin, which is often curative for this disease (2). Impaired DNA repair mechanisms help make HPV⁺ OPSCCs sensitive as a group to these therapies (3), which kill tumors by increasing reactive oxygen species (ROS) generation beyond cellular antioxidant capacity, leading to uncontrolled damage to DNA, protein, and lipids (4). High cure rates for HPV⁺ OPSCCs treated with standard therapy have sustained interest in de-escalating cisplatin and radiation therapy, whose toxic sequelae leave lifelong disabilities in survivors (5). However, attempting to supplant cisplatin with anti-EGFR therapy was unsuccessful (6, 7), and radiation de-escalation efforts remain a work in progress. The ability to personalize therapy based on more accurate predictive biomarkers would greatly enhance treatment de-escalation efforts for this disease; however, such innovation remains hindered by limited understanding of the mechanisms underlying divergent responses among HPV⁺ OPSCCs to cisplatin and radiation.

Variation among HPV⁺ OPSCCs in ability to neutralize oxidative stress could contribute to differences in treatment response and survival. Mitochondria contain a large portion of a cell’s antioxidant machinery and not only neutralize ROS byproducts of oxidative metabolism but also allow cancers to survive extrinsic stressors such as hypoxia and nutrient deprivation (4). NADPH generation within mitochondria fuels the glutathione peroxidase and peroxiredoxin antioxidant systems, which are also partly intrinsic to mitochondria (8). Mitochondrial biogenesis and related antioxidant capacity are upregulated by oncogenic drivers in the RAS pathway (9) and contribute to the aggressive phenotype of RAS-driven cancers by mitigating oxidative stress (10). Less is known about mitochondrial metabolism’s contribution to poor outcomes among HPV⁺ OPSCCs, which are primarily driven by viral factors and show greater vulnerability as a group to oxidative damage. However, there is emerging evidence that HPV⁺ OPSCC progression may be enhanced by upregulation of oxidative metabolism by diverse mechanisms (11, 12).

A direct role for HPV oncoproteins in regulating mitochondrial function has not been demonstrated to our knowledge. Many viruses have evolved to evade the type I interferon–mediated innate immune responses initiated on the mitochondrial membrane through RIG-I/MDA5/MAVS signaling (13). Viruses have developed mechanisms to inhibit these signals, disrupt mitochondrial membrane potential, and/or impair mitochondrial function by altering fission-fusion dynamics (13). HPV 16’s E6 and E7 oncoproteins jointly suppress innate immune gene expression (14), but it is unknown whether this phenomenon is mediated through effects upon mitochondria. Any suppressive effects of E6 and E7 on mitochondrial function to serve HPV replication might secondarily limit growth of HPV⁺ OPSCCs in hostile microenvironments and sensitize to cisplatin and radiation. Under this hypothesis, the variable viral oncoprotein levels among HPV⁺ OPSCCs are predicted to impact treatment responses.

This study evaluated whether variations in HPV oncoprotein levels among HPV⁺ OPSCCs contribute to diversity in their responses to cisplatin and radiation. We found that high oxidative metabolic gene expression is associated with low HPV E6 levels and reduced survival, whereas high E6 levels correlate with decreased mitochondrial function and sensitivity to therapy. Genetic manipulation of E6 expression in HPV⁺ cancer cell lines revealed that E6 can suppress mitochondrial biogenesis by depletion of p53, a transcription factor for expression of transcriptional co-activator peroxisome proliferator–activated receptor gamma co-activator 1α (PGC-1α) (15). PGC-1α acts with its DNA-binding nuclear partner, estrogen-related receptor α (ERRα), as a master transcriptional regulator for mitochondrial biogenesis (16). Our findings uncover clinically relevant variations in E6 levels and PGC-1α/ERRα pathway activation that provide a framework for personalizing therapy.

High oxidative metabolic gene expression is associated with decreased HPV⁺ OPSCC survival. Three cohorts were identified that offer publicly available RNA-Seq data from treatment-naive HPV⁺ OPSCCs with survival outcomes annotation. In The Cancer Genome Atlas (TCGA), RNA-Seq data are available from 53 OPSCCs that express high-risk HPV transcripts (17, 18). We used this data to assess the relationship between oxidative metabolic gene expression and overall survival (OS). Expression was quantified for the Molecular Signatures Database’s Hallmark_Oxidative_Phosphorylation gene set (19), which contains 200 genes involved in oxidative metabolism and/or the TCA cycle. Upregulation for each gene was defined as expression ≥ 1 standard deviation above the mean log₁₀-transformed level for all cases. The number of upregulated transcripts was used to divide the cases into tertiles. The lowest tertile had better 3-year OS than the highest (P = 0.019) and intermediate tertiles (P = 0.043) (Figure 1). This finding was tested in the JHU HPV⁺ OPSCC cohort (20). Dividing this cohort (n = 47) into tertiles by the identical methodology used for TCGA cases (Figure 1) showed better 3-year OS for the lowest versus the highest tertile (P = 0.048). Wide variation in therapy was apparent from characteristics of the JHU and TCGA cohorts (Supplemental Table 1; supplemental material available online with this article; https://doi.org/10.1172/jci.insight.159600DS1), making the relationship between survival and response to radiation plus chemotherapy unclear. Thus, we examined a subcohort of patients within the broader VU HPV⁺ OPSCC cohort (21) who received only radiation plus cisplatin as primary therapy (n = 37). Despite small sample size, dividing these cases into tertiles revealed better 3-year OS (P = 0.041) in the lowest versus the highest tertile (Figure 1). This association of high expression of mitochondrial components with decreased survival across 3 cohorts led to the hypothesis that high mitochondrial mass contributes to the resistance to cisplatin and radiation seen in a minority of HPV⁺ OPSCCs.

Figure 1

Association of high tumor oxidative metabolic gene expression with reduced patient survival. Three HPV⁺ OPSCC patient cohorts were divided into tertiles based on the number of upregulated Hallmark_Oxidative_Phosphorylation Gene transcripts in each cohort (top). Upregulation was defined as log₁₀-transformed expression of 1 standard deviation above the mean for each cohort. When a dividing line between tertiles spanned multiple patients with equally upregulated transcripts, those patients were grouped into the higher tertile. Overall survival was estimated by the Kaplan-Meier method (bottom), and the log-rank test was used for pairwise comparisons among tertiles. JHU, Johns Hopkins University; VU, Vanderbilt University.

Linking features of poor prognosis patients to high mitochondrial function and cisplatin resistance. The oxidative metabolic gene expression seen in patients with reduced survival was evaluated in patient-derived xenograft (PDX) and cell line models for prediction of increased mitochondrial mass, mitochondrial function, and cisplatin response. RNA-Seq was performed for the 7 HPV⁺ OPSCC PDX models previously established by us from treatment-naive patients (22). Upregulated oxidative metabolic transcripts in the PDXs were quantified by the methodology used in Figure 1. Mitochondrial mass in the PDXs was estimated using levels of mitochondrial DNA measured by a quantitative PCR (qPCR) assay that defines the ratio of a mitochondrial gene, MT-CO1, to a nuclear gene, B2M, and then infers mitochondrial DNA levels from a standard curve (23) (Supplemental Figure 1A). A positive linear correlation between the number of upregulated transcripts in the PDXs and MT-CO1/B2M (Figure 2A) supported that our method of quantifying oxidative metabolic gene expression effectively estimates mitochondrial mass. Cisplatin treatment of the PDXs was performed in vivo to define the relationship between mitochondrial mass and therapy response. Cisplatin resistance was quantified using rate-based tumor/control (T/C) values, and this measures growth of treated tumors relative to controls using an exponential model that incorporates all data points in both growth curves (24) (Supplemental Figure 1B). Rate-based T/C values in the PDXs grown in NOD/SCID/IL-2Rγ^−/− (NSG) mice showed a positive linear correlation with mitochondrial mass defined by MT-CO1/B2M ratio (Figure 2A). Relationships among mitochondrial mass, mitochondrial function, and cisplatin resistance were further evaluated in a panel of 7 HPV⁺ cancer cell lines of head and neck origin. Oxidative phosphorylation was quantified by measuring basal oxygen consumption rate (OCR) using the Cell Mito Stress Test on the Seahorse XF platform. A strong positive correlation was shown between mitochondrial mass and basal OCR in the cell lines (Figure 2B and Supplemental Figure 1C). In addition, NADPH/NADP⁺ ratio was used to quantify antioxidant capacity and showed strong positive correlation with mitochondrial mass (Figure 2B). Last, cisplatin IC₅₀ in the cell lines showed a linear positive correlation with MT-CO1/B2M ratio, basal OCR, and NADPH/NADP⁺ ratio (Figure 2C). These findings in preclinical models show that the elevated oxidative metabolic gene expression in poor prognosis patients correlates with high mitochondrial mass, oxidative metabolism, and antioxidant capacity. This analysis provided models for subsequent experiments to pursue mechanistic determinants of these metabolic features and assess their impact on therapy response.

Figure 2

Correlation analyses among oxidative metabolic gene expression, mitochondrial mass, oxidative phosphorylation, and cisplatin resistance in HPV⁺ OPSCC models. (A) Mitochondrial mass was measured in a panel of HPV⁺ PDXs (n = 7) by DNA qPCR to define MT-CO1/B2M. Correlation analysis was performed for MT-CO1/B2M versus number of upregulated Hallmark_Oxidative_Phosphorylation genes by RNA-Seq (left) and versus cisplatin resistance in vivo (right) measured using rate-based T/C value. (B) Correlation analysis was performed in cell lines (n = 7) for MT-CO1/B2M versus basal OCR by Seahorse assay (left) and versus NADPH/NADP⁺ measured by enzyme cycling-based colorimetric assay (right). (C) Correlation analysis was performed in cell lines for cisplatin IC₅₀ by WST-1 assay versus MT-CO1/B2M, basal OCR, and NADPH/NADP⁺. Pearson correlation coefficients were used to calculate r values for scatterplots with confidence intervals, and P values were determined by t distribution.

Among HPV oncogenes, high full-length E6 is associated with low mitochondrial mass and function. Diversity in viral transcript levels among HPV⁺ OPSCCs led us to evaluate potential associations between viral gene expression and mitochondrial mass. Among the HPV⁺ OPSCCs in TCGA, wide variation was apparent among levels of the early viral transcripts with known roles in malignant transformation (E2, E4, E5, E6, E7) (Figure 3A). The highest coefficient of variance was seen for the full-length, nonspliced form of E6 (fl-E6), a major cancer driver that is the only form of E6 that targets p53. Fl-E6 levels varied by more than 200-fold and were reduced in the highest tertile for oxidative metabolic gene expression relative to the lowest tertile (Figure 3B). This association was lost when total E6 transcripts or E6* were considered (Figure 3B). No other oncogenic viral transcripts showed association with oxidative metabolic gene expression in TCGA (Supplemental Figure 2A). Although E6 splicing data were not obtainable for the JHU and VU patient cohorts, associations similar to those observed for fl-E6 in TCGA were also present in experimental models of HPV⁺ OPSCC. In the PDXs, fl-E6 showed a negative linear correlation with upregulation of oxidative metabolic transcripts and mitochondrial mass (Figure 3C), whereas no correlations were detected for other oncogenic viral transcripts (Supplemental Figure 2B). Similarly, fl-E6 was negatively correlated with MT-CO1/B2M ratio, basal OCR, and NADPH/NADP⁺ in the cell lines (Figure 3D), whereas such relationships were not apparent for other oncogenic HPV transcripts (Supplemental Figure 2C). These findings reveal a reciprocal relationship between viral oncogene expression and both mitochondrial mass and mitochondrial function that is specific to fl-E6.

Figure 3

Fl-E6 levels are associated with low mitochondrial mass and function in patient tumors and models. (A) RNA-Seq data for the HPV⁺ OPSCCs in TCGA (n = 53) showing coefficients of variation among the mRNAs for viral oncogenes (left) and the range of fl-E6 expression across the cohort (right). Data shown as single value. (B) Comparison of fl-E6, total E6, and spliced E6 (E6*) RNA-Seq levels (right) for the intermediate and highest tertiles compared with the lowest tertile of the Hallmark_Oxidative_Phosphorylation gene expression in TCGA. Adjusted P values were calculated using a 1-way Welch’s ANOVA with Dunnett’s multiple comparisons procedure. Mean ± SD. (C) PDX fl-E6 mRNA levels by RNA-Seq versus number of upregulated Hallmark_Oxidative_Phosphorylation Genes (left) and versus mitochondrial mass (right). (D) Cell line fl-E6 levels normalized to 18S by quantitative real-time PCR (qRT-PCR) versus mitochondrial mass (MT-CO1/B2M), basal OCR (Seahorse assay), and antioxidant capacity (NADPH/NADP⁺). Pearson correlation coefficients for scatterplots were used to calculate r values with confidence intervals, and P values were determined by t distribution.

Increasing fl-E6 levels represses mitochondrial biogenesis and depletes antioxidant capacity. Fl-E6’s association with reduced mitochondrial mass and function prompted the question of whether manipulating E6 levels would alter mitochondrial biogenesis and oxidative phosphorylation. Because selective genetic silencing of E6 is hindered by E6 and E7 arising from the same transcript, we tested our hypothesis by increasing fl-E6 expression in VU147T and SCC154, the 2 HPV⁺ head and neck cancer lines with highest mitochondrial mass and lowest fl-E6 (Figure 3D). Lentiviral transduction of HPV 16 fl-E6 was performed with a nonspliceable fl-E6 form (V42L) expressing a hemagglutinin (HA) tag (25) and used to generate 3 fl-E6–overexpressing clones and 3 vector control clones for each cell line. The 3 clones expressed fl-E6 mRNA at 60- to 300-fold (Supplemental Figure 3A) over endogenous levels, thus modeling the range of differences between the human tumors in TCGA with highest versus lowest fl-E6 (Figure 3A). Fl-E6 overexpression did not alter in vitro growth (Supplemental Figure 3B) but produced the anticipated decrease in levels of E6’s canonical p53 target (Figure 4A). Increasing fl-E6 broadly repressed mRNA levels for 6 genes involved in oxidative metabolism: SDHB (complex II), UQCRC2 (complex III), MTCO2P12 (complex IV), ATP5F1A (complex V), ACADM (acyl-CoA dehydrogenase medium chain-1), and MDH1 (malate dehydrogenase-1) (Supplemental Figure 3C), and the same effects were shown to be present in the initially transfected pools prior to clone generation (Supplemental Figure 3D). Similarly, 5 mitochondrial components examined by Western blot using a commercial antibody cocktail (NDUFB8, SDHB, UQCRC2, MTCO2P12, and ATP5F1A) were reduced by fl-E6 overexpression (Figure 4A and Supplemental Figure 3E). To determine whether these effects were generalizable beyond malignancy, fl-E6 was transduced by lentivirus into N-tert/1 keratinocytes expressing E7 (N-tert/E7 keratinocytes) (26). Three biologic replicates expressing fl-E6 did not show altered in vitro growth (Supplemental Figure 3, F and G), but reduction in expression of mitochondrial components was observed at the mRNA and protein level (Supplemental Figure 3, H and I). Further characterization of VU147T and SCC154 cells overexpressing fl-E6 supported reductions in MT-CO1/B2M (Figure 4B) and basal OCR (Figure 4C). Increased fl-E6 levels also depleted mitochondrial antioxidant capacity, as reflected by reduced NADPH/NADP⁺ ratio (Figure 4D). These results demonstrate fl-E6’s ability to attenuate mitochondrial mass, oxidative phosphorylation, and antioxidant capacity in HPV⁺ cancer cells and nontransformed cells expressing E7.

Figure 4

Increasing fl-E6 levels represses mitochondrial biogenesis and antioxidant capacity. (A) Western blot of parental cells, 3 lentiviral fl-E6 transduced clones, and 3 vector control clones showing levels of HA-tagged E6, p53, and 5 electron transport chain components (top) and their densitometric quantification (bottom). Adjusted P values were calculated using a 1-way ANOVA with Holm-Šidák multiple comparisons procedure. The effect of lentiviral fl-E6 expression is shown for (B) mitochondrial mass (MT-CO1/B2M ratio by DNA qPCR), (C) basal OCR (Seahorse assay), and (D) antioxidant capacity (NADPH/NADP⁺ by enzyme cycling-based colorimetric assay). Super-plots display mean ± SEM of 3 technical replicates for each of 3 biologic replicates. Different colors represent biologic replicates. Light colors represent technical replicate values, and corresponding dark colors represent their means. P values are based on unpaired t tests.

Depletion of mitochondrial antioxidant capacity by fl-E6 potentiates therapy responses. The efficacy of radiation and cisplatin for many HPV⁺ OPSCCs indicates their limited capacity to maintain redox balance upon increasing ROS. We hypothesized that depletion of mitochondrial antioxidant capacity by fl-E6 contributes to therapy sensitivity in this setting. To test this possibility, radiation dose responses were measured in colony-forming assays for VU147T and SCC154 cells, which contain low endogenous fl-E6. Lentiviral fl-E6 overexpression enhanced radiation dose responses (Supplemental Figure 4A) and decreased LD₅₀ (Figure 5A). The same effect occurred with a chemical ROS inducer, 2,3-dimethoxy-1,4-naphthoquinone (DMNQ) (Figure 5B and Supplemental Figure 4B). Comparable effects of fl-E6 on cisplatin dose responses (Supplemental Figure 4B) and IC₅₀ were shown by WST-1 assay (Figure 5C) and confirmed using a BrdU incorporation assay (Supplemental Figure 4C), which avoids the potential for mitochondrial mass differences to confound cell number measurement in WST-1 assays. Because cisplatin may act by multiple mechanisms beyond ROS induction (27), we tested whether cisplatin sensitization by fl-E6 depended on mitochondrial antioxidant capacity depletion. MitoTEMPO, a reagent that increases antioxidant capacity selectively in mitochondria at a 10 μM dose (28), fully reversed the cisplatin sensitivity created in VU147T and SCC154 cells by fl-E6 (Figure 5D and Supplemental Figure 4D). A partial reversal of cisplatin sensitivity was seen using the weaker antioxidant MnTMPyP, a cell-permeable mimetic of manganese superoxide dismutase, at its standard 100 μM dose (29) (Figure 5D and Supplemental Figure 4D). Effects of increased fl-E6 on in vivo growth and cisplatin responses were evaluated in xenografts of SCC154 and VU147T. In vivo growth was maintained by both cell lines after lentiviral fl-E6 transduction, but, in contrast to in vitro (Supplemental Figure 3B), slowed to variable extents (Figure 5E). Solid tumors from fl-E6-transduced cells showed fl-E6 mRNA expression at 15- to 60-fold the endogenous levels and maintained repression of oxidative metabolic transcript expression (Supplemental Figure 4E). Control tumors continued growth during treatment, confirming cisplatin resistance of these 2 low fl-E6 cell lines in vivo, and stable Fl-E6 expression markedly sensitized both to cisplatin treatment (Figure 5F and Supplemental Figure 4F). These results provide evidence that antioxidant capacity depletion by high fl-E6 levels contributes to the treatment sensitivity in HPV⁺ cancer; they further suggest that attenuated fl-E6 expression confers therapy resistance to certain tumors and enhances in vivo growth by mitigating oxidative stress.

Figure 5

Depletion of mitochondrial antioxidant capacity by fl-E6 potentiates therapy responses. In vitro LD₅₀ and IC₅₀ were measured by simple logistic regression. (A) LD₅₀ for radiation response is derived from number of cells surviving 10 days after 2, 4, and 6 Gy. IC₅₀ is quantified by WST-1 assay upon treatment with (B) DMNQ (C) or cisplatin and by BrDu assay upon treatment with (D) cisplatin ± 10 μM MitoTEMPO or ± 100 μM MnTMPyP. Values represent mean ± SEM for 3 biologic replicates. P values are from unpaired t tests. (E) In vivo subcutaneous tumor growth is shown in presence and absence of lentiviral fl-E6 expression. (F) Treatment responses for both conditions are shown during intraperitoneal treatment with 1 mg/kg cisplatin 3 times/w for 2 weeks. P values of interaction were based on 2-way ANOVA.

Fl-E6 represses the PGC-1α/ERRα axis by reducing p53-dependent PGC-1α promoter activation. The central role of transcriptional co-activator PGC-1α in mitochondrial biogenesis led us to evaluate its potential role in fl-E6–mediated depletion of mitochondrial mass. Silencing of PGC-1α with siRNA in VU147T and SCC154 cells reduced mRNA levels of mitochondrial components to a similar degree as transient overexpression of fl-E6, and overexpressing PGC-1α in addition to fl-E6 fully reversed the effect of fl-E6 alone (Supplemental Figure 5). Similarly, stable fl-E6 expression in VU147T and SCC154 cells reduced mRNA levels for PGC-1α (Figure 6A) but not its homolog PGC-1β (Supplemental Figure 6A) and concomitantly decreased PGC-1α protein (Figure 6B and Supplemental Figure 6B). Levels of ERRα, PGC-1α’s obligate binding partner whose expression is positively regulated by PGC-1α in a feed-forward loop (16), also decreased at the mRNA and protein level (Supplemental Figure 6C). PGC-1α levels similarly decreased in the N-tert/E7 keratinocytes with fl-E6 expression (Supplemental Figure 6, D and E), showing E6’s effect on PGC-1α to be generalizable to other cell contexts. Therefore, effects of fl-E6 on PGC-1α gene transcription were assessed using a luciferase reporter linked to the PGC-1α promoter. The reporter indicated a 50%–75% decrease in PGC-1α promoter activity in VU147T and SCC154 cells upon fl-E6 overexpression (Figure 6C), consistent with p53’s direct positive regulation of the PGC-1α promoter in cancer cells (15). To discriminate among the multiple targets of E6 that might affect PGC-1α transcription, the same reporter was transfected into HEK293 cells along with mutant forms of E6 (Figure 6D). As anticipated, the nonspliceable fl-E6 construct and HPV 16 wild-type E6 comparably repressed reporter function, whereas the most common spliced form (E6*I) did not. PGC-1α promoter repression was maintained upon deleting the C-terminal PDZ-binding motif (E6Δ146-151) (30), showing E6’s effect to be independent of the multiple E6 target interactions mediated by PDZ domains. By contrast, E6 mutants with reduced capacity to bind p53 (8S9A10T) or reduced E6AP interaction (I128T) had reduced ability to repress the reporter. Because 8S9A10T and I128T only partly prevent p53 degradation and may not be p53-selective among E6’s targets (31), p53’s role was further evaluated using congenic HCT116 cells that are p53 wild-type versus null. The p53^–/– background reduced reporter activity relative to that in p53^+/+ HCT116 cells, and fl-E6 repressed reporter function in p53^+/+ cells but not in p53^–/– cells (Figure 6E). These findings demonstrate fl-E6’s capacity to antagonize PGC-1α expression and function and indicate a large portion of this activity is mediated by reduced p53.

Figure 6

Fl-E6 represses PGC-1α levels by reducing p53-dependent PGC-1α promoter activity. (A) qRT-PCR for PGC-1α normalized to 18S after lentiviral fl-E6 expression in HPV⁺ cancer cell lines. Super-plots represent mean ± SEM for 3 biologic replicates with 4 technical replicates. P values are from unpaired t tests. (B) Western blot for PGC-1α in 3 lentiviral fl-E6 transfected clones and 3 vector control clones. (C) PGC-1α promoter-driven luciferase reporter activity in HPV⁺ cancer cell lines in presence and absence of lentiviral fl-E6 expression. P values are based on paired t tests. (D) HEK293 cells cotransfected with the PGC-1α promoter-driven luciferase reporter plus the shown mutants of HPV 16 E6 (top) before measuring luciferase activity (bottom). (E) TP53-WT and -null HCT116 cells were cotransfected with the PGC-1α promoter-driven luciferase reporter plus fl-E6 or empty vector. Reporter activity is normalized using background from a vector control reporter and expressed as mean ± SEM for 3 technical replicates. All results are representative of at least 3 independent experiments. Adjusted P values were calculated using a 1-way ANOVA with Holm-Šidák multiple comparisons procedure.

PGC-1α/ERRα axis activation in HPV⁺ OPSCCs with high mitochondrial mass and cisplatin resistance. Our observation that increasing E6 represses mitochondrial biogenesis by reducing p53-dependent PGC-1α expression led to us to evaluate this phenomenon for physiologic relevance in PDXs and patients, where fl-E6 levels vary widely. PGC-1α mRNA in the patient cohorts and PDX panel was captured by RNA-Seq at very low levels, often not meeting the detection threshold, and thus was instead quantified in the PDXs by qRT-PCR. As predicted, PGC-1α mRNA showed a positive linear correlation with the number of upregulated Hallmark_Oxidative_Phosphorylation transcripts and MT-CO1/B2M ratio (Figure 7A). Likewise, PGC-1α mRNA in the PDXs correlated negatively with fl-E6 mRNA and positively with protein levels of E6’s degradation target p53 (Figure 7B and Supplemental Figure 6F). The predicted negative correlation between PGC-1α expression and cisplatin resistance measured by rate-based T/C was confirmed in the PDXs (Figure 7B). To overcome technical inability to detect PGC-1α in the clinical cohorts by RNA-Seq, we examined as a surrogate marker its obligate binding partner ERRα, whose role as transcriptional target of PGC-1α/ERRα creates a feed-forward loop that amplifies PGC-1α’s activity (16). ERRα was quantifiable in the cohorts in Figure 1 where RNA-Seq was performed using frozen tissue (TCGA and JHU) but not FFPE (VU), thus providing a readout for PGC-1α/ERRα axis activation. A positive linear correlation was observed between ERRα levels and number of upregulated Hallmark_Oxidative_Phosphorylation transcripts in the PDX panel and both patient cohorts (Figure 7C). Moreover, ERRα mRNA correlated negatively with fl-E6 in TCGA, the only cohort with fl-E6 annotation (Figure 7D). Last, ERRα mRNA levels positively correlated with the number of upregulated transcripts in the Hallmark_P53_Pathway gene set in both cohorts (Figure 7E), with upregulation defined by the same criteria used for Hallmark_Oxidative_Phosphorylation genes in Figure 1. In the context of our mechanistic findings in cell lines, these concordant observations in PDXs and patients support the concept that variable fl-E6 levels among HPV⁺ OPSCCs lead to differential PGC-1α/ERRα pathway activity by a p53-mediated mechanism, in turn causing variability in mitochondrial mass that affects therapy responses.

Figure 7

PGC-1α/ERRα axis activation in tumors with high mitochondrial mass and therapy resistance. PGC-1α was measured in HPV⁺ OPSCC PDXs by qRT-qPCR normalized to 18S. All other transcripts were measured by RNA-Seq. Pearson correlation coefficients for scatterplots were used to calculate r values with confidence intervals, and P values were determined by t distribution. (A) PDX PGC-1α mRNA versus the number of upregulated Hallmark_Oxidative_Phosphorylation transcripts (left) and versus mitochondrial mass (right) (MT-CO1/B2M by DNA qPCR). (B) PDX PGC-1α versus fl-E6 mRNA (left), p53 levels from Western blot densitometry (middle), and in vivo cisplatin response by rate-based T/C (right). (C) ERRα mRNA versus number of upregulated Hallmark_Oxidative_Phosphorylation transcripts in the PDXs (left) and the HPV⁺ OPSCCs in TCGA cohort (middle) and the JHU cohort (right). (D) Fl-E6 versus ERRα mRNA in TCGA cohort. (E) ERRα mRNA versus number of up-regulated Hallmark_P53_Pathway transcripts in TCGA cohort (left) and the JHU cohort (right).

The impact of variation in viral oncoprotein levels on clinical behavior of HPV⁺ cancers has been unclear to date. Furthermore, unknown molecular mechanisms distinguishing favorable versus unfavorable responses to therapy have hindered personalizing management for HPV⁺ OPSCCs. Our results reveal an ability of E6-mediated p53 downregulation to deplete mitochondrial antioxidant capacity and thus sensitize HPV⁺ OPSCCs to standard therapy, which relies upon induction of oxidative stress. We provide evidence that differing fl-E6 levels across HPV⁺ cancers variably repress the PGC-1α/ERRα pathway, leading to diversity in mitochondrial mass that likely affects treatment response. This finding provides a potentially novel mechanism by which levels of a viral oncoprotein may govern treatment outcomes in HPV⁺ OPSCC. Therapy resistance may prove prospectively identifiable by measuring fl-E6 levels and downstream components of the PGC-1α/ERRα axis and be mitigated through therapeutically targetable nodes in that pathway.

Because both E6 and E7 are viewed as essential to HPV-mediated immortalization (32), it was surprising to find wide variation in E6 levels among HPV⁺ OPSCCs. Specifically, unspliced fl-E6, the only E6 form that degrades its canonical target p53 (33), was downregulated in tumors with high mitochondrial mass, whereas E6* forms were not reduced (Figure 3B). Downregulation of fl-E6 by splicing might serve to uncouple its expression level from that of E7, which arises from the same polycistronic mRNA but is not spliced. Our data do not address whether the resulting E6* forms have functions that further contribute to the treatment-sensitive phenotype of the tumors with low fl-E6. This possibility is raised by the observation that HPV 18 E6* protein interferes with E6-directed p53 degradation (34) and thus could enhance the effect of splicing on reducing fl-E6 function. In addition, overexpressing E6* in cervical cancer cell lines was shown to increase ROS, DNA damage, and apoptosis (35). These effects of E6* may explain why expressing wild-type E6 similarly increases ROS and DNA damage (36). Of note, we did not detect significant increase in ROS upon stably expressing fl-E6 in SCC154 and VU147T cells in vitro (data not shown). Thus, increasing fl-E6 alone does not appear sufficient to upset the redox balance of untreated cells in culture despite reducing mitochondrial mass, and the E6* forms deriving from wild-type E6 expression might explain the previously observed increases in ROS (36).

Our findings add to understanding of how variable HPV oncogene expression affects oropharyngeal cancer clinical behavior. HPV E2 disruption by viral genome integration in a subset of cancers is traditionally thought to derepress E6/E7 expression, and there is some evidence suggesting that loss of E2 expression portends worse prognosis (37). Increased E1^E4 mRNA (38) and decreased E2F target gene dysregulation downstream of E7 (22) have also been linked to worse survival in small clinical cohorts, albeit without mechanistic basis. In addition, E2/E4/E5 expression was observed to predominate over E6/E7 in some OPSCCs (39), although this alternative viral gene expression pattern was not prognostic. Increased E6 splicing is observed in cancer relative to normal viral replication (40) and may indicate that the fl-E6 levels best serving malignant progression are lower than those that serve the viral life cycle. Similarly, chronic in vitro cisplatin exposure is reported to reduce fl-E6 expression (41) and may represent adaptation to mitigate oxidative stress. Taken together with our findings, these prior observations suggest that certain HPV⁺ cancers achieve competitive advantage through partially preserving p53 function by downregulating fl-E6.

The p53-mediated transcriptional upregulation of PGC-1α observed in other cancer types (15) agrees with our observed mechanism of E6 repressing mitochondrial biogenesis but may not fully account for E6’s effects in this context. Our observations notably diverge from findings in mouse embryonic fibroblasts, where p53 was activated by telomere shortening and directly bound and repressed PGC-1α, thus contributing to growth arrest, senescence, and apoptosis (42). By contrast, both wild-type and mutant p53 can act as positive regulatory cofactors for PGC-1α at promoters of its target genes in multiple cancer types (43, 44). Increased p53 also drives expression of sestrins (45), which mediate PGC-1α activation by stimulating its phosphorylation by AMPK (46). In addition, p53 supports mitochondrial DNA replication via interaction with mitochondrial transcription factor A (47). Beyond degrading p53, fl-E6 may regulate PGC-1α by inhibiting certain histone-modifying enzymes that also have nonhistone substrates. Prominent among them are the p300/CPB acetyltransferases, which acetylate p53 to enhance its stability (48) and interact with PGC-1α (49, 50) to remodel chromatin at its target promoters. PGC-1α is also activated by certain methyltransferases that are inhibited by fl-E6, including PRMT1 and set7/9 (51–54). E6 also targets TIP60, another acetyltransferase whose substrate BRD4 activates a super-enhancer of PGC-1α expression (55). These mechanisms highlight AMPK, PRMT1, set7/9, p300/CBP, TIP60, and BRD4 as targetable nodes upstream of PGC-1α/ERRα with potential utility for phenocopying E6’s desirable effects on therapy response.

The PGC-1α/ERRα complex and its myriad downstream effects on mitochondrial metabolism offer additional potential targets for treatment-resistant HPV⁺ OPSCCs. Pharmacologic ERRα inhibition was recently achieved using diaryl ether-based thiazolidinediones, which interfere with its DNA binding (56). Targeting p300/CBP’s co-activator function for PGC-1α (57) offers another approach to PGC-1α/ERRα axis repression and may bring an additional benefit of repressing DNA repair in certain HPV⁺ tumors with p300-activating mutations (58). Furthermore, inhibitors of electron transport chain complexes, TCA cycle enzymes, and glutamine metabolism in clinical trials for various treatment-refractory cancers (59) may have value for a subset of HPV⁺ OPSCCs of similar metabolic phenotype. Although aggressively targeting oxidative metabolism in the clinic remains hampered by toxicity, the severely impaired DNA damage repair phenotype in HPV⁺ OPSCCs (3) may widen the therapeutic window by allowing responses at lower drug doses than those needed for other cancers.

To date, the multimarker transcriptional profiles with prognostic signals for HPV⁺ OPSCC (21, 60) have not had a clear biologic basis or been generalizable enough for clinical application. Whereas E6 is unlikely to achieve utility as a single biomarker, mitochondrial mass is related to other features of HPV⁺ OPSCCs that hold promise as prognostic and predictive biomarkers. Image-based quantification of hypoxia as a source of radioresistance is in clinical trials to guide therapy de-escalation for HPV⁺ OPSCC (61). Because hypoxia induces oxidative stress, tumors with the antioxidant capacity to survive hypoxia may exploit the same antioxidant mechanisms to resist treatment. The benefit of increased mitochondrial function in this context is consistent with decreases in E6/E7 expression and p53 degradation (62) observed under hypoxia in HPV⁺ cervical cancer cell lines. There is also evidence that genetic silencing of negative regulators of NRF2, a driver of antioxidant gene transcription with downstream and upstream relationships to PGC-1α, confers poor prognosis in HPV⁺ cancers (12). Therefore, fl-E6 levels may merit joint consideration with NRF2 regulation, hypoxia, and other features linked to antioxidant capacity in developing multimarker predictors of therapy response.

Patient cohorts. Data for the Head and Neck Squamous Cell Carcinoma TCGA cohort were downloaded from the Genomic Data Commons (project TCGA-HNSCC) and Broad GDAC Firebrowse repository. A total of 53 HPV⁺ cases were selected from this broader TCGA cohort based on oropharyngeal primary site and previous detection of high-risk HPV type E6/E7 transcript expression (17). RNA-Seq and clinical information for the 47 HPV⁺ OPSCCs in a JHU cohort are previously described (20). RNA-Seq and clinical data were obtained for the 37 HPV⁺ OPSCCs treated with primary radiation plus chemotherapy within a published VU cohort (21).

Bioinformatic analysis of TCGA RNA-Seq data. RNA-Seq reads within FastQ files for each TCGA OPSCC sample (n = 53) were aligned to HPV genomes using Bowtie2 and samples with alignment selected and identified as containing HPV 16, HPV 33, HPV 35, or HPV 56 genomes. Gene expression values were represented by the normalized reads of reads per kilobase of transcript per million mapped reads. HPV genome sequences and annotations with gene coordinates were obtained from RefSeq (https://www.ncbi.nlm.nih.gov/refseq/). Reads with mapped coordinates corresponding to regions specific for each gene were counted. Fl-E6 expression was quantitated in a region nonoverlapping with E6*. The probe region specific to E6* overlaps fl-E6 and thus represents total E6. E6* is calculated by subtracting fl-E6 from total E6. HPV 16 gene expression analysis used GenBank K02718.1, with the following gene coordinates for quantitation: E6: 227-408, total E6: 83-226, E7: 562-858, E2: 2814-3331, E4: 3332-3619, and E5: 3863-4099. HPV 33 gene expression analysis used GenBank M12732.1, with the following gene coordinates for quantitation: E6: 232-413, total E6: 109-230, E7: 573-866, E2: 2814-3325, E4: 3326-3577, and E5: 3854-4081. HPV 35 gene expression analysis used GenBank X74477.1, with the following gene coordinates for quantitation: E6: 233-414, total E6: 110-231, E7: 562-861, E1: 868-2713, E2: 2882-3293, E4: 3294-3584, and E5: 3814-4065. HPV 56 gene expression analysis used GenBank: X74483.1, with the following gene coordinates for quantitation: E6: 234-415, total E6: 102-232, E7: 572-889, E2: 2806-3221, E4: 3222-3577, and E5: 3854-4081.

Growth and RNA-Seq of PDXs. The 7 treatment-naive HPV⁺ OPSCC PDXs used in this study were established and passaged as described by us in the subcutaneous flank of NSG mice (22). The PDXs were harvested for RNA-Seq in triplicate upon reaching 1 cm³. Volumes were calculated as width² × length/2. RNA isolation from 7 PDXs was performed with the Qiagen RNA Kit. Sample quality was assessed using Agilent RNA 6000 Nano Reagent on Agilent 2100 Bioanalyzer and quantified by Qubit RNA HS assay (Thermo Fisher Scientific). Ribosomal RNA depletion was performed with Ribo-zero rRNA Removal Kit, HMR (Illumina Inc.). Library preparation was performed with the NEBNext Ultra II Non-directional Synthesis Module (New England Biolabs). Samples were sequenced on Illumina NextSeq High Output sequencer and aligned to combined genome from hg38 (human), mm10 (mouse), and the high-risk HPV genomes — 16, 18, 31, 33, 35, 39, 45, 51, 52, 56, 58, 59, 66 and 68 — using STAR with default parameters. HPV sequences in all samples mapped to HPV 16 exclusively. After STAR first-pass alignment, reads mapped to mm10 were removed, and reads aligned to human genome (hg38) and HPV genomes were retained. Raw counts for each viral transcript (E2, E4, E5, L2, L1, total-E6, E7, and E1^E4) were obtained from reads mapped to each gene based on the alignment file. Fl-E6 was computed as the ratio of the average coverage level in the first intron (approximate full-length E6 level) divided by the average coverage level in the first exon (approximate all E6 transcript level). RnaSTAR was used to identify the first introns for E6*, which are, in the format of donor-acceptor, 227-408 for HPV 16. The average coverage of read pairs within the intron (i.e., 227-408 for HPV 16) was divided by the average coverage of read pairs within exon1 (i.e., 82-227 for HPV 16), where coverage was computed by dividing the number of reads by the feature width. RNA-Seq data reported here were submitted to National Center for Biotechnology Information Gene Expression Omnibus (GSE193388).

Cell lines and culture conditions. SCC090 (RRID:CVCL_1899), SCC152 (RRID:CVCL_0113), SCC154 (RRID:CVCL_2230), HEK293(RRID:CVCL_0045), HCT116 (RRID:CVCL_0291), and HCT116 p53^–/– (RRID:CVCL_HD97) cells are from ATCC. UM-SCC047 (RRID:CVCL_7759) and UMSCC104 (RRID:CVCL_7712) are from Sigma-Aldrich. UDSCC2 (RRID:CVCL_E325) and 93VU147T (RRID:CVCL_L895) cell lines were gifts from Silvio Gutkind (UCSD, San Diego, California, USA) and Hans Joenje (VU Medical Center, Amsterdam, the Netherlands), respectively. Cell lines were authenticated using the Identify Mapping Kit (Coriell). They were cultured in DMEM-F12 (Gibco), 10% heat-inactivated FBS, 1% penicillin-streptomycin, and 1% nonessential amino acids. N/tert-1/E7 cells were cultured in keratinocyte serum-free medium (Thermo Fisher Scientific).

qRT-PCR and Western blotting. Total RNA was extracted using the RNeasy Plus Micro Kit (QIAGEN). DNase I (AMPD1, Sigma-Aldrich) was used to remove genomic DNA from RNA samples. Reverse transcription was performed with oligo-dT plus random decamer primers using Superscript II (Thermo Fisher Scientific). QRT-PCR was performed with SYBR green master mix (Thermo Fisher Scientific) using Step-one Plus Real-Time PCR. Primers are in Supplemental Table 2. Western blots were performed with luminol-based chemiluminescent substrate (Bio-Rad), and bands in x-ray films were quantified with ImageJ software (RRID:SCR_003070; NIH). Antibodies are detailed in Supplemental Table 3.

Transient overexpression and silencing. Fl-E6 was expressed transiently using the pLentiN-16E6no* plasmid (RRID:Addgene_37445) and its pLentiN vector control (RRID:Addgene_37444). Transient PGC-1α overexpression was achieved with the pcDNA4 PGC-1α plasmid (RRID:Addgene_10974). A total of 4 μg of each plasmid was delivered per well in 6-well plates using Lipofectamine 2000 (Thermo Fisher Scientific) and incubated for 72 hours. SiRNA silencing of PGC-1α was performed using the mix of the 3 validated targeting siRNAs and scrambled negative control included in the PPARGC1A TriFECTa DsiRNA Kit (IDT) per supplier instructions. SiRNA was delivered at 50 nM concentration per well in 6-well plates using Lipofectamine RNAiMax (Thermo Fisher Scientific) and incubated for 72 hours.

Lentiviral transduction. The pLentiN-16E6no* plasmid containing full-length E6 and its pLentiN vector backbone control were used for stable fl-E6 expression. Lentiviral preparation and transduction were performed as described (25). Briefly, VU147T and SCC154 cells grown on 10 cm plates were incubated overnight with 1 mL of retroviral supernatant plus 4 mL of media containing 5 μg polybrene/mL (MilliporeSigma). Blasticidin selection at 2–5 μg/mL was performed for 3 to 4 weeks before bulk cultures were harvested and plated at limiting dilution in 96-well plates. Three empty vector and fl-E6-transduced single-cell clones were randomly chosen after 2–3 weeks of growth. N/tert-E7 cells were similarly infected with the control and E6 viruses before blasticidin selection, propagation, and analysis of bulk cultures.

Luciferase assays. HPV 16 E6 mutant constructs 16E6, 16E6 (V42L), 16E6-8S9A10T, 16E6-I128T, and 16E6* are previously described (30). The pGL3 luciferase reporter plasmid was used as the backbone for 2,423 bp human PGC-1α promoter. Plasmids were cotransfected using Lipofectamine 2000. Fold changes were determined by measuring the ratio of firefly to Renilla luciferase activity and comparing it to the ratio for control vector-transfected cells. Plasmid details are in Supplemental Table 4.

Evaluating mitochondrial mass and function. Mitochondrial DNA content was determined by qPCR as described (23). OCR was measured on the Seahorse Bioscience XFe96 Extracellular Flux Analyzer per manufacturer protocol using 0.001 mg/mL oligomycin, 2.5 μM carbonyl cyanide p-trifluoro-methoxyphenyl hydrazone, and 2 μM each of antimycin and rotenone (Sigma-Aldrich). Total protein per well was determined afterward with the Pierce BCA Protein Assay Kit (Thermo Fisher Scientific). Results were analyzed using XFe Wave software (Seahorse Bioscience), and 1 × 10⁵ cells were used to measure NADPH and NADP⁺ levels using an NADP⁺/NADPH colorimetric quantitation kit (Sigma-Aldrich).

Pharmacologic and radiation therapy experiments. Cell viability after 72-hour treatment with cisplatin, DMNQ, MnTMPyP, and/or MitoTEMPO (Sigma-Aldrich) was determined by WST-1 and/or BrdU assays (Sigma-Aldrich). Cisplatin dissolved in saline was administered intraperitoneally in vivo at 2 mg/kg. For irradiation, cells were harvested during exponential growth and plated at 2,000 cells/mL in 60 mm dishes containing 2 mL media. After 12 hours, they received 2, 4, or 6 Gy using a Gammacell 40 irradiator (Best Theratronics). Media were changed postirradiation. Colonies were stained with crystal violet (Thermo Fisher Scientific) and counted 10 days later.

Statistics. When comparing 2 groups for cell viability or gene expression, significance was calculated using 2-tailed unpaired t tests. Pearson correlation coefficients were based on analyzing the means. Two-way ANOVA was used when comparing xenograft growth curves. Log-rank tests were used for Kaplan-Meier analyses. P values were calculated for multiple comparisons using a 1-way ANOVA corrected with Dunnett’s or Holm-Šidák procedure. Tests used are indicated in figure legends. Analyses were performed using Prism (GraphPad Software). A P value less than 0.05 was considered significant.

Study approval. Mouse experiments were performed under Wistar Institute IACUC protocols 201166 and 201178 (Philadelphia, Pennsylvania, USA). PDXs were previously established under University of Pennsylvania IRB-approved protocol 417200 “Head and Neck Cancer Specimen Bank” (principal investigator: DB) by signing a combined informed consent and HIPAA form for use of tissue for research. All clinical data shown are obtained from deidentified, publicly available data sets.