Development of a Hypoxia Gene Expression Classifier with Predictive Impact for Hypoxic Modification of Radiotherapy in Head and Neck Cancer

The specific issue of hypoxia-induced radiation resistance has been acknowledged for more than 100 years and methods to describe and overcome the obstacle have occupied researchers in the field of radiotherapy ever since (1–3). Therapeutic initiatives such as hyperbaric oxygen (4), hypoxic radiosensitizers (5, 6), and hypoxic cytotoxins (7) to control tumor hypoxia have been promising, but the question still remains whether we are evaluating and directing the hypoxia-targeted therapy to the right patients.

Techniques to detect relevant tumor hypoxia have continuously evolved and the approaches have included both direct measurements in the tumor by an oxygen electrode (8, 9), infusion and detection of exogenous hypoxia tracers [pimonidazole (10) and various PET-tracers (11–16)], or quantification of endogenous hypoxia markers expressed by the tumor cells under hypoxic conditions (17–19). Despite progression in these fields, no common treatment strategy and identification approach about hypoxia in tumors has yet been implemented in a clinical setting.

In this study, we aimed at generating a method for characterizing the hypoxic status of a tumor based on quantifying the gene expressions of hypoxia-responsive and hypoxia-specific genes within the tumor biopsy, and furthermore generate a method that could improve the ability of individualizing treatment in accordance to such characterization. Thus, the experimental plan involved identifying specific hypoxia-induced genes responding with a significant increase in expression correlating to a radiobiological relevant oxygen level (20) that could additionally identify patients having benefit from hypoxic modification of radiotherapy. Such “hypoxia-regulated genes” have previously been suggested in the literature (21–23), but to our knowledge none of the developed hypoxia gene expression signatures (24–26) have yet shown to be predictive for beneficial treatment strategies and consequently implemented in the clinic.

To identify relevant genes for hypoxia classification, we have previously focused especially on the influence of hypoxia and pH on gene expression in vitro (27). From these studies, 29 genes characterized by being upregulated under hypoxic conditions and furthermore being independent of pH fluctuations were identified as potential classifier genes across a range of human squamous cell carcinoma cell lines. The identified genes were upregulated more than 2-fold under hypoxia in at least 4 of 5 cell lines. In addition, we chose to investigate and quantify the gene CA9, which has frequently been associated to hypoxia in the literature (refs. 28–30; Table 1).

We first validated the in vivo relevance of the in vitro identified genes in a xenograft study by applying the hypoxia radiotracer ¹⁸F-flouroazomycin arabinoside (FAZA) in the detection and isolation of hypoxic and nonhypoxic tumor tissue. Following gene expression quantification, we compared the genetic response in the 2 areas to detect whether an upregulation in the hypoxic area could be registered. To proceed into the clinical setting, we then developed a hypoxia classifier on gene expressions quantified from hypoxia-evaluated human tumor biopsies. This classifier contained 15 of the in vitro identified hypoxia-responsive genes that could best discriminate between “more” and “less” hypoxic human head and neck squamous cell carcinomas (HNSCC; evaluated by oxygen electrode measurements). In an independent data set, where patients were randomized for either placebo or hypoxic modification in combination with radiotherapy, we finally evaluated whether the developed 15-gene hypoxia classifier could identify patients with hypoxic tumors and thereby identify patients that would potentially benefit from hypoxic modification. Our findings indicate that with this strategy, it has been possible to develop a 15-gene hypoxia classifier, which attain both prognostic and predictive impact.

In vitro experiments were based on human squamous cell carcinoma cell lines UTSCC5, UTSCC14, UTSCC15 [oral carcinoma, unknown human papilloma virus (HPV) status, established by Dr. Reidar Grenman, University of Turku, Finland, obtained from Dr. Michael Baumann], FADU_DD (a subline of FaDu, an undifferentiated hypopharyngeal carcinoma, HPV-negative, obtained from Dr. Michael Baumann), and SiHa (uterine cervix carcinoma, HPV-positive, obtained from the American Type Culture Collection). Xenograft tumors for in vivo validation were generated with cell lines of UTSCC33 (oral carcinoma, unknown HPV status, established by Dr. Reidar Grenman, University of Turku, Finland, obtained from Dr. Michael Baumann), FADU_DD, and SiHa. None of the cell lines have been tested and authenticated. The hypoxia classifier was generated from 58 head and neck cancer biopsies archived as formalin-fixed and paraffin-embedded samples (FFPE). The oxygen status of these tumors had previously been evaluated in accordance to the relative number of oxygen electrode measurements less than 2.5 mm Hg in their metastatic lymph nodes (9, 31, 32). To test the prognostic and predictive impact of the developed classifier, we isolated mRNA from 326 accessible archival supraglottic larynx or pharynx tumor samples (FFPE) from the randomized, double-blinded DAHANCA 5 trial (6). Gene expression could be quantified in 323 samples, and these were individually classified as either “more” or “less” hypoxic. In the original study, 414 tumors had been randomized to either placebo or hypoxic modification with nimorazole [4-[2-(5-nitro-1H-imidazol-1-yl)ethyl]morpholine] in conjunction with conventional radiotherapy (62–68 Gy, 2 Gy/fx, 5 fx/wk). Those treated with nimorazole experienced an improved outcome in comparison with those treated with placebo. There was no significant differences in terms of patient and tumor characteristics (N- and T-stage, tumor site, tumor differentiation, age, gender, HPV status) and locoregional tumor failure between the 91 patients where RNA was not accessible and the 323 patients included in the present analysis (data not shown).

The subjected genes were selected from the data from a previously published in vitro study, where the above-mentioned cell lines were exposed to different oxygen concentrations and pH levels (7.5 or 6.3; ref. 27). Gene expression was analyzed with microarray (Affymetrix–Human Genome U133 Plus 2.0 Array). In addition, the expression of the pH-dependent gene CA9 was evaluated. (Fig. 1A).

In autoradiographic studies of xenograft tumors, FAZA was used as exogenous tracer revealing hypoxic tumor areas (<10 mm Hg; ref. 33). Immediately after cryosection of the excised tumor, autoradiography was carried out. This procedure was carried out at −20°C to minimize degradation of mRNA. By demarcating FAZA-positive areas (H) and FAZA-negative areas (N) from the autoradiography a computer-assisted (ImageGauge) 1:1 template was made to describe the hypoxic status inside the tumor. Guided by hematoxylin and eosin (H&E) staining of every fifth section (Fig. 1B), necrotic areas were avoided if possible and the demarcated areas were dissected. Tissue from multiple corresponding sections, but representing each area, was pooled to achieve tissue for quantitative PCR (qPCR) quantification. Every fifth tumor section was left in total (M) before RNA extraction, cDNA preparation, and qPCR. All sections were preserved in RNAlater-ICE at −80°C prior to dissection.

RNA from fresh-frozen tissue was extracted by using RNeasy-kit (Qiagen) according to the manufacturer's instructions. RNA from FFPE tumor biopsies was extracted by using a fully automated, bead-based RNA isolation procedure (34). cDNA was generated by using the High Capacity cDNA Archive kit (Applied Biosystems; ABI) and gene expression was quantified by using qPCR. cDNA based on FFPE samples was preamplified according to the manufacturer's details (TaqMan PreAmp, ABI) before real time qPCR. To detect transcripts of interest, TaqMan Gene Expression assay (ABI) was used for all potential classifier and reference genes (Supplementary Table S1). Reactions were carried out on an ABI Prism 7900 Sequence Detector (ABI) and in duplicate. Seven potential reference genes were analyzed: ACTR3, CHCHD1, NDFIP1, SPCS2, and TMEM85 were selected on the basis of our previous in vitro results (27); RPL37A and TFRC have previously been described (35, 36). The geNORM application available in RealTime Statminer (Intergromics) was used on the training set of 58 HNSCC samples, and 3 optimal reference genes (RPL37A, ACTR3, and NDFIP1) were selected for normalization. Data were analyzed by using RealTime Statminer (Intergromics) and ΔCt values were generated by normalizing to the geometric mean of the 3 reference genes. When analyzing FFPE samples, Ct values more than 35 or with an SD more than 0.3 were dismissed and interpreted as empty wells. Gene expression levels were quantified as 2^−ΔCt and log2-transformed before building the classifier.

Gene expression in the H, N, and M areas of the xenograft tumors were compared with Wilcoxon's signed rank analysis.

The prognostic impact of the 15-gene hypoxia classifier was evaluated on the placebo group of the DAHANCA 5 data set. This subgroup had been treated with radiotherapy alone (6). The predictive impact for hypoxic modification with nimorazole in conjunction with radiotherapy was evaluated by comparing the response to placebo or nimorazole among the groups classified as having “more” or “less” hypoxic tumors, respectively. A multivariable Cox proportional hazard analysis was used to evaluate the importance of independent tumor characteristics and treatment, with respect to the risk of locoregional failure. All patients were observed for at least 5 years or until death. Follow-up was completed as part of the original study. End point in the evaluation of the classifier was cumulative incidence of locoregional tumor failure defined as persistence or recurrence of the disease in the primary tumor (T site) or regional lymph nodes (N site) after radiotherapy and compared at 5 years. Death within 3 years after finalized treatment and without failure was considered a competing risk in the analysis (37). Statistical analyses were carried out with STATA 10 software. All P values are 2-sided with a level of significance at 5%.

To confirm the in vivo hypoxia specificity of the 30 genes, subcutaneous xenograft tumors in nude mice were established from 3 cell lines of human squamous cell carcinoma, of which 2 were used in the in vitro experiments. Before excision of the evolved tumor, mice were administered with the hypoxia tracer FAZA (33) which allowed ex vivo visualization of hypoxic tumor areas by autoradiography. A computer-assisted 1:1 template of demarcated FAZA-positive and FAZA-negative areas was generated from the autoradiography to guide the dissection of hypoxic and nonhypoxic tumor tissue (Fig. 1B and C). Tissue was isolated and gene expression was quantified from each area, respectively. A Wilcoxon's signed rank analysis of the intratumor variability of gene expression measured with quantitative real-time reverse transcriptase PCR showed that all the suggested genes were significantly upregulated in hypoxic tumor areas (H) compared with nonhypoxic tumor areas (N; Fig. 1D). To mimic the clinical scenario where biopsies may include both well-oxygenated and hypoxic areas, we also compared the gene expression from nonhypoxic tumor areas with the gene expression in total neighbor sections (M) that represented more heterogeneous tissue material. With this, we expected to verify whether the hypoxia-induced upregulation would be traceable in terms of gene expression, irrespective of being quantified from a mix of both hypoxic and nonhypoxic tumor areas. We found all but 3 genes to be significantly upregulated (Fig. 1E), which supported their potential role as hypoxic markers.

Identification of the most informative genes for relevant hypoxia classification in human biopsies was carried out by quantifying gene expression in 58 head and neck cancer biopsies that had previously been hypoxia evaluated and ranked in accordance to oxygen electrode measurements in their metastatic neck nodes (9, 31, 32). In the effort of building the hypoxia classifier we separated and categorized the 58 tumors constituting the training set into a “more” hypoxic group containing the 10 tumors with the highest frequency of low oxygen electrode measurements and a “less” hypoxic group containing the remaining 48 tumors (Fig. 2A). With this split of the hypoxia-ranked tumors, we obtained 2 groups with the largest possible distance between mean gene expression levels among the groups and thereby the greatest discrimination in terms of gene expression. In short, the separation of the tumors into these 2 groups was determined by the ratio of between to within variations (B/W) in expression levels (see Supplementary Data on classification). We obtained the highest B/W ratios when the “more” hypoxic group consisted of the 10 most hypoxic tumors (Fig. 2B). Subsequently, an independent tumor would be classified as belonging to the predefined group (“more” or “less” hypoxic), where the distance from the gene expression level of the independent tumor and to the mean gene expression level of the predefined group was lowest (Supplementary Data on classification; Supplementary Table S2).

We then determined the number of genes to constitute the most optimal classifier. By use of a “leave one out” cross-validation approach (38), each sample from the training set was classified, one by one, as belonging to either of the 2 groups. The classification was based on all samples from the training set, excluding the sample being classified. Each of the 58 samples was then classified with combinations of the 30, 29, 28, and so on best separating genes. We found that a combination of 15 genes was optimal to classify the highest number of tumors into the same group as they were predefined into (Fig. 2C). As the final classifier, we chose the 15 candidate genes (Table 1) that were present with the highest frequency in the 58 “leave one out” classifications with a 15-gene classifier (Fig. 2D). To verify the prognostic relevance of the developed 15-gene hypoxia classifier, we categorized the 156 HNSCC biopsies from patients in the placebo group of the DAHANCA 5 trial as either “more” or “less” hypoxic. These patients had been treated with conventional radiotherapy and without hypoxic modification. The prognostic value is illustrated in Figure 2E, where patients categorized as having “more” hypoxic tumors suffered a significantly higher cumulative incidence of locoregional failure at 5 years as compared with patients categorized as “less” hypoxic (79% vs. 54%, P = 0.001).

To test the classifiers potential for predicting which tumors could benefit from hypoxic modification of radiotherapy with nitroimidazoles, we applied the 15-gene classifier on 323 HNSCC biopsies from patients enrolled in the DAHANCA 5 trial. These patients had been randomized to either placebo or hypoxic modification with nimorazole in conjunction with conventional radiotherapy (6). Overall, the hypoxia classifier categorized 114 tumors (35%) as “more” hypoxic and 209 tumors (65%) as “less” hypoxic. In the group classified as “more” hypoxic (Fig. 3A), those treated with nimorazole in conjunction with radiotherapy carried a significantly reduced cumulative incidence of locoregional tumor failure after 5 years when compared with those treated with placebo and radiotherapy (46% vs. 79%, P = 0.0001). In the group classified as “less” hypoxic (Fig. 3B), we observed a uniform prognosis in terms of cumulative incidence of locoregional tumor failure (46% vs. 54%, P = 0.28), irrespective of whether the patient had been treated with nimorazole or placebo in conjunction with the radiotherapy. This indicated a beneficial effect of additive treatment with nimorazole that was restricted to those tumors classified as “more” hypoxic. A Cox proportional hazards analysis for risk of locoregional failure was carried out. The following parameters were included in the model: 15-gene classification (“more” vs. “less” hypoxic), nodal status (N-positive vs. N-negative), tumor stage (T3-4 vs. T1-2), hypoxic modification (nimorazole vs. placebo), HPV status (p16-positive vs. p16-negative), gender (male vs. female), and age (<60 vs. ≥60). The 15-gene classifier was found to have independent prognostic importance with an HR of 1.41 (95% CI, 1.03–1.94). A test for interaction showed that the response to nimorazole was significantly different in the “more” hypoxic tumors compared with the response in the “less” hypoxic tumors (P = 0.003). In a separate multivariate analysis, nimorazole was an independent factor associated with good prognosis in “more” hypoxic tumors (HR = 0.42; 95% CI, 0.25–0.68), whereas there was no significant effect of hypoxic modification with nimorazole in “less” hypoxic tumors (HR = 0.98; 95% CI, 0.67–1.44).

In principle, there are at least 3 different methods to measure hypoxic status of a tumor: the physical approach, with an oxygen sensing electrode inserted into the tumor (39); a metabolic approach with preinjection and visualization of metabolically bound exogenous hypoxia tracers (10, 14, 33); or the more biology-driven approach characterized by the measurement of one or more endogenous genes being expressed under hypoxic conditions in the tumor (17, 28). Despite negative influence of hypoxia on the outcome of HNSCC after radiotherapy and despite evidence that modification of hypoxia can benefit the patients with hypoxic tumors (3, 6), such methods for pretheraputic hypoxic tumor evaluation have not yet been implemented in the clinical setting. This is probably due to an understandable combination of disadvantages either in the form of physical inconveniencies, need for preinjection of tracer, or a lack of sensitivity or specificity of the hypoxia markers (17). In this project we tried to apply the advantages of 2 of these previously established methods to measure hypoxia to refine and improve the third approach into a more clinically relevant and applicable form. Thus, we developed a hypoxia gene expression classifier enabling us to characterize an independent tumor as either “more” or “less” hypoxic, based on the quantification of hypoxia-responsive genes in the tumor biopsy.

The initially suggested potential classifier genes were characterized by being upregulated under hypoxic conditions and furthermore independent of pH fluctuations in the previous in vitro studies (27). The pH parameter was included as studies have shown that hypoxia-induced genes can also be regulated by other factors such as inflammation and pH fluctuations (17, 36, 40). Therefore, independency to one or more of these factors could hypothetically improve the hypoxia specificity of the gene in question. To improve the clinical applicability and the choice of potentially predictive genes contained in a hypoxia-specific signature/classifier, we aimed at identifying hypoxia-induced and pH-independent genes. In addition to the pH-independent genes, we also chose to evaluate the pH-dependent gene CA9 due to its frequent linkage to hypoxia in the literature (41–43).

At first, we verified that the genes being upregulated under hypoxic in vitro conditions were also upregulated in hypoxic areas of in vivo xenograft tumors (Fig. 1D). Furthermore the hypoxia-induced gene expression could also be measured from total xenograft sections which resembled tumor biopsies (Fig. 1E). We therefore believe that these genes can be used as indicators of clinically relevant hypoxia. We also registered that although 27 of the genes were significantly upregulated in the total sections compared with the nonhypoxic tumor tissue, the magnitude of upregulation varied considerably. Thus, a fraction of the suggested genes probably carried more hypoxia-specific information than all of the genes together. In a recent study by Marotta and colleagues (22), a comparable xenograft approach was applied in the description of in vitro versus in vivo hypoxia-responsive genes in a single gliosarcoma model. Few genes were common to the 2 growth conditions, and in general they found striking differences between the global in vitro and in vivo hypoxic mRNA profiles. In contrast to Marotta and colleagues, we did not aim at identifying a new set of in vivo generated hypoxia-responsive genes through an explorative supervised cDNA microarray analysis in a xenograft study, but used it as an in vivo verification of the hypoxia-induced upregulation already observed in the previous in vitro studies. Furthermore, the genes from our in vitro studies were selected on the basis of several different cell lines, and this probably explains why we in general observe a good correlation between our in vitro and in vivo settings. Interestingly, 13 of the 30 suggested hypoxia-responsive genes included in our study were also identified either in the in vitro or the in vivo studies from the Marotta and colleagues. This supports the idea that common hypoxia-induced genes can be used to evaluate hypoxia in different tumor sites (26).

In contrast to previous studies on hypoxia gene expression classifiers, we wanted to build a classifier based on the hypoxic status of the tumors. We therefore developed the final 15-gene classifier in a training set of 58 HNSCC tumors that were evaluated and ranked according to the number of low oxygen-electrode measurements of their metastatic neck nodes (31, 32, 44). We also wished to develop a robust classifier that can be used routinely to classify individual FFPE tumor biopsies. Compared with the previously published hypoxia signatures (22, 24–26) we identified some genes that are in common. In the most comprehensive study, Buffa and colleagues used a network-based approach involving both coexpression networks and gene set enrichment analysis to identify a gene expression signature that correlates with prognosis across different tumor sites (26). That there is some overlap between the genes identified by these very different approaches supports the use of gene expression as a global marker for hypoxia. However, because of the difference between the approaches, we would also expect to find genes that are not in common. First, we did not develop the classifier on the basis of its prognostic value, but solely on the association with hypoxic status in the tumor. Second, some of the coexpression networks in the network-based approach include CA9 (25, 26). Although CA9 failed in our in vitro analysis which generally excluded pH-dependent genes, we chose to include it in the analysis. In the xenograft study, CA9 was indeed upregulated in hypoxic areas consistent with the findings by Marotta and colleagues (22). However, CA9 was not correlated to the hypoxic status of the tumors in the training set, and is not among the genes in the final 15-gene classifier. Finally, we have specifically focused on test genes and reference genes that can be measured robustly in FFPE samples and where gene expression measurements in an unknown sample can be used directly to classify the tumor in a clinical setting. Although the samples in the validation set were 20 to 24 years old, we could successfully classify 323 of 326 tumors. Whether these exact classification levels can be directly applied to tumors of other cellular origin still remains to be examined.

The relevance of the classification levels in HNSCC was tested when classifying the 323 patients from the DAHANCA 5 trial as either “more” or “less” hypoxic. This validation set was unique as both prognostic and predictive evaluation was possible. We verified the prognostic relevance of the classifier by focusing on the outcome of the classified tumors in the placebo group. These tumors were treated with radiotherapy alone. We observed a significantly poorer outcome of those tumors classified as “more” hypoxic, when compared with the “less” hypoxic tumors (Fig. 2E). We believe this difference can be explained by a more pronounced resistance to radiotherapy in the “more” hypoxic tumors, thus a less efficient treatment response.

By including those randomized for hypoxic modification in the analysis, we also showed the predictive impact of the classifier. Those tumors classified as “more” hypoxic responded significantly better to the radiosensitizer nimorazole than compared with those classified as “less” hypoxic (Fig. 3A and B). This was confirmed in a multivariate analysis. Thus, we conclude that when quantifying the expression of hypoxia-induced genes in the tumor biopsy, the classifier identifies a subgroup of candidate HNSCC patients for hypoxic modification of radiotherapy. Therefore, the 15-gene hypoxia classifier attains both prognostic and predictive potential and suggests that hypoxic modification of radiotherapy should only be tailored to a subgroup of patients with gene expression classified “more” hypoxic tumors.

No potential conflicts of interest were disclosed.

Cell lines were obtained courtesy of Dr. Reidar Grenman and Dr. Michael Baumann. We thank I.M. Thuesen, B. Kierkegaard, A. Baden, and M. Johansen for technical assistance.

The work was supported by grants from the Danish Cancer Society, the Danish Medical Research Council, the Danish Council for Strategic Research, the EC FP7 project METOXIA, and CIRRO—the Lundbeck Foundation Centre for Interventional Research in Radiation Oncology.

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.