Integrative Analysis of Head and Neck Cancer Identifies Two Biologically Distinct HPV and Three Non-HPV Subtypes

Head and neck squamous cell carcinoma (HNSCC) is the sixth most common nonskin cancer worldwide with an annual incidence of approximately 600,000 cases and a mortality rate of 40% to 50% (1, 2). The major known risk factors are environmental exposures to tobacco products, alcohol, and infection with high-risk human papillomaviruses (HPV). The incidence of HPV⁺ tumors is rising rapidly and HPV status is now the strongest prognostic marker (3).

Although HNSCC is now widely viewed as composed of two distinct clinical entities, HPV⁺ and HPV⁻, a therapeutically relevant molecular classification system remains elusive. Prior attempts at profiling this clinically heterogeneous disease have been hampered by lack of information about HPV status and/or absence of significant numbers of HPV⁺ cases, small sample size, and lack of annotation (4, 5). Currently, neither molecular classification nor validated biomarkers are used in clinical practice. Instead all patients with HNSCC are treated independent of the underlying biology based on stage and anatomic location, typically using a combination of surgery, radiation, and chemotherapy (6). Cetuximab, an anti-EGFR antibody, is the only approved targeted therapy for HNSCC, but predictive biomarkers remain to be identified and overall cetuximab has a response rate of 13% (7).

In the current study, we investigate a fully clinically annotated, very homogenous patient cohort of locoregionally advanced HNSCC (including 44% HPV⁺ tumors), all treated uniformly with organ-preserving chemoradiotherapy together with all available other HNC cohorts, including the cancer genome atlas head and neck cancer cohort (TCGA), to identify HNSCC subtypes and their biologic and prognostic characteristics. We integrate gene expression-based consensus clustering, cross-platform approaches for validation, copy number profiling, and HPV status to discover subtype-specific, translationally relevant biology.

The resulting new taxonomy of HNSCC closely correlates with multiple HNC-relevant oncogenic processes and a biology-based taxonomy of HNSCC. In particular, we identify two subtypes of HPV⁺ head and neck cancer and two subtypes (one HPV⁺ and one HPV⁻) with a prominent immune phenotype characterized by CD8⁺ T-cell infiltration, as well as a subtype with near exclusive HER ligand/EGFR copy number-driven signaling. As such, our classification is broadly applicable to biomarker and therapy development for head and neck cancer and will inform future clinical and preclinical studies, biomarker development, and personalized HNC care.

OCT blocks of frozen tissue samples were obtained from the University of Chicago Head and Neck Cancer tissue bank (Chicago, IL; Institutional Review Board-approved protocol UCCCC#8980; Supplementary Experimental Procedures and Supplementary Fig. S1).

A total of 171 locoregionally advanced HNSCC specimens obtained before treatment with concurrent chemoradiotherapy (all patients received organ-preserving chemoradiotherapy with curative intent) were selected for this study. Our aim was to select at least 100 cases to be able to perform consensus clustering later. Samples were collected and banked between 1997 and 2010, but our focus was on recent cases. A section was cut from OCT-frozen blocks and stained with hematoxylin and eosin (H&E). An expert HNC pathologist reviewed slides to determine and circle the area with the highest tumor content. Samples with microscopic tumor or tumor content <60% were excluded from further analysis. Guided by the H&E-stained slides, the region with the highest tumor content was cut from the OCT blocks, pulverized using CryoPrep (Covaris) and homogenized in lysis buffer from an AllPrepRNA/DNA/Protein Mini kit (Qiagen) or from an RNA/DNA/Protein Purification Kit (Norgen Biotek) using an Ultrasonicator (Covaris). DNA, RNA, and protein were isolated from each sample using the respective kit and following the manufacturer's protocol.

Gene expression profiling was done using Agilent 4 × 44 Kv2 expression arrays on 134 samples (Supplementary Fig. S1 and Supplementary Experimental Procedures). To evaluate copy number aberration (CNA) in HNSCC, we selected a panel of 86 and subsequently 147 candidate genes known to have frequent CNA in cancer and analyzed these using the Nanostring nCounter assay (NanoString Technologies; ref. 8; Supplementary Experimental Procedures). Data and materials availability: Gene Expression Omnibus Accession number: GSE40774 (Embargoed until publication; http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?token=jjsnbqwecyoyirg&acc=GSE40774).

To discover gene expression-based subtypes, we used newly generated gene expression data (Agilent, n = 134), as well as a data based on Illumina BeadArrays (n = 131; Yordy and colleagues; in preparation), and publicly available HNSCC microarray data from two studies based on Affymetrix HG133 arrays (n = 106; refs. 9, 10; see Supplementary Experimental Procedures for details).

IHC and immunofluorescence were performed using routine methods (Supplementary Experimental Procedures).

Microarray gene expression profiles of 371 HNSCC samples were derived from several datasets (Supplementary Table S1), including our cohort (n = 134, Agilent; Supplementary Fig. S1), an Illumina cohort (Yordy and colleagues; unpublished results), and two Affymetrix cohorts (9, 11). These cross-platform datasets were used as the discovery set for the classification of HNSCC samples. All datasets were analyzed as described in the Supplementary Experimental Procedures and outlined in Supplementary Fig. S2A. This resulted in 821 reliable and informative genes (Supplementary Table S2), based on which five subtypes were derived separately in each dataset (Supplementary Fig. S2B). These five subtypes slightly differ across the three datasets, but fall into three very stable cross-platform groups (termed supergroups) on the basis of hierarchical clustering of their centroids (Fig. 1A, B, and C). All centroids in the three cross-platform supergroups show a strongly positive silhouette width (average = 0.47), indicating that centroids fit their respective cluster (Fig. 1B; ref. 12). The high correlation of our three supergroups and previously identified subtypes of lung SCC (13) and HNSCC (14) confirms the validity of our classification (Fig. 1C). We applied the nearest centroid approach, a standard classification method, to assign new samples from the validation datasets (14–16) to the supergroup with closest matching centroid. Supergroups were validated in two independent validation datasets from different platforms, including Affymetrix and RNA-Seq, suggesting broad applicability (Fig. 1E). According to prior nomenclature in other cancer types and molecular characteristics (12, 17, 18), we named the three supergroups inflamed/mesenchymal (IMS), basal (BA), and classical (CL), respectively.

HNSCC tumors are composed of two distinct disease entities: HPV⁻ and HPV⁺ tumors. We investigated the relationship between HPV status and the three supergroups [basal (BA), classical (CL), and inflamed/mesenchymal (IMS)]. First, we determined the HPV status on our cohort using HPV DNA-based PCR and RNA-based qPCR (see Supplementary Experimental Procedure); for the TCGA dataset, we used exome-sequencing data and RNA-Seq data to determine the DNA and RNA of HPV viruses (see Supplementary Experimental Procedure). The HPV status of other public datasets was determined using our unpublished HPV gene expression signature. In all the datasets, HPV⁺ tumors were not gathered into one group, but fall into two distinct groups: (i) the inflamed/mesenchymal (IMS) supergroup and (ii) the classical (CL) supergroup. No HPV⁺ tumors fell in the basal supergroup (Fig. 2A). It is evident that HPV⁺ HNSCCs are composed of two distinct gene expression subtypes, namely, IMS-HPV and CL-HPV. Interestingly, there is biologic overlap of these HPV⁺ samples in a supergroup and HPV⁻ tumors in the same supergroup (e.g., immune or mesenchymal differentiation). The non-HPV subtypes were named as basal, CL-nonHPV, and IMS-nonHPV.

One distinctive feature of the basal subtype is the significant enrichment for hypoxia signaling, represented by hypoxia-responsive genes such as HIF1A, CA9, and VEGF (Fig. 2B and Supplementary Fig. S3A). In addition, the gene expression profile of the basal subtype closely matches a published hypoxia gene signature (Fig. 2C) which validates the hypoxic phenotype (19). The second key characteristic of the basal subtype is strong enrichment for neuregulin signaling, including EGFR, AREG (amphiregulin), and NRG1 (neuregulin 1/heregulin; Fig. 2B and Supplementary Fig. S3A). Overexpression of epithelial markers such as P-cadherin (CDH3) and cytokeratins (KRT1, KRT9) is another distinctive characteristic of the basal subtype (Fig. 2B and Supplementary Fig. S3A) similar to basal breast cancer (20). Moreover, a published EMT signature also demonstrates the elevated epithelial pathway activities in the basal subtype (Fig. 2C; ref. 21). Consistent with the high expression of cytokeratins in the basal subtype, the morphology data also demonstrated that basal tumors are highly keratinizing and well differentiated (Supplementary Fig. S4).

The most distinctive feature of the classical supergroup is the significant enrichment for putrescine (polyamine) degradation pathway (Supplementary Fig. S3B), which is relevant for detoxification, for example, related to tobacco use. Increased polyamine levels are associated with increased cell proliferation (22). Consistently, a published proliferation signature also indicates that the classical supergroup has a higher proliferation rate compared with the other groups (Supplementary Fig. S4; ref. 23).

Although the CL-HPV and CL-nonHPV subtypes share similarities, they are still two distinct disease entities, reflected in many biologic pathways. Cell-cycle genes, such as mini-chromosome maintenance proteins (MCM2 and MCM10), cell division cycle protein kinase (CDC7), and -related genes (CDKN2A, E2F2, and RPA2) are overexpressed in the CL-HPV subtype. The two subtypes also show significant difference in tobacco use with 74% heavy smokers in CL-nonHPV compared with 42% heavy smokers in CL-HPV (Fig. 2B and Table 1, P = 0.01, Fisher exact test). Consistent with smoking status, xenobiotic metabolism pathway genes AKR1C1, AKR1C3, and ALDH3A1, which are known to be associated with smoking, are enriched in CL-nonHPV (Fig. 2B; refs. 24, 25).

The distinguishing features of the IMS group is expression of immune response genes like CD8, ICOS, LAG3, and HLA-DRA (Fig. 2B) related to the infiltration of CD8⁺ T lymphocytes in tumors.

Mesenchymal genes such as vimentin (VIM), matrix metalloproteinases (MMP9), and S100A4 also show increased expression in the IMS group (Fig. 2B), which typically associates with increased metastatic risk (8, 26). Epithelial markers such as P-cadherin (CDH3) and cytokeratins (KRT1, KRT9) are downregulated, suggesting epithelial-to-mesenchymal transition (EMT; Fig. 2B; ref. 21). Consistently, a published EMT signature (21) is upregulated in the IMS group and there is downregulation of epithelial differentiation markers (Fig. 2C).

The two subtypes in the IMS group show a significant difference in cell-cycle pathways and smoking-associated pathways (Fig. 2B). Similar to the CL-HPV subtype, the IMS-HPV subtype has significantly higher cell-cycle pathway activities, in which HPV is known to play a critical role. The IMS-HPV subtype shows a higher proliferation rate indicated by a published proliferation signature (Fig. 2C; ref. 23) and IMS-HPV subtype tumors are nonkeratinizing and poorly differentiated according to the morphology review by light microscopy (Supplementary Fig. S3).

Our cohort of 134 HNSCC samples with 44% HPV⁺ is representative for the patient population presenting to a large referral center for advanced HNC. There is an increasing proportion of HPV⁺ cases (Table 1). All patients in this new cohort had locoregionally advanced disease and were treated with concurrent chemoradiotherapy. Consistent with previous reports, patients with oropharyngeal tumors are mainly in the two HPV subtypes: IMS-HPV and CL-HPV (71%, P = 4.19 × 10⁻¹¹). In contrast, oral cavity tumors are overrepresented in the basal group (72%, P = 1.04 × 10⁻⁵). Seventy-four percent of the CL-nonHPV tumors and 85% of IMS-nonHPV are heavy smokers, which is significantly higher than other subtypes.

Five-year survival rate was assessed for our Agilent, TCGA, and Illumina cohorts, respectively. In all the cohorts, the HPV subtypes have a significantly higher 5-year survival than non-HPV subtypes. Furthermore, between the two HPV subtypes, the IMS-HPV subtype shows a trend toward higher 5-year survival than CL-HPV (Fig. 3A) consistent across all cohorts. We used Kaplan–Meier estimator to measure the overall survival (OS for the subtypes in the above three cohorts, separately. The log-rank test was used to evaluate the Kaplan–Meier curve difference among the subtypes. We found that OS differs significantly among the five subtypes in our Agilent cohort (P = 0.037) and the TCGA cohort (P = 0.003) but not in the Illumina cohort (P = 0.589; Fig. 3B). Furthermore, the Kaplan–Meier analysis on the combined dataset of the three cohorts showed a clearer difference in OS of the five subtypes (P = 8e-06), which was more significant than any cohort independently (Fig. 3B). Importantly, the IMS-HPV subtype is different from the CL-HPV subtype by better survival (Fig 3A and B) suggesting a potential positive impact of immune response.

Copy number analysis was executed in all samples. Initially, 101 samples were analyzed targeting 75 genes previously described to have CNA (27–29). A GISTIC-derived strategy was used to identify significant CNA. Significant copy number gains were detected in 12 genes, including PIK3CA (3q26.3), VEGFA (6p21), EGFR (7p11), MYC (8q24), and CCND1 (11q13; Fig. 4A and B and Supplementary Table S3). Significant copy number losses were detected in 16 genes, including CDKN2A (9p21) and RB1 (13q14; Fig. 4B and Supplementary Tables S3). Subsequently, we aimed to further investigate these aberrant regions by scanning 72 additional genes within copy number altered regions in 55 of the samples. We identified that TP63, SOX2 (3q26-28), and PIK3CA (3q26.3; Supplementary Table S4) are coamplified in most tumors.

Furthermore, CNAs for a lot of oncogenes and tumor suppressors are found to have strong associations with subtypes (Supplementary Table S5). Amplification of EGFR (7p12), CCND1 (11q13), and FADD (11q13), as well as deletion of FHIT (3p14) and CDKN2A (9p2) with few exceptions are only observed in the nonHPV subtypes, particularly in the basal subtype (Supplementary Table S5). Amplification of MYC (8q24) and ITGB4 (17q25) is significantly enriched in the basal subtype (Supplementary Table S5). ITGB4 is also overexpressed in the basal subtype, and it is reported that ITGB4 mediates cell adhesion and plays a crucial role on the initiation, progression, and metastasis of solid tumors (16, 30).

The copy number gains of genes (PRKCI, PIK3CA, and DCUN1D1) in the 3q26-28 region can be seen in all subtypes, but are significantly higher in the classical subtypes (Supplementary Table S5). Although the frequency of E2F3 (6p22) amplification in our cohort is only 5%, most of the amplification events occur in the classical subtypes, including both CL-HPV and CL-nonHPV subtypes (Supplementary Table S5). Consistently, E2F3 is overexpressed in classical subtypes. E2F3 encodes a transcription factor important for cell-cycle regulation and DNA replication, and its amplification and overexpression are known to associate with invasive tumor growth and rapid tumor cell proliferation in urinary bladder cancer (31). Importantly, the two HPV subtypes differ significantly in these regions, suggesting different biology of two HPV subtypes. The amplification of 3q26-28 was validated by different markers such as TP63 and SOX2 (Fig. 4C).

Using IHC, we validated the ability of key markers (SOX2, Cyclin D1/CCND1, p16/CDKN2A, and TP63) to differentiate subtypes (Supplementary Fig. S5).

One of the key findings based on gene expression is the discovery of immune-related marker expression (CD8A/B) in IMS tumors (Fig. 5A). We performed multicolor immunofluorescence to identify tumor-infiltrating CD8⁺ lymphocytes in these tumors (Fig. 5B). The enrichment of cytotoxic T-cell infiltration in mesenchymal tumors is present even when considering oropharynx tumors only, suggesting independence of anatomic location and lack of contamination from normal lymphoid tissue (Fig. 5C).

Although previous attempts at expression profiling of HNSCC identified distinct tumor subtypes (11, 14, 19), a therapeutically relevant molecular classification remains elusive due to lack of information about HPV status and/or absence of significant numbers of HPV⁺ cases, small sample size, and lack of predictive biomarkers. The HNSCC subtypes proposed in our study are identified in an unsupervised, cross-cohort and cross-platform way, which is supported by four discovery cohorts and four independent validation cohorts, which together total 938 patients. We showed that our subtypes are remarkably similar to those found in LUSC (17) and also have a strong correlation with the previous identified HNSCC subtypes (14). Moreover, our expression profiling includes a large and representative number of HPV⁺ tumors, allowing us to identify two biologically distinct HPV⁺ subtypes.

Consequently, we now propose that HNSCC can be classified into five distinct subtypes—two HPV and three non-HPV subtypes. The close correlation of the five HNSCC subtypes with morphologic characteristics, molecular processes, survival, and copy number changes supports a biologic and clinical basis for this classification, making it unlikely that subtypes are attributable to chance, artifact, or bias.

The most important finding of this report is the identification of two distinct HPV subtypes. Previous clinical observations suggest that HPV⁺ tumors are diverse, for example, the subgroup of patients with HPV⁺ tumors do not respond well to therapy (3, 32). Our classification into two distinct HPV subtypes provides a biologic basis for clinical heterogeneity and suggests that differential treatment approaches might be required for HPV⁺ tumors. CL-HPV and IMS-HPV subtypes exhibit significant difference in many aspects such as morphology, molecular processes, CNAs, and clinical features. More than 40% of the CL-HPV tumors show keratinization, whereas none of the IMS-HPV tumors are keratinizing. Accordingly, IMS-HPV subtype tumors are more poorly differentiated compared with CL-HPV subtype tumors. The IMS-HPV subtype exhibits significantly elevated expression of mesenchymal markers, whereas the CL-HPV subtype shows higher proliferation. Compared with the IMS-HPV subtype, the CL-HPV subtype presents significantly more canonical genomic aberrations associated with squamous cell carcinoma, such as amplification of 3q26-27. For example, E2F3 amplification, which is associated with cell proliferation, is only found in the CL-HPV subtype but not in the IMS-HPV subtype, which correlates with the proliferation signature seen by expression in the CL-HPV subtype. The differences in the biologic patterns found between the two HPV subtypes are supported by a trend toward better survival in the IMS-HPV subtype compared with the CL-HPV subtype, and this observation holds true across all cohorts, although it does not correlate as strongly with prognosis as, for example, tobacco use.

Another key finding of this report is the recognition of the differential expression of immune markers in the subtypes. Cancer immune surveillance is considered to be a factor in the body's ability to prevent cancer (33). Evading immune surveillance has been considered as an emerging hallmark of cancer (34). Importantly, activity of checkpoint blockade was recently reported in HNC using anti–PD-1 therapies, and PD-L1 expression correlates with an inflamed phenotype that is consistent with and validates our inflamed/mesenchymal (IMS) subtypes (35, 36).

Furthermore, accumulating evidence indicates that cancer immune suppression can be driven by hypoxia (37, 38). We detected that the basal subtype is highly enriched for HIF1A signaling. Immunosuppressive factors such as VEGF (39) are upregulated in the basal subtype, and a lack of immune-related markers is evident. Alternative to hypoxia-induced activation, HIF1A can be activated under normoxic conditions by EGFR signaling (40, 41), which is also found to be significantly activated in the basal subtype tumors. EGFR amplification and MYC amplification are found in the basal subtypes. Taken together, the immunosuppression of the basal subtype could be driven by either hypoxia or potentially also EGFR signaling, and further validation will be important as immune-targeted therapies are being developed for HNC.

Currently, no clinically relevant EGFR biomarkers exist and anti-EGFR therapy is administered indiscriminately despite a low single-agent response rate of 13% (7, 42). Our results suggest that the basal subtype could serve as a candidate predictive biomarker in evaluating the benefit of an anti-EGFR therapy, as well as a hypoxia-targeting therapy.

In contrast with the absence of immune markers of the basal subtype, the IMS subtypes, including both the IMS-HPV and the IMS-nonHPV subtypes, exhibit a strongly activated immune phenotype, validated by CD8⁺ T-cell infiltration on IHC. Infiltration of tumors with CD8⁺ cytotoxic T lymphocytes has been associated with a favorable prognosis in several tumor types and may be a predictive biomarker for cancer immunotherapy (43). Furthermore, novel immunomodulatory therapies (e.g., PD-1 checkpoint targeting drugs) have shown activity in squamous cell carcinomas (44) and evidence of activity in HNSCC was recently reported. We demonstrate that both IMS subtypes independent of HPV status are characterized by CD8⁺ T-cell infiltration (Fig. 5A and B) irrespective of anatomic location (Fig. 5C). Immune response against HPV⁺ is of particular interest and evidence of immune escape was reported recently (45, 46). This is the first report that HPV⁺ tumors have two subtypes that exhibit a different immune phenotype. CD8⁺ T-cell infiltration is therefore one possible explanation for the trend toward a more favorable prognosis of IMS-HPV subtype tumors, and consistent with findings in other cancer types. Although further investigations are needed to explain the difference in immune response between the two HPV subtypes, one can hypothesize that differences in the host, tumor, and viral genetic background, as well as their interplay, are etiologic.

Given the widespread use of IHC for diagnostic purposes in routine pathology practice (e.g., p16, p63), we evaluated the ability of protein expression to differentiate subtypes (Supplementary Fig. S5). Although additional optimization and validation are necessary, development of either an IHC or a qPCR-based assay could potentially also allow classification of HNSCC similar to what is used for breast cancer (47).

In conclusion, we propose a new taxonomy of HNSCC based on five subtypes with a comprehensive overview of HPV⁺ as well as HPV⁻ HNSCC biology validated across 938 HNSCC patient tumors. There are significant translational implications with respect to immunotherapy, presence of two biologically distinct HPV subtypes, and HER-targeted therapies that will inform biomarker development and personalized care efforts that are already being pursued.

T.Y. Seiwert is a consultant/advisory board member for Merck/MSD and Bayer/Onyx. No potential conflicts of interest were disclosed by the other authors.

Conception and design: M.K. Keck, Z. Zuo, T.P. Stricker, M. El Dinali, R. Salgia, E.E.W. Cohen, K.P. White, E.E. Vokes, T.Y. Seiwert

Development of methodology: M.K. Keck, Z. Zuo, A. Khattri, T.P. Stricker, M. El Dinali, S. Aktolga, Z. Lei, S.G. Rozen, T.Y. Seiwert

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): M.K. Keck, T.P. Stricker, M. Imanguli, P. Fang, R.R. Weichselbaum, M.W. Lingen, M.D. Story, K.K. Ang, E.E.W. Cohen, K.P. White, E.E. Vokes, T.Y. Seiwert

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): M.K. Keck, Z. Zuo, A. Khattri, T.P. Stricker, C. Brown, D. Rieke, J. Bragelmann, R. DeBoer, M. El Dinali, S. Aktolga, S.G. Rozen, R. Salgia, R.R. Weichselbaum, M.W. Lingen, E.E.W. Cohen, T.Y. Seiwert

Writing, review, and/or revision of the manuscript: M.K. Keck, Z. Zuo, D. Rieke, J. Bragelmann, R. DeBoer, S. Aktolga, P. Tan, R. Salgia, R.R. Weichselbaum, M.W. Lingen, M.D. Story, K.K. Ang, E.E.W. Cohen, K.P. White, E.E. Vokes, T.Y. Seiwert

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): M.K. Keck, D. Rieke, K. Endhardt, R. DeBoer, M. El Dinali, S. Aktolga, T.Y. Seiwert

Study supervision: E.E. Vokes, T.Y. Seiwert

Other (funding): T.Y. Seiwert

The authors thank Drs. Janet Rowley and Yusuke Nakamura for their critical feedback, suggestions, and guidance, the Gleason Family, and the Achatz Research Fund for their support of HNC research (to T.Y. Seiwert).

This work was supported by an ASCO Translational Professorship (to E.E. Vokes). Dr. T.Y. Seiwert was supported by The Flight Attendant Medical Research Institute (FAMRI) Young Investigator Award.

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.