Plasma Metabolic Phenotypes of HPV-Associated versus Smoking-Associated Head and Neck Cancer and Patient Survival

Head and neck squamous cell carcinoma (HNSCC) can be divided into two groups that are clinically, genetically, and possibly metabolically distinct: human papillomavirus (HPV)-associated HNSCC and smoking-associated HNSCC (sometimes viewed as HPV-unassociated). The former consists of HPV-positive oropharyngeal tumors with little to no patient smoking history, while the latter consists of HPV-negative oropharyngeal and all other HNSCC sites from patients with a substantial smoking history (i.e., >10 pack-years); HPV-positive oropharyngeal cancer patients with a substantial smoking history are likely a mix of HPV-associated and smoking-associated HNSCC. Clinically, HPV-associated HNSCC is less aggressive, responds better to treatment, and on average patients live two-to-seven times longer than smoking-associated HNSCC (1). Genetically, “hotspot” locations in the PIK3CA oncogene are mutated in HPV-associated tumors, while TP53 tumor suppressor loss mutations are in most smoking-associated tumors (2, 3). These genes are known to regulate tumor metabolism. For instance, wild-type p53 dampens glycolysis and promotes mitochondrial oxidative phosphorylation (OXPHOS) by regulating glucose transporters (GLUT) and key enzymes in the glycolysis pathway (4).

It has been hypothesized that loss or inactivation of p53 may promote the Warburg effect (5), that is, the reprogramming of tumor metabolism to favor glycolysis at the expense of OXPHOS (6). Noticed in the 1920s by Otto Warburg, this atypical metabolic state provides the energy and macromolecules for rapid proliferation and metastasis (7–10). It links to nearly all hallmarks of cancer (11). The Warburg effect can be seen, literally, via 18F-Fluorodeoxyglucose PET (18F-FDG-PET; ref. 12). Higher tumor uptake of the fluorescent glucose marker means higher glycolytic activity, which is directly associated with poor prognosis and treatment resistance in HNSCC (13, 14). In fact, the Warburg effect may be a factor in HNSCC treatment resistance. Increased glycolysis and decreased OXPHOS promotes resistance to ionizing radiotherapy (15) and studies suggest modulating glycolysis may resensitize a tumor (16, 17). Recent studies in HNSCC cell lines suggest that a Warburg-like metabolic phenotype differentiates HPV-positive from HPV-negative cells (18, 19). The Warburg effect may thus explain, in part, the differences in presentation, progression, and survival observed between these HNSCC groups. However, it has yet to be investigated as a potential blood biomarker.

Our goal was to investigate the metabolic differences of HPV-associated HNSCC relative to smoking-associated HNSCC in pretreatment patient plasma using untargeted high-resolution metabolomics (HRM) and test whether those differences were prognostic of overall survival independent of HPV status. HRM is the comprehensive measurement of small molecules arising from metabolic reactions occurring in a biospecimen—that is, metabolites (20). It uses gas or liquid chromatography paired with high-resolution mass spectrometry to quantify thousands of metabolites via their mass-to-charge ratio (m/z; ref. 21). The metabolites are matched to reference databases (22) or identified with authentic standards (23). While single metabolite investigations used to be the norm, recent advances in high-throughput measurement, bioinformatics, and pathway analysis have shifted the field toward a systems-biology approach, gaining further insight into disease via dysregulated metabolic pathways and profiles (24–26).

This was a prospective cohort of patients with HNSCC recruited at radiation oncology clinics affiliated with Emory University Hospital (Atlanta, GA) from 2013 to 2016 before undergoing intensity-modulated radiotherapy (IMRT) with or without chemotherapy (27). In this analysis, all patients gave written informed consent to participate in the parent study and all ancillary studies. Recruitment and data collection for the parent study and the current ancillary study were guided by the ethical principles of the Belmont Report and reviewed and approved by the Emory University Institutional Review Board. Inclusion criteria were histological HNSCC with no distant metastasis, ≥21 years of age, and no evidence of uncontrolled metabolic, hematologic, cardiovascular, renal, hepatic or neurologic disease. Patients with simultaneous primaries or major psychiatric disorders were excluded. Demographic and clinical variables were collected through chart review and standardized questionnaires at study entry. Prior to IMRT, blood was collected into chilled EDTA tubes for the immediate isolation of plasma. Plasma was separated by centrifugation at 1,000 × g for 10 minutes at 4°C then aliquoted into siliconized polypropylene tubes and stored at −80°C.

HPV status was ascertained via review of the tumor pathology report. We deemed a subject HPV positive if p16 was positive or if HPV DNA was detected. Smoking status was ascertained via questionnaire of tobacco habits; the questionnaire did not allow calculation of pack-years. We defined our smoking-associated HNSCC group as all current or former smokers with a non-oropharyngeal or HPV-negative oropharyngeal tumor; non-oropharyngeal tumors are rarely attributed to HPV (28) and thus not uniformly tested. We defined our HPV-associated HNSCC group as all subjects who were HPV-positive never smokers. We endeavored to mitigate group misclassification thus we deliberately excluded HPV-positive oropharyngeal cancer patients with a smoking history. We believe these patients are a mix of smoking-associated and HPV-associated HNSCC and by including them would introduce misclassification. In a 2010 seminal article, Ang and colleagues demonstrated that HPV-positive oropharyngeal cancer patients who are smokers with >10 pack-years—deemed as having “Intermediate risk” for survival—are clinically distinct from the “Low risk” group (i.e., HPV-positive, ≤10 pack-years) and the “High risk” group (i.e., HPV-negative, >10 pack-years; ref. 29).

Vital status was ascertained via linkage to the Georgia Comprehensive Cancer Registry maintained at Emory University (Atlanta, GA). Overall survival was defined as the time from study entry to death by any cause, or to the last contact up to August 2020.

We used LC/MS protocols developed at the Emory Clinical Biomarkers Laboratory (30, 31). Frozen aliquots from randomized samples were removed from −80°C, thawed, and extracted with ice-cold acetonitrile. Supernatants were added to a 4°C autosampler. Each sample was analyzed in triplicate (10 μL) using hydrophilic interaction liquid chromatography (HILIC) with positive electrospray ionization and Fourier transform Orbitrap high-resolution mass spectrometry (Dionex Ultimate 3000, Q-Exactive HF, Thermo Fisher Scientific). Analyte separation was achieved using a 2.1 mm × 100 mm × 2.6 μm Accucore HILIC column (Thermo Fisher Scientific) and a gradient elution of 2% formic acid, water, and acetonitrile starting at 10%, 10%, 80%, for 1.5 minutes with a linear increase to 10%, 80%, 10% at 6 minutes and held for 4 minutes per injection.

The high-resolution mass spectrometer operated at 120,000 resolution and an m/z range of 85.0000–1275.0000. Probe temperature, capillary temperature, sweep gas, and S-Lens RF levels were 200°C, 300°C, 1 arbitrary units (AU), and 45 AU, respectively. Positive tune for sheath gas, auxiliary gas, sweep gas, and spray voltage settings were 45 AU, 25 AU, and 3.5 kV, respectively. Spectral peaks were extracted and aligned using apLCMS (32) and xMSanalyzer (33). Mass, retention time, and ion intensity were recorded for unique m/z spectral features. Features with a coefficient of variation ≥75% were removed. Feature intensities were median summarized across triplicates and corrected for any batch effects (34). Features were log₂ transformed and quantile normalized. We restricted our investigation to m/z features annotated via the Kyoto Encyclopedia of Genes and Genomes (KEGG) and Human Metabolome (HMDB) reference databases (22).

We identified differences in metabolic pathways of HPV-associated HNSCC relative to smoking-associated HNSCC via feature selection (step 1) followed by metabolic pathway enrichment analysis (step 2). We reduced each significant pathway via principal component analysis (PCA; step 3) and estimated the first principal components (PC) with overall patient survival via Cox models (step 4). On the basis of the pathway results, we classified each patient, post hoc, into a binary glycolytic or OXPHOS metabolic phenotype and estimated that classification with overall survival (step 5).

Step 1: We selected m/z features that differentiated HPV-associated versus smoking-associated HNSCC in any of three models routinely used in untargeted metabolomics analyses: age at study entry, sex, and race adjusted logistic regression using a Benjamini–Hochberg FDR P < 0.10; Partial least squares discriminant analysis using a variable importance in projection >2; Support vector machine using recursive feature elimination. Step 2: We identified enriched metabolic pathways from the features using Mummichog, a pathway analysis tool that maps m/z to KEGG metabolic networks and uses resampling to obtain a Fisher exact test P value (24)—a similar methodology to gene set enrichment analysis. Pathways with ≥3 metabolite hits and P < 0.05 met our threshold for statistical significance. We note the identification level (35) of the pathway metabolites and plot a heatmap. Step 3: We performed PCA on the metabolites that comprised each significant pathway. Step 4: The first PC for each pathway was z-standardized and regressed on overall patient survival using Cox proportional hazards adjusted for HPV and smoking status, age at study entry, sex, race, tumor stage and patient Eastern Cooperative Oncology Group (ECOG) status; we found no statistically significant interaction by these variables. Further adjustment for body mass index (BMI), alcohol consumption, marriage status, treatment, feeding tube, and time of blood draw did not change our results and were not included in final models. The HRs estimated the increase in mortality per 1 SD change in the pathway-specific PC. We conducted a subanalysis restricting to oropharyngeal and non-oropharyngeal carcinoma subjects to investigate whether results hold across tumor site. Step 5: We dichotomized the PCs at their study medians and used the three most significant pathways to classify each subject as having a glycolytic or an OXPHOS phenotype (i.e., if ≥2 of the 3 pathways indicated glycolysis the patient was classified that way and vice versa with OXPHOS). We regressed the glycolytic versus OXPHOS phenotype on overall survival with a Cox model and plotted a Kaplan–Meier survival curve.

We used publicly available R packages for feature selection (xmsPANDA 1.1.76; ref. 36), pathway analysis (Mummichog 2.0; ref. 24), and metabolite annotation (xmsAnnotator 1.3.2; ref. 22). PCA used the “prcomp” R package. Cox regression used SAS version 9.4 (37).

Compared with smoking-associated, the HPV-associated HNSCC subjects on average were younger (56.8 vs. 60.1 years), more obese (29.0 vs. 26.1 kg/m²), White (87% vs. 72%), and male (89% vs. 62%; Table 1). As expected, the HPV-associated tumors were predominantly in the oropharynx (96%) and less advanced (29% were stage IV vs. 64% for smoking-associated). Also expected, the smoking-associated subjects had higher mortality than HPV-associated group [unadjusted HR = 7.6; 95% confidence interval (CI) = 2.8–25.5; P < 0.001; Kaplan–Meier curve in Supplementary Fig. S1].

From 16,342 extracted features, 3,075 matched to KEGG or HMDB reference databases. After filtering for missing values, 2,410 were left for analysis. The three feature selection models identified 191 metabolites that were differentially expressed between HPV-associated and smoking-associated HNSCC (Supplementary Table S1; Supplementary Fig. S2).

From 91 pathways, pathway enrichment analysis identified six statistically significant metabolic pathways from the 191 metabolites (Fig. 1; Supplementary Table S2). Notably, all six pathways are relevant to mitochondrial OXPHOS or glycolysis—demonstrating a potential Warburg-like metabolic phenotype. The pathways related to OXPHOS, that is, bile acid biosynthesis (P < 0.0001), 21-carbon steroid metabolism (P = 0.007), and octadecatrienoic acid beta-oxidation (P = 0.012), were upregulated in HPV-associated HNSCC and downregulated in smoking-associated HNSCC. The pathways related to glycolysis, that is, galactose metabolism (P = 0.001), vitamin B6 metabolism (P = 0.014), and starch and sucrose metabolism (P = 0.02), were the opposite. Four pathways implicated in HNSCC did not meet significance but followed the Warburg-like pattern: selenoamino acid metabolism (P = 0.065), phosphatidylinositol phosphate (P = 0.076), and hexose phosphorylation (P = 0.091) were upregulated in smoking-associated HNSCC, while omega-3 fatty acid metabolism (P = 0.065) was upregulated in HPV-associated HNSCC.

For subsequent analyses, we reduced the six pathways to four because of metabolite overlap (five of the six 21-carbon steroid metabolites overlapped with bile acid biosynthesis, and all three of the starch and sucrose metabolites overlapped with galactose). In Table 2, we annotate the 23 metabolites that comprised the bile acid biosynthesis (n = 9), galactose metabolism (n = 6), octadecatrienoic acid beta-oxidation (n = 5), and vitamin B6 metabolism (n = 3) pathways. Twelve features were unique (* = identified via an authentic standard): adenosine 5′-monophosphate*, dihydroxy-5beta-cholestane, cholic acid, galactosylglycerol, 3-ketolactose, lactose*, inosine diphosphate, raffinose*, stachyose*, pyridoxamine*, pyridoxic acid*, and pyridoxamine 5-phosphate. The other 11 had multiple matches. A heatmap shows the relative intensity of each metabolite and the pathways between HPV-associated and smoking-associated HNSCC (Fig. 2). In general, HPV-associated patients had higher levels of bile acid biosynthesis and octadecatrienoic acid beta-oxidation and lower levels of galactose metabolism and vitamin B6 metabolism; the smoking-associated patients with HNSCC had the opposite.

We reduced each pathway to its first PC and regressed that on overall mortality using adjusted Cox regression (PCA plots in Supplementary Fig. S3; standardized loading factors in Table 2). Three of the four pathways were significantly associated with overall mortality in unadjusted models (Table 3). The bile acid biosynthesis and the octadecatrienoic acid PCs were inversely associated with mortality (HR = 0.52 per SD; 95% CI, 0.38–0.72; P < 0.001; and HR = 0.54 per SD; 95% CI, 0.38–0.78; P < 0.001, respectively). The galactose PC was positively associated mortality (HR = 1.72 per SD; 95% CI, 1.09–2.72; P = 0.02); Kaplan–Meier curves in Supplementary Fig. S4). To clarify, a subject at the lowest 20th percentile in the bile acid PC would have one-third the mortality of a subject at the highest 20th percentile (HR = 0.52⁁1.68 SD = 0.33). Most importantly, the associations for bile acid and octadecatrienoic acid were independent of HPV status, smoking history, age at study entry, and sex (Table 3). In fact, bile acid biosynthesis remained statistically significant after further adjustment for race, tumor stage, and patient ECOG status (HR = 0.58 per SD; 95% CI, 0.37–0.92; P = 0.02). Restricting to tumor site had little effect on the bile acid biosynthesis pathway estimate but modest effects on the other pathways' estimates; fewer deaths among the oropharyngeal subjects may have precluded statistical significance (Table 3).

We used median dichotomized PCs from bile acid biosynthesis, galactose metabolism and octadecatrienoic acid beta-oxidation, to classify each subject with a glycolytic or mitochondrial OXPHOS phenotype (Fig. 3B). The Kaplan–Meier curve shows stark separation in survival (P < 0.0001) with the OXPHOS phenotype having a 93% estimated 3-year survival compared with 59% for the glycolytic phenotype (Fig. 3A). That survival difference was larger than what we observed for HPV-associated versus smoking-associated HNSCC, at 88% and 69%, respectively (Supplementary Fig. S1). The glycolytic versus OXPHOS Cox model showed the metabolic phenotype was significantly associated with mortality independent of HPV status and smoking history, age at study entry, sex, race, tumor stage, and patient ECOG status (HR = 3.17; 95% CI, 1.07–9.35; P = 0.04). Among HPV-associated subjects, 24% were classified as having the glycolytic phenotype (74% OXPHOS). Among the smoking-associated subjects, 34% of were classified as having the OXPHOS phenotype (66% glycolytic). This suggests that a glycolytic/OXPHOS phenotype may indeed be an independent prognostic biomarker and not merely an extension of HPV positivity or smoking status.

In this prospective study, we observed a glycolytic versus mitochondrial OXPHOS metabolomic phenotype in HNSCC patient plasma that was associated with patient survival independent of key clinical factors including HPV and smoking. Relative to smoking-associated HNSCC, the metabolic pathways of bile acid biosynthesis and octadecatrienoic acid beta-oxidation were upregulated in HPV-associated patient plasma, while galactose metabolism and vitamin B6 metabolism were downregulated. These pathways help regulate the balance between OXPHOS and glycolysis—a balance frequently dysregulated in tumors, termed the Warburg effect.

Four prior studies have investigated serum or plasma metabolomics in patients with HNSCC (38–41), but to our knowledge, this is the first study to investigate metabolic differences of HPV-associated versus smoking-associated HNSCC. Perhaps the best comparison with our study is one by Jung and colleagues that studied HPV-positive versus HPV-negative HNSCC cell lines (18). Jung reported that HPV-positive cells rely on mitochondrial OXPHOS for energy metabolism while HPV-negative cells rely on glycolysis—results nearly identical to ours. Jung also found HPV-negative cells had increased levels of hypoxia-inducible factor-1α (HIF1α) and hexokinase II (HK-II), regulating enzymes of the hexose phosphorylation pathway, the first step in glycolysis. We found the same pathway upregulated in smoking-associated HNSCC though the result was not statistically significant. Importantly, Jung found that pharmacologic reversal of the glycolytic phenotype increased a cell's sensitivity to radiation therapy, suggesting a potential drug target. The similarity of our findings in patient blood to those in cell cultures suggests that circulating metabolites may be putative biomarkers of HNSCC tumor metabolism.

The Warburg-like plasma phenotype we observed has two parts. The first part was the dysregulated pathways related to glycolysis (i.e., galactose, starch, and sucrose, phosphatidylinositol phosphate, and hexokinase metabolism). Of our findings, glycolysis has the strongest prior evidence. It is a hallmark of HNSCC proliferation and survival, driven by oncogenic signaling (e.g., p53 loss and PI3KCA mutation), altered metabolic enzymes (e.g., GLUT-1, HK-II), and a hypoxic tumor environment (e.g., HIF1α; refs. 42, 43). Tumor hypoxia activates HIF1α, a transcription factor which increases a tumor's reliance on glycolysis while reducing mitochondrial OXPHOS (43). Using gene expression profiles, Keck and colleagues found hypoxia to be the discriminating feature of the basal HNSCC phenotype—a phenotype made up of HPV-negative smokers (44). Higher levels of glycolytic metabolites found in the plasma of smoking-associated patients with HNSCC may indeed be the downstream effect of tumor hypoxia. Hypoxia, and its glycolytic switch, is associated with aggressive disease (45), treatment failure (16), and poor survival (46) in HNSCC. Glycolytic metabolites in serum are significantly associated with HNSCC relapse (40), and in our study, upregulated glycolytic pathways in plasma were associated with a more aggressive tumor type and poor survival.

The second part of the Warburg-like phenotype we observed was the dysregulated pathways related to mitochondrial OXPHOS (i.e., bile acid biosynthesis, octadecatrienoic acid beta-oxidation, and omega-3 fatty acid metabolism). Bile acid biosynthesis was our most important finding. It was significantly upregulated in HPV-associated HNSCC and was positively associated with longer survival independent of HPV, smoking, and other clinical factors. It is understudied in HNSCC but deserves closer examination. Bile acids are synthesized in the liver from cholesterol to aid in the digestion of dietary fat. That process, however, is regulated via a host of nuclear receptor signaling pathways, notably FXR, TGR5, and PPAR, that are now recognized for suppressing cancer progression (47–50). Bile acids may also be related to tumor hypoxia as they have been shown to destabilize HIF1α in prostate, lung, and mammary cancer cell lines (51). Conversely, tumor hypoxia may represses CYP7A1, a key bile acid biosynthesis enzyme though many factors may affect this pathway (52).

Like bile acids, we found octadecatrienoic acid beta-oxidation (i.e., oxidation of 18-carbon lipids such as omega-3 and omega-6 fatty acids) upregulated in HPV-associated HNSCC and positively associated with longer survival independent of HPV and smoking. Pro-oncogenic pathways generally promote synthesis of fatty acids over oxidation because proliferating cells require lipids for biological membranes (53). Apart from lipoprotein(a), serum lipids appear to decrease as HNSCC progresses (54–59) and short-chain fatty acids may suppress HNSCC growth (60). Our findings agree that higher levels of lipid oxidation, in general, is a beneficial marker in HNSCC.

The Warburg effect is likely just one aspect of tumor metabolism linked to HNSCC progression. Molecular changes underlying HNSCC progression, such as p53 loss, EGFR mutation, and mutations in the signaling of PI3K-AKT, NOTCH1, and ALK1, parallel metabolic changes in the tumor (10, 61–64). Metabolomics may thus provide a rich source of HNSCC biomarkers particularly if metabolic reprogramming drives treatment resistance. Resistance to radiotherapy, cisplatin chemotherapy, and cetuximab immunotherapy may be directly related to reprogrammed glucose metabolism and mitochondrial-generated oxidative stress (15, 65, 66). Thus, metabolic biomarkers may give an early indication of HNSCC treatment efficacy and may be targets for novel metabolic therapies (10, 64).

The recruitment of the parent study limited our metabolome-wide analysis to HNSCC patient plasma and not tumor specimens. In addition, we did not have healthy controls as a comparison group. However, we analyzed plasma metabolic differences in two important clinical subtypes and assessed those differences with survival—something no prior study has done. The exclusion of HPV-positive patients with a smoking history gave us distinct HPV-associated and smoking-associated groups but may limit our generalizability. The extracted metabolites are likely influenced by recent diet and the subjects did not fast before blood draw. Information on diet was not available, but we did adjust the models for factors that influence diet such as smoking, alcohol, age, sex, BMI, and the time of the blood draw. Finally, untargeted metabolomics is a “shotgun” approach that measures the relative intensity of each mass-to-charge feature and we cannot be certain of the annotation, outside of laboratory standards (23). To improve confidence in our annotation, we restricted our analysis to metabolites matched to the KEGG and HMDB reference databases. However, a feature may match to multiple reference metabolites and we cannot distinguish between them.

In conclusion, we observed a Warburg-like metabolomic phenotype in HNSCC patient plasma that may aid patient prognosis. We expound on the metabolic pathways that differentiate HPV-associated from smoking-associated HNSCC, two groups who differ in their presentation, treatment sensitivity, and survival. We are working to independently validate the results we observed with bile acids, lipid oxidation, and galactose metabolism so that they warrant future investigations as predictors of treatment efficacy and as therapeutic targets for treatment resistant tumors.

R.C. Eldridge reports grants from Georgia Clinical and Translational Science Alliance during the conduct of the study. A.H. Miller reports non-financial support from Boehringer Ingelheim outside the submitted work. K.A. Higgins reports other support from AstraZeneca and grants and other support from RefleXion Medical during the conduct of the study. N.F. Saba reports receiving compensation for consult services to GSK, Merck, CUE, Kura, Vaccinex, and Pfizer. No disclosures were reported by the other authors.

R.C. Eldridge: Conceptualization, data curation, software, formal analysis, funding acquisition, methodology, writing–original draft, project administration, writing–review and editing. K. Uppal: Conceptualization, formal analysis, supervision, methodology, writing–review and editing. D.N. Hayes: Supervision, writing–review and editing. M.R. Smith: Writing–review and editing. X. Hu: Writing–review and editing. Z.S. Qin: Methodology, writing–review and editing. J.J. Beitler: Resources, writing–review and editing. A.H. Miller: Resources, writing–review and editing. E.C. Wommack: Data curation, project administration, writing–review and editing. K.A. Higgins: Resources, writing–review and editing. D.M. Shin: Resources, writing–review and editing. B. Ulrich: Data curation, writing–review and editing. D.C. Qian: Data curation, writing–review and editing. N.F. Saba: Resources, supervision, writing–review and editing. D.W. Bruner: Conceptualization, supervision, writing–review and editing. D.P. Jones: Conceptualization, resources, methodology, writing–review and editing. C. Xiao: Conceptualization, data curation, supervision, funding acquisition, project administration, writing–review and editing.

R.C. Eldridge would like to thank the Georgia Clinical and Translational Science Alliance (UL1TR002378 and KL2TR002381), and C. Xiao would like to thank the National Institute of Nursing Research (K99/R00NR014587 and R01NR015783) and the Winship Cancer Institute via the NCI (P30CA138292) for funding the conducted research. The content is solely the responsibility of the authors and does not necessarily represent the official views of the Georgia CTSA or the NIH.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.