Mesenchymal gene expression subtyping analysis for early-stage human papillomavirus-negative head and neck squamous cell carcinoma reveals prognostic and predictive applications

Patients with oral cavity squamous cell carcinoma (OCSCC) are predominantly human papillomavirus (HPV)(−), and treatment typically involves surgical resection ± neck dissection, followed by radiation ± chemotherapy. We previously described four mRNA expression patterns (classical, atypical, basal, and mesenchymal), each with unique genomic features and prognosis. Here, we examine the clinical utility of gene expression subtyping in head and neck squamous cell carcinoma (HNSCC) and introduce potentially predictive applications in HPV(−) OCSCC. A retrospective genomic database analysis was performed including 562 HNSCC patients from MD Anderson (MDA-GSE41116) and The Cancer Genome Atlas (TCGA). Samples were assigned molecular subtypes (classical, atypical, basal, and mesenchymal) using an 88-gene classifier. HPV status was determined by gene expression. The clinical endpoint was overall survival censured at 36 months. The Kaplan–Meier plots and log-rank tests were used to investigate associations between clinical variables and survival. Of the 418 TCGA training patients who met analysis criteria, nearly 20% presented as stage I/II. Among node(−) OCSCC patients, the mesenchymal subtype is associated with worse survival (hazard ratio (HR) = 2.4, p = 0.021), offering a potentially actionable biomarker in otherwise early-stage, low-risk disease. This was confirmed in the MDA validation cohort. Node(−) non-mesenchymal OCSCC patients had far better survival compared to node(−) mesenchymal, and all node(+) patients had similarly poor survival. These findings suggest that the mesenchymal subtype is associated with poor survival in surgically resected, early-stage, node(−) OCSCC otherwise expected to have favorable outcomes. These findings highlight the potential value of gene expression subtyping as a pathology adjunct for prognostication and treatment decision-making in OCSCC patients.

Patients with oral cavity squamous cell carcinoma (OCSCC) are predominantly human papillomavirus (HPV)(−), and treatment typically involves surgical resection ± neck dissection, followed by radiation ± chemotherapy. We previously described four mRNA expression patterns (classical, atypical, basal, and mesenchymal), each with unique genomic features and prognosis. Here, we examine the clinical utility of gene expression subtyping in head and neck squamous cell carcinoma (HNSCC) and introduce potentially predictive applications in HPV(−) OCSCC. A retrospective genomic database analysis was performed including 562 HNSCC patients from MD Anderson (MDA-GSE41116) and The Cancer Genome Atlas (TCGA). Samples were assigned molecular subtypes (classical, atypical, basal, and mesenchymal) using an 88gene classifier. HPV status was determined by gene expression. The clinical endpoint was overall survival censured at 36 months. The Kaplan-Meier plots and log-rank tests were used to investigate associations between clinical variables and survival. Of the 418 TCGA training patients who met analysis criteria, nearly 20% presented as stage I/II. Among node(−) OCSCC patients, the mesenchymal subtype is associated with worse survival (hazard ratio (HR) = 2.4, p = 0.021), offering a potentially actionable biomarker in otherwise early-stage, low-risk disease. This was confirmed in the MDA validation cohort. Node(−) non-mesenchymal OCSCC patients had far better survival compared to node (−) mesenchymal, and all node(+) patients had similarly poor survival. These

Introduction
Head and neck squamous cell carcinoma (HNSCC), including cancers of the oral cavity, oropharynx, nasopharynx, hypopharynx, and larynx, is one of the most common cancers worldwide (1). In the United States, it is estimated that there were approximately 66,000 new cases and 14,00 deaths in 2021 (1). The majority of HNSCC cases are associated with heavy tobacco and alcohol use, although over the last 30 years, there has been an increase in the incidence of human papillomavirus (HPV)-related cancers, primarily in the oropharynx. While treatment of HNSCC depends on multiple tumor and patientrelated factors, the three main modalities used in HNSCC management are surgical resection, radiation therapy, and chemotherapy. Patients with early-stage tumors are generally treated with a single modality therapy, while those with advanced-stage tumors often require multiple modalities. Oncologic outcomes in HNSCC are driven largely by stage at presentation: the 5-year overall survival for stage I-II and III-IV HNSCC cases is approximately 70%-90% and 40%-60%, respectively.
While the majority of early-stage HNSCC cases are curable with surgery or radiation, it is notable that 10%-30% of HPVnegative HNSCC cases without pathologically aggressive features still experience a relapse event (2). Oral cavity squamous cell carcinoma (OCSCC or OC) is the most common head and neck cancer, comprising 1/3 of all cases, with the vast majority of cases presenting as HPV-negative and associated with tobacco use. Dependent on clinical staging, OC treatment involves surgical excision of the primary tumor with or without neck dissection, followed by radiation with or without chemotherapy.
There have been significant advances in our understanding of HNSCC molecular biology and genomic tumor heterogeneity. Based on earlier work in lung cancer (3), our group and others described four mRNA expression patterns (classical, atypical, basal, and mesenchymal) demonstrating unique genomic features and prognostic significance (4,5). These subtypes show varied biology and may be helpful in prognostic assessments complementing risk stratification based on HPV status, stage, anatomic site, and other characteristics (4,5). The basal subtype is characterized by over-expression of genes functioning in cell adhesion including COL17A1, and growth factor and receptor TGFA and EGFR (5). The mesenchymal subtype displays over-expression of genes involved in immune response (6,7) and is characterized by the expression of genes associated with epithelial-to-mesenchymal transition (EMT) including VIM, DES, TWIST1, and HGF (5). It has been suggested previously that EMT pathways are important in the initiation of nodal metastasis (8). The classical subtype is characterized by over-expression of genes related to oxidative stress response and xenobiotic metabolism and is most strongly associated with tobacco exposure (9)(10)(11)(12). The atypical subtype, which includes both HPV and non-HPV tumors, is characterized by elevated expression of CDKN2A, LIG1, and RPA2 and is associated with low EGFR expression (5). These four gene-expression-based head and neck cancer subtypes have been validated in other studies, including The Cancer Genome Atlas (TCGA) head and neck cancer cohort (4)(5)(6).
In this study, we examine the potential clinical utility of gene expression subtyping in HNSCC, with an emphasis on evaluating this biomarker among early-stage HPV-negative cancers. Our findings provide further evidence to support the clinical utility of gene expression subtyping in HNSCC within the context of clinical site, stage, and treatment as well as to introduce the potential for predictive applications of gene expression subtyping analysis in HPV-negative HNSCC.

Patients and datasets
The study was conducted in accordance with the Declaration of Helsinki and the International Conference on Harmonization Good Clinical Practice guidelines and was approved by the Institutional Review Boards of Washington University in St. Louis (IRB#201706088) and the University of Tennessee Health Science Center (IRB# 17-05549-XP).
Two datasets were identified in the public record: 1) TCGA HNSCC (n = 520) (6) and 2) large institutional cohort (n = 42) (13). For statistical analyses, cases were considered if they had clinical parameters of N stage and overall survival. Model fits and analyses used all available patients having complete relevant data. Patient statistics and demographics are shown in Table 1. TCGA data were sourced from the Broad Institute Genomic Data Commons (GDC) (14), and the institutional cohort was obtained from the gene expression omnibus (GEO) (GSE41116) (15).

Gene expression analysis
For TCGA, the upper quantile normalized RNA-seq expression values by expected maximization (RSEM) (16) were downloaded from GDC (gdac.broadinstitute.org/, accessed 12/4/ 2015) and log2-transformed. All samples were assigned a molecular subtype using the nearest centroid classification method as previously reported by Dabney (17), describing each sample as belonging to one of four molecular subtypes (basal, mesenchymal, atypical, or classical). Briefly, the HNSCC centroid predictor is a set of vectors, each comprised of typical gene expression values for one of the subtypes, and uses a set of 838 feature genes selected to distinguish the four molecular subtypes (5). By calculating the distance (1 − Pearson correlation coefficient) between each sample and each centroid, the algorithm determines the subtype to which the sample is most similar based on the predictor gene set. Each sample is then uniquely assigned to the subtype to which the distance was the shortest. For the purpose of developing a more parsimonious and potentially clinically relevant predictor, a reduced, 88-gene centroid predictor was developed (Supplementary Table 1). Fivefold cross-validation using all 520 samples and the ClaNC software package was used to identify the number of genes required for strong separation of the subtypes and sufficient agreement with the previously developed gold standard. Candidates for the reduced set were all genes in the gold standard classifier, and an additional set of genes (348) was chosen for high observed mean and variance in the entire data set. Here the standard ClaNC approach was modified by requiring an equal proportion of high and low genes per subtype in the final model rather than selecting genes based on extreme absolute values of the ClaNC t-statistic. Calculation of the coefficients in the final nearest centroid classifier excluded samples with low gold standard classifier call strength (20% per subtype were excluded), where call strength was the commonly used silhouette score, and the coefficients themselves were within-subtype gene medians after each gene had been centered by its overall median. Heat maps displaying expression profile patterns by subtype calls were generated using the Complexheatmap package in R.

Clinical data and statistical analysis
Paired clinic data were obtained from GDC (gdac..broad. instituteorg/accessed 12/4/2015) (14). To account for limitations in median follow-up, and for comparison to prior work, all survival times longer than 36 months were truncated and censored at 36 months. In general, clinical parameters were represented as presented in downloaded clinical datasets. HPV positivity was assessed by RNA-seq evaluation of HPV-aligned sequences in HPV types 16,18,33, and 35 at levels >1,000 counts. HPV reference sequence data were based on the PaVe website: https://pave.niaid.nih.gov/. Read counts >1,000 for HPV RNA-seq (TCGA) or HPV E6 gene expression (18) were used as the criterion for ongoing HPV replication and an HPVpositive tumor designation. Other parameters of interest included gender, age, smoking, T stage, N stage, and overall stage. Associations between two categorical variables were evaluated using Fisher's exact test and the chi-square test. Associations between categorical and continuous variables were evaluated using the Kruskal-Wallis test. The Kaplan-Meier plots and the log-rank test were used to assess univariate associations between survival and study parameters. Cox models were used to check associations with adjustment for potential confounders. The R survival package was used for all statistical analyses.

Clinical and molecular groups
We first considered the clinical characteristics of 418 patients from TCGA dataset meeting our eligibility criteria to understand the generalizability of our results to the greater population of head and neck cancers. The median age of TCGA HNSCC cohort was 60 years, only slightly younger than that reported in the American population for this disease, which is 63 years (1). Twenty-five percent of patients were female compared to 27% of patients in the American HNSCC population. Seventy-eight percent of patients in the dataset admitted some degree of smoking, which is consistent with prior reports (19). Eighty-one percent of patients presented with at least stage III disease, consistent with most datasets studied by molecular profiling. Consistent with the head and neck cancer disease course, only one patient was known to have metastatic disease at presentation, although this data field was missing for more than half of the patients. Larger, more advanced tumors tend to yield more sufficient tissue for molecular profiling. Nonetheless, nearly 20% of patients presented as stage I and II, offering a unique opportunity to assess risk profiles in earlystage patients. Somewhat unexpectedly, only 63% of patients had a record of radiation in the dataset, which seems low considering    that most of the 81% of advanced stage patients might be expected to receive radiation as part of multidisciplinary treatment. Considering subgroups, we noted that 71% of nodepositive oral cavity patients reported radiation versus 85% of the node-positive non-oral cavity. Forty-three percent of nodenegative oral cavity patients were radiated, compared to 50% of the node-negative non-oral cavity. Whether the low percentage represents underutilization of the standard of care based on patient-specific factors or under-reporting of radiation in the database is not known. That said, the trends were as expected in which the highest rates of radiation were found in node-positive non-OC patients, for whom nearly all patients would have a recommendation for radiation based on National Comprehensive Cancer Network (NCCN) guidelines, either as part of concurrent chemoradiation or surgery followed by radiation. The lowest rates were among node-negative oral cavity patients, many of whom could be treated with single modality surgery, according to NCCN guidelines. We then interrogated the distribution of molecular subtypes as a function of the anatomic site, where subtypes were determined by applying the centroid subtype predictor to all samples. We document that in the 418 TCGA samples meeting eligibility criteria, the distribution of molecular subtypes was nearly identical to that of the original TCGA HNSCC report of 279 cases published in 2015: 30% basal, 26% mesenchymal, 18% classical, and 26% atypical versus 31% basal, 27% mesenchymal, 16% classical, and 26% atypical (6) ( Figure 1A and Table 1). Consistent with other reports, there is a strong association of subtype with the anatomic site. We observed enrichment in oral cavity tumors among the mesenchymal and basal subtypes, atypical samples primarily in the oropharynx and to a lesser extent the larynx, and classical subtype in the larynx. In an unexpected finding, we found that although lymph node involvement is observed in all molecular and anatomic sites as expected, we observed a statistically significant association with lymph node positivity in the mesenchymal tumors (p = 0.03), where the proportion of positive nodes in mesenchymal was 67% compared to 55% in non-mesenchymal subtypes. Overall, nodenegative OC patients demonstrated significantly better survival than node positive OC patients (Figures 2A), however we also found a significant association between Mesenchymal-and nonmesenchymal subtype and overall survival in OC, which our group and others had previously reported in smaller cohorts as a statistically significant association for overall survival and mesenchymal subtype (4-6, 13). The finding was again observed in this cohort ( Table 2) (Figures 2B, 3). Since nodal status is a component of the overall cancer stage, itself defining patient prognosis, we considered that lymph node involvement might either be confounding for the worse prognosis for mesenchymal tumors or, more interestingly, be in the causal pathway of poor prognosis.

Oral cavity cohort
For the purposes of defining a cohort in which questions of clinical management and prognosis might be more specifically considered, we isolated patients with oral cavity squamous cancers ( Figure 1B), a group generally treated by a more explicitly clinical pathway. In general, patients with OC are treated primarily with surgery in all cases where a tumor is expected to be resected with negative margins. Early-stage A B FIGURE 1 Gene expression heat maps including 838 gene classifier genes as described previously (5)    a higher clinical stage and the overall worse prognosis observed for mesenchymal patients in multiple datasets. Interestingly, although the mesenchymal patients were overall of higher nodal status, the association with the T stage was less clear. Among node-positive OC patients, mesenchymal and nonmesenchymal patients had nearly identical nodal stage distribution of N1, N2, and N3. By contrast, among nodenegative patients, 69% (18 of 26) of mesenchymal patients were T1 or T2, whereas only 43% (40 of 93) of nonmesenchymal patients were T1 or T2. Only 19% (5 of 26) of node-negative mesenchymal patients were T4 compared to 38% (35 of 93) of non-mesenchymal patients. Of T1-T2 mesenchymal tumors, 49% (17 of 35) were node-positive, compared to 44% (31 of 71) of T1-T2 non-mesenchymal patients. By contrast, of mesenchymal T3-T4 tumors, 82% (36 of 44) were node-positive compared to 55% (66 of 119) of nonmesenchymal T3-T4 tumors. Additionally, in mesenchymal patients, T3-T4 tumors were associated with node positivity (p = 0.003), whereas the T stage was not associated with node status in non-mesenchymal patients (p = 0.13). Summarizing the associations for OC patients, we conclude that mesenchymal patients are both more likely to develop nodal metastasis, and they are more likely to do this at an earlier T stage. At the higher T stage, mesenchymal patients were much more likely to be node-positive.
Given the association between three clinical prognostic factors (anatomic site, and T and N stage) with a validated molecular marker (mesenchymal subtype), we considered both stratified and multivariate models of prognosis to investigate potential subgroups as well as the independent contribution of each factor (Table 2). We demonstrate, as expected, that clinical outcomes differ as a function of the anatomic site, T stage, and N stage. We then considered substrata as a function of nodal status, observing a striking finding in which those patients who were node-negative mesenchymal subtype demonstrated identical risk of death to patients who were node-positive mesenchymal or node-positive non-mesenchymal (hazard ratio (HR) = 2.4, p = 0.021). In other words, the mesenchymal molecular subtype conveyed all the risk of positive nodes whether nodes were clinically present or not. Among node-positive patients, the added risk of the mesenchymal subtype was no longer observed (HR = 1.3, p = 0.3). Given the unexpected nature of this finding, we sought independent datasets of oral cavity cancer for the purposes of validation, noting perhaps the largest being a set of well-characterized tumors from MD Anderson. Quite strikingly, the results were nearly identical, with patients of the non-mesenchymal OC group showing overall excellent survival and node-negative mesenchymal patients, node-positive mesenchymal patients, and node-positive non-mesenchymal patients all with similarly poor survival (Figure 3).
We then considered the remaining mesenchymal patients from non-OC sites, of which there were only 30 out of 418 TCGA patients, divided roughly equally between the larynx and oropharynx. The non-OC mesenchymal patients were 1/3 nodenegative and 2/3 node-positive. Unlike OC patients, the nodenegative patients did extremely well overall. Although the sample number was only 10 patients, there was no suggestion that non-OC mesenchymal patients did worse than nonmesenchymal patients in the node-negative state. By contrast, node-positive mesenchymal patients did poorly overall compared to node-positive non-OC patients. We then excluded HPV(+) patients, where the results were somewhat attenuated, but mesenchymal patients still fared worse than nonmesenchymal patients.

Discussion
This study confirms several important previous findings regarding gene expression subtypes in OC and HNSCC more broadly. First, we demonstrate that the mesenchymal and basal subtypes comprise the majority of OC cases. Second, we confirm previous reports suggesting that the mesenchymal subtype is associated with worse outcomes across all HNSCC sites. Importantly, we also identify novel and more nuanced associations between gene expression subtypes, stage, site, and treatment with important implications for future treatment stratification. Remarkably, we demonstrate that the mesenchymal subtype is associated with poor survival even in the setting of early-stage, node-negative OC treated with surgical resection. In contrast, our data demonstrate that mesenchymal subtype cases have favorable outcomes compared to other gene expression subtypes in early stage, non-OC cases, the majority of which were treated with definitive radiation therapy. These findings highlight the potential value of gene expression subtyping as an adjunct to pathology for treatment decision-making. Gene expression subtypes provide an objective method of molecular classification of HNSCC based on unsupervised clustering and are reflective of important differences in tumor biology. The four gene expression subtypes in HNSCC have been validated in multiple datasets, and similar classifications have been developed for lung cancer (3)(4)(5)20). In the present study, we demonstrate differences in gene expression subtype distribution by anatomic site. As previously reported, OC is comprised primarily of mesenchymal and basal subtypes, while classical is the predominant subtype in laryngeal squamous cell carcinoma (LSCC). Our group and others (13) have previously demonstrated the prognostic value of the mesenchymal subtype in HNSCC. The mesenchymal subtype is associated with EMT and predisposes to increased tumor invasiveness and lymph node metastases (5,(21)(22)(23). Recently, a partial EMT signature has shown further evidence of the importance of a mesenchymal phenotype in OC, suggesting that the transition from epithelial to mesenchymal phenotype represents a spectrum.
This study provides a more refined examination of the prognostic value of gene expression subtype in HNSCC that is specific to early-stage HNSCC. While the mesenchymal subtype is prognostic of worse survival in early-stage OC, there is no significant difference in outcomes between mesenchymal and other subtypes in non-OC early-stage tumors. Therapeutic decision-making and treatment dilemmas in HNSCC are anatomic site specific, and our data suggest that the potential clinical application of molecular subtyping should be considered within this context. These data also highlight the potential predictive application of gene expression subtype analysis in HNSCC. OC is generally treated surgically, with adjuvant radiation and chemotherapy reserved for advanced-stage tumors or adverse pathologic features such as positive margins and extranodal extension. Nevertheless, there is a subset of OC patients who recur even with early-stage disease and in the absence of adverse pathologic features (2). We have previously shown that the mesenchymal subtype is associated with an increased risk of occult nodal metastasis in the setting of clinically node-negative disease and suggest that a gene expression classifier applied to early stage HNSCC could potentially be used to assist in therapeutic decision-making (24).
We considered early on that risk conveyed by the mesenchymal subtype may simply be replaced in risk models by a gene expression-based lymph node positivity predictor. To investigate, we constructed a predictor of lymph node positivity similar to those reported by others (4,25). Briefly, our predictor had an accuracy of 66% in the OC training fit and only 57% in the independent test set, which was comparable to accuracies obtained by Chung et al. (2004) (4) (53%-60%) when OC tumors were included in classifier training. They found that removing OC tumors from training improved classifier performance; however, OC being a focus of our work, we decided to pursue the mesenchymal subtype as a biomarker of outcomes.
In the context of radiation, we considered possible explanations for the poor survival experienced by mesenchymal patients in some strata but not others. The overall inferior survival of patients with EMT signature, the most prominent component of the mesenchymal subtype, has been suggested by many prior reports (8,21,22,26). Mesenchymal tumors are characterized by EMT programs of gene expression, as well as inflammatory signatures that might be associated with worse outcomes. However, such programs might not be expected to have differential outcomes with respect to node-positive versus node-negative disease. One possible explanation would be differential treatment, especially radiation. In node-negative patients, radiation was administered at overall similar rates between mesenchymal and non-mesenchymal patients, 40% and 45% respectively, suggesting that differential radiation alone would not explain differential outcomes. The role of chemotherapy in nodenegative patients would be in conjunction with radiation and would only be limited to patients with positive margins, and as such, differential chemotherapy usage is also an unlikely explanation for differential outcomes. As expected, radiation is more common in node-positive patients. Among node-positive patients, mesenchymal patients were radiated at slightly lower rates overall, 62% versus 76%. Similarly, in the non-OC sites (larynx and oropharynx), node-negative mesenchymal and nonmesenchymal patients were radiated at somewhat higher rates, 56% and 49%, respectively, consistent with higher rates of radiation-based treatment of these disease sites. Node-positive non-OC mesenchymal and non-OC, non-mesenchymal patients were radiated at the highest rates of 73% and 88%, respectively, likely a combination of primary chemoradiation and adjuvant radiation cases. Although speculative, it is at least possible that part of the difference between patient groups might be due to the use of radiation, in which the poor prognosis of early-stage OC mesenchymal patients can be at least somewhat attenuated in higher stage compared to non-mesenchymal patients when they are radiated. This would argue that increased radiation of nodenegative OC mesenchymal patients might be beneficial.
We believe that while the conclusions of this study are well supported by the data presented, this study was limited by its retrospective nature and unavailability of patient treatment status from TCGA. Furthermore, the study would benefit from the availability of a larger validation dataset.
Predicting recurrence or relapse events in early-stage, HPVnegative HNSCC remains a significant challenge for clinicians. Despite significant improvements in our understanding of HNSCC molecular biology and prognostication, there is a paucity of biomarkers that have been developed to address this specific issue. Our data suggest that a gene expression classifier applied to early-stage HNSCC could potentially be used to assist in therapeutic decision-making.

Data availability statement
The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding author.

Ethics statement
The studies involving human participants were reviewed and approved by Institutional Review Boards of Washington University in St. Louis (IRB#201706088) and the University of Tennessee Health Science Center (IRB# 17-05549-XP). The ethics committee waived the requirement of written informed consent for participation.