Abstract
Affymetrix U133plus2 GeneChips were used to profile 59 head and neck squamous cell cancers. A hypoxia metagene was obtained by analysis of genes whose in vivo expression clustered with the expression of 10 well-known hypoxia-regulated genes (e.g., CA9, GLUT1, and VEGF). To minimize random aggregation, strongly correlated up-regulated genes appearing in >50% of clusters defined a signature comprising 99 genes, of which 27% were previously known to be hypoxia associated. The median RNA expression of the 99 genes in the signature was an independent prognostic factor for recurrence-free survival in a publicly available head and neck cancer data set, outdoing the original intrinsic classifier. In a published breast cancer series, the hypoxia signature was a significant prognostic factor for overall survival independent of clinicopathologic risk factors and a trained profile. The work highlights the validity and potential of using data from analysis of in vitro stress pathways for deriving a biological metagene/gene signature in vivo. [Cancer Res 2007;67(7):3441–9]
- hypoxia
- head and neck cancer
- gene profiling
- breast cancer
- microarrays
Introduction
Head and neck cancer represents the fifth most common cancer in men and the eighth in women worldwide, with ∼600,000 new cases each year ( 1). Surgery, radiotherapy, and chemotherapy play a role in the management of the disease, and 5-year survival rates for patients with advanced cancers are ∼50% ( 2, 3). Many factors contribute to this poor prognosis, including late presentation of disease, nodal metastases, and the failure of advanced cancers to respond to conventional treatments ( 4).
Hypoxia, a key factor defining tumor behavior, is a characteristic of head and neck squamous cell carcinomas (HNSCC; refs. 5, 6). Although hypoxia results in cell death if it is severe or prolonged, cancer cells can adapt to this hostile environment, allowing them to survive and proliferate. It is, in part, this ability to adapt to a hostile environment that defines their malignant potential and characterizes a more aggressive phenotype ( 7). Up to 1.5% of the human genome is estimated to be transcriptionally responsive to hypoxia ( 8).
One way cells respond to reduced oxygen levels is through the hypoxia-inducible factor (HIF)-1 pathway. The HIF DNA-binding complex is a heterodimer composed of α and β subunits. In normoxia, HIF-1α undergoes rapid hydroxylation and degradation. In hypoxia, hydroxylation is prevented, stabilized HIF-1α binds to HIF-1β, and the heterodimer binds to hypoxia response elements in genes involved in the regulation of cellular processes, such as cell proliferation, angiogenesis, cell metabolism, apoptosis, and migration ( 7).
Tumor hypoxia is an independent adverse prognostic factor in many tumors, including HNSCC ( 6, 7). Evidence showing that hypoxia is important in tumor progression ( 9) and prognosis ( 7) has spurred research into developing therapies that target hypoxic cells. Therapeutic strategies include modification of the hypoxic environment or targeting components of the HIF-1 signaling pathway ( 10, 11). Although these approaches have shown some promising results, it remains difficult to identify hypoxic tumors and those patients most likely to benefit from hypoxia modification therapy.
Currently, the level of tumor oxygenation is assessed by direct or indirect methods. The main direct approach is to measure intratumoral pO2 with polarographic electrodes ( 6). Indirect techniques being explored include measuring the immunohistochemical expression of hypoxia-regulated proteins, such as carbonic anhydrase 9 (CA9) and HIF-1α ( 12, 13). The latter approach is attractive for routine clinical use. It is limited, however, by their variability of expression within a tumor, the lack of hypoxia specificity of individual proteins, and the complex interrelationships between molecular pathways in different types of tumors. Establishing tumor type–specific or tumor type–independent hypoxia gene signatures would be a major advance in this area of research.
Therefore, we assessed the mRNA profile of HNSCC samples and defined an in vivo hypoxia metagene by clustering around the RNA expression of a set of known hypoxia-regulated genes. The validity of the method was confirmed by showing that the metagene contained many previously described in vitro–derived hypoxia response genes and was prognostic for treatment outcome in independent data sets.
Materials and Methods
Sample procurement and RNA extraction. Oxford and South Manchester Ethics Committees approved the work, and all patients gave informed consent. The study comprised patients with histologically proven, previously untreated HNSCC who were undergoing surgical resection with curative intent. Tumor samples taken at operation were placed in RNAlater (Applied Biosystems, Warrington, United Kingdom) for up to 12 h before cryopreservation in liquid nitrogen. Where possible, biopsies of healthy mucosa were taken from the contralateral side and used as “normal” samples. Subsequently, the samples were divided and half paraffin embedded for histologic analysis. RNA was extracted from the remaining tissue using Tri Reagent (Sigma-Aldrich, Poole, United Kingdom) and washed using the Qiagen spin column (Qiagen, Crawley, United Kingdom). RNA quality and quantity were confirmed using the NanoDrop ND-1000 spectrophotometer and the Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA).
Immunohistochemistry. Immunohistochemistry for CA9 and HIF-1α was done on sections from formalin-fixed, paraffin-embedded tumor biopsies (one section per protein). CA9 immunohistochemistry has been described previously and involved a mouse monoclonal anti-human antibody (M75) raised to the external domain of CA9, a gift from Profs. S. Pastorekova and J. Pastorek ( 14). The primary HIF-1α antibody was ESEE122, a mouse monoclonal IgG1 antibody to human HIF-1α (University of Oxford, Oxford, United Kingdom). Tumor membrane CA9 expression was estimated over the whole section as a percentage. Tumor nuclear HIF-1α expression was assessed in 10 separate fields under intermediate (×200) magnification. The percentage of tumor cells with nuclear reactivity in each field was recorded as a continuous variable, and the overall score for each tumor was the mean percentage score across the 10 fields.
Microarray data. Affymetrix U133plus2 GeneChips were used. Biotinylated target RNA was prepared with minor modifications from the manufacturer's recommendations. A description of procedures is available online. 9 Target RNA generated from each sample was then processed as per the manufacturer's recommendation using an Affymetrix GeneChip Instrument System. 10 Briefly, spike controls were added to 10 μg of fragmented cRNA before overnight hybridization to HGU133plus2 arrays. Arrays were then washed and stained with streptavidin-phycoerythrin before imaging on an Affymetrix GeneChip (3000) scanner. 11 Following scanning, data files were processed using Affymetrix GeneChip Operating Software version 1.1.1.052. Subsequent data processing was done using affy ( 15), simpleaffy ( 16), and GCRMA ( 17) packages from BioConductor and R. mRNA was estimated using GCRMA, and log2 values were used.
After scanning, array images were assessed visually to confirm scanner alignment and the absence of significant bubbles or scratches. The 3′/5′ ratios for glyceraldehyde-3-phosphate dehydrogenase (GAPDH) and β-actin fell within acceptable limits (0.9–1.45; 1.14–2.87), and BioB spike controls were present on all arrays, with BioC, BioD, and CreX also present in increasing intensity. When scaled to a target intensity of 100 (using the MAS5 algorithm, as implemented in simpleaffy), scaling factors for all arrays were within acceptable limits, as were background, Q values, and mean intensities. Quality control data are provided as Supplementary Material.
The hypoxia metagene. A flow chart summarizing the method used to derive the hypoxia metagene is in the Supplementary Material. The in vivo expression of genes that correlated strongly with 10 hypoxia-regulated genes was investigated. Genes representing a wide range of known hypoxia-regulated pathways were chosen, including those involved in angiogenesis (VEGF and ADM), glucose transport and glycolysis (GLUT1, PDK-1, and ENO1), linking mitochondrial metabolism and glycolysis (HK2), pH regulation (CA9), reuse of ATP breakdown products (AK3), and inhibition of cell cycle progression (CCNG2). PFKB3 was also used, which regulates the overall rate of glycolysis and is the most active isozyme of the four that regulate the concentration of fructose-2,6-bisphosphate, which activates the rate-limiting enzyme 6-phosphofructo-1-kinase. Genes selected have been widely investigated justifying their classification as hypoxia-responsive genes, particularly in terms of promoter analysis, and have been frequently reported in gene array analysis of hypoxia-induced cell lines. Pearson correlation was used within the framework of significance analysis of microarrays (SAM; refs. 17– 19). To account for multiple testing, the local false discovery rate (FDR) was calculated for each probe set using SAM 2.21 within BioConductor. For each gene, a local FDR of <0.05 was used as cutoff to determine cluster membership. This tells us which probe sets are associated with the hypoxia proband genes in HNSCC in vivo. The probe sets chosen comprised up-regulated and down-regulated genes; the up-regulated genes made up the hypoxia metagene. However, the presence of a probe set in a single cluster does not guarantee the involvement of its target gene in hypoxia regulation in vivo. Some probe sets might have binding artifacts or defective probes; the in vitro and in vivo behavior of the genes could differ or probe sets might be aggregated from genes unrelated to hypoxia pathways. To minimize random aggregation, the number of clusters in which a probe set appeared was calculated. Probe sets appearing in >50% of the clusters were chosen. This conservative cutoff was selected after a Monte-Carlo simulation was done computing cluster aggregation around randomly selected sets of genes; this simulation showed that the probability of having an overlap >50% between 10 randomly selected genes was P < 0.001.
Calculation of hypoxia scores. The pattern of expression of the up-regulated and down-regulated genes was investigated in the 59 HNSCC using both a hierarchical clustering algorithm with Bayesian optimization of the number of clusters and a continuous hypoxia score (HS), which ranked patients according to expression of genes in the hypoxia metagene. The clustering method has been described previously ( 20) but used a two-way approach by clustering on the genes first and then using the median centroids of the gene clusters to cluster on the samples. HSs were calculated from the distribution of the level of RNA expression of the up-regulated and down-regulated genes (HS-up and HS-down) for each patient. Median values were chosen to avoid effects of outliers and nonnormality of the distribution. Furthermore, to make HS-up and HS-down general and transportable to any data set, the fractional rank was computed, where patients were ranked based on their score, and rank was normalized between 0 and 1.
Statistical analyses. Two published data sets were analyzed: series of 60 HNSCC ( 21) and 295 breast cancers ( 22– 24). Details on the treatment and follow-up are in the cited publications. Briefly, HNSCC patients underwent a variety of treatments, including various combinations of surgery, radiotherapy, and/or chemotherapy, and had a median follow-up of 16 months. Breast cancer patients underwent mastectomy or breast-conserving surgery followed by radiotherapy if indicated, and the median duration of follow-up was 7.8 years. HSs were calculated from the publicly available microarray data set as described above. Univariate and multivariate Cox survival analyses were carried out. HSs were introduced together with relevant clinical variables, and stepwise backward likelihood was used to select covariates that correlated with outcome. A flow chart in the Supplementary Material summarizes the approach used to analyze the published data sets.
Results
Samples and RNA extraction. Sixty-eight patients with primary HNSCC were recruited, and samples from 9 were found to contain <10% tumor cells. Patient demographics for the remaining 59 patients are presented in Table 1 . Normal tissue samples were obtained from 11 patients. All RNA was of high quality and met the manufacturer's quality control criteria.
Summary of HNSCC patient data
Evaluation of HIF-1α and CA9 protein expression to obtain a hypoxia metagene. There was a significant correlation between HIF-1α and CA9 protein expression in the 59 tumors (Spearman ρ = 0.53; P < 0.001; Supplementary Fig. S1). Supervised clustering on the percentage expression of the two proteins showed that CA9 mRNA correlated strongest with CA9 protein expression (Spearman ρ = 0.60), but the overall FDR was high for both lists and a heterogeneous pattern of gene expression was produced, no obvious groupings and only 1% overlap between the top 50 genes in each list. This was true both when using Pearson and Spearman correlation within SAM. Comparison of the genes with a list of 245 published hypoxia-regulated genes taken from relevant reviews (Supplementary Table S1; refs. 7, 8, 11, 25, 26) showed that none of the top 20 HIF-1α–associated genes appeared in the literature list and only 3 of the top 20 CA9-associated genes. Therefore, clustering on percentage of protein expression of a hypoxia-regulated gene was unreliable for deriving a hypoxia metagene.
Generation of the hypoxia metagene by clustering HNSCC around the RNA expression of known in vitro hypoxia-regulated genes. Clustering around the CA9 RNA expression identified several genes that correlated strongly with CA9 and were known to be hypoxia associated (10 of the top 20 were in the literature list). This suggested that the approach might identify other hypoxia-related genes in vivo. Furthermore, the cohort of patients were divided broadly into two groups (Supplementary Fig. S2), with one group showing high expression of some known hypoxia-regulated genes.
A single gene would not reflect the molecular diversity of the response of a cell to hypoxia and was considered a high-risk approach due to possible binding artifacts, differences between in vitro and in vivo expression for single probes or genetic variability in response. Therefore, 10 widely investigated hypoxia-regulated genes were selected for the clustering process (VEGF, ADM, GLUT1, PDK-1, ENO1, HK2, PFKB3, CA9, AK3, and CCNG2). Genes clustering around the 10 feeder genes were compared with 245 hypoxia-associated genes taken from the literature. Two genes showed an unusual clustering pattern: CCNG2 clustered genes correlated inversely with the literature genes and PFKFB3-associated genes correlated with none of the literature genes. The remaining eight clusters contained a varying percentage of literature genes; this percentage increased as the list of highest ranked genes was shortened (14–30% when the 50 most correlated genes were considered). This improved when genes appearing in more than one cluster were considered; 67% of the up-regulated genes found in >70% of clusters were in the literature list (Supplementary Fig. S3). Obviously, the number of genes selected using this method decreased as the amount of overlap increased. A final level of 50% overlap was chosen (i.e., transcripts appearing in >50% of clusters were considered). This was based on a Monte-Carlo simulation and offered a balance between identifying new genes and including previously reported genes.
There were 210 probe sets in the final gene list, 129 up-regulated and 81 down-regulated, representing 99 (126 probe sets with known gene identities) and 57 (68 probe sets with known gene identities) genes, respectively. Of these, 27% ( 27) and 2% ( 1), respectively, were present in the literature list (i.e., they were known to be hypoxia regulated). Table 2 lists the 99 up-regulated genes that comprised the hypoxia metagene. Supplementary Tables S2 and S3 list the up-regulated and down-regulated genes with correlation coefficients.
The 99 genes in the hypoxia metagene
Composition of the hypoxia metagene. The up-regulated gene list ( Table 2) contained several genes previously described ( 7, 11) as induced by hypoxia, including those involved in glucose metabolism (e.g., ALDOA and GAPDH), angiogenesis (ANGPTL4 and PGF), cell migration (TPBG), and the regulation of apoptosis (BNIP3). The known genes represent a broad spectrum of molecular pathways involved in cellular response to hypoxia. Some of the genes have also been described in RNA profiles of HNSCC [e.g., IL8 ( 27) and PLAU ( 28)]. Genes not previously reported to be up-regulated under hypoxia include MTX1, BCAR1, HES2, PSMA7, and SLCO1B3.
Pattern of expression of the hypoxia metagene in HNSCC. Investigation of the pattern of gene expression in the 59 HNSCC showed that the genes were clustered into two main groups, whereas tumors were clustered into three groups. The two gene clusters are a consequence of how the genes were selected (i.e., positively and negatively correlated with known hypoxia-regulated genes). The three tumor clusters are characterized by differing patterns of expression, with one showing high expression of the 99 up-regulated genes and low expression of the 57 down-regulated genes. However, when bootstrapping was applied to the tumor clusters (random sampling of patients with replacement), there was a high level of variability in the classification of tumors into clusters. This suggested that the hypoxia metagene was better described as a continuous or graded characteristic rather than a defined binary status, and it might be difficult to assign reliably future patients to a specific cluster. Using a continuous score also enables patient risk adjustment (based on hypoxia) on an individual basis.
Two HSs were calculated for each patient from the distribution of the level of RNA expression of the list of up-regulated and down-regulated genes (HS-up and HS-down). The mean and SD of HS-up were 8.76 ± 0.69 (median, 8.80; range, 5.68–10.04). In Fig. 1A , patients were ranked based on HS-up and hierarchical clustering was carried out on all the up-regulated and down-regulated genes. On the X axis, tumors are ordered by increasing HS-up so those on the far right have high expression of hypoxia-associated genes. Guidelines indicate median and quartiles of HS-up for the 59 patients. The Y axis shows the results of hierarchical clustering on the whole gene list to illustrate which transcripts cluster together (e.g., different transcripts of the same gene cluster together and one cluster contains VEGF, CA9, and ANGPTL4). Using HS-down to rank patients produced a different, but qualitatively similar, heatmap as HS-up and HS-down were inversely correlated (ρ = −0.75; P < 0.001).
A, dendogram of 11 normal tissue samples and 59 HNSCC. Two HSs were calculated for each patient by looking at the distribution of the level of RNA expression of up-regulated and down-regulated hypoxia-associated gene (HS-up and HS-down). The median of this distribution was then calculated for each patient. X axis, patients are ordered by increasing HS-up; Y axis, results of hierarchical clustering on all the genes. The guidelines indicate the median and quartiles of the HS-up for the 59 patients. B, Kaplan-Meier plot of recurrence-free survival in 60 patients with HNSCC stratified according to HS-up (the hypoxia metagene). The patients are from a published study involving the molecular classification of HNSCC using patterns of gene expression ( 21). C, patients were first stratified according to quartiles, which suggested the separation into the two groups shown [i.e., the highest HS-up quartile (0.75−1) versus the remaining three quartiles (0−0.75)]. The number of events/patients in the two arms was 10/45 and 8/15, respectively.
The pattern of expression of the hypoxia metagene in normal tissue. Normal tissue samples showed a pattern of expression similar to the normoxic tumors (i.e., low median score of the up-regulated cluster when a continuous HS-up was considered; Fig. 1A). For example, the median expression of VEGF was similar in normal tissue (9.7) and “oxygenated” tumors (9.5) but significantly lower than in “hypoxic” tumors (11.0; P = 0.005). There was also no correlation between tumor and normal tissue expression levels in the same patient. In a paired comparison using SAM with a FDR of <0.05, 55 of the 210 transcripts (representing 44 of the genes) in the gene list were expressed significantly higher in tumor versus normal tissue.
Hypoxia metagene as a prognostic factor in an independent HNSCC data set. As no treatment outcome data were available for our HNSCC data set due to a short follow-up, it was not possible to train our hypoxia metagene on outcome or test its significance on the present data set. Therefore, to explore whether our signature could provide prognostic information, it was applied to a publicly available data set ( 21). The characteristics of our 59 metagene-generating and the 60 literature HNSCC were compared using χ2. There were no significant differences in T, N, pT, and pN status, but our series contained more oropharynx and fewer larynx cancers (P = 0.010) and more poorly differentiated tumors (P = 0.045) than the published data set. The data set described a four-class expression profile, one of which was associated with a poor recurrence-free survival. Patients in the Chung data set were stratified according to HS-up quartiles or two groups (upper quartile versus rest). High HS-up expression was an adverse prognostic factor for recurrence-free survival ( Fig. 1); the hazard ratio (HR) for the upper quartile versus the rest was 3.64 [95% confidence interval (95% CI), 1.43–9.31]. A multivariate Cox regression analysis was done, which included HS-up entered as a continuous fractional rank, Chung's intrinsic profile, age, gender, smoking history, alcohol, clinical stage, lymph node status, tumor grade, tumor subsite (hypopharyngeal versus others), and array batch ( Table 3 ). These factors were chosen as they were used in the multivariate analysis in the original study. Significant adverse prognostic factors were Chung's intrinsic profile (P = 0.022), poor tumor differentiation (P = 0.003), and HS-up (P = 0.012; Fig. 2A ). There was no significant difference in survival when HS-down was considered. When the individual genes were analyzed separately, VEGF, GLUT1, and HK2 predicted a worse prognosis but were associated with large confidence intervals ( Fig. 2B). A HS for the whole literature list had no prognostic significance in either univariate or multivariate analysis.
Cox multivariate regression analyses of survival in published data sets
Cox multivariate analysis derived HRs with 95% CIs. A, variables used in the original publication of an independent HNSCC patient data set were entered in the model along with HS-up. HS-up was entered as a continuous variable (fractional rank varying from 0 to 1), and the HR reported represents the risk for an increasing HS-up quartile (the HR for the two ends of the spectrum was 14.83; 95% CI, 1.80−122.35). B, HR for the RNA expression of the individual genes used to derive the hypoxia metagene for a Cox multivariate analysis of recurrence-free survival in the independent HNSCC data set. No data were available for ENO1. C, variables used in the original publication in an independent breast cancer data set were entered into a Cox multivariate regression model along with HS-up. The HR for two ends of the HS-up spectrum was 3.58 (95% CI, 1.53−8.39).
Hypoxia metagene as a prognostic factor in another cancer type. The ability of HS-up to predict outcome was investigated in a breast cancer series ( 22). Increasing HS-up was a significant adverse prognostic factor for metastasis-free ( Fig. 3A and B ) and overall survival when patients were stratified by either quartiles or two groups. The HRs for metastasis-free survival were 1.62 (95% CI, 1.33−1.98) for quartiles and 2.83 (95% CI, 0.78−3.38) for two groups. In a multivariate analysis, including HS-up with the clinical variables included in the original study ( Table 3), HS-up retained significance for metastasis-free (P = 0.003; Fig. 2C) and overall (P = 0.002) survival. In a multivariate analysis, including the 70-gene trained intrinsic signature, although HS-up was not significant for metastasis-free survival, it retained significance for overall survival (P = 0.003; Supplementary Table S4). When a cell line–derived wound response signature ( 22) was included, HS-up retained independent prognostic significance for metastasis-free (HR, 2.80; 95% CI, 0.09−7.22; P = 0.033) and overall (HR, 3.01; 95% CI, 1.06−8.58; P = 0.038) survival. There was only 2.2% overlap between the wound response and our hypoxia signature (9 of 401 known genes: LDLR, MNAT1, NME1, SLC16A1, CDCA4, CORO1C, NUDT15, PSMA7, and PSMD2). When both the 70-gene trained and wound response signatures were included in the analysis, HS-up retained prognostic significance for overall survival but the wound signature did not.
Kaplan-Meier plots of metastases-free survival in 295 patients with breast cancer. The data were taken from Chang et al. ( 22). A, stratified according to HS-up quartiles. The numbers of events/patients for increasing quartiles were 12/73, 18/73, 20/74, and 38/73. B, stratified by highest HS-up quartile (0.75−1) versus the remaining three quartiles (0−0.75). The number of events/patients for the two arms was 50/220 and 38/73, respectively. C, stratified using the Chi gene expression signature of cellular response to hypoxia ( 29). Two arms are for patients with high (mean score, >0; 35 events/77 patients) and low (mean score, <0; 53 events/218 patients) mean RNA expression of hypoxia-associated genes. D, stratified using the hypoxia gene signature derived in this article but using the Chi et al. ( 29) approach to compute the hypoxia score. The two arms are for patients with high (mean score, >0; 42 events/87 patients) and low (mean score, <0; 46 events/208 patients) mean RNA expression of hypoxia-associated genes.
Comparison with a published cell line–derived hypoxia signature. A gene expression signature of cellular response to hypoxia in vitro was recently published ( 29). There was little overlap between the Chi and our signature, with only eight genes in both lists. Whereas our signature contained 27% of the literature genes, the Chi profile contained 17%. To compare the different approaches used to derive hypoxia signatures, the above analyses on the Chang data set were repeated, including the Chi signature. In the Chi et al. article, mean RNA expression was used as a hypoxia signature metric rather than the median used here and patients were stratified according to high versus low gene expression. Analyses were carried out, therefore, using both methods, and in both cases, our approach proved superior. In the breast cancer series of van't Veer, the Chi signature quantified as a median score was a significant prognostic factor for metastasis-free survival in univariate (HR for quartile increase, 1.39; 95% CI, 1.14−1.70; P = 0.004, log-rank test) and multivariate (HR, 2.74; 95% CI, 1.29−5.83; P = 0.0090) analyses. These data compare less favorably with our signature in both univariate (HR for quartile increase, 1.62; 95% CI, 1.33−1.98; P < 0.0001, log-rank test; Fig. 3A) and multivariate (HR, 3.58; 95% CI, 1.53−8.39; P = 0.0030) analyses. When both signatures were included in the multivariate analysis, the Chi signature lost prognostic significance. Similar results were obtained when considering overall survival. The analyses were repeated using mean RNA expression of hypoxia-associated genes and stratifying by high versus low hypoxia response groups as done by Chi et al. ( 29). Similar results were seen in univariate (Chi: HR, 2.15; 95% CI, 1.40−3.30; P < 0.001, log-rank test; ours: HR, 2.67; 95% CI, 1.76−4.06; P < 0.001, log-rank test; Fig. 3C and D) and multivariate (Chi: HR, 1.91; 95% CI, 1.22−3.00; P = 0.0047; ours: HR, 2.06; 95% CI, 1.29−3.31; P = 0.0019) analyses. Again, when both signatures were included, only ours retained prognostic significance. Likewise, repeating the analyses for overall survival produced similar results.
Discussion
The limitations of using protein markers to reflect tumor hypoxia. High expression of HIF-1α and CA9 is associated with adverse prognosis in HNSCC ( 12, 30). Although high expression of HIF-1α and CA9 was thought to reflect the hypoxic nature of a tumor and activation of the HIF pathway, other studies reported no association with survival ( 31, 32) or association for only one factor ( 33). Some of these anomalous findings have been explained by the different half-lives for CA9 (days) and HIF-1α (minutes) proteins ( 34). It is more probable that, because hypoxia influences many biological pathways, a single factor is incapable of adequately describing this complex response.
In this study, a cluster analysis based on the protein expression of HIF-1α or CA9 produced gene profiles with minimal overlap. The CA9 protein and RNA expression correlated, indicating that mRNA is readily translated into the protein in HNSCC. In contrast, the post-translational modification of HIF-1α ( 11) explains the lack of correlation for HIF. A cluster analysis around the mRNA expression of CA9 identified several genes recognized as hypoxia inducible involving a variety of molecular signaling pathways, suggesting the potential to identify novel hypoxia-regulated genes. Furthermore, the cluster analysis divided the cohort of patients into two groups, with one showing high expression of known hypoxia up-regulated genes and thus possibly reflecting hypoxic versus oxygenated tumors.
The hypoxia metagene based on clustering HNSCC around the RNA expression of known in vitro hypoxia-regulated genes. Ten genes shown to be hypoxia regulated in vitro were used to cluster other genes in vivo. To test the validity of this approach in identifying hypoxia-regulated genes, clustered genes were compared with a set of 245 genes reported in literature reviews as hypoxia regulated. The approach clustered other well-recognized hypoxia-regulated genes, suggesting that the method might identify novel hypoxia-associated genes. Our approach highlighted an interesting finding about differences in hypoxia response in vivo versus in vitro. Two of our selected genes showed unexpected clustering patterns. Genes clustering around CCNG2, involved in cell cycle regulation, correlated inversely with genes in the literature list. None of the genes clustering around PFKFB3, important in the metabolism of fructose, correlated with any of the genes in the literature list. Interestingly, both are up-regulated in vitro in response to hypoxia. Tumor type–specific differences in hypoxia response might also be important, as CCNG2 was shown to be down-regulated in another gene profiling study in HNSCC ( 27).
In addition, some studies advocate the use of laser microdissection to ensure that only the tumor RNA is extracted ( 35), but this ignores the effects the tumor may have on the surrounding stroma. Therefore, in this study, we have used tumor and surrounding stroma to extract RNA while quantitatively confirming the presence of tumor. Several genes were identified that clustered with the hypoxia-associated genes but have not previously been shown to be hypoxia regulated [e.g., MTX1, a component of a transport complex in the outer mitochondrial membrane ( 36); BCAR1, which encodes a docking protein critical in controlling integrin-dependent signaling ( 37); PSMA7, involved in ubiquitin-dependent protein catabolism and plays a part in the regulation of HIF-1α ( 38); and SLCO1B3, a member of the organic anion transporting polypeptide family of transporters ( 39)].
Deriving a metric for the clinical application of the hypoxia metagene. The hypoxia metagene clustered the 59 HNSCCs studied into three tumor groups, which were not robust to bootstrapping. Thus, a continuous score was used to rank the tumors from low to high hypoxia metagene expression ( Fig. 1). While this work was in progress, Chi et al. ( 29) published a hypoxia gene profile derived in a different way. They analyzed the expression of genes induced in vitro by hypoxia in a panel of normal cell lines and showed a dependency on both cell type and the intensity/duration of hypoxia ( 29). They used a set of common genes expressed in the normal cell types studied for their profile, which predicted relapse in ovarian and breast cancers.
In both studies, HSs were derived to provide a metric applicable to other data sets. Chi et al. ( 29) used the mean gene expression levels of the genes in their signature, whereas in our study median values were calculated as the median is a more robust statistic not affected by outliers. Furthermore, to make HS general and transportable, the ranked score was used rather than raw values. It is of interest that, in our profile, 27% of the genes were previously recognized as hypoxia regulated in the literature, suggesting the identification of new pathways regulated by hypoxia. There were nine transcripts/eight genes that appeared in both signatures: P4HA1, PGK1, CA9, ANGPTL4, NDRG1, LDHA, SLC2A1, and SLC6A8. These eight genes reflect a variety of pathways (glycolysis, angiogenesis, extracellular matrix deposition, and pH regulation) and may represent a core set for future classifications. It is of interest that VEGF is not present in the Chi hypoxia gene profile.
The hypoxia metagene as a prognostic factor. In both univariate and multivariate analyses, the hypoxia metagene was a significant prognostic factor for recurrence-free survival in an independent HNSCC data set. In multivariate analysis, HS-up (P = 0.012) compared favorably with the original study four-class intrinsic classifier (P = 0.022). Interestingly, the down-regulated genes did not correlate with outcome, suggesting that they play a less significant part in tumor treatment response. Another finding of interest was that HS-up correlated with the intrinsic profile in the Chung study (Kruskal-Wallis χ2 = 22.3; P < 0.001); specifically, HS-up was higher in cluster 1 than in the other three Chung clusters. Cluster 1 correlated with a poor recurrence-free survival in the Chung study and had a high expression of genes that were associated with activation of the epidermal growth factor receptor pathway. Hypoxia up-regulated genes in Chung cluster 1 were SLC16A2, TGF-α, and KRT14.
In the breast cancer data set, our hypoxia metagene was superior to both the wound response and Chi signatures. When the trained intrinsic classifier was included in the multivariate analysis, only our hypoxia metagene retained independent prognostic significance for overall survival. The difference in efficacy of the two hypoxia signatures may reflect the way they were developed: from normal cells in vitro (Chi signature) versus a combination of tumor cells in vitro plus cancers in vivo (our signature).
It is of interest to note that both Chi ( 29) and our studies initially investigated a specific tumor type, renal cancer or HNSCC, respectively, yet found applicability across tumor types. Ein-Dor et al. ( 40) highlighted the lack of overlap between expression profiles, which predict cancer treatment outcome and showed that many equally predictive gene lists could be produced from the van't Veer breast cancer signature. It was suggested that this is due in part to the many genes that correlate with survival. However, Shen et al. ( 41) analyzed four independent microarray studies to derive an interstudy validated metasignature associated with breast cancer prognosis, which was comparable or better at providing prognostic information compared with the intrinsic signatures. With the publication of further hypoxia metagenes ( 42), therefore, our study and that of Chi et al. ( 29) are working toward the derivation of a future hypoxia metasignature with potential application across a range of tumor types. A rationale behind the derivation of phenotype-specific profiles is that metagenes derived from defined pathways may be more powerful markers of prognosis but also particularly of value in selecting patients for specific interventions.
In summary, hypoxia results in molecular changes that promote an aggressive phenotype and reduce the efficacy of conventional treatments, resulting in a significant therapeutic challenge. Using a combined in vitro and in vivo approach, a hypoxia metagene was derived that outperformed the original intrinsic profile of a published HNSCC data set and was an independent prognostic factor in a widely studied breast cancer microarray series. Our combined in vitro and in vivo strategy has potential for the future development of metasignatures, which are not based on the extensive training of intrinsic profiles. In conclusion, gene signatures can be derived in vivo that reflect biological phenotypes relevant in determining cancer patient prognosis.
Acknowledgments
Grant support: Cancer Research UK, National Translational Cancer Research Network of the United Kingdom, and EU Integrated Project ACGT (FP6-IST-026996).
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.
We thank Dr. Jo Cresswell for administrative assistance in running the project, Stuart Pepper for advice on the microarray work, and the Molecular Biology Core facility of the Paterson Institute (Manchester, United Kingdom), which houses the Cancer Research UK Affymetrix facility.
Footnotes
-
Note: Supplementary data for this article are available at Cancer Research Online (http://cancerres.aacrjournals.org/).
-
S.C. Winter and F.M. Buffa are joint first authors. C.M.L. West and A.L. Harris are joint senior authors.
-
Supplemental and microarray data will be made available at: http://bioinformatics.picr.man.ac.uk/experiments.
-
↵9 http://bioinf.picr.man.ac.uk/mbcf/Downloads/GeneChip_Target_Prep_Protocol-CR-UK_v4.pdf
-
↵10 http://www.affymetrix.com/support/technical/manual/expression_manual.affx
-
↵11 http://bioinf.picr.man.ac.uk/mbcf/Downloads/GeneChip_Hyb_Wash_Scan_Protocol-CR-UK_v4.pdf
- Received September 7, 2006.
- Revision received December 13, 2006.
- Accepted January 25, 2007.
- ©2007 American Association for Cancer Research.