Development and Validation of a Combined Hypoxia and Immune Prognostic Classifier for Head and Neck Cancer

Head and neck cancer (HNC) is the sixth most common cancer worldwide, with an annual incidence of around 600,000 cases. It has a poor outcome with a 5-year mortality rate of around 50% (1). The main risk factors include smoking, alcohol, and infection with high-risk human papillomaviruses (HPV), the latter conferring considerably better survival outcomes than the former two (2).

While many prognostic biomarkers have been described previously (2, 3), all include a combination of clinical factors and HPV status, and currently there is no widely accepted molecular classification. Importantly, there are no treatment response classifiers. As a result, patients receive treatment based on disease stage and patient/clinician preference, rather than on the biology of the tumor.

Recently, there have been advances in our understanding of the biologic and molecular characteristics of HNC. A comprehensive multi-platform genomic characterization by The Cancer Genome Atlas (TCGA) Network confirmed four previously defined subtypes: classical, mesenchymal, basal, and atypical (4–6). A more recent study identified five subgroups: HPV-positive mesenchymal, non-HPV mesenchymal, HPV classical, non-HPV classical, and basal (non-HPV) subtypes (7). When combined, the findings of these studies suggest three main subgroups: an inflamed/mesenchymal subtype demonstrates high expression of immune response genes and mesenchymal markers. A classical subtype (regardless of HPV status) has increased genomic aberrations associated with tobacco use. Finally, a basal subtype is highly enriched for hypoxia signaling with a lack of immune-related markers. Hypoxia is known to induce immunosuppression, both directly (e.g., via inhibition of T-cell proliferation and effector cytokine production) and indirectly (e.g., through metabolic competition, upregulation of coinhibitory receptors, or recruitment/conversion of immunosuppressive cell populations; refs. 8–14). Furthermore, HIF-1α signaling can be switched on by EGFR signaling (15, 16), which was also increased in some basal cancers.

Using this as a basis for a molecular classification, we then researched the availability of potential therapeutic paradigms for each of the subgroups. Hypoxic modifiers have been shown to have significant effects on survival in HNC (17, 18). More recently, evidence has emerged for the efficacy of immune checkpoint inhibitors in HNC (19, 20).

Combining these together, we sought to develop and validate a prognostic molecular classifier, based on immune response and hypoxia status, first in silico from fresh-frozen tissues, and then using targeted RNA sequencing of formalin-fixed, paraffin-embedded (FFPE) tissue samples to facilitate ease of clinical application. We then evaluated whether an IHC-based signature could substitute for, and simplify, the developed molecular signature, to further facilitate clinical implementation and cost reduction.

The in silico discovery cohort (characteristics and treatment detailed in Table 1 and Supplementary Table S1, respectively: TCGA) comprised samples from the HNC dataset in the TCGA, which contains whole-transcriptome gene expression data from fresh-frozen samples. Because HPV positivity is an independent prognostic factor for HNC (2), only cases with known HPV status were considered; this yielded a cohort of 275 patients, the majority being HPV negative.

The in silico validation cohort comprised microarray gene expression data for fresh-frozen samples from 134 patients with locoregionally advanced HNC [ref. 7; detailed in Table 1; Supplementary Table S1 (Chicago)], including both HPV-positive and -negative cases.

We then evaluated FFPE diagnostic biopsy or surgical resection samples from a retrospective cohort of 163 patients with oropharyngeal cancer, recruited to the PET-NECK or Predictr clinical studies (Table 1; Supplementary Table S1, Correlate). A cohort of 12 patients who underwent tonsillectomy for management of a nonmalignant process (snoring), recruited via the Oromouth study, served as controls. Ethical approval for use of tissue samples in translational research was granted by North West - Preston Research Ethics Committee (Reference: 16/NW/0265). p16 expression was utilized as a surrogate marker for HPV status. IHC staining for p16 was performed using the CINtec Histology Kit (Roche Laboratories); samples with ≥70% strong diffuse nuclear and cytoplasmic staining of tumor cells were considered positive (21).

Studies were conducted in accordance with the Declaration of Helsinki. Patients gave informed written consent to participate or anonymized samples were used in accordance with Human Tissue Act rules of the United Kingdom, and the above Research Ethics Committee approval.

For in silico development, expression data for genes comprising the hypoxia (Eustace and colleagues; ref. 22) and immune signatures (CIRC; ref. 23), plus additional genes of interest (including other immune-related genes and genes frequently mutated in HNC; Supplementary Table S2) were filtered from HT-seq gene count files downloaded from the TCGA HNC dataset (ref. 6, https://portal.gdc.cancer.gov/projects/TCGA-HNSC). Each dataset was normalized by dividing the expression values by the sum of expression values of the analyzed genes for each sample, then log₂ transformed; z-scores were calculated for these values using R v3.3.2 (https://www.r-project.org/).

Cluster analyses for genes/samples were obtained from the log₂ values matrix for each dataset using Spearman distance and Ward criterion in R. These clusters were then used to plot heatmaps using the z-score values matrix for color intensities. Survival was assessed by Cox regression analyses using Kaplan–Meier curves to compare the three highest hierarchical sample groups on each heatmap. In addition, a Cox regression multivariate was used to calculate the combined effect of HPV status and heatmap groups, using the Survival package in R (https://cran.r-project.org/web/packages/survival/index.html). Outcome measure was overall survival (OS) from treatment end date to death or last follow-up and censure.

For in silico validation, data from a previously published cohort (7) were interrogated for expression of the above gene signatures. Data for four genes COL4A6 (Eustace and colleagues), CD80, IFNG, and PDCD1LG2 (CIRC) were unavailable.

Using the data from the development TCGA cohort, to identify a difference of 25% in 3-year survival (from 75% immune^high to 50% hypoxia^high and mixed groups) with a power of 80% and a two-sided alpha of 0.05, we would require 138 cases; at 85% or 90% power we would require 158 or 185 cases, respectively (using the Kelsey calculation method and OpenEpi v3.0 software).

Gene expression was quantified using Illumina TruSeq Targeted RNA technology. RNA was extracted from 3–8 FFPE sections of 10 μm thickness using the PureLink FFPE RNA Isolation Kit (Thermo Fisher Scientific), DNase treated, and then quantified using the Qubit RNA BR Assay Kit (Thermo Fisher Scientific). Libraries were prepared using a TruSeq Targeted RNA Custom Panel Kit (Illumina), according to the manufacturer's instructions. The panel comprised 83 genes (see Supplementary Table S2 for full details). The resultant libraries were pooled and sequenced using the Illumina MiSeq platform, with a MiSeq reagent kit v3 150 cycle and paired end reads.

Quality assessment and filtering of the reads were carried out using FastQC v0.11.2 (http://www.bioinformatics.babraham.ac.uk/projects/fastqc) and PrinSeq Lite v0.20.4 (http://prinseq.sourceforge.net/), respectively. Reads were aligned to appropriate reference gene sequences obtained from GenBank (see Supplementary Table S2 for reference numbers), using BWA samse v0.7.15 (http://bio-bwa.sourceforge.net/). The obtained SAM files were then sorted and indexed using SAMtools v1.2 (http://samtools.sourceforge.net/) and read count by gene files was obtained using subread featureCounts v1.5.0-p1 (http://subread.sourceforge.net/). The expression matrix comprising genes within the analyzed signatures was then submitted to a similar normalization and analysis as described for the discovery dataset.

All staining was performed on FFPE sections of 4- to 5-μm thickness. Hematoxylin and eosin (H & E) staining was carried out using a Leica Autostainer xl on program 1. Automated IHC staining was carried out using a Leica Bond Max with standard F protocol and the following primary antibodies: Novocastra Liquid Mouse Monoclonal CA-IX Antibody clone TH22 (Leica Biosystems), 1:100 dilution; PD-L1 (E1L3N) XP Rabbit mAb (Cell Signaling Technology), 1:200 dilution.

For manual IHC staining, sections were deparaffinized in xylene then rehydrated in graded concentrations of alcohol. Following heat-induced epitope retrieval, slides were stained using the Novolink Polymer Detection System (Leica Biosystems) and the following primary antibodies: FLEX ready-to-use polyclonal rabbit anti-human CD3 (Dako) and mouse monoclonal anti-human LAG3 Antibody clone 17B4 (LSBio), 1:250 dilution.

IHC slides were scored independently by two trained calibrated pathologists; when scoring was discrepant, slides were rescored and a consensus reached. H & E slides were scored for tumor-infiltrating lymphocytes (TIL) as described by Ward and colleagues (3). CD3 was semiquantitatively estimated (range 1–4) using the mean score from three representative high power fields, where 1 = no, or sporadic CD3⁺ cells, 2 = moderate numbers of CD3⁺ cells, 3 = abundant occurrence of CD3⁺ cells, and 4 = highly abundant occurrence of CD3⁺ cells; and was evaluated separately for three regions: intratumoral, peritumoral stroma, and tumor advancing margin (adapted from ref. 24). Median values were used for survival analysis dichotomization. Cases were also assigned an immune profile (desert, excluded, or inflamed; ref. 25) based on CD3 expression using the following cutoffs: immune desert, score of <2 in all regions; immune excluded, difference of >1 between stromal versus intratumoral score; immune inflamed, score of ≥2 in two or more regions.

PD-L1 was evaluated on both tumor cells and immune-infiltrating cells (morphologically identified as lymphocytes and macrophages/dendritic cells; adapted from refs. 26, 27). For tumor cells, the percentage of positive cells was estimated for whole sections; four different cutoffs defining positivity were considered: ≥1%, ≥5%, ≥25%, and ≥50% (28). The pattern of staining was also recorded, as either constitutive (most/all of tumor cells positive) or inductive (peripheral staining on the interface between tumor and stroma or discrete regions of staining within tumor nests). Immune cells within tumor nests or within the peritumoral stroma were semiquantitatively scored as follows: 0 = negative, 1 = low abundance, and 2 = high abundance. Tumor and inflammatory cell expression were also pooled to give a combined positive score (CPS). The presence/absence of immune cells surrounding the tumor and forming a PD-L1–positive cordon was also noted.

LAG3 expression on tumor-infiltrating immune cells was scored as the mean percentage of positive cells from 10 high power fields for two regions: intratumoral and peritumoral. CA-IX was evaluated for tumor cells only; the percentage of positive cells was estimated on whole sections. For both LAG3 and CA-IX, four different cutoffs defining positivity were considered: ≥1%, ≥5%, ≥25%, and ≥50%.

Univariable Cox regression was used to identify IHC factors associated with OS. Follow-up was censored at 60 months in all analyses. Variables were analyzed as linear continuous measures and dichotomized according to prestated thresholds. Survival by category for each marker is displayed in Kaplan–Meier plots. Multivariable Cox regression models assessed whether prognostic value was altered with adjustments for important clinical factors (age, gender, T stage, N stage, p16 status, smoking status, and alcohol consumption). Analysis was conducted as complete case (123/161 patients), and using multiple imputation with chained equations (24 imputations) to impute missing data. All analyses were performed using Stata 15.

Multi-color immunofluorescent staining was carried out on the Leica Bond Rx using an Opal 7-color fIHC kit and the following primary antibody fluorochrome combinations: CD3/Opal 520, CA-IX/Opal 570, PD-L1/Opal 650, and CK [monoclonal mouse anti-human cytokeratin (concentrate) antibody clone AE1/AE3]/Opal 650. Staining was imaged with the Vectra Automated Quantitative Pathology Imaging System, and 7–11 representative fields per case were quantitatively analyzed using Inform Advanced Image Analysis Software v 2.3 (PerkinElmer).

TCGA data from 275 HNC cases (detailed in Table 1; Supplementary Table S1) were interrogated for expression of hypoxia- and immune response–related genes using a combination of two previously developed signatures: a hypoxia signature composed of 26 genes (Eustace and colleagues; ref. 22) and an immune response signature containing 28 genes (CIRC; ref. 23); individual genes are detailed in Supplementary Table S2. Bioinformatics analysis, using unsupervised two-dimensional hierarchical clustering, identified distinct HNC patient subgroups with coordinate expression of hypoxia- and immune response–related genes (Fig. 1A). Three subgroups were classified as follows: (i) hypoxia^high/immune^low, comprised almost exclusively of HPV-negative cases, (ii) hypoxia^low/immune^high, enriched for HPV-positive cases, and (iii) a mixed subgroup, containing both HPV-positive and -negative cases (Fig. 1A; gene ID is provided in Supplementary Fig. S1). Kaplan–Meier analysis revealed that the three subgroups have different survival profiles, with the hypoxia^low/immune^high subgroup having the best prognosis, corresponding to a 5-year OS of 71.4% (Fig. 1B). The hypoxia^high/immune^low and mixed subgroups have similar profiles, with 5-year OS figures of 50.8% and 48.6%, respectively (P = 0.0015, HPV-adjusted P = 0.0005; adjustment for other standard prognostic variables is considered in Supplementary Table S3). Of note, the hypoxia and immune gene signatures were not independently prognostic (Supplementary Fig. S2). Furthermore, consideration of HPV-negative cases alone, within the predefined hypoxia-immune subgroups, maintained comparable survival profiles (P = 0.001).

We then sought to validate the hypoxia-immune signature, first in silico using data for 134 HNC cases comprising the previously published Chicago cohort (ref. 7; detailed in Table 1; Supplementary Table S1). Bioinformatics analysis identified comparable patient subgroups with those described above (hypoxia^high/immune^low, hypoxia^low/immune^high, and mixed; Fig. 1C). As described previously, the hypoxia^high/immune^low subgroup was predominantly HPV negative. Although HPV-positive cases were again overrepresented in the hypoxia^low/immune^high subgroup, approximately one third of cases in this group were HPV negative. Likewise, the mixed class comprised both HPV-positive and -negative cases, suggesting that the hypoxia-immune signature classifies beyond HPV status. Correlation of gene expression profiles with survival confirmed that subgroups have distinct outcomes. Here again, the hypoxia^low/immune^high subgroup have the best prognosis, 5-year OS 76.2%, whereas the hypoxia^high/immune^low and mixed classes have inferior outcomes with 5-year OS of 65.1% and 60.0%, respectively, (P = 0.012; HPV-adjusted P = 0.0025).

The second validation cohort comprised 163 cases of oropharyngeal squamous cell carcinoma; consistent with their anatomic location a higher percentage of cases were HPV positive (Table 1; Supplementary Table S1: Correlate; flow of patient samples through study is summarized in Supplementary Fig. S3). FFPE tissue samples were available from this cohort enabling comparative gene and protein expression analysis. Gene expression was quantified using Illumina TruSeq Targeted RNA technology. Quality control measures, based on a minimum number of mapped reads of 70,000, reduced the number of informative cases to 123; the heatmap for these and the 12 normal control samples is shown in Fig. 1E. As described previously, hierarchical clustering identified three patient subgroups, two with coordinate expression of hypoxia and immune signature genes. Normal controls clustered with the hypoxia^low/immune^high subgroup, where p16-positive cases were also overrepresented. Conversely, the small hypoxia^high/immune^low subgroup was composed of a majority of p16-negative cases. The mixed subgroup again included both HPV-positive and -negative cases, supporting the proposal that the hypoxia-immune signature classifies beyond HPV status (Fig. 1E). Despite the smaller subgroup size, Kaplan–Meier analysis revealed that the three subgroups have distinct outcome profiles, mirroring the TCGA results (Fig. 1F), with 5-year OS rates of 79.5%, 64.4%, and 55% for the hypoxia^low/immune^high, mixed and hypoxia^high/immune^low subgroups, respectively (P = 0.097; HPV-adjusted P = 0.0003).

Hypoxia may result from rapid tumor growth with inadequate neovascularization. To address the possibility that hypoxia-immune subgroup classification was significantly influenced by tumor volume the distribution of T stage within the three subgroups was examined. A correlation with T stage was not consistently observed (Supplementary Fig. S4).

The analyses described above identified three distinct HNC patient subgroups with coordinate expression of hypoxia- and immune-related gene signatures. However, within each group individual genes may differ in their expression and contribution to the endotype. Subsequent exploratory analyses therefore mapped expression of individual genes within the three subgroups in an attempt to understand their distinct contributions to each endotype. This analysis focused on the TCGA dataset where whole transcriptome data are available, enabling interrogation of genes of interest not included in the Chicago microarray or Correlate targeted sequencing datasets. As illustrated in Fig. 2 and Supplementary Fig. S5A, the pattern of expression for individual hypoxia-related genes was relatively broad with considerable overlap between subgroups, particularly the hypoxia^high/immune^low and mixed classes. However, a subset of genes showed greater differential expression (hypoxia^high/immune^low > mixed > hypoxia^low/immune^high) including CA9, SLC2A1 (GLUT1), and SLC16A1 (MCT1), the products of which are involved in pH regulation and glucose metabolism (Fig. 2, top; Supplementary Fig. S5A). Twenty-three of the 26 genes represented in the hypoxia signature showed significant differential expression between the hypoxia^low/immune^high and hypoxia^high/immune^low or mixed subgroups.

Immune response–related genes described by the CIRC signature were all differentially expressed, albeit to varying extents, according to the expected order: hypoxia^low/immune^high > mixed > hypoxia^high/immune^low (Fig. 2; Supplementary Fig. S5B). The majority of immune genes of interest outside of the CIRC signature also demonstrated a similar pattern of differential expression (Fig. 2; Supplementary Fig. S5C). This was true of genes whose products would be expected to mediate anti- or protumoral roles; for example, markers of cytotoxic and regulatory T-cell function. Genes associated with the myeloid lineage were also expressed at significantly higher levels within the hypoxia^low/immune^high subgroup; although markers of effector/suppressor function [e.g., ARG1 and NOS2 (iNOS)] were not uniformly coexpressed. Other genes of particular interest included those encoding activating receptors (e.g., CD40, CD137, and OX40) or inhibitory receptors (e.g., CTLA-4, LAG3, and PD-1) the ectonucleotidases CD39/CD73 and indoleamine 2,3-dioxygenase (IDO); as therapies targeting these markers are currently in clinical practice, trials, or in development. IDO1 is one of six genes comprising an “IFNγ signature” that has shown predictive value for response to anti-PD-1 mAb therapy in HNC and other cancers (19). The other five genes, all part of the CIRC signature, display a similar expression pattern: hypoxia^low/immune^high >> mixed > hypoxia^high/immune^low; these include IFNG itself and the chemokines CXCL9/CXCL10, which may be important for T-cell homing to tumors. Finally, EGFR was significantly overexpressed in the hypoxia^high/immune^low and mixed subgroups relative to the hypoxia^low/immune^high group.

While gene signatures potentially better reflect the complex interactions between immune cells and the hypoxic tumor microenvironment, simplified IHC signatures may be more clinically applicable. Furthermore, gene expression may not correlate with end protein expression due to posttranscriptional or translational control and protein degradation. Finally, gene signatures neither discriminate between expression within tumor or stromal compartments, nor identify cell type–specific expression. Therefore, we next examined whether an IHC signature comprising TILs, CD3, PD-L1, and LAG3 as immune-response markers and CA-IX as a marker of hypoxia could stratify patients into the corresponding subgroups with different prognoses.

Results for individual markers are summarized in Supplementary Table S4; representative images showing the different patterns of expression and their association with OS are shown in Fig. 3 and Supplementary Fig. S6, with the corresponding statistical analyses in Supplementary Table S5. In agreement with previous studies, high TIL status (assessed on H & E staining) correlated with p16 status (3), and patients with the highest TIL numbers had a significantly better OS [HR 0.27; 95% confidence interval (CI), 0.12–0.66; P = 0.004 on univariate analysis]. Similarly, high numbers of T cells (measured by CD3 staining) correlated with improved OS [combined CD3: HR, 0.39; 95% CI, 0.2–0.77; P = 0.007 on univariate analysis). Because the precise location of T cells within the tumor microenvironment may impact on their prognostic significance (29), T-cell frequencies were independently evaluated for three locations: within tumor nests, within the peritumoral stroma, and at the tumor margin. An elevated frequency of intratumoral T cells showed the greatest prognostic value (HR, 0.33; 95% CI, 0.16–0.67; P = 0.002) and this was independent of clinical factors and p16 status [adjusted HR (aHR), 0.43; 95% CI, 0.19–0.97; P = 0.043; Fig. 3; Supplementary Table S5).

The prognostic impact of PD-L1 expression was dependent on both the cell types considered and the cutoff selected (Fig. 3C and H; Supplementary Fig. S6). Higher expression on tumor cells correlated with improved survival using one of the four tested cutoffs (25%; aHR, 0.33; 95% CI, 0.13–0.89; P = 0.025). Similarly, a trend for improved survival was observed for higher expression on intra-/peritumoral immune cells (HR, 0.45; 95% CI, 0.19–1.05; P = 0.066). In addition, the presence of a PD-L1–positive immune cordon around tumor nests was associated with better OS (aHR, 0.19; 95% CI, 0.05–0.68; P = 0.011). Using a CPS, which considers PD-L1 expression on both tumor and inflammatory cells, strengthened prediction of OS (Supplementary Table S5). Thus CPS based on ≥1%, ≥5%, and ≥25% cutoffs all predicted OS on multivariate analysis; only the ≥50% cutoff did not.

Results for the checkpoint receptor, LAG3, expressed on inflammatory cells showed a trend for improved survival with higher expression (Fig. 3D and I), particularly when these cells were located within tumor nests (25% cutoff: HR, 0.34; 95% CI, 0.12–0.96; P = 0.042). However, these were not significant when adjusted for other tumor factors.

Two-thirds of tumors contained hypoxic regions as measured by CA-IX expression (≥10% cutoff). There are conflicting reports in the literature regarding the prognostic value of CA-IX expression in HNC (reviewed in ref. 30). Here, CA-IX expression was not prognostic in the complete case analysis (Fig. 3J), although there was a trend for improved survival in the subset of patients with very low or absent expression when multiple imputation analysis was used (HR 3.89, 95% CI, 0.95–16.00, P = 0.060; aHR 9.02, 95% CI, 1.96–41.59, P = 0.005; Supplementary Fig. S6).

Since both RNASeq and IHC data were available for the majority of samples from the Correlate cohort, we compared normalized gene versus protein expression levels and found a significant positive correlation for all examined hypoxia or immune markers (Supplementary Fig. S7). Subsequent analyses evaluated levels of IHC marker expression within the three heatmap gene cluster–defined subgroups (hypoxia^high/immune^low, hypoxia^low/immune^high, and mixed). All immune-related markers showed significantly higher levels of protein expression in the hypoxia^low/immune^high versus hypoxia^high/immune^low subgroups; expression levels for the mixed subgroup were intermediate. As expected, CA-IX protein expression was significantly decreased in the hypoxia^low/immune^high versus hypoxia^high/immune^low subgroups and intermediate in the mixed group (Fig. 4A).

It has been suggested that cancers, including HNC, can be assigned one of three immune profiles, based on the frequency and distribution of T cells within the tumor and microenvironment (25). Thus “immune deserts” are characterized by minimal T-cell presence, “immune-excluded” tumors are surrounded by abundant peritumoral T cells, but permit limited intratumoral penetration, and “inflamed tumors” contain relatively high numbers of T cells in both compartments—here T-cell function is apparently compromised. To explore a possible relationship between tumor immune profile and the three gene cluster–defined subgroups as above, all cases within the Correlate cohort were assigned an immune phenotype (desert/excluded/inflamed) based on the magnitude and pattern of CD3 staining. As illustrated in Fig. 4B, hypoxia^high/immune^low tumors were overrepresented in the “immune desert” category and significantly reduced in the inflamed subgroup. The reverse pattern was observed for hypoxia^low/immune^high tumors; while tumors classified as mixed on the basis of hypoxia-immune gene expression were more evenly distributed across all three immune profiles.

In the final series of experiments, we examined the spatial distribution of different cell populations within the tumor microenvironment, in particular the relationship between regions of tumor hypoxia (as measured by CA-IX expression), immune cell localization and PD-L1 expression. Fifteen cases, representative of the three hypoxia-immune subgroups identified in the gene expression analyses, were stained using Opal Multiplex IHC assays. Representative images are shown in Fig. 5 and Supplementary Fig. S8. Staining for cytokeratin (as a marker for tumor cells) confirmed that CA-IX expression is restricted to tumor cells (Supplementary Fig. S8A). A subset of tumors with high levels of CA-IX expression displayed an “immune desert” phenotype, with few T cells present in either the stromal or tumoral compartments (Fig. 5A; hypoxia^high/immune^low). In other cases, abundant T cells were observed within the peritumoral stroma, but minimal numbers appeared to have infiltrated tumor nests, consistent with an “immune excluded” phenotype (Fig. 5B; mixed; Supplementary Fig. S8B and S8C). The final category comprised “inflamed” tumors with abundant T cells present in the peritumoral stroma and significant intratumoral infiltration, particularly in the absence of hypoxia as indicated by no/minimal CA-IX expression (Fig. 5C, hypoxia^low/immune^high). Overall, we observed a significant inverse correlation between CA-IX positivity and the prevalence of CD3⁺ T cells within tumor nests (Fig. 5D, left; r = −0.5464; P = 0.0377,) but not within the stromal compartment alone (Fig. 5D, middle; r = −0.2393; P = 0.3982) or the combined tumor environment (Fig. 5D, right; r = −0.4571; P = 0.0889), consistent with hypoxia-mediated inhibition of T-cell migration and/or survival within tumor nests. No significant correlation was observed between tumoral expression of CA-IX and PD-L1 (Supplementary Fig. S8D), which is unsurprising given that hypoxia (via HIF1α) is only one of multiple potential pathways leading to PD-L1 expression. Excepting “constitutive” cases, strong tumoral PD-L1 expression was restricted to cells at the stromal interface, found in close proximity to CD3-positive T cells (Fig. 5B and C, yellow triangles). The positive correlation between stromal PD-L1 and CD3 T cells (Fig. 5E) is consistent with a dominant role for IFNγ-mediated expression.

Recent studies have advanced our understanding of the biological and molecular characteristics of HNC (4–7). However, there remains an unmet need for clinically relevant classifications to guide therapy selection. Key features associated with the biologically distinct subtypes include immune landscape, expression of EGFR/HER, and hypoxia (7, 31). Given the strong evidential basis for efficacy of hypoxia modification (reviewed in ref. 18) and immune checkpoint inhibitors (19, 20) in the treatment of HNC, we explored a combination of hypoxia and immune status classifiers. Use of the hypoxia-immune signature in an unsupervised cross-cohort manner identified three distinct HNC subtypes. Subgroup classification correlated with survival, with the hypoxia^low/immune^high subgroup having the best prognosis. The hypoxia^high/immune^low and mixed subgroups have similar survival profiles; however, they have distinct transcriptional patterns and immune profiles, consistent with activated pathways that could be targeted by different therapeutic interventions, hence the importance of separating them.

Given the widespread use of IHC for diagnostic purposes in routine pathology practice, we evaluated whether a simplified IHC-based signature could substitute for the developed hypoxia-immune gene expression signature. Expression of individual immune response markers, notably intratumoral CD3⁺ T cells and PD-L1 expression (assessed using a CPS) correlated with superior outcome, but did not improve the prognostic value of clinical factors alone. CA-IX, as a marker of hypoxia, did not have prognostic power.

Finally, consistent with hypoxia-driven immune suppression, multiplex IHC staining identified an inverse relationship between CA-IX and infiltration of CD3 T cells into tumor regions. The observed spatial distribution of hypoxia and immune response markers within the tumor microenvironment correlated strongly with our classification based on gene expression, confirming its validity.

An inverse correlation between hypoxia-related gene expression and antitumoral immune responses is consistent with previous studies in HNC (7, 31). As mentioned above, hypoxia drives immune suppression via multiple mechanisms. It may be achieved through HIF1α-mediated activation of inhibitory pathways, including CD47 (32), the adenosine-generating enzymes CD73 and CD39 (33), and immune checkpoint receptors such as PD-L1 (11, 12). Perhaps paradoxically, many inhibitory markers are more highly expressed in “inflamed” tumors (hypoxia^low/immune^high), where they may be indicative of an ongoing “exhausted” antitumor immune response (34), with the potential for reactivation via immune checkpoint blockade. Tumor cell metabolic adaptations to hypoxia, such as increased glucose uptake and lactate production, also act to promote and perpetuate an immunosuppressive microenvironment (9, 10). Because of differences in their metabolic programs these factors would be expected to have greater impact on effector versus regulatory T cells (35). Alternatively, because most hypoxic tumor areas contained minimal T cells (immune desert or excluded profile), immune suppression may be mediated mainly by mechanisms that inhibit T-cell migration into tumors (e.g., VEGF and CXCL12) or promote T-cell death (e.g., FasL; reviewed in 36).

As described here and reported previously (31), “inflamed” (hypoxia^low/immune^high) tumors demonstrate robust myeloid gene expression profiles. The balance of antitumoral versus immunosuppressive function of myeloid cells, including myeloid-derived suppressor cells (MDSCs) and tumor-associated macrophages (TAMs), may be modulated by tumor hypoxia, which favors suppression (8, 11, 13). Mechanistically, hypoxia (predominantly via HIF signaling) can tip the balance of cytokine/chemokine/effector molecule production toward those with protumoral functions (e.g., IL10, VEGF, arginase, and iNOS) and also upregulate expression of T-cell inhibitory ligands (e.g., PD-L1).

Hypoxia can upregulate EGFR expression (37, 38) and promote ligand-independent EGFR signaling (39), both mechanisms potentially augmenting tumor glycolytic metabolism and consequent metabolic competition. EGFR overexpression is observed in a high percentage of HNC, particularly HPV-negatives cases, and is associated with poor prognosis; other mechanisms of upregulation include mutations, gene duplications, and protein stabilization (reviewed in ref. 40). As previously mentioned, EGFR signaling may lead to hypoxia-independent stabilization of HIF1α and consequent upregulation of glycolytic metabolism (15, 16); providing an alternative (not mutually exclusive) explanation for low intratumoral T-cell infiltration (31). However, it is noteworthy that hypoxic modification appears most beneficial for patients with HPV-negative tumors (41, 42), consistent with an important role for hypoxia-mediated HIF1α upregulation. Furthermore, it has been reported that oropharyngeal tumors, especially HPV-positive tumors are less hypoxic (43).

T-cell infiltration of tumors (especially by CD8⁺ cytotoxic T lymphocytes) has been associated with a favorable prognosis in several tumor types, including HNC (3, 24, 29, 44). Likewise, PD-L1 expression, particularly within immune cells in the tumor microenvironment, can have both prognostic and predictive importance (26, 27). In our study, higher numbers of TILs and CD3⁺ T cells, particularly within tumor nests, correlated with improved survival. PD-L1 expression was also linked with better prognosis (dependent on cutoff), especially when both tumor and immune cell expression were considered as a CPS. The differential prognostic impact of intermediate (CPS > 1% and <50%) versus high (≥50%) PD-L1 expression might reflect alternative mechanisms controlling expression. Intermediate expression (found on immune cells and tumor cells at the stromal interface) may be indicative of adaptive immune resistance in the face of an ongoing antitumoral immune response (45, 46); whereas high levels are associated with tumor intrinsic expression (e.g., loss of PTEN or EGFR activation; refs. 47, 48).

Previous studies in diverse tumor settings have correlated response to immune checkpoint blockade (particularly involving the PD-1/PD-L1 axis) with “cancer immune phenotype” (25).Thus “inflamed” tumors exhibit higher response rates compared with those having “immune desert” or “immune excluded” profiles, although even in the former cases a response is not assured. Our combined analysis of gene and protein expression in HNC suggests that hypoxia may be an important factor distinguishing hot (inflamed) and cold (immune desert or excluded) tumors. In agreement with this, it has recently been reported that targeted hypoxia reduction restores intratumoral T-cell infiltration in a mouse model of prostate cancer (49).

We recognize certain limitations in our study, including the potential for bias due to different tumor:stroma ratios within samples; for example, TCGA samples comprise >80% tumor cells, therefore reducing analysis of the stromal compartment. Minor intercohort variations in gene clustering and prognostic value of the classifier (as measured by 5-year OS) may reflect (i) different RNA quantification platforms (microarray vs. whole-transcriptome sequencing vs. targeted sequencing), (ii) fresh-frozen versus FFPE tissue, and (iii) differences in case mix (anatomic site, TNM/clinical stage, smoking, alcohol consumption, and HPV status) or treatment. Of note, the consistent performance of the classifier despite the heterogeneous nature of the study cohorts and assays is an indication of the strength of the signature (50). Currently, subgroup classification is based on gene clustering (three highest order clusters) rather than defined cutoffs. For clinical application, future work will need to define appropriate cutoffs, for example using median values (22) and/or develop a clinically applicable continuous scoring system.

The identification of our hypoxia-immune prognostic classifier for HNC suggests that differential treatment approaches might be required for patient subgroups. Although hypoxic modification (nimorazole) and immune checkpoint therapies (PD-1/PD-L1 inhibitors) have shown single-agent activity in HNC, the response rate is relatively low; for example, only 13%–18% for PD-1 in the recurrent/metastatic setting (19, 20). Our data suggest that combinations or sequential use of hypoxia-modifying and/or immunomodulatory drugs may be beneficial. For example, in patients with a hypoxia^high/immune^low or mixed signature correlating with poor prognosis, treatment with hypoxia modifiers may sensitize to chemo/radiotherapy, and EGFR inhibitors or agents targeting MDSCs/TAMs may reverse the immunosuppressive environment. Sequential treatment with immune checkpoint inhibitors may then be required to prevent inhibition or dampening of the emerging immune response. Cases that are hypoxia^high/immune^high might benefit from concurrent treatment with hypoxia modifiers and immune oncology agents. Importantly, our data on expression of individual genes also indicate potential for immunotherapy treatment strategies beyond PD-L1, for example combinations of PD-L1 and anti-LAG3 treatments for the hypoxia^low/immune^high subgroup.

In conclusion, we developed and validated a prognostic molecular classifier based on hypoxia and immune status. This classifier may have clinical application to guide the use of hypoxia modification and targeted immunotherapies such as checkpoint inhibitors. We would recommend validation of its prognostic value and assessment of any potential predictive power (e.g., using treatment with nimorazole, EGFR inhibitors, or anti-PD-1/PD-L1) in a prospective setting.

C.J. Bagnall is an employee of AstraZeneca. A.D. Beggs is a consultant/advisory board member for Illumina Inc, Bristol Myers Squibb, Sharp Life Sciences, and Ono Pharmaceuticals Ltd. T.Y. Seiwert reports receiving other commercial research support from Bristol-Myers Squibb, Jounce Therapeutics, and Merck/MSD and is a consultant/advisory board member for Merck/MSD, Bristol-Myers Squibb, AstraZeneca, Aduro, Bayer, Celgene, Innate, Loxo Oncology, and Nanobiotix. H. Mehanna is an employee of and holds ownership interest (including patents) in Warwickshire Head and Neck Clinic, reports receiving commercial research grants from AstraZeneca and GlaxoSmithKline, and speakers bureau honoraria from AstraZeneca, MSD, Sanofi Pasteur, and Merck. No potential conflicts of interest were disclosed by the other authors.

Conception and design: J.M. Brooks, N. Lal, R.J. Spruce, N. Batis, S.P. Lee, G. Middleton, D.A. Tennant, B.E. Willcox, H. Mehanna

Development of methodology: J.M. Brooks, A.N. Menezes, N. Lal, C.J. Bagnall, S.V. von Zeidler, D.A. Tennant, A.D. Beggs, J.-B. Cazier, B.E. Willcox, H. Mehanna

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): J.M. Brooks, N. Lal, H.R. Valentine, J.L. Bryant, M. Hartley, A. Khattri, A.D. Beggs, C.M.L. West, H. Mehanna

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): J.M. Brooks, A.N. Menezes, M. Ibrahim, L. Archer, N. Lal, R. Bao, A. Khattri, K.U.E. Ogbureke, G. Middleton, A.D. Beggs, J. Deeks, J.-B. Cazier, B.E. Willcox, T.Y. Seiwert, H. Mehanna

Writing, review, and/or revision of the manuscript: J.M. Brooks, A.N. Menezes, L. Archer, N. Lal, R.J. Spruce, N. Batis, J.L. Bryant, R. Bao, S.P. Lee, K.U.E. Ogbureke, G. Middleton, A.D. Beggs, J. Deeks, C.M.L. West, J.-B. Cazier, B.E. Willcox, T.Y. Seiwert, H. Mehanna

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): J.M. Brooks, J.L. Bryant, M. Hartley, B. Kaul, G.B. Ryan, H. Mehanna

Study supervision: J. Deeks, J.-B. Cazier, H. Mehanna

This work was funded by Queen Elizabeth Hospital Birmingham Charity and The Get A-head Charitable Trust - 16/5/741; with support from South Egypt Cancer Institute, Assiut (to M. Ibrahim); NIHR Research Methods Fellowship NIHR-RMFI-2016-07-28 (to L. Archer); CRUK Clinical PhD studentship (to N. Lal); the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior - Brasil (CAPES) - Finance Code 001 (to S.V. von Zeidler); Major Centre (to H.R. Valentine); US-UK Fulbright Commission, Distinguished Chair Award (to K.U.E. Ogbureke); The Center for Research Informatics of The University of Chicago Biological Science Division, The Institute for Translational Medicine/Clinical and Translational Award (NIH 5UL1TR002389-02), and The University of Chicago Comprehensive Cancer Center Support Grant (NIH P30CA014599; to R. Bao); the NIHR Birmingham Biomedical Research Centre at the University Hospitals Birmingham NHS Foundation Trust (to J. Deeks); and NIHR Manchester Biomedical Research Centre (to C.M.L. West) and the University of Birmingham.

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.