DNA-Methylome–Based Tumor Hypoxia Classifier Identifies HPV-Negative Head and Neck Cancer Patients at Risk for Locoregional Recurrence after Primary Radiochemotherapy Open Access

Hypoxia elicits numerous negative effects in the tumor-microenvironment interface and plays a pivotal role in the development of therapy resistance (1). Tumor hypoxia is considered a negative prognostic and predictive biomarker for numerous cancer entities (2–5). Head and neck squamous cell carcinoma (HNSCC) is a prototypic cancer for which the negative impact of hypoxia on radiochemotherapy (RCHT) has been well established (6, 7).

In particular, the presence of hypoxia in human papillomavirus (HPV)–negative patients with HNSCC confers poor prognosis. Patients with hypoxic tumors may benefit from stratified therapy intensification (8, 9). However, the robust, feasible, and reliable detection of tumor hypoxia remains a major challenge. From the earlier invasive implantation of probes (10) to non-invasive PET tracers such as FAZA/FMISO/HX4 (11, 12) to more recent MRI-based approaches (13), none of these methods for tumor hypoxia assessment have received regulatory approval. This lack of biomarkers limits the proper stratification of patients for therapy adaptation and may provide a plausible explanation for the failure of most prospective randomized trials to demonstrate a significant benefit for innovative drugs targeting tumor hypoxia. For instance, in the MAESTRO phase III trial for patients with locally advanced unresectable or metastatic pancreatic cancer, the addition of the hypoxia-activated prodrug Evofosfamide to Gemcitabine improved median progression-free survival [hazard ratio (HR), 0.77; P = 0.004] but narrowly missed its primary outcome of overall survival (OS; HR, 0.84; P = 0.059; ref. 4).

A growing body of data indicates that alternative molecular surrogates, such as transcriptional response deciphered by gene expression analysis, can assist clinicians in selecting patients with hypoxic tumors that might benefit most from therapy intensification (5, 14–16). However, transcriptional signatures are vulnerable to several parameters such as surgical clamping time, rapid mRNA degradation by RNAses, and a substantial drop in RNA quality after standard formalin-based fixation and paraffin embedding (FFPE) of the tissue (Fig. 1A). We evaluated the prognostic ability of three previously established gene expression–based signatures (GES) of hypoxia in RNAseq from frozen tissue of The Cancer Genome Atlas (TCGA) HNSCC cohort (ref. 16; Fig. 1B). However, these signatures could not be validated in patients treated with primary RCHT in a multi-centric retrospective cohort from the German Cancer Consortium Radiation Oncology Group (DKTK-ROG; ref. 15). Their prognostic power was limited to smaller tumors (gross tumor volume, GTV<19 cm), where tumor hypoxia was associated with significantly decreased rates of locoregional control (LC) and OS (15). Therefore, we sought to investigate alternative approaches for a more robust detection of tumor hypoxia status.

DNA methylation is an epigenetic marker of a tissue that conserves the cell of origin signatures (17) and has also been shown to reflect cellular perturbations induced by a spectrum of stimuli such as aging (18) or cancer disease states (19, 20). Chronic hypoxia is associated with global changes in DNA methylation in prostate cancer cell lines (21). Genome-wide DNA methylation analysis by hybridization against small oligonucleotide probes immobilized on microarrays is readily available from DNA isolated from FFPE biomaterial (22). Consequently, we sought to develop a DNA methylation–based classifier for tumor hypoxia (hypoxia-M; Fig. 1C). The training set was the TCGA-HNSCC cohort, in which two hypoxic GES were prognostic for outcomes in HPV(⁻) HNSCC. Hypoxia-M was validated in a multi-centric, retrospective primary RCHT cohort of DKTK-ROG. In addition, the prognostic ability of the Hypoxia M classifier was explored in TCGA pan-cancer cohort.

The TCGA-HNSCC cohort was a retrospective cohort consisting of 278 patients with oral cavity, oropharynx, hypopharynx, and larynx tumors (23). HPV status was determined on the basis of genome-wide RNA-sequencing data (RNA-seq, HPVRNA readcounts >1,000), as described previously (23) and was available for all patients. 36 patients with HPV RNA-positive tumors were excluded from this analysis, leaving 242 patient samples. Gene expression data using genome-wide RNA-seq were generated by the TCGA consortium, as previously documented (16, 24). Level 3 RNA-seq data were downloaded, and rlog transformation was performed using the DESeq2 package in R (RRID:SCR_000154; ref. 25). DNA methylation profiling was performed on fresh-frozen tissues using Infinium HumanMethylation 450 (450K) array (23). Functional normalization of DNA methylation was performed using the Minfi package in R (26, 27). The TCGA-HNSCC cohort was used as the training cohort.

The DKTK-primary HNSCC RCHT trial is a multicentric, retrospective primary arm of the DKTK-ROG, consisting of 158 patients with histologically proven HNSCC of the hypopharynx, oropharynx or oral cavity (15). All patients underwent definitive RCHT. Radiotherapy (RT) was delivered to the tumor region and regional lymph nodes (LN) according to the standard protocol, with a boost to the tumor and involved regional LNs. HPV status was determined as previously described (28, 29). Briefly, genomic DNA was extracted from 5-μm FFPE tumor sections using the QIAmp DNA FFPE tissue kit (Qiagen). The LCD-Array HPV 3.5 Kit (CHIPRON, GmbH, Berlin, Germany) was used for analysis of HPV DNA status and genotyping according to the manufacturer's instructions. p16 IHC staining (p16-IHC) was assessed using the CINtec Histology kit (Roche mtm laboratories AG; refs. 28–30). Tumors with strong and diffuse nuclear and cytoplasmic staining in more than 70% of the tumor cells were considered as p16 positive (28–30). Patients with HPV DNA positivity were excluded, leaving 88 samples for validation. DNA methylation profiling was performed using Illumina 450 K microarrays (Illumina) as previously described (8). Functional normalization of DNA methylation was conducted using the minfi package in R (26, 27) and used as a validation cohort. This study was conducted in accordance with recognized ethical guidelines (Declaration of Helsinik). Local ethics Committees at all eight DKTK partner sites granted IRB approval for retrospective data collection and analysis of clinical and biological data.

We previously investigated the impact of two GES in the training cohort from Toustrup and Lendahl (Fig. 1B; ref. 16). Using k-means clustering (KMC), patients were assigned to “GES high” versus “GES low” (16). In addition, three GES from Winter and colleagues (31), Ragnum and colleagues (32), and Buffa and colleagues (33) were used to assign patients continuous hypoxia scores as previously described (5). A cutoff value of >0 was used to classify tumors as having “GES high” versus “GES low” levels of hypoxia.

Hypoxia assignments were also determined for the validation cohort, as described previously (34). Briefly, gene expression for the hypoxia signatures and endogenous control genes (ACTR3, NDFIP1, RPL37A, B2M, GNB2L1, RPL11, and POLR2A) were generated in this cohort using NanoString design code sets and Element Reagents (NanoString Technologies), and two-class KMC was used to assign tumors to high versus low hypoxia classes (34).

OS data were available for download in the training cohort. Data on distant metastasis (DM) or local recurrence (LR) were not available. In the validation cohort, all clinical outcomes were defined as the time (in months) from the first day of RT treatment to the time of the first event. The primary endpoint of this study was LR. The time to progression was defined as the time to LR or DM, whichever occurred first. Data for patients without progression were censored at the time of death or last follow-up. Time to death, LR, and DM were calculated from the start of RT until documented death, LR, and DM, respectively. Review of RT treatment plans and clinical imaging data were performed at the treatment sites to confirm locoregional recurrences arising from the irradiated volume (28, 30).

All statistical analyses were performed using R (35) in Rstudio environment (36).

First, we selected the most hypervariable methylation probe. The methylation values (beta) were used (distribution of values from 0 to 1). Probes with a mean value less than 0.05 or more than 0.95 were filtered. The standard deviation (SD) for each probe was then calculated. The difference (delta, Δ) between the mean and SD was calculated [Δ = absolute value (mean-SD)]. Probes with a high Δ>0.6 or low Δ<0.2 were filtered out. Finally, the 50th percentile for the SD was calculated for the remaining probes. Probes with an SD larger than or equal to the 50th percentile were kept, leaving 55,552 probes for analysis.

The Toustrup and Lendahl signatures had the strongest association with OS in the training cohort (16). We used the consensus of these two signatures to assign the hypoxia status in a binary fashion (“consensus GES high” vs. “else”), as shown in Fig. 1C. Logistic regression was then performed to identify methylation probes associated with hypoxia in the training cohort, TCGA-HNSCC cohort (n = 242) using “glm” with family = “binomial” from the stats package in R (35). Probes with p-FDR ≤0.05 were included in the random forest (RF) classification (Fig. 1C).

Differentially methylated probes (DMP) that were significant in logistic regression were then included in an RF classification in the training set using “randomForest” (37).

In the training cohort, the forest was trained using GES assignment based on the Toustrup/Lendahl combination of signatures for 500 trees (quarternary classification: “consensus high,” “consensus low,” “Toustrup high,” and “toustrup low”). The seed was set to 1 and the out-of-bag error was calculated. An RF classifier (Hypoxia-M) was built using the top 5% DMPs (ranked by decreasing order of importance based on the mean decrease in the Gini Index).

The forest was subsequently used to independently predict the risk groups for patient samples in the validation cohort (DKTK-ROG). In addition, the impact of using hypervariable probes as a starting point (based on the decreasing order of SD, using the 85th percentile of SD of probes as a cutoff value) to train an RF classifier (Hypoxia-M V2) was investigated in Supplementary Methods (Supplementary Figs. S1, S2, and S4).

Kaplan–Meier (KM) survival curves were compared with estimate the difference in clinical outcomes between Hypoxia-M assigned groups using a P value of <0.05 for significance in “survival” (38) and “ggsurv” (39) R packages. For Hypoxia-M–high versus low-risk groups, median time to clinical event and 95% confidence intervals (CI) were calculated in “survival” package.

Multivariate Cox regression analysis was performed to adjust for relevant parameters. A permuted Fisher's test (two-sided) or a Student t test (two-sided) was applied for categorical versus continuous data to detect significant differences in clinical parameters across cohorts (training vs. validation) and intra-cohorts across risk groups (Hypoxia-M high vs. Hypoxia-M low-risk groups). The correlation between the Hypoxia-M assignment and all the GES assignments (in all cohorts) was calculated using the stats package in R. For all correlations, a P value was calculated.

Enumeration of 22 immune cell type subsets based on RNA-seq data for the training cohort was performed using CIBERSORT (40). In the validation cohort, IHC for Cytotoxic T Cells (CD8) and immune checkpoints (PD1 and PD-L1) was conducted as described previously (refs. 15, 41, 42; see Supplementary Methods). Immune cell enrichment was evaluated in conjunction with the methylation classifier using a Fisher's permuted test (two-sided).

Four distinct molecular subtypes based on mRNA expression (classical, atypical, basal, and mesenchymal) have been previously described in the TCGA-HNSCC cohort (23, 43). A Fisher's exact test was used to evaluate the association between hypoxia-M and the RNA subtype values.

A prognostic classifier of hypoxia/immune cell enrichment has been previously described in the TCGA-HNSCC cohort (44). Assignment of patients with TCGA-HNSCC into three subgroups (low hypoxia/high immune group, high hypoxia/low immune group, and intermediate group) was kindly provided by the authors. Fisher's exact test and multivariate regression analysis were used to evaluate the association between this classifier and Hypoxia-M (Supplementary Table S1).

To integrate DMPs with gene expression analysis, the significant DMPs identified on logistic regression were mapped to various gene features (transcription start sites TSS1500 and TSS200, 1st exon, gene body, or Untranslated Regions (UTRs, 5′ or 3′, according to the Illumina Manifest; ref. 45). The methylation values of all probes mapped to the same feature (e.g., all probes mapped to MARCH11 TSS1500) were averaged to obtain the corresponding feature (MARCH11_TSS1500). Pearson's correlation analysis was performed between methylation features and the expression level of corresponding genes. All genes where the level of expression showed a correlation of >0.2 or ←0.2 to their methylation features were retained, and differential gene expression analysis was performed for these genes between Hypoxia-M–high versus Hypoxia-M–low tumors. At an FDR<0.05, we identified genes that were differentially methylated and differentially regulated as functions of hypoxia-M. KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway analysis was performed as previously described to identify the top regulated pathways.

The IGCG/TCGA Pan-Cancer Analysis of Whole Genomes Consortium (in this article abbreviated as TCGA-pancancer) recently reported on divergent mutational processes in hypoxic versus non-hypoxic tumors (5). For 1,188 tumor samples, hypoxia scores were calculated using the GES from Winter, Buffa, and Ragnum, as previously described (5). OS (vital status, days to death, days to last alive) are available in the Supplementary Materials from the corresponding publication (5). Publically available methylation data (Illumina 450 K Human Infinium Methylation) were downloaded from the Genomic Data Commons data portal for 619 tumor samples from this study. We excluded duplicate, locally recurrent, and metastatic tumors from the analysis. Matching survival information and hypoxia scores were available for 521 primary tumors from the publication Supplementary Materials and were used as an independent validation cohort.

Bhandari and colleagues (5) analyzed the genomic correlates of the GES signatures of hypoxia (Ragnum), including mutational signatures, copy-number events, percentage genomic aberrations, and mitochondrial mutations. Integration with these genomic correlates was performed using Pearson's correlation based on the Supplementary Data available from the publication (see Appendix).

DNA-Methylome data were deposited in the ArrayExpress database at EMBL-EBI under accession number (E-MTAB-12431, https://www.ebi.ac.uk/biostudies/studies/E-MTAB-12431).

We previously validated the prognostic impact of the two previously discovered hypoxia-GESs in the TCGA-HNSCC cohort (Fig. 1B; ref. 16). The Toustrup GES, comprising 15 pH-independent genes that were differentially expressed under hypoxic conditions, was developed in human HNSCC cell lines (14). This GES proved to be both prognostic and predictive in a prospective randomized phase III clinical trial from the Danish Head and Neck Cancer Study (DAHANCA-5), whereby patients with lower GES hypoxia scores had better LC. For patients with higher GES hypoxia scores, the administration of nimorazole (a hypoxia modifier) improved LC (46). Lendahl-GES (n = 30 genes) was developed based on a meta-analysis of publicly available microarray gene expression data for seven different cancer cell lines under hypoxia experimental conditions (47).

Although these GES signatures were prognostic in TCGA-HNSCC, they could not be validated in a DKTK-ROG cohort of patients treated with primary RCHT (15). In comparison with RNA biomaterials, DNA methylation is less vulnerable to degradation and cross-linking owing to standard processing and fixation methods. Consequently, we sought to translate the tumor hypoxia assignment from the transcriptome to the methylome.

Clinicopathological characteristics of HPV-negative patients with HNSCC included in the training cohort (TCGA-HNSCC, n = 242) and validation cohort (DKTK-ROG cohort, n = 88) are shown in Fig. 2 (Table 1; Supplementary Table S2). Compared with the validation cohort, patients in the training cohort were older at diagnosis (62 vs. 58 years, P < 0.005), more often female (30% vs. 18%, P < 0.005), predominantly oral cavity tumors (66% vs. 11%, P < 0.005), and had lower TNM stage at diagnosis [T1-T2 tumors (33% vs. 10%, P < 0.005); N0–N1 tumors (52% vs. 13%, P < 0.005), and stage IV tumors (50% vs. 92%, P < 0.005)]. There was no significant difference in the rate of smoking (83% vs. 69%, P = 0.25). The median time to the last follow-up was shorter in the training cohort (29.8 vs. 54.4 months). Finally, there was a trend toward increased mortality in the validation cohort (56% vs. 44%, P < 0.08), although the time to death was similar in both cohorts (30 vs. 27.5 months, respectively, P < 0.1).

A total of 5,378 probes were differentially methylated between tumors with high hypoxic GES (n = 96) and intermediate–low hypoxic GES (n = 146; FDR≤0.05) in the training cohort. To improve future data reproducibility, we dropped probes that were not represented in the newer Infinium 850 Methylation EPIC Beadchip (n = 249 probes) leaving 5129 probes for analysis.

An RF classification was trained on these 5129 methylation probes using a quaternary assignment of hypoxia [“Consensus high” n = 96, “Consensus low” (n = 94), “Toustrup GES high” (n = 25), and “Lendahl GES high” (n = 27)]. The Hypoxia-M classifier was built (Fig. 1C) using the 95th percentile of the most important DMPs (n = 257, based on the mean decrease in the Gini Index). The out-of-bag error was 39.7% in the forest. A total of 128/242 patients were classified as “Consensus High” (53%), 112/242 patients as ”Consensus Low” (46%), and 2 patients (1%) as “Lendahl High.” No patients were classified into “Toustrup High.” Patients in the predicted “Consensus High” group had a significantly increased incidence of death 66/128 (52%) versus 1/2 (50%) in the predicted “Lendahl high group” versus 40/112 (36%) in the predicted “Consensus low” group; (P < 0.015). Given the predominantly binary classification of patients, we retained a binary classification based on RF, and patients predicted to be “Consensus High” were grouped as “Hypoxia-M High” (n = 128), and the rest were “Hypoxia-M Low” (n = 114; Fig. 3A). The KM curves for OS for the hypoxia-M and GES assignments are shown in Fig. 3B.

In multivariate analysis, hypoxia-M was an independent prognosticator of poor OS (HR, 1.66; 95% CI, 1.1–2.48; P = 0.015). There was no difference in OS depending on tumor localization” [Advanced age and nodal status (pN2) were also inversely associated with OS (HR, 1.02; P = 0.047 and HR, 1.6; P = 0.053, respectively; Fig. 3C]. Patients with Hypoxia-M–low tumors were more likely to be female versus male (39% vs. 22%, P < 0.05; Supplementary Table S3). Lendahl-GES and consensus of Toustrup/Lendahl GES were also prognostic in the training cohort (P < 0.015; Fig. 3D). Otherwise, Hypoxia-M–high versus -low tumors were balanced for clinicopathological factors. Similarly, Hypoxia-M V2 classifier was prognostic for OS (P = 0.01) and after multivariate analysis (HR, 1.59; P = 0.026; Supplementary Fig. S1).

The hypoxia-M signature was significantly associated with previously described RNA molecular subtypes (classical, mesenchymal, atypical, and basal; ref. 23) whereby 75% of tumors predicted to be hypoxia-M high belonged either to the classical subtype (45%) or basal subtype (30%), respectively (P < 10⁻¹¹). In comparison, 69% of the tumors predicted to be Hypoxia-M low belonged to either the mesenchymal (45%) or atypical subtypes (24%; Fig. 3E; Supplementary Appendix S4).

In addition, hypoxia-M was inversely correlated with the presence of T cells CD8 (r = −0.21), T cells CD4 activated (r = −0.18), and M1 macrophages (r = −0.18) based on CIBERSORT deconvolution, and positively correlated with the presence of eosinophils and mast cells (P < 0.05; Fig. 3E; Supplementary Table S4).

A total of 5,129 methylation probes were differentially methylated between tumors with high and low hypoxia by GES. Collapsing probes that mapped to the same gene and methylation features resulted in 4,313 unique methylation features (Fig. 3F). RNA-seq data for 2,329 corresponding genes were available for correlation.

After filtering gene expression–methylation feature correlations equal to or greater than 0.2 in absolute value, 1,180 genes remained. We performed a differential gene expression analysis between Hypoxia-M–high and Hypoxia-M–Low tumors using “limma” package in R. At an FDR<0.05, 619 genes were differentially expressed between the two groups (Supplementary Appendix S6).

Five out of the 41 genes in the Toustrup and Lendahl GES were among the differentially expressed genes, representing a significant enrichment for GES genes on a two-sided Fisher's exact test (P < 0.015). Eleven of 140 genes from all GES used in this study (Toustrup, Lendahl, Ragnum, Winter, and Buffa) were significantly upregulated in hypoxic tumors, representing a significant enrichment (P < 0.02; NDRG1, BNIP3, TRIP13, TFAP2C, P4HA2, S100A10, SDC1, CDKN3, ENO1, PAWR, and TPD52L2)

In addition, the top 10 upregulated genes in hypoxia-M were GXP2 (direction of methylation in Hypoxia: hypomethylated), NDRG1 (hypomethylated), BNIP3 (hypomethylated), NEFL (hypermethylated), STC2 (hypermethylated), PTHLH (hypomethylated), SHANK2, IGFL3, SLC7A11, EN1 (hypermethylated). The most downregulated genes in hypoxia-M were SELP, CD27, GVIN1, IRF8, MMRN1, FAM56C, PRKCB, and PIK3CG.

KEGG pathway analysis mapped 105 genes to 30 significantly regulated pathways (Fig. 3G; Supplementary Appendix S7). The top downregulated pathways were cell adhesion molecules, with 19 genes mapped to that pathway, followed by Th1 and Th2 cell differentiation. HIF1a pathway was also significantly altered with 10 genes mapped to the pathway, including four genes with increased expression (CDKN1A, ENO1, HK1, and IGFR) and six genes with decreased expression (FLT1, IFNG, IGF1, PIK3CD, PRKCB, and TLR4).

Applying Hypoxia-M to stratify HPV(⁻) patients with HNSCC in the validation cohort from the DKTK-ROG, 36/88 (41%) patients were predicted to be in the Hypoxia-M–high group and 52/88 (59%) patients in the Hypoxia-M–low group (Fig. 4A; Supplementary Table S5). Compared with the hypoxia-low group, patients in the hypoxia-high group had higher rates of disease progression [25/36 (69%) vs. 21/52 (40%), P = 0.002], death [25/36 (69%) vs. 24/52 (46%), P = 0.02], and local recurrence [LR; 22/36 (61%) vs. 15/52(29%), P = 0.0006]. However, there was no difference in the rate of DM [7/36 (19%) vs. 8/52 (15%), P = 0.4]. On multivariate regression, the Hypoxia-M classifier remained prognostic for OS, LR, and disease progression after adjusting for age, smoking status, GTV, and anatomical site (HR range, 2.34–4.31; P < 0.03). GTV was independently associated with poor OS, LR, and DM (HR, 1.01; P < 0.02; Fig. 4B and C). In addition, Hypoxia-M V2 was prognostic for OS, LR, and disease progression but not DM (Supplementary Fig. S2).

IHC data available for DKTK-ROG showed enrichment for CD8 infiltrating lymphocytes in Hypoxia-M–Low tumors (73% vs. 45% of Hypoxia-M+, P = 0.029; Supplementary Table S6). This effect was particularly significant at the periphery of tumors, with 71% of Hypoxia-M–low tumors staining for CD8 versus 35% of Hypoxia-M–high tumors (P = 0.0035). There was no difference in PD-L1 staining between the two groups (32% vs. 29%, P < 0.9; Fig. 4D; Supplementary Appendix S8). Total CD8-score was prognostic for LR, OS and progression but not DM (Supplementary Fig. S3).

The demographic characteristics of the TCGA-Pancancer cohort are shown in Supplementary Table S7. Within the TCGA-Pancancer cohort (n = 521), 60 patients (12%) were classified in the hypoxia-M–high group and 461 (88%) in the Hypoxia-M–low group (Fig. 5A). Patients with high hypoxia-M tumors had a significantly higher probability of death [27/60 (45%) vs. 117/461 (25%), P < 0.00002; Fig. 5B]. On multivariate regression analysis, hypoxia-M was associated with worse OS (HR, 1.83; P = 0.04; Fig. 5C). Hypoxia-M V2 was also prognostic (HR, 1.81; P = 0.03; Supplementary Fig. S4).

In addition, the tumor stage was independently associated with worse OS (HR, 1.91; P = 0.03). After stratification by tumor histology, renal cell carcinomas were associated with better outcomes (HR, 0.35; P < 0.01).

There was a significant enrichment for hypoxic tumors in bladder cancer (14/23; 61%), cervical cancer (10/20; 50%), head and neck cancers (16/42; 38%), and lung cancer (10/48; 21%). The Hypoxia-M classifier was prognostic globally (HR, 1.83; P = 0.04) and stratified by subsite for squamous carcinoma and bladder cancer. In this subgroup analysis, 46/104 (44%) patients were classified as having Hypoxia-M tumors with a worsened HR for OS (HR, 2.23; 95% CI, 1.02–4.86; P = 0.04). Furthermore, the utility of patient assignments by hypoxia signatures across tumor histology in the TCGA-Pancancer cohort was compared (Fig. 5D). GES has limited utility for patient stratification in squamous carcinomas and bladder cancer, with >90% of patients classified as having high levels of hypoxia according to Buffa and colleagues (33). In addition, Winter, Buffa, and Ragnum GES were strongly prognostic in CNS tumors (glioblastoma and lower-grade gliomas), as well as renal cell carcinoma (clear cell, chromophobe, and papillary).

Hypoxia-M was prognostic in squamous carcinomas and bladder cancer, with 0% of CNS and renal cell tumors being hypoxic, limiting its usefulness in these tumor-subtypes. Furthermore, no GES was a prognostic factor for adenocarcinomas.

Next, we evaluated the correlation between the assignment by hypoxia-M and the assignment by hypoxia GES (Fig. 5E).

In the TCGA-HNSCC cohort, Hypoxia-M correlated with the different GES (the strongest correlation with Toustrup/Lendahl GES (0.47, P < 0.00001), followed by Lendahl GES (0.46, P < 0.00001), and Toustrup GES (0.38, P < 0.00001; Fig. 5E). Hypoxia-M also correlated significantly with a higher T stage (T3) and was inversely correlated with male sex. Lendahl and Toustrup GES showed different predilections for anatomical sites, with assignment according to Lendahl GES being significantly correlated with laryngeal localization and anticorrelated with oral cavity location, and Toustrup GES was correlated with oropharyngeal location. There was no correlation between hypoxia-M and anatomical site or age.

In the validation cohort, hypoxia-M correlated with the Toustrup/Lendal GES (0.36, P < 0.0001), Lendahl GES alone (0.43, P < 0.00001), and Toustrup GES (0.24, P = 0.04). There were no correlations between T, N or other clinical covariates. The Lendahl and Toustrup GES correlated significantly with stage (0.49) and N status (0.46; P < 0.05).

Within the TCGA-pancancer cohort, hypoxia-M was significantly correlated with Winter (0.48), Buffa (0.39), and Ragnum (0.35) GES. Hypoxia-M, Buffa, Ragnum, and Winter GES were also significantly correlated with T stage (r = 0.11–0.15), N-stage (r = 0.17–0.36; P < 0.05), and tumor stage (r = 0.19–0.28; P < 0.05). All signatures were significantly correlated with bladder tumors, HNSCC tumors, and cervical cancers. Renal cell carcinomas, thyroid tumors, and hepatocellular carcinomas were inversely correlated with hypoxia signatures.

This study reports a novel approach for the establishment of an epigenetic DNA-methylome–based classifier for tumor hypoxia, hypoxia-M. Assignment by two prognostic GES of hypoxia in patients with HNSCC with HPV-negative tumors served as the basis for training hypoxia-M as an independent prognostic biomarker of OS. Hypoxia-M was validated in HPV-negative patients with HNSCC treated with primary RCHT in a multicentric DKTK-ROG trial. Patients in the hypoxia-M high group had significantly higher risks of LR (HR, 4.13), disease progression, and death. GES-based biomarkers were shown to be prognostic in smaller-sized tumors in the DKTK-ROG trial, that is, patients with GTV <19 cm (15). In contrast, hypoxia-M was prognostically independent of tumor size. There was no difference in GTV size between Hypoxia-M–high and –low tumors.

The discrepancy between GES and hypoxia-M may be partly explained by the differences in the biomaterials. The training cohort was selected for Hypoxia-M classifier discovery owing to the availability of matched transcriptome (RNA-seq) and DNA methylome (Illumina HumanMethylation arrays) data derived from fresh-frozen tissue in TCGA. In contrast with the global normalization algorithms applied to transcriptome studies, gene expression analysis of the selected hypoxia GES panels in the validation cohort required normalization against a limited set of endogenous control genes. In patients treated with primary RCHT, such as those in the DKTK validation cohort, tumor material stems primarily from biopsies and may often be limited in quantity. Moreover, RNA was extracted from FFPE tissues with known biomaterial quality limitations, such as degradation and cross-linking after fixation, which may negatively affect the biological signal of the GES approach. FFPE-extracted DNA biomaterial was robustly assayable for methylome analysis in relatively low amounts (100–300 ng) as demonstrated in this study. In addition, DNA methylation reprogramming takes a longer time to occur, and may therefore be reminiscent of more chronic tumor hypoxia. By contrast, transient alterations in tissue oxygen levels due to ischemia/reperfusion during surgery/biopsy may lead to rapid perturbations in gene expression.

Consistent with emerging evidence about the interplay between hypoxia and immune cell depletion in cancer (44, 48), we found that Hypoxia-M high tumors had lower T-cell infiltration (CD8 and CD4) and reduced M1 Macrophages in the TCGA-HNSCC cohort. This finding was validated by IHC staining for CD8 T cells in the DKTK-ROG cohort, whereby hypoxia-M high tumors had lower CD8 T-cell staining. This is consistent with the reported characterization of cellular responses to hypoxia in solid tumors (49, 50), whereby hypoxic niches are infiltrated by immunosuppressive cells such as T regulatory cells (51), myeloid-derived suppressor cells (52), and tumor-associated M2 macrophages (50). In addition, T-cell activation is altered under hypoxic conditions (53). Further corroborating this interplay is the recent development of a GES of both hypoxia and immune cell infiltration in the TCGA-HNSCC cohort (44). This GES assigned patients into hypoxia low/immune cell high, hypoxia high/immune cell low, and mixed groups, and was prognostic of outcome.

Finally, we explored the prognostic ability of this classifier in a TCGA pan-cancer cohort in which gene signatures of hypoxias were used to identify hypoxia-associated mutational processes. In the TCGA-pancancer cohort, the classifier was robust and prognostic for squamous cell carcinoma and bladder cancer. Bladder cancer tumors were identified as the most hypoxic tumors (14/23, 61%), followed by cervical cancer (10/20, 50%), and HNSCC (36%). These results are promising because hypoxia modification is associated with a 13% benefit in OS in bladder cancer (54) and warrants further exploration.

Hypoxia-M consisted of 257 methylation probes and mapped to 218 genes/transcripts. There were prominent genes among the DMPs that were known to be associated with hypoxia (STC2, SLC6A3, HIGD1B, VGLL4, PTPRJ, ANKDD1A, and NCOA2), and treatment resistance (EN1, SPINK7, HOXA11-AS1, and TP53TG1), and a number of immune-related genes, such as IL19, ZAP70, DSE, PRIMA1, CLEC12A and CLEC16A, NOTCH3 were differentially methylated. There were also many genes involved in metabolic signaling, including ACADL, OAT, PDX1, TBC1D4, ACLY, ADPGK, ALDH6A1, CBR4, genes involved in epithelial–mesenchymal transition (EMT) pathways (ELMO1, ITGA8, STK32A, ZAK, PALLD, MNX1, NR2F2, UBE29, FNDC3B, and NR2F1), and genes associated with glutathione synthesis (GCLM, CHAC2, and GLRX3).

Notably, Stanniocalcin-2 (STC2, methylation probe: cg03846076) was the most important DMP in the random forest classifier (by mean decrease in Gini). STC2 was also among the top 10 upregulated genes under hypoxia in differential gene expression analysis. Although the function of STC2 has not been precisely elucidated, HIF-1alpha binds to the promoter region of STC2 to increase its expression in the hypoxia (55, 56). Stanniocalcin-2 is also a glycolytic gene that is postulated to be associated with the Warburg effect in osteosarcoma (57). It has also been shown to be a mediator of metastasis in HNSCC (58), which is resistant to EGFR-TKI inhibition in non–small lung cancer (59).

Engrailed-1 EN1 (cg09601912) was differentially methylated and upregulated under hypoxic conditions. EN1 is a transcription factor associated with the proliferation and migration of triple-negative breast cancer (60). EN1 amplifies TGFβ signaling and mediates its profibrotic effects in myofibroblasts (60). Glutoredoxin 3 (GLRX3) or PICOT is an oxireductase that protects cells from DNA damage by attenuating replication-induced instability (61, 62).

Among the silenced genes, immune-related genes such as PIK3CG (cg15881332, methylated), which is involved in leukocyte chemotaxis (63), ZAP70 (cg17913347, demethylated), a tyrosine kinase that initiates T-cell response (64), whose methylation is a prognostic factor in chronic lymphocytic leukemia (65) and ELMO1 (methylated, cg06159352), a docking protein that is involved in cell migration but also in promoting inflammation (66). The apoptosis genes CASP7 (cg07462448, methylated), cell adhesion molecules ITGA8 (cg19944763, cg07073391, methylated), and SLIT3 (cg21122370, methylated). Silencing SLIT3 has been shown to increase EMT and cell migration in hepatocellular carcinoma (67, 68).

In addition to evaluating Hypoxia-M classifier probes (n = 257), we evaluated all differentially expressed genes between Hypoxia-M–High versus hypoxia-M–low tumors (n = 5129 probes). The top upregulated genes included GPX2, also a glutathione peroxidase that plays a role in protecting cells against oxidative damage. In one study, the expression of GPX2 abrogated hypoxia-induced DNA damage in colorectal cancer and the use of tirpazamine–re-sensitized DNA to reactive oxygen species (69). NDRG1, BNIP3 are hypoxia signature genes (Toustrup and Lendahl GES; refs. 3, 47). Parathyroid hormone–related Peptide (PTHrP) has been shown to promote EMT and upregulate levels of HIF1-alpha, VEGF, and MMP-9 in cell lines of colorectal cancer (70).

The top downregulated genes in hypoxia-M were mostly immune-related. This includes SELP (activates inflammatory cytokines and leukocytes), CD27 (a regulator of B-cell activation), GVIN1 (induced by IFN), IRF8 (expression of tumor-associated macrophages), MMRN1 (activation of platelets), CECR1 or ADA2 (adenosine deaminase), which are also involved in immune system activation, and PRKCB (involved in B-cell activation). Together, these alterations in gene methylation and gene expression suggest that Hypoxia-M–high tumors have a more aggressive phenotype, which is in line with their worsened clinical outcomes.

Methylation-based classifiers are attractive. In this study, Hypoxia-M was consistent in identifying hypoxic tumors in all cohorts, irrespective of the tissue material (fresh-frozen or FFPE). The RF showed a stable correlation with the GES (range, 0.34–0.46) that was measured in different tissues, different gene expression platforms (RNAseq vs. Nanostring), and in different tumor types.

Novel robust and reliable biomarkers for stratification and therapy adaptation are urgently needed, especially for HPV-negative patients with HNSCC. Increased tumor hypoxia has been associated with resistance to RT, decreased local control, and poorer prognosis in HNSCC (71) and the use of hypoxia-modifying strategies has been postulated to improve OS in patients with HNSCC (71, 72). However, in clinical practice, tumor hypoxia has no influence on decision-making or patient stratification (73). This is mainly due to the lack of approved biomarkers that are robust, reliable, affordable, reliable, and easy to use. The discovery of hypoxia-M has the potential to change this situation. It may serve clinicians to identify patients who might benefit most from particle RT with high linear energy transfer carbon or heavier ions that are less dependent on the oxygen enhancement ratio and/or assist in stratifying suitable patients for treatment with hypoxia modifiers (71, 74, 75). Validation of Hypoxia-M is envisioned in the currently recruiting prospective HNprädBio study of DKTK-ROG in a cohort of patients treated with primary RCHT (www.clinicaltrials.gov, NCT02059668).

C. Rödel reports grants from DKTK (German Cancer Consortium) during the conduct of the study. A. Linge reports involvement in a publicly funded (German Federal Ministry of Education and Research) project with the companies Medipan (2019–2022), Attomol GmbH (2019–2022), GA Generic Assays GmbH (2019–2022), Gesellschaft für medizinische und wissenschaftliche genetische Analysen (2019–2022), Lipotype GmbH (2019–2022), and PolyAn GmbH (2019–2022), none of which, to the best of A. Linge’s knowledge, were involved in the preparation of this article. I. Tinhofer reports personal fees from Merck Serono outside the submitted work. M. Krause reports other support from Attomol GmbH, GA Generic Assay, PolyAn GmbH, Gesellschaft für medizinische und wissenschaftliche genetische Analysen Dresden mbH (Mewigen), Lipotype GmbH, Medipan GmbH, Cybertron GmbH, KDS Radeberger Werkzeugbau GmbH, and Resintec GmbH and grants from NIH and Merck outside the submitted work. M. Stuschke reports grants and personal fees from AstraZeneca and personal fees from Janssen-Cilag, Bristol-Myers Squibb, and Sanofi-Aventis outside the submitted work. S.E. Combs reports personal fees and nonfinancial support from Roche, AstraZeneca, Medac, Dr. Sennewald Medizintechnik, Elekta, Accuray, BMS, Brainlab, Daiichi Sankyo, Icotec AG, Carl Zeiss Meditec AG, and HMG Systems Engineering outside the submitted work. C. Belka reports grants from Government Agency HMGU and Government Agency BMBF during the conduct of the study as well as personal fees from Merck Darmstadt outside the submitted work. A. Stenzinger reports grants and personal fees from Amgen, Bayer, Incyte, BMS, and MSD and personal fees from Astra Zeneca, Takeda, Illumina, and Eli Lilly during the conduct of the study. M. Baumann reports grants from DKFZ during the conduct of the study; in addition, M. Baumann is a CEO and Scientific Chair of the German Cancer Research Center (DKFZ, Heidelberg) and is responsible for collaborations with a large number of companies and institutions worldwide. In this capacity, M. Baumann has signed contracts for research funding and/or collaborations, including commercial transfers, with industry and academia on behalf of M. Baumann’s institute(s) and staff; M. Baumann reports membership on several supervisory boards, advisory boards, and boards of trustees. J. Debus reports grants from Accuray, Elekta, Siemens, Raysearch, and Viewray during the conduct of the study as well as grants from Merck Serono outside the submitted work. A. Abdollahi reports grants and other support from Merck KGaA, FibroGen Inc., Bayer, Roche, and Merck Serono and other support from BMS and BioMedX outside the submitted work. No disclosures were reported by the other authors.

B. Tawk: Conceptualization, resources, data curation, software, formal analysis, funding acquisition, validation, investigation, visualization, methodology, writing–original draft, writing–review and editing. K. Rein: Resources, visualization, project administration. C. Schwager: Resources, data curation, software, formal analysis, validation, investigation, methodology. M. Knoll: Resources, software, formal analysis, methodology. U. Wirkner: Resources, data curation, data aquisition. J. Hörner-Rieber: Resources, data curation, data acquisition. J. Liermann: Resources, data curation, data acquisition. I. Kurth: Resources, data curation, data acquisition. P. Balermpas: Resources, data curation, data acquisition. C. Rödel: Resources, data curation, data acquisition. A. Linge: Resources, data curation. S. Löck: Resources, data curation. F. Lohaus: Resources, data curation. I. Tinhofer: Resources, data curation, data acquisition. M. Krause: Resources, data curation. M. Stuschke: Resources, data curation. A.L. Grosu: Resources, data curation. D. Zips: Resources, data curation. S.E. Combs: Resources, data curation. C. Belka: Resources, data curation. A. Stenzinger: Resources, data curation. C. Herold-Mende: Resources, data curation. M. Baumann: Resources, data curation, data acquisition. P. Schirmacher: Resources, data curation. J. Debus: Resources, data curation, formal analysis, methodology, project administration. A. Abdollahi: Conceptualization, resources, data curation, formal analysis, supervision, funding acquisition, methodology, writing–review and editing.

This study was supported by Zentrum für Personalisierte Medizin (ZPM-Network BW, Project MAX-PRO), the German Cancer Consortium (DKTK), the Helmholtz Cross-Program Initiative Personalized Medicine (iMed) project on “Multi-Scale Integrative Biology of HNSCC,” and intramural funds from the National Center for Tumor Diseases (NCT) Heidelberg Radiation Oncology Program. B. Tawk was supported by a stipend of the German Cancer Research Center (DKFZ) Clinician Scientist Program, supported by the Dieter Morszeck Foundation (https://www.dkfz.de/en/clinicianscientist/index.html). We would like to thank Barbara Schwager and Claudia Rittmüller for their excellent technical assistance, Dr. M. Bewerunge-Hudler at the DKFZ Genomics and Proteomics Core Facility team for providing excellent support for DNA methylation microarray experiments, as well as all members of the DKTK-ROG and related multidisciplinary departments at Heidelberg Core and all DKTK-partner sites who contributed to the success of this trial.

The publication costs of this article were defrayed in part by the payment of publication fees. Therefore, and solely to indicate this fact, this article is hereby marked “advertisement” in accordance with 18 USC section 1734.