Clinical Validity of a Prognostic Gene Expression Cluster-Based Model in Human Papillomavirus–Positive Oropharyngeal Carcinoma

PURPOSE Under common therapeutic regimens, the prognosis of human papillomavirus (HPV)–positive squamous oropharyngeal carcinomas (OPCs) is more favorable than HPV-negative OPCs. However, the prognosis of some tumors is dismal, and validated prognostic factors are missing in clinical practice. The present work aimed to validate the prognostic significance of our published three-cluster model and to compare its prognostic value with those of the 8th edition of the tumor-node-metastasis staging system (TNM8) and published signatures and clustering models. METHODS Patients with HPV DNA-positive OPCs with locoregionally advanced nonmetastatic disease treated with curative intent (BD2Decide observational study, NCT02832102) were considered as validation cohort. Patients were treated in seven European centers, with expertise in the multidisciplinary management of patients with head and neck cancer. The median follow-up was 46.2 months (95% CI, 41.2 to 50), and data collection was concluded in September 2019. The primary end point of this study was overall survival (OS). Three-clustering models and seven prognostic signatures were compared with our three-cluster model. RESULTS The study population consisted of 235 patients. The three-cluster model confirmed its prognostic value. Two-year OS in each cluster was 100% in the low-risk cluster, 96.6% in the intermediate-risk cluster, and 86.3% in the high-risk cluster (P = .00074). For the high-risk cluster, we observed an area under the curve = 0.832 for 2-year OS, significantly outperforming TNM 8th edition (area under the curve = 0.596), and functional and biological differences were identified for each cluster. CONCLUSION The rigorous clinical selection of the cases included in this study confirmed the robustness of our three-cluster model in HPV-positive OPCs. The prognostic value was found to be independent and superior compared with TNM8. The next step includes the translation of the three-cluster model in clinical practice. This could open the way to future exploration of already available therapies in HPV-positive OPCs tailoring de-escalation or intensification according to the three-cluster model.


INTRODUCTION
Among head and neck squamous cell carcinomas (HNSCCs), the incidence of oropharyngeal cancer (OPC) has been increasing in the past few decades. This is mainly due to the increasing number of human papillomavirus (HPV)-related OPCs. 1 The first findings of HPV infection in HNSCC have been reported by European research groups in the 1980s, 2 followed by a first causal association by expression analysis and a link to the OPC in the 1990s. 3 The association between HPV-positive OPC and the positive prognostic impact has been described approximately 10 years later. 4 The favorable prognosis of patients with HPV-positive OPC led to explore treatment de-escalation in recent clinical trials. However, expectations were in contrast with the results because this approach was associated with poorer outcomes. Thus, standard of care has not been changed. 5,6 Therefore, there is still a strong need for trustable prognostic factors. 7,8 Different preclinical and clinical studies have been performed leading to the identification of several genomic, epigenomic, transcriptomic, and proteomic characteristics of HPVpositive OPCs. [9][10][11][12][13][14] At present, the most relevant clinical prognostic factors are stage, extranodal extension (ENE), and smoking exposure. 15 ENE is currently part of the AJCC and UICC staging system for HPV-negative HNSCCs, but not for HPV-positive ones. In particular, in p16-negative OPC and non-OPC HNSCCs, clinical or radiologic ENE is considered for cN staging, and histopathologic ENE for pN one. Nonetheless, its negative prognostic impact was found in HPV-positive OPC as well, so some authors suggested its ASSOCIATED CONTENT

Data Supplement
Author affiliations and support information (if applicable) appear at the end of this article.
Accepted on September 20, 2021 and published at ascopubs.org/journal/ po on October 27, 2021: DOI https://doi. org/10.1200/PO. 21. 00094 implementation in refining the cN of the 8 th edition of the tumor-node-metastasis staging system (TNM8). 16 In the field of gene expression (GE) profiling, the pioneering studies of Slebos 9 and Pyeon 17 (referred to the initial characterization of potential biomarkers able to differentiate HPV-positive from HPV-negative HNSCC) and numerous omics characteristics of HPV-positive OPCs were identified in preclinical and clinical studies. [9][10][11][12][13][14] Additionally, the analysis of HNSCC transcriptomic data highlighted the existence of HPV-positive distinctive prognostic clusters. In particular, we published a meta-analysis of publicly available transcriptomic data sets of HPV-positive HNSCC, in which three different clusters were characterized. Our three-cluster model was commented and compared with previous HPV-positive HNSCC subtypes in a recent comprehensive review. 18 Moreover, Dhawan 19 performed a meta-analysis of OPC data sets with available survival data to determine whether currently available GE signatures could prognosticate outcome more accurately than established clinicopathologic predictors.
Exploiting a well-characterized homogenous HPV-positive OPC cohort with locoregionally advanced disease treated with curative intent, 20 the present study aims to validate the prognostic significance of our previously identified threecluster model and to compare its prognostic value with those of TNM8 and the most relevant already published signatures and clustering models.

Patient Selection
Patients with HNSCC with locoregionally advanced disease (stages III and IV according to TNM 7th edition 21 ) and treated with curative intent had been included in the project Big Data and Models for Personalized Head and Neck Cancer Decision Support (BD2Decide). 20 The study was approved by institutional review boards and ethical committee in March 2016, and patients consented to enrollment. The study population included 1,537 patients diagnosed and treated between 2008 and 2017, and the follow-up closed in September 2019. In this work, data and analyses are focused on patients with HPV-positive OPC. To reduce confounding variables affecting oncologic outcome and consequentially prognostic accuracy of the proposed model, we excluded (1) patients with p16-positive but HPV DNA-negative OPC 22 and (2) patients receiving unimodal treatment (surgery without postoperative radiation or exclusive radiotherapy without concomitant systemic treatments) for clinical stage T1N . 1, T2N . 1, and T3-T4 any N disease.

Knowledge Generated
The stratification of the proposed three-cluster model confirmed its prognostic value, significantly outperforming the 8 th edition of the tumor-node-metastasis staging system and other available transcriptomic signatures. The three clusters also differed for their biological peculiarities. The biological characterization highlighted the differential expression of key components of the tumor microenvironment as potentially targetable pathways. Relevance This study confirmed the presence of three different prognostic and biological clusters for HPV-positive oropharyngeal squamous cell carcinoma. Moreover, it supported the potential of introducing a gene expression model in HPV-related oropharyngeal carcinoma to provide a valuable biological basis for the design of clinical studies in the context of treatment tailoring. generated using Affymetrix human Clariom D microarrays (Affymetrix, Santa Clara, CA). 20 MIAME-compliant data are deposited in the GEO repository (GSE163173).

Bioinformatics Analysis
After microarray data normalization, the HPV-positive OPC selected cohort was classified according to the threecluster model. 24 Tumor microenvironment, the immune and stroma scores, and other cell components were inferred by xCell. 25 Pairwise gene set enrichment analysis (GSEA) 26 comparisons were performed to discern which Hallmark gene sets 27 were differentially enriched. In the final data set, we assessed the overall survival (OS) prognostic discrimination ability of our three-cluster model by area under the curve (AUC) at 2 years, and its performance was challenged against TNM8 and publicly available prognostic models and signatures: Regarding the other models, a metagene approach was applied: A score was computed in our data set by ssGSEA using the GSVA R package 33 providing an unsupervised method for estimating variation of gene set enrichment through the samples. AUC at 2 years was assessed through timeROC R package. 34 To enable a better understanding of strengths and limits of the available signatures, studies and their characteristics are summarized in the Data Supplement.

Statistical Analysis
Differences between the survival curves of the three clusters were assessed using the log-rank test and R package survival. Primary end point was OS, and secondary end point was disease-free survival (DFS). For survival analyses, we set 2 years of observation. Further details are provided in the Data Supplement.

Definition of the Final Data Set
In the BD2Decide project, 624 patients with stage III-IV (TNM 7th edition) OPC were included, and 377 patients were p16-positive. For validation analysis, we analyzed the cohort composed of 286 p16-positive cases with informative transcriptomics. Among them, 235 p16-positive and HPV-positive cases (14 excluded because of HPV DNAnegative and 37 excluded because of single-modality treatments) were suitable as a final data set for clinical correlations (Fig 1). No significant differences were recorded between the cohorts (Data Supplement). In the final data set, the median follow-up was 46.2 months (95% CI, 41.2 to 50), 79% of patients were men (185), and the median age was 60 years (Table 1). Of the 225 cases with regional lymph node involvement, 18% had a radiologic and/or clinical ENE, defined as the radiologic and/or clinical evidence of cancer invasion through nodal capsule or beyond. 35 Table 1).

Prognostic Analysis
The prognostic stratification according to the three-cluster model was tested in both the validation (Data Supplement) and final data sets (Fig 2A) and confirmed the rank order of the three-cluster model obtained in the original metaanalysis (Data Supplement). 24  Survival and recurrence probability were predicted with statistical significance by TNM8 (Fig 2B), not by TNM7 (Data Supplement). Through an alluvial plot, the distribution of patients according to the TNM7 and TNM8 classification and the three clusters was explored (Fig 2C; Data Supplement).

Tumor Biological Landscape
Three signatures reported molecular patterns associated with HPV status; 9,17,37 these signatures and their performance in stratifying HPV tumors were investigated and compared with our three-cluster model (Fig 3). The stratification ability of pioneering signatures 9,17 was not confirmed because of the limited sample size (Fig 3A and  Fig 3B). On the other hand, a significant difference was observed with Zhang clustering (Fig 3C), and the distribution in the three clusters confirmed the similarity, as already highlighted in the review by Qin et al. 18 The samples of the final data set, subdivided according to the threecluster model, were analyzed for difference in cell composition by xCell bioinformatics tool (Data Supplement), essentially confirming the results obtained in the original study. 24 The highest level of immune score was recorded in low risk (P = 2.8 × 10 −6 ; Fig 4A). When the immune populations were explored in detail, we observed a significant enrichment of total or memory CD4 T cells, effector memory CD8 T cells, and B cells in low risk, followed by  Fig 4B); the detailed analysis of the microenvironment revealed an abundance in fibroblasts, keratinocytes, and endothelial cells in high risk, while epithelial cells characterized both high and intermediate risks (Data Supplement). To gain a deeper insight into the biological and functional differences among the three clusters, analysis of the 50 hallmark gene sets in GSEA was performed in the final data set (Fig 4C). According to GSEA,

DISCUSSION
This is an independent validation of a three-cluster model, previously developed through meta-analysis of publicly available data sets. 24 The main difference between the original and present data sets relies on the current study   cohort. Our original meta-analysis 24 had a limitation since the cases were HNSCC HPV-positive, but 32% of them were non-OPC. The present analysis considered only HPVpositive locoregionally advanced OPC patients treated with curative intent. All the HPV-positive not OPCs together with HPV DNA-negative or suboptimally treated patients were excluded from the final data set. This homogeneous patient series and the availability of controlled clinical and pathologic data led to the refinement of the three-cluster model on the largest data set of OPC presently available. Unsupervised clustering analysis is a popular approach to discover novel biological-oriented disease subtypes. This approach has deeply been applied in breast cancer starting from the early 2000s; the signature was refined and its last version, PAM50, added prognostic value to the traditional pathologic, histologic, and biological parameters. 38 Clustering models of HNSCC associated with HPV have been already developed essentially defining two clusters associated with HPV-positive tumors. However, these studies 28,37,39 contained non-OPC HPV-positive HNSCCs and HPV-negative OPCs. In the original meta-analysis, 24 we observed that most of the HPV-positive non-OPC cases belonged to the high-risk cluster. It is noteworthy that the majority of patients (47%, 24 of 51) excluded from the validation set, because of HPV DNA-negative or suboptimally treated, were part of the high-risk cluster. Thus, the results obtained on OPC final data set indicated that, after removal of potential confounding factors, high risk maintained the worst prognosis. The present work definitely confirmed the prognostic value of the three-cluster model in HPV-positive OPCs. The prognostic accuracy of the three-cluster model (AUC = 0.832 for 2-year OS) outperformed the TNM8, the current gold standard for prognostic stratification in oncology, and at bivariate   analysis, its prognostic role was independent. In addition to the three-cluster model, three prognostic signatures, even if developed in all sites HNSCC or multiple cancers or OPC HPV-negative cases, were also able to outperform TNM8. They significantly separate the OPC HPV-positive cases with an AUC . 0.76. Interestingly, these signatures were obtained by meta-analysis of more than a 1,000 samples and/or were validated in different published studies. On the contrary, the other three signatures and models, built on smaller or heterogeneous HNSCC cohorts, did not show a prognostic role, suggesting the importance of larger or homogeneous study design in the context of the omics studies. The prognostic performance of the three-cluster model was not influenced by underlying clinical factors. The only difference in terms of clinical characteristics among the three clusters was ENE prevalence. Indeed, high-risk patients had a higher frequency of ENE plus disease, reflecting a more aggressive disease. 35 This consistent association corroborates the robust independent prognostic ability of the proposed three-cluster model. Low risk was characterized by the best prognosis and included 46% of patients. These tumors expressed a high immune score and significant enrichment of immune active cells (total or memory CD4 T cells, effector memory CD8 T cells, and B cells). They were characterized by a low RSI score and the lowest hypoxia scores. Our low-risk cluster is probably similar to the HPV-IMU identified by Zhang et al, 37 also confirmed by the experimental comparison of Qin et al. 18 The worst prognosis was observed in high-risk patients (13% of the study cohort). This group, Molecular Subtyping in HPV-Positive Oropharyngeal Carcinoma characterized by the highest keratinocytic differentiation status, is similar to the HPV-KRT cluster of Zhang et al 37 and was confirmed in the comparison of Qin et al. 18 Moreover, the high-risk cluster was characterized by the lowest immune score and an enrichment in GSEA pathways associated with the immune-exhausted HNSCC tumors, 38 which is also defined as cold immune. 18 The enrichment in hypoxia hallmark, consistent with the immune suppression and radioresistance (high RSI), was confirmed by the highest expression of hypoxia signatures.
The novelty introduced by our work was the identification of the intermediate-risk cluster (41% of the study cohort), never characterized before by other clustering models. 14-16 Despite our efforts, this was the less biologically and functionally characterized cluster. Nevertheless, biologically we were able to separate this group from the patients with good (low risk) or poor (high risk) prognosis. In detail, it expressed an intermediate immune score not associated with exhausted immune cells, a marginal radioresistance and low expression of De Cecco hypoxia signature. Notably, we observed an evident enrichment in metabolic pathways in both high and intermediate risks. Further investigation and deeper comments about possible treatments applicable to patients with HPV-positive OPC are reported in the review of Qin et al 18 on the basis of the experimental comparison of the three-cluster model 24 data with the other published HPV-positive twocluster models. 14,17,28,37,39 In addition, we observed that in our study population, the majority of current or previous smokers had a relevant tobacco exposure (. 10 pack-years in 71% of cases). Nevertheless, the frequency of current or previous smokers in our study population (66%) was lower than the one reported in previously published studies (82% in the RTOG0129 trial 36 ). This is consistent with the observed smoking reduction in Western countries, where cigarette smoking rates have been reducing in the past 50 years.
Moreover, future projections are forecasting a further reduction in smoking prevalence. 40 These epidemiologic variations imply that the observations found in our study population might need further updates in the next decades, to study whether behavioral changes will be reflected by different biological alterations, with possible clinical impact. Potential limitations of the three-cluster model are the observed dropout rate (25%) in our GE case material compared with the original clinical set and the lack of specification of the viral genotype for each case. The first limitation could be due to the fact that for most OPC cases, where chemoradiation is usually preferred over surgery, the analyzed specimen was an initial small biopsy. However, independent studies show that an increase in sample size, corresponding to our dropout rate, had limited effect on the clustering performance. 37 The second drawback could be attributed to the HPV DNA tests used. However, this limitation is mitigated by the fact that the vast majority of HPV-positive OPCs are known to be related to the HPV16 genotype. 41 Indeed, we do not anticipate major negative impacts of the cited drawbacks on our analysis.
In conclusion, in patients with HPV-positive OPC with locoregionally advanced disease treated with curative intent, there is still a strong need for defining trustable prognostic factors, and our results seem to fulfill this necessity. The rigorous clinical selection of the cases included in the BD2Decide study corroborates the robustness of the proposed three-cluster model, which resulted in an independent and superior prognostic factor compared with TNM8, the current gold standard for prognostic forecasting in clinical practice. The next step includes the required translation of this GE clustering in clinical practice, following a framework applied in other tumors. This could open the way to future exploration of already available therapies in HPV-positive OPCs tailoring de-escalation or intensification according to the three-cluster model.

AUTHORS' DISCLOSURES OF POTENTIAL CONFLICTS OF INTEREST
The following represents disclosure information provided by authors of this manuscript. All relationships are considered compensated unless otherwise noted. Relationships are self-held unless noted. I =Immediate Family Member, Inst = My Institution. Relationships may not relate to the subject matter of this manuscript. For more information about ASCO's conflict of interest policy, please refer to www.asco.org/rwc or ascopubs.org/po/author-center. Open Payments is a public database containing information reported by companies about payments made to US-licensed physicians (Open Payments).

Frank Hoebers
Consulting or Advisory Role: Bristol Myers Squibb