Molecular Classification of Lymph Node Metastases Subtypes Predict for Survival in Head and Neck Cancer

Lymph node metastases (LNM) represent an aggressive yet curable state for many solid tumors. The clinical significance of LNMs reflects, in part, an intermediate metastatic step during the sequential spread of cancers from the primary site to the lymph nodes and then to distant organs (1–3). Alternatively, LNMs may not directly lead to distant metastases but reflect the biological selection of more aggressive subpopulations that are associated with increased local and/or regional recurrence (4–6). The presence of LNMs is associated with decreased survival in gastric cancers, melanomas, and breast cancers among others (7–9). Similarly, LNMs are one of the most important prognostic features in head and neck squamous cell cancers (HNSCCs) (; ref. 10) where the survival of patients with LNMs is approximately half the survival of patients without LNMs (11, 12). However, despite these worse outcomes, LNMs represent a heterogeneous disease state as evidenced by the inability of existing lymph node staging systems to accurately predict locoregional control, distant metastasis, and/or survival (13). Previous investigators have tried to stratify patients with LNMs using clinicopathologic features including lymph node density, tumor, and/or lymph node location and extracapsular extension, among others (14, 15). Although several groups have sought to biologically classify LNMs in HNSCCs, and in other cancers, these efforts have primarily focused on differences between LNMs and the primary tumor rather than defining distinct subtypes within LNMs (16, 17). Given that approximately 50%–60% of patients with HNSCC present with LNMs (18), identifying LNM subtypes based on intrinsic molecular differences would likely provide novel biological and prognostic information.

Here, we used integrative clustering of 72 LNMs, to identify three distinct LNM-specific molecular subtypes that were associated with different rates of locoregional control and survival. The most aggressive LNM subtype displayed a lymph node–specific gene expression signature that was prognostic across multiple cancer types and was enriched for genes associated with an invasive subtype.

From a database of nonmetastatic HNSCC (6), we identified 34 patients with pathologic lymph node–positive HNSCC treated with postoperative radiotherapy at University of Illinois (Chicago, IL) between March 26, 2001 and March 31, 2013 (UIC HNSCC dataset; Supplementary Fig. S1). The institutional review board (IRB) waived informed consent due to the difficulty in obtaining consent of patients with archived pathologic specimens. Patients were excluded due to stage I–II disease (n = 131), treatment with definitive radiation and/or did not have a therapeutic lymph node dissection prior to radiotherapy (n = 352), no pathologically involved lymph nodes on dissection (n = 89), no pathologic accession numbers and/or tissue blocks were not available (n = 79). The remaining 9 patients were excluded due to poor sequencing quality (n = 8) or the presence of HPV16 transcripts (n = 1). The resulting dataset is referred to as UIC HNSCC dataset. Biospecimens and deidentified linked patient data were collected as approved by University of Illinois at Chicago IRB guidelines (IRB 2011-1075) and data analysis was performed as approved by University of Illinois at Chicago and University of Chicago (IRB 2015-1008) IRB guidelines in accordance with the ethical standards of the responsible committee on human experimentation and with the Helsinki Declaration of 1975, as revised in 2000.

Follow-up was calculated from the completion of radiotherapy to the date of first failure or last follow-up. Margin status, perineural invasion, lymphovascular space invasion, extracapsular extension, number of positive lymph nodes, and number of lymph nodes dissected were based on the final pathology report. Treatment package time was calculated as the time from ablative surgery to the last day of radiotherapy. Events were determined from the last day of radiotherapy. Patterns of local (LF), regional (RF), or distant failure (DF) were documented as sites of first failure. Locoregional control (LRC) was calculated as the time to local or regional failure. Relapse-free survival (RFS) was calculated as the time to any failure or death from any cause. Overall survival (OS) was calculated as the time to death from any cause.

Statistical analysis was performed using JMP version 9 (SAS Institute). All tests of statistical significance were two-sided, and significance was defined as P < 0.05. The χ² test was used to compare between discrete variables and t test between continuous variables. Differences between medians were assessed using the Wilcoxon test. Survival analysis was performed for all patients. Censoring was noninformative. All events were calculated using Kaplan–Meier estimator with log-rank test, and the differences were compared using Cox regression models. For univariate analysis (UVA), we selected factors known to impact oncologic outcomes as well as patient and treatment characteristics. We used Cox proportional hazards model to examine the effects of these different risk factors on event outcomes for LRC. Cox multivariable analysis (MVA) was performed to adjust for explanatory confounding variables prognostic on UVA. Patient characteristics that were not recorded were not included during statistical analysis. Detailed description of patient and clinical variables are contained in Supplementary Information.

Archived formalin-fixed, paraffin-embedded (FFPE) tissue blocks of LNMs and matched primary tumors with the corresponding H&E slide were analyzed by expert head and neck pathologists (O. David and R.J. Cabay). The area of greatest density of tumor cells with at least 70% cancer cells were identified on the hematoxylin and eosin (H&E) slide and mapped to the corresponding FFPE block. Punch biopsies (2 mm) were obtained for the indicated region. RNA was extracted from isolated specimens using the RecoverAll Total Nucleic Acid Isolation Kit and stored at −80°C until further analysis. mRNA and miRNA libraries were sequenced on a HiSEQ2500 machine using standard reagents and protocols provided by Illumina. Sequencing data have been deposited at the European Genome-phenome Archive (EGA), hosted by the EBI and the CRG, under accession number EGAS00001003233. Read alignment, quantification, normalization are described in the Supplementary Information.

We applied an Integrative Clustering of Multiple Genomic Data Types (iCluster+; ref. 19) to the matched mRNA and miRNA expression data of the lymph node samples alone as well as primary tumor and lymph node samples combined. The algorithm was run for N = 2–6 with the number of lambda as 144. For each N, we obtained a model that contains a deviance ratio metric. The metric was selected on the basis of minimum Bayesian information criterion (BIC) for each lambda value, and can be interpreted as the percentage of total variation explained by the current model. The optimal number of clusters N was determined at the point when increasing N, the percent explained variation no longer significantly increases. Consensus clustering is described in the Supplementary Information.

To identify differentially expressed mRNAs among samples grouped by iCluster method, we removed heteroscedascity from the normalized log₂-transformed CPM data using the voomWithQualityWeights function from the R/Bioconductor package limma. We then fit a linear model for each gene using the limma algorithm, adjusted for patient covariate, and ranked the genes for differential expression using the empirical Bayes method with robust option enabled. For miRNAs, we applied the limma method to identify differentially expressed miRNAs among the samples grouped by iCluster method. We first estimated the relative quality weights for each array using the arrayWeightsSimple function, and then fit a linear model for each probeset adjusted for batch effect, followed by ranking probesets for differential expression using empirical Bayes method. The differentially expressed genes/miRNAs were identified with the Benjamini–Hochberg procedure for multiple testing adjustment and fold change. The false discovery rate (FDR) controlling adjusted P value and fold change threshold were set at 0.05 and 1.5 for mRNAs and 0.1 and 1.5 for miRNAs, respectively.

Raw gene feature counts were mapped to Entrez ID using the R/Bioconductor package org.Hs.eg.db v3.5.0. Low/nonexpressed genes with less than 1 log₂-transformed CPM across the minimum number of samples in any iCluster group were excluded from subsequent analysis using the R/Bioconductor package edgeR. Quality weighted, TMM normalized, and log₂-transformed CPM were calculated using limma-voom (verson 3.34.5). Gene-set enrichment analysis was performed with planned contrasts of one iCluster group against the remaining groups using GSEA (20) for comparison within the UIC HNSCC dataset or the R/Bioconductor package EGSEA (version 1.6.1; ref. 21) for comparison across other publicly available databases including HNSCC (ref. 16; Array Express E-MEXP-44), breast cancer (ref. 22; GSE56493), and melanoma (ref. 23; Array Express E-GEOD-65904) databases. Independent GSEA/EGSEA analyses were performed for gene lists provided by MSigDB v5.2 (http://software.broadinstitute.org/gsea/msigdb) containing hallmark gene set and KEGG pathway gene set. Enriched pathways from EGSEA were selected according to the median rank across 12 enrichment methods implemented in the package. Over-representation analysis was performed on differentially expressed genes that were the predicted targets of differentially expressed miRNAs. miRNA target prediction was carried out using miRWalk 2.0.

To build a classifier to distinguish samples between Group 2 and Groups 1 and 3 combined identified by iCluster method, the normalized mRNA expression data of the metastatic lymph nodes from 33 patients was split into a training set and test set, where 75% samples were in training set and 25% samples were in test set. The group ratio remained unchanged during the partition. For the training set, we first filtered genes with near zero-variance. We then identified highly correlated genes with a pair-wise absolute correlation coefficient greater than 0.75, and removed those with the largest mean absolute correlation. We further removed potential linear dependencies of the data using the findLinearCombos function from the R package caret (version 6.0). We applied the preProcess function to center and scale the training and test data by mean and SD, followed by rescaling data to −1 and 1. We applied Prediction Analysis of Microarrays (PAMR, version 1.55), a nearest shrunken centroid classification algorithm, on the training set (24). A 10-fold cross-validation was performed to obtain cross-validated misclassification error rate for a range of threshold values. The threshold associated with the minimum cross-validated misclassification error rate was retained. We evaluated the classifier with this threshold using the held-out test data. Performance metrics such as accuracy, balanced accuracy, sensitivity, specificity, positive prediction value (PPV), negative prediction value (NPV), Cohen Kappa, Matthew correlation coefficient (MCC), and area under the curve (AUC) were calculated using the confusionMatrix function from the caret package and an in-house R script. We ran this procedure 100 times randomly and selected the model based on the classifier's performance measure described above. For models with same performance measure, we arbitrarily selected the less complex one with the number of genes (features) less than 100 and the overall cross-validation error rate less than 10%. The final classification model, containing 73 genes, was trained on 22 Group 2 samples and 33 Group 1 and 3 samples; evaluated using the held-out test data of 7 Group 2 samples and 11 Group 1 and 3 samples. The accuracy, no information rate, sensitivity, specificity, positive predictive value, negative predictive value, and misclassification errors are described in Supplementary Fig. S2.

To evaluate the LNM classifier performance in other tumor types, we collected LNM tumor expression data from three published studies. The studies contain microarray expression data from breast cancer (ref. 22; GSE56493), melanoma (ref. 23; GSE65904, 10,000 most variable genes), and HNSCCs (ref. 16; MDACC HNSCC dataset; E-TABM-114) from Gene Expression Omnibus (GEO) and ArrayExpress database. For each dataset, we downloaded the normalized gene expression data from the original studies. The probesets (Affymetrix) or probes (Illumina) were converted to Ensembl gene identifiers. The expression data of the 73 genes were scaled to −1 and 1. Values for the missing genes were set to −1. The LNM classifier was then applied to the data to predict the test samples as Group 2 or Group 1 and 3 combined.

We evaluated the LNM classifier's ability of distinguishing different survival group in our patients with HNSCC with the comparison of four published gene signatures (25–28). For each signature, we collected the gene signatures from the original publication and converted the gene symbol or alias to Ensembl gene identifier, then computed the score based the gene expression in our dataset. We used the following criteria to group our patients into two risk groups and performed the Kaplan–Meier survival estimates for the risk groups:

•  (i) Enokida classifier: PI>75 quartile – PI high, otherwise – PI low

•  (ii) RSI classifier: RSI < median - Responder, otherwise – Non-responder

•  (iii) 13 OSCC classifier: RS > median – High risk, otherwise – Low risk

•  (iv) 172-gene classifier: RS > mean – High risk, otherwise – Low risk

We identified 34 patients with human papillomavirus (HPV)-negative HNSCC who had accessible primary tumors and matched LNMs (6). These patients were treated with surgical resection of the primary tumor and lymph node dissection followed by radiotherapy with or without chemotherapy (Fig. 1A). Because differences in coding and noncoding gene expression may distinguish clinically relevant subtypes in LNMs, we performed RNAseq and miRNAseq on 29 tumors and 72 matched LNM (termed UIC HNSCC dataset). Supplementary Table S1 lists the clinicopathologic characteristics of this group where the majority of patients had N2–N3 disease (81.8%) and received postoperative chemoradiation (78.8%). The 2 year LRC, progression-free survival, and OS in the entire cohort was 66.7%, 55.1%, and 52.7%, respectively, similar to previous reported outcomes in HPV-negative HNSCCs (refs. 11, 12; Supplementary Fig. S3).

To identify distinct metastatic subtypes, we performed integrative clustering on the 72 LNMs to identify potential lymph node–specific subtypes. Integrative clustering identified three distinct molecular subtypes: 27 LNMs were classified as Group 1, 29 LNMs were classified as Group 2, and 16 LNMs were classified as Group 3. Of the 16 patients with multiple LNMs analyzed, 5 patients had LNM that all belonged to the same molecular subtype and 11 patients had LNMs belonging to multiple subtypes (Supplementary Table S2).

We assessed the clinical relevance of the Group 1, Group 2, and Group 3 LNM subtypes by comparing clinicopathologic variables, disease recurrence, and survival. Because 11 of 17 patients with multiple LNMs displayed more than one LNM subtype, we classified patients with multiple LNMs as Group 2 if patients had ≥2 Group 2 LNMs. If patients with multiple LNMs had <2 Group 2 LNM, these patients were classified according the most prevalent LNM subtype. On univariate analysis, patients with Group 3 LNM was associated with greater lymphovascular space invasion (P = 0.02) and patients with Group 2 LNM was associated with fewer numbers of lymph nodes removed at the time of surgery (P = 0.03; Supplementary Table S3). Patients with Group 2 LNMs displayed significantly worse LRC (1y LRC Group 1: 91.7% vs. Group 2: 0.0% vs. Group 3: 72.9%; P = 0.002), RFS (P = 0.002), and OS (1y OS Group 1: 92.9% vs. Group 2: 44.4% vs. Group 3: 100.0%; P = 0.02; Fig. 1B–E; Supplementary Fig. S4). On multivariate analyses, Group 2 subtype remained the only predictor for LRC compared with either the Group 1 subtype [HR: 0.1; 95% confidence interval (CI), 0.01–0.71; P = 0.02] or the Group 3 subtype (HR: 0.2; 95% CI, 0.03–0.99; P = 0.04; Supplementary Table S4). There was no difference in LRC between Group 1 and Group 3 subtypes (HR: 2.04; 95% CI, 0.27–17.82; P = 0.48). The median number of lymph nodes removed was not a significant predictor of LRC on univariate analysis (Supplementary Table S4).

The classification of these three molecular subtypes was primarily due to differences in mRNA expression and not miRNA expression. Using pair-wise comparisons, we detected 2,039 differentially expressed mRNAs (DEmRNA) and 24 differentially expressed miRNAs (DEmiRNA; Fig. 1F; Supplementary Fig. S5A and S5B; Supplementary Table S5). Consensus clustering, consensus cumulative distribution function, and proportion of ambiguous clustering of DEmRNAs confirmed that three molecular subtypes were the optimal classification of LNMs (Fig. 1G–I). In contrast, consensus clustering of DEmRNAs using other predefined cluster numbers of K = 2, 4, 5, or 6 exhibited less optimal cluster groupings (Supplementary Fig. S6). Furthermore, more than three predefined clusters failed to further stratify differences in LRC. Predefined clustering of DEmiRNAs did not identify any subtypes having differences in LRC (Supplementary Fig. S5C–S5G). Therefore, molecular classification of LNMs in HPV-negative HNSCCs using DEmRNAs defined three distinct subtypes where one subtype, Group 2, had significantly worse LRC, RFS, and OS.

We performed gene-set enrichment analysis (GSEA) to better understand the distinct biological pathways characterizing the different LNMs subtypes (Fig. 2; Supplementary Fig. S7). The Group 2 LNM subtype was associated with pathways characterizing an invasive subtype with enriched pathways including epithelial–mesenchymal transition, apical junction, TGFβ signaling, angiogenesis, hypoxia, extracellular matrix receptor interactions, regulation of the actin cytoskeleton, and focal adhesion (Fig. 2A and D). Group 1 was associated with pathways characterizing an immune subtype with enriched pathways including allograft rejection, T-cell receptor signaling pathway, and chemokine receptor signaling pathway among others (Fig. 2B). Group 3 was associated with pathways characterizing a proliferative/metabolic subtype with enriched pathways including MYC targets, basal transcription factors, mismatch repair (Fig. 2C). Consistent with its enrichment of immunologic pathways, Group 1 had a higher estimated proportion of B cells, CD4⁺ T cells, and CD8⁺ T cells (Fig. 2E; Supplementary Fig. S8). Of note, Group 1 and Group 3 often displayed an inverse pattern of enriched gene sets. In contrast, Group 2 differed from Group 1 and Group 3 primarily in gene sets involved in the cell membrane and the extracellular matrix. Furthermore, pathway analysis using differentially expressed miRNA–mRNA pairs identified immune-regulated pathways specific to the Group 2 LNMs (Supplementary Table S6).

Previous reports in HNSCCs as well as in other cancers have demonstrated that the transcriptional profiles of LNMs were similar to the primary tumor (16, 17, 29, 30). Consequently, we sought to determine whether unbiased classification methods would identify similar molecular subtypes in matched primary tumors. We performed RNAseq and miRNAseq on 29 matched primary tumors from patients with LNMs. Integrative clustering identified two molecular subtypes in primary tumors. However, these primary tumor subtypes, Tumor Group 1 and Tumor Group 2, did not differ in LRC (1y LRC: 63.5% for Tumor Group 1 vs. 64.0% for Tumor Group 2; P = 0.92; Fig. 3A). Unlike LNMs, primary tumors did not have any DEmRNAs between groups, but there were 14 DEmiRNAs that confirmed the two subtypes as optimal stratification based on consensus clustering and proportion of ambiguous clusters (Fig. 3B and C). Given that unbiased approaches did not identify similar molecular subtypes in primary tumors, we performed clustering of primary tumors based on the DEmRNAs identified in the LNMs subtypes. Clustering based on the LNMs DEmRNAs identified 6 clusters that displayed similar rates of LRC (Fig. 3D–F). To exclude the possibility that the unique expression profiles of LNMs, which were not seen in primary tumors was due to normal lymph node rather than unique metastatic subtypes, we performed RNAseq and miRNAseq of 8 nonmetastatic lymph nodes isolated from 8 different patients in our cohort. Hierarchical clustering of either miRNA or mRNA expression data from these uninvolved lymph node specimens demonstrated that normal lymph nodes clustered separately from the LNMs cohort. (Supplementary Fig. S9). Furthermore, hierarchical clustering of mRNA expression data from two publically available microarray datasets of normal lymph nodes demonstrated that normal lymph nodes clustered separately from the LNMs in our cohort (Supplementary Fig. S10). Consequently, the molecular subtypes identified in LNMs were not observed in the matched primary tumors.

The ability to identify clinically relevant molecular subtypes in LNMs but not in primary tumors may reflect the biological selection of more aggressive subpopulations at the metastatic site. To determine whether LNMs represented more aggressive variants, we assessed the ability of 4 prognostic HNSCC classifiers to predict locoregional recurrence and RFS using the primary tumor or LNMs gene expression profiles. We used (i) a 172-gene signature derived from semisupervised clustering and validated in a meta-analysis of 646 primary HNSCCs (25), (ii) a 30-gene prognostic index (PI) validated in 123 primary oral cavity carcinomas (26), (iii) a 10-gene radiosensitivity index (RSI) derived from a systems biology model and validated in 92 primary HNSCCs (27), and (iv) 13-gene risk score derived from an L1-penalized Cox proportional hazard regression model and validated in 103 primary oral cavity cancers (28). When the primary tumor gene expression profiles were used in analyses, all four prognostic classifiers failed to stratify patients with worse LRC and 2 of 4 classifiers failed to stratify patients with worse RFS (Fig. 4; Supplementary Fig. S11). In contrast, when the LNM gene expression profiles were used in analysis, 3 of the 4 classifiers stratified patients with worse LRC (173-gene signature P = 0.04; 10-gene RSI P = 0.02; 13-gene RS P = 0.001) and all four classifiers stratified patients with worse RFS. Therefore, the gene expression profiles of LNMs better predicted outcomes than the gene expression profiles of primary tumors.

Previously published HNSCC prognostic signatures (25–28) were able to distinguish the aggressive Group 2 variant with 61.1% to 80.6% accuracy (Supplementary Table S7). Consequently, classifiers developed from primary tumors to identify high-risk subtypes may not adequately identify all LNMs that are high-risk subtype. We used a discovery set of 52 randomly selected LNMs to generate a nearest shrunken centroid-based 73-gene binary classifier that distinguished Group 2 LNMs from Group 1 and Group 3 LNMs (Groups 1 and 3; Supplementary Table S8). Using the remaining 18 LNM specimens as a testing cohort, this 73-gene classifier had an accuracy of 0.89 (95% CI: 0.65–0.99; P = 0.01) with a sensitivity of 1.00, specificity of 0.82, positive predictive value of 0.78, and negative predictive value of 1.00. The ROC plot of this 73-gene classifier had an AUC of 0.987 (95% CI, 0.951–1.000; Supplementary Fig. S12). As expected, the 73-gene LNMs classifier identified a clinically significant Group 2 subtype having worse LRC when applied to LNMs, but not primary tumor specimens in the UIC HNSCC cohort (Fig. 5A). We next assessed the extent to which this 73-gene classifier predicted outcomes using the gene expression profiles from LNMs or non-LNMs in other cancers. We utilized a prospectively obtained breast cancer dataset (22) isolated from initial recurrences in LNMs (n = 44), primary tumors (n = 19), liver metastases (n = 27), or other sites (n = 30; Supplementary Table S9). In addition, we utilized a dataset of retrospectively obtained melanoma specimens (23) isolated from the primary tumor (n = 16), LNMs (n = 139), distant metastases (n = 23), in transit metastases (n = 15), or local recurrences (n = 11; Supplementary Table S10). Group 2 subtypes displayed worse overall survival and disease-specific survival when applied to the transcriptomes of breast and melanoma LNMs, respectively (Fig. 5B and C). In contrast, when other non-LNM sites were analyzed, Group 2 subtypes were not associated with differences in outcomes consistent with our observations in HNSCC primary tumors.

To confirm that this 73-gene classifier identified a Group 2 LNM subtype, we performed ensemble GSEA (EGSEA) to identify pathways upregulated in Group 2 subtypes in breast cancer LNM (22), melanoma LNM (23), as well as in a small retrospective series of LNMs from the MDACC HNSCC cohort (16). Group 2 subtypes identified from UIC HNSCC, MDACC HNSCC, breast, and melanoma cohorts had a significant degree of concordance (Kendall concordance coefficient; Hallmark W = 0.67; P = 3.2 × 10⁻⁹; KEGG W = 0.52; P = 6.13 × 10⁻¹⁶). Of the 49 shared Hallmark gene sets, the top three median ranked pathways were epithelial–mesenchymal transition (median rank: UIC HNSCC 5, Breast 3.5, Melanoma 1 and MDACC HNSCC 2), coagulation (median rank: UIC HNSCC 2, Breast 5.5, Melanoma 6 and MDACC HNSCC 7.5), and hypoxia (median rank: UIC HNSCC 14.5, Breast 5, Melanoma 4 and MDACC HNSCC 4.5; Fig. 6A). In the 187 shared gene pathways of the KEGG database, the top three median ranked pathways were ECM–receptor interactions (median rank: UIC HNSCC 16, Breast 5.5, Melanoma 5, and MDACC HNSCC 8), protein digestion and absorption (median rank: UIC HNSCC 16, Breast 10.5, Melanoma 3, and MDACC HNSCC 23), and focal adhesion (median rank: UIC HNSCC 15.5, Breast 10.5, Melanoma 5.5, and MDACC HNSCC 23; Fig. 6B and C). Therefore, in four different datasets encompassing HNSCCs, breast cancers, and melanomas, this LNM-specific classifier identified a Group 2-like subtype displaying similar biological subtypes and clinical outcomes.

Here, we identified three distinct molecular subtypes of LNMs that were associated with differences in clinical outcomes after surgery and postoperative radiotherapy. These subtypes, termed Group 1, Group 2, or Group 3, displayed an immune, an invasive and a metabolic/proliferative phenotype, respectively. In contrast to Group 1 and Group 3 LNM subtypes, the Group 2 LNM subtype was associated with significantly worse locoregional recurrence, RFS, and OS. These subtypes are likely specific to LNMs for several reasons. First, LNMs subtypes were not observed in the matched primary tumors. Second, a similar Group 2 subtype was also identified in breast, melanoma, and MDACC HNSCC LNM datasets and displayed similar gene expression patterns as well as increased risk for recurrence and death. Similar to the absence of these LNMs subtypes in primary HNSCCs, this Group 2 subtype was not associated with recurrence and survival differences when classifying non-LNMs sites. Furthermore, the lack of nonmetastatic lymph nodes clustering separately from the LNMs support the observation LNMs represented distinct molecular subtypes that were not observed in primary tumors. Even though we identified these LNMs subtypes using a single institutional retrospective dataset encompassing 34 patients, we validated similar observations in existing prospective and retrospective datasets totaling 373 patients across three different cancer types. Thus, we describe a novel molecular classification of LNMs metastasis, which has both biologic and prognostic significance in HNSCCs and likely extends to LNMs in other cancers.

We defined a 73-gene classifier to identify a Group 2-like subtype across multiple cancer types. Although a minority of samples in training and testing sets included LNMs from the same patient, this classifier was not dependent on patient-specific characteristics because it distinguished a similar Group 2 subtype in external validation sets. Furthermore, a separate analysis using samples from the same patient that were restricted to either the training or testing sets maintained high sensitivity (100%), specificity (92.9%), and accuracy (95.2%) as the original classifier. This classifier also predicted outcomes in the breast cancer dataset when applied to LNMs (P = 0.04) but not when applied to other sites (P = 0.45; data not shown). Still, the overall patient classification may also be dependent upon the number of LNMs tested as fewer LNMs may underestimate the classification of high-risk patients.

Previous attempts to biologically classify HNSCCs have mostly relied upon analyses incorporating known patient outcomes to identify mutations and/or differentially expressed genes that correlated with worse survival (25, 28, 31–33). In contrast, fewer studies have used unbiased strategies to classify HNSCCs into distinct biological subtypes that would predict different biological outcomes (26, 34, 35). In either case, previous reports have almost always focused on the primary tumor which likely masks small subpopulations of aggressive variants due tumor heterogeneity (25, 26, 34–39). To this end, LNMs subtypes identified here were not detected in matched primary tumors using unbiased integrative clustering, clustering of DEmRNAs, or a 73-gene LNM-specific classifier. Our findings are consistent with a few reports that demonstrate lymph node and/or distant metastases are molecularly distinct from the primary tumor (40). Unlike our results describing LNM-specific subtypes, most previous studies of LNMs in HNSCCs and other cancers demonstrated mutational and transcriptional profiles that were most similar to the matched primary tumors (16, 17, 29, 41–43). However, these studies were not designed to identify LNM-specific molecular subtypes because analyses combined both the primary tumors and LNMs from multiple patients. By comparing LNMs to primary tumors, the variance between patients likely confounded the identification of differences between primary tumors and LNMs and, therefore, was not optimized to identify LNM-specific differences. Furthermore, compared with previous studies (16, 17, 44–46), we used 5- to 7-fold more LNMs samples that provided greater power to identify distinct LNM subtypes. Rather, our data are consistent with the premise that LNMs represent biologically selected subpopulations that display unique characteristics that are distinct from the genetically heterogeneous primary tumor. Furthermore, this selection of distinct subpopulations in LNMs likely reflects the observation that a single patient may have multiple LNMs assigned to different molecular subtypes. Finally, these LNMs may reflect the selection of cells that were not more invasive than other but that proliferated in nodal tissues and evaded the immune system.

Although the number of patients in our study was relatively small, the focus of LNMs enabled us to identify unique LNMs subtypes that were predictive of clinical outcomes. Previous studies using larger sample sizes of only primary tumors have identified gene expression patterns that correlated with survival. However, many of these previously published gene expression patterns could not predict outcomes when using the gene expression patterns of primary tumors in our cohort. When applied to the gene expression patterns of LNMs, these previously published gene expression classifiers better predicted outcomes suggesting that LNMs rather than primary tumors better reflected the biology that dictates patient outcomes. Furthermore, the distinct microenvironment in lymph nodes may explain the inability to detect similar LNMs subtypes in the primary tumor. Therefore, classifying the molecular differences present in LNMs will provide novel biological and prognostic information.

Unlike the miRNA expression patterns in primary tumors, few, if any studies have reported miRNA expression patterns in LNMs of HNSCCs. Although the inclusion of miRNA may have facilitated the differentiation of molecular subtypes observed in LNMs, mRNA expression was sufficient to differentiate the molecular subtypes of LNMs indicating that the miRNA features were less likely to have an impact. Furthermore, we focused on mRNA expression data in the 73-gene classifier because mRNA expression was sufficient to predict LNM-specific subtypes and could be validated using external datasets. In contrast, the lack of miRNA expression data in external datasets limited us from using this data in such a classifier. Furthermore, the miRNA data, alone, was not able to stratify the LNM-stratify subtypes highlighting the importance of mRNA expression in LNMs.

These lymph node–specific molecular subtypes represented an invasive (Group 2), inflammatory (Group 1), and metabolic/proliferative subtype (Group 3). Previous groups have described 3 to 6 subtypes in primary HNSCCs (35, 38, 47, 48). Pathways enriched in the Group 2 LNM subtype were similar to pathways enriched in a basal hypoxia or angiogenesis subtype in primary HNSCCs (34, 35, 38, 47, 48). Similarly, pathways enriched in the Group 2 or Group 3 LNMs subtype were similar to pathways enriched in the inflamed/mesenchymal subtype (35, 38, 47, 48) or classical-atypical subtype (34, 35, 47, 48) in primary HNSCCs, respectively. However, in contrast to our observations, differences in molecular subtypes of primary HNSCCs were not often associated with differences in clinical outcomes. Although Chung and colleagues did find that intrinsic classification system in addition to other clincopathologic features were associated with worse recurrence-free survival (34, 37), we observed that the LNM-subtype was the sole factor associated with cancer recurrence. While Keck and colleagues (35) observed survival differences between HNSCC subtypes, these differences were mainly due to differences in survival between HPV-positive and HPV-negative HNSCC subtypes. In contrast, our study identified subtypes of HPV-negative HNSCCs that were associated with significant survival differences. One explanation for these differences between our findings and others is that LNMs represented the biological selection of aggressive variants compared with a more genetically heterogeneous primary tumor.

LNMs likely reflect biological selection of clinically relevant subpopulations as well as display LN-specific gene expression patterns. The biological selection of LNMs was supported by the observation that existing prognostic classifiers in HNSCCs (25–28) identified patients at increased risk of locoregional failure and RFS primarily when applied to the transcriptional profiles of LNMs but not to the transcriptional profiles of primary tumors.

However, these existing classifiers only predicted the Group 2 LNM subtype with 60%–80% accuracy. In contrast, a novel 73-gene classifier identified Group 2-like subtypes in multiple cancer types only when applied to the gene expression patterns of LNMs, but not when applied to non-lymph node sites indicating that LNMs also reflect lymph node–specific changes in the gene expression patterns of metastases. Similarly, in liver metastases, a selective gene signature was associated with adverse outcomes in patients with a luminal A breast cancer subtype (49). Thus, we report a novel LNM-specific classifier that predicts outcomes using both retrospective and prospective datasets encompassing multiple cancer types.

In conclusion, using unbiased approaches, we have identified a unique subtype of LNMs in HNSCCs that displayed an invasive subtype and was associated with high rates of rapid disease recurrence. LNMs likely reflected both the selection of more aggressive variants as well as specific signatures shared across multiple cancer types. These results provide the basis for the molecular classification of LNMs in HNSCCs, and likely in other cancers, that may be used to predict outcomes and to identify pathways that can be targeted for therapeutic gain.

No potential conflicts of interest were disclosed.

Conception and design: L. Huang, O. David, R. Zhong, M.T. Spiotto

Development of methodology: L. Huang, O. David, V. Macias, R. Zhong, M.T. Spiotto

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): O. David, R.J. Cabay, K. Valyi-Nagy, R. Zhong, B. Wenig, L. Feldman, M.T. Spiotto

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): L. Huang, R.J. Cabay, R. Zhong, R. Weichselbaum, M.T. Spiotto

Writing, review, and/or revision of the manuscript: L. Huang, O. David, R.J. Cabay, K. Valyi-Nagy, V. Macias, R. Zhong, L. Feldman, R. Weichselbaum, M.T. Spiotto

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): R. Zhong, R. Weichselbaum, M.T. Spiotto

Study supervision: R. Zhong, R. Weichselbaum, M.T. Spiotto

Other (review of histologic slides and identification and marking of tumor and lymphoid tissue to be used for preparation of tissue microarray): O. David

Other (pathology verification of regions of interest and tissue microarray construction): V. Macias

The authors thank Drs. Nicholas Tobin and Thomas Hatschek (Karolinska Institute and University Hospital) for access to the TEX dataset. This work was supported by Burroughs Wellcome Career Award for Medical Scientists (1010964; to M.T. Spiotto) and a grant from the NIH (R01DE027445-01; to M.T. Spiotto). This work was supported in part by the Virginia and D.K. Ludwig Fund for Cancer Research (R. Weichselbaum and M.T. Spiotto). The Center for Research Informatics is funded by the Biological Sciences Division at the University of Chicago with additional funding provided by The Institute for Translational Medicine/Clinical and Translational Award (NIH 5UL1TR002389-02; Center for Research Informatics), and The University of Chicago Comprehensive Cancer Center Support Grant (NIH P30CA014599; Center for Research Informatics); Burroughs Wellcome Career Award for Medical Scientists 1010964 (to M.T. Spiotto); NIH/NIDCR R01DE027445-01 (to M.T. Spiotto), Institute for Translational Medicine, NIH/CTSA 5UL1TR002389-02 (Center for Research Informatics, CRI); The University of Chicago Comprehensive Cancer Center Support Grant (NIH P30CA014599) (CRI); and Virginia and D.K. Ludwig Fund for Cancer Research (R. Weichselbaum and M.T. Spiotto).

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.