Purpose: Hepatocellular carcinomas (HCC) have an unpredictable clinical course, and molecular classification could provide better insights into prognosis and patient-directed therapy. We hypothesized that in HCC, certain microenvironmental regions exist with a characteristic gene expression related to chronic hypoxia which would induce aggressive behavior.
Experimental Design: We determined the gene expression pattern for human HepG2 liver cells under chronic hypoxia by microarray analysis. Differentially expressed genes were selected and their clinical values were assessed. In our hypothesis-driven analysis, we included available independent microarray studies of patients with HCC in one single analysis. Three microarray studies encompassing 272 patients were used as training sets to determine a minimal prognostic gene set, and one recent study of 91 patients was used for validation.
Results: Using computational methods, we identified seven genes (out of 3,592 differentially expressed under chronic hypoxia) that showed correlation with poor prognostic indicators in all three training sets (65/139/73 patients) and this was validated in a fourth data set (91 patients). Retrospectively, the seven-gene set was associated with poor survival (hazard ratio, 1.39; P = 0.007) and early recurrence (hazard ratio, 2.92; P = 0.007) in 135 patients. Moreover, using a hypoxia score based on this seven-gene set, we found that patients with a score of >0.35 (n = 42) had a median survival of 307 days, whereas patients with a score of ≤0.35 (n = 93) had a median survival of 1,602 days (P = 0.005).
Conclusions: We identified a unique, liver-specific, seven-gene signature associated with chronic hypoxia that correlates with poor prognosis in HCCs. Clin Cancer Res; 16(16); 4278–88. ©2010 AACR.
Hepatocellular carcinoma (HCC) has an increasing incidence and high mortality. Prognosis is difficult to predict for individual patients. We started with an in vitro model of chronic hypoxia and determined differentially expressed genes. The clinical value of these chronic hypoxia–related genes was assessed using bioinformatics in four data sets with HCC patients. We identified a very small gene set with important prognostic value in HCC patients. We developed a score based on these seven genes which correlates with survival and early recurrence. The score should make it feasible to better predict prognosis for an individual patient. This method requires the determination of only seven genes in contrast with most of the previous microarray studies that needed large cohorts for clustering and determining relative gene expression. Furthermore, we suggest that gene expression under chronic hypoxia might lead to new therapeutic targets, and the hypothesis-driven experimental methods we used could be applied to other malignancies as well.
Hepatocellular carcinoma (HCC) is the sixth most common malignancy in the world and the third most common cause of cancer-related deaths (1). Every year, 600,000 new cases are diagnosed and almost just as many patients die annually of this disease. The most important risk factor for the development of HCC is cirrhosis, which is present in 80% of patients.
HCC generally develops over many years and unless the patients at risk are frequently screened, the disease is often diagnosed at a later stage (2). Only a small proportion of patients are eligible for surgical resection or transplantation; the majority have to rely on systemic chemotherapy, embolization, or supportive regimens (3). Patients with HCC have a heterogeneous course of disease. Therefore, therapeutic decisions for a specific patient would greatly benefit from a reliable prognostic indicator considering the treatment burden, duration, expected success rate, and costs.
Several studies have tried to identify gene sets with prognostic or diagnostic relevance by using microarray analysis. Each study resulted in its own classification with a specific separation into groups (4). All these microarray studies show remarkably little overlap and it is difficult to find a clear correlation between the molecular classes and prognosis. The results of these studies seem to be center-dependent because of the different microarray techniques used, the small heterogeneous cohorts that were studied, and the different clinical variables used for the evaluations (5).
Although it has been recognized that the microenvironment plays a role in tumor biology, it has never been used as a rational concept to study HCC. One of the factors that seems to affect cancer cell behavior and patient prognosis is hypoxia (6). Hypoxia is associated with poor prognosis in solid tumors and with the development of resistance to chemotherapeutic agents and radiation (7, 8). Although HCC is a hypervascular malignancy, there are regions with hypoxia (9). Hypoxic regions are already present in the early stage, when the vasculature is not sufficiently extended, and in more advanced stages, when rapid cell proliferation induces hypoxia (10). Furthermore, the neovascularization is unorganized and the structural and functional defects lead to both acute hypoxia and chronic hypoxia (11, 12). Moreover, liver cancer usually develops in a cirrhotic environment in which the blood flow is already impaired.
We hypothesized that in HCC, there are regions with sustained hypoxia which induces a characteristic gene expression pattern. In addition, the extent of hypoxic gene expression determines the aggressiveness, or in general, the prognosis.
Materials and Methods
Chronic hypoxic gene expression
We used the human hepatoblastoma cell line HepG2 as an in vitro model. Cell culture experiments were done within 6 months after the cell line was obtained from American Type Culture Collection.
HepG2 cells were grown in a humidified incubator (Sanyo MCO-18M O2/CO2 incubator; 20% O2, 5% CO2 at 37°C) in Williams medium E (WEM; InVitrogen) supplemented with 10% FCS, 2 mmol/L of l-glutamine, 20 mU/mL of insulin, 50 nmol/L of dexamethasone, 100 units/mL of penicillin, 100 μg/mL of streptomycin, 2.5 μg of fungizone, 50 μg/mL of gentamicin, and 100 μg/mL of vancomycin (WEM-C).
Dynamic behavior of hypoxia-related gene expression was studied over a period of 72 hours at different oxygen concentrations using real-time PCR and immunohistochemistry. From these results (Supplementary Figs. S1 and S2; Table 4), we concluded that 72 hours and 2% of O2 represents a chronic hypoxia–like condition distinct from the acute situation during the first 24 hours.
Differentially expressed genes between HepG2 cells cultured for 72 hours at 2% O2 (hypoxia) and cells cultured in parallel in 20% O2 (normoxic) were determined with microarray analysis. For the microarray analyses, two experiments were executed in parallel. Cells were seeded at 3 × 106 in 75 cm2 tissue culture flasks (n = 4) at 20% O2 and were grown until 70% confluence. After reaching near-confluence, two flasks were placed in a humidified incubator under hypoxic conditions (2% O2, 5% CO2 at 37°C), whereas the other flasks (n = 2) remained in normoxic conditions (20% O2). Cells were cultured for 72 hours in these different oxygen conditions, and after 3 days, cells were harvested. RNA samples were labeled and hybridized on dual-color Agilent Human Whole Genome Oligo Microarray (Agilent). The Agilent technology uses one glass array for the simultaneous hybridization of two populations of labeled, antisense cRNAs obtained from two samples (hypoxic and normoxic control) leading to a ratio for expression. Probes with corrected P values of <0.01 and a fold change of >2 were defined as highly significant and designated as the in vitro hypoxia gene set (Supplementary Table S2).
Clinical data sets and patient characteristics
To determine the clinical relevance of in vitro gene expression, we compared our findings in the largest microarray data sets with corresponding clinical information that is available in public databases. We used four publicly available HCC studies, which are summarized in Table 1 (13–17). By testing only the differentially expressed genes in the four data sets, the reliability of genome-wide research can be increased because the number of genes tested is in proportion with the number of patients studied (18).
Development of a prognostic gene set
We used the first three data sets as training sets to optimize our in vitro hypoxia gene set and to investigate the prognostic correlation. The fourth data set (Chiang et al.), was subsequently used to independently validate the signature. In the first step to defining a robust score from the different data sets, we used a global test (19) to investigate whether the hypoxia genes were associated with prognosis under a Q2-null hypothesis (20). This approach should be advantageous in that it is less dependent on the array platform used in different laboratories (Affymetrix, Agilent, Stanford, etc.). Moreover, by starting from a small subset of in vitro–determined hypoxia genes, this method provides more insight into the degree of relationship between the different genes found to be upregulated or downregulated. This enrichment-based method was then used to investigate whether our hypoxia set separates the good and poor prognostic characteristics in the three data sets individually. Because in all four data sets, another prognostic factor was reported, we also had to use a different prognostic factor in every data set. From Boyault et al. (13), the FAL index was used (21), which is a measure for chromosomal instability and a high score (>0.128) was associated with poor prognosis. In Wurmbach et al. (16), vascular invasion was used (22, 23). Lee et al. (14, 15) used the different prognostic clusters that correlate with survival (cluster A with poor prognosis and cluster B with good prognosis), and Chiang et al. (17) used the Barcelona Staging Classification (BCLC; ref. 24). Finally, based on the z scores, we determined which genes in our hypoxia set were most influential in differentiating between good and poor prognosis (19).
Chronic hypoxia gene signature
We started with a cell culture model and determined the differentially expressed genes in HepG2 cells that were cultured for 72 hours at either 20% oxygen or in hypoxic conditions at 2% oxygen. We used Agilent technology with color flip on two independent experiments in duplicate, resulting in eight ratio values. A total of 37,707 spots showed a representative signal, of which 3,592 (8%; representing 3,259 genes) had a fold change of >2 and a P value of <0.05 (1,879 upregulated and 1,713 downregulated). To control the false discovery rate, multiple testing corrections were done and probes with a corrected P value of <0.01 and a fold change of >2 were selected (25). This resulted in a list of 265 highly significant genes (207 upregulated and 58 downregulated; Supplementary Table S2), designated as the in vitro hypoxic gene set (the microarray data has been deposited at the National Center for Biotechnology Information under no. GSE 15366).
Minimal prognostic gene set
We used presently available published data sets to investigate the prognostic correlation with our in vitro–derived hypoxia gene set and to develop an optimized, reduced prognostic gene list (Fig. 1). To test whether the overall expression pattern of these 265 hypoxia genes was significantly related to the prognostic factor considered for each of the three training data sets, the global test of Goeman et al. was used (19). This resulted in a significant enrichment of the hypoxia gene set for all three training sets (P = 0.03595 for Boyault, P < 0.00001 for Lee, and P = 0.0064 for Wurmbach).
Next, when keeping only the significant genes with a z score of >1, 130 genes remained for the data set of Lee et al., 43 genes for Boyault et al., and 58 genes for Wurmbach et al. Finally, genes for which the direction of altered expression did not correspond to the direction observed in vivo were removed. With this approach, we were able to downsize our hypoxia gene set to seven genes, the hypoxia signature of which was found to overlap between the three training data sets (Fig. 1).
Of the genes in our hypoxia signature, four genes were upregulated and three downregulated under hypoxia (Table 2). When we compared the expression of these seven genes in different normal tissues, we found that they have a unique expression pattern in liver not found in other tissues (Supplementary Fig. S3). Furthermore, we tested our seven-gene set related to HCC in other malignancies, such as colon carcinoma and colorectal liver metastasis (GSE3294, GSE 2630, and GSE6988). The seven-gene set had a poor prognostic performance in these data sets (data not shown).
Chronic hypoxia as a prognostic marker
The seven genes were used to define the hypoxia score:
UP, the in vivo upregulated genes (n = 4), and DOWN, the in vivo downregulated genes (n = 3; Table 2).
This hypoxia score was used to classify patients in the three training sets and to calculate the area under the ROC curve (AUC), to assess the predictive performance in all data sets. The hypoxia score based on these seven genes could significantly divide patients into those with or without vascular invasion (Wurmbach; AUC, 88.9%), with a FAL index of >0.128 and ≤0.128 (Boyault; AUC, 72.8%), and with cluster A and cluster B gene expression (Lee; AUC, 84.9%; Fig. 2A). To independently test the performance of our hypoxia score, we used the Chiang data set as validation set, with the BCLC classification as prognostic characteristic. The seven genes significantly separate the BCLC group 0/A/B from group C (AUC, 91%) as well as the group 0/A from B/C (AUC, 71.5%; Fig. 2B).
Survival and early recurrence
With the development of the hypoxia score, we were able to test whether the score correlates with survival and recurrence. We conducted a retrospective survival analysis on 135 patients of the study by Lee et al. (MedCalc Software, version 11.0.1). The patient data is summarized in Table 3. We first determined the Cox proportional hazards ratio for survival because our hypoxia score is a continuous variable. Indeed, the hypoxia score significantly increased the risk of death (hazard ratio, 1.39; 95% confidence interval, 1.09-1.76; P = 0.007). If we used a cutoff value of 0.35 for the hypoxia score (log rank test; P = 0.0018), we were able to show significant differences in survival in 135 patients with a Kaplan-Meier survival curve (Fig. 3A). The median survival for patients with a hypoxia score of >0.35 (n = 42) was 307 days, whereas the median survival for patients with a hypoxia score of ≤0.35 (n = 93) was 1,602 days (P = 0.002). For recurrence in HCC patients, it has been suggested to make a differentiation between early recurrence (<2 years) and late recurrence (>2 years; refs. 26, 27). Early recurrence is the result of dissemination of the primary tumor and tumor characteristics determine the risk of recurrence. On the other hand, recurrence after 2 years is usually a second primary tumor that arises in a cirrhotic liver and has no relation with the first tumor. Because our hypoxia score is determined on the tumor tissue itself, we tested if it could predict early recurrence. We calculated a significant Cox proportional hazard ratio of 2.92 (95% confidence interval, 1.33-6.39; P = 0.007), which means that with an elevation of the hypoxia score with 0.1 point, the risk of developing a recurrence is 29.2% higher. Again, when we use a cutoff of 0.35 for the hypoxia score, the Kaplan-Meier curve shows a significant difference in early recurrence (P = 0.005; Fig. 3B).
We also conducted a multivariate statistical analysis to determine the performance of the hypoxia score compared with other variables, such as AFP level (≤300 versus >300), tumor size (≤5 versus >5 cm), and differentiation grade (1-2 versus 3-4; Table 4). For the 135 samples, it was already shown that the clustering (A versus B and hepatocyte versus hepatoblast) has a very strong correlation with survival and recurrence (both P < 0.001; ref. 15). Multivariate analysis with Cox proportional hazards regression, including the hypoxia score, again showed that clustering based on A/B/hepatoblast is the only significant predictor of both poor survival and early recurrence (P = 0.009 and P < 0.001). However, the differentiation between A and B is based on 406 genes and the differentiation between a hepatoblast or hepatocyte phenotype is based on 941 genes. Furthermore, different samples are needed to perform clustering based on these genes. In contrast, our hypoxia score is based on seven genes and can be calculated for an individual patient. We therefore explored the performance of the hypoxia score without integration of clustering. The hypoxia score is the only significant variable for both survival (P = 0.02) and early recurrence (P = 0.008); this correlation was, however, not as strong as that found for the original clustering (A/B/hepatoblast), but required a much larger set of genes and is only valid in that single study.
Both in the study of Boyault et al. as well as in the study by Chiang et al., the authors divided their patients into different molecular subgroups. Using their classifications, we found that the hypoxia score corresponded with the subgroups that had a poor prognosis (Fig. 4A and B). This subgroup is associated with proliferation and mTOR signaling and chromosomal instability and is seen more frequently in patients with advanced HCC.
The diagnosis for HCC is currently based on a combination of imaging techniques, serum markers, and histology (28, 29); and for a reliable diagnosis, a multidisciplinary approach is required. There are, however, no clear criteria available that could accurately predict the prognosis and several approaches have been attempted to develop rational criteria. Great progress has been made developing clinically relevant cancer signatures especially for breast cancer (30, 31) and, to a lesser extent, for colon cancer (32, 33). For patients with HCC, a similar method has not been developed until now. Recently, several studies have tried to identify gene sets with prognostic or diagnostic relevance by microarray analysis (34–36). Prognosis is a key factor in the selection of a therapy and therefore molecular classification has become increasingly important for the development of targeted therapy such as the multikinase inhibitor sorafenib (37). Although a lot of progress has been made in the molecular classification of HCC over recent years, it is difficult to transpose the results from one study to another study or patients.
The rationale behind this study was to start from the hypothesis that during the development of a solid tumor such as HCC, there is a chronic hypoxic condition which is different from acute hypoxia. Hypoxia is a well-known characteristic of solid tumors and has an effect on the aggressiveness of tumors, induces neoangiogenesis, anaerobic metabolism, and promotes invasiveness.
Central to this study is the working hypothesis that chronic exposure to hypoxia leads to an adaptive gene expression profile which influences the aggressive behavior of the tumor cells. To test our hypothesis independently of patient selection and variability, we decided to start from cell cultures. Here, we did not want a concomitant role for mutations in the p53 gene or liver cell lines that express fragments of hepatitis B or C (38). Therefore, we selected the human HepG2 cell line, although it has hepatoblastoma characteristics (39). The second selection was made to study hypoxia at 2% O2, which is well below the oxygen tension found in normal livers (40). Additional experiments at different oxygen tensions support the selection of 2% O2 for 72 hours in vitro (Supplementary Figs. S1 and S2; Table 4). Subsequently, the gene expression in this model was compared with gene expression determined by microarrays in different groups of HCC patients to identify the genes relevant for HCC and not solely for hepatoblastoma.
In contrast with some of the previous studies on HCC, we did not limit the number of genes we wanted to investigate by a priori selection, and used the Agilent 44 K microarray, which covers all the known genes. Although the dynamics of gene expression indicate that after an adaptation period of 72 hours, the gene expression is not as strongly altered as during the first 24 hours, we still found that 8% of the genes were significantly changed at 72 hours.
Starting with the group of 265 highly significant genes extracted from the microarray study of HepG2 cells, we followed a sequence of analytic steps and developed a very robust seven-gene prognostic signature using the method of Goeman et al. (19). This seven gene prognostic set was applied to four data sets of HCC patients. It could significantly divide patients with different known risk factors into the three studies (training sets) used to develop the seven-gene prognostic signature. Furthermore, in the fourth study (17), this set could significantly separate the BCLC group 0/A/B from C or BCLC group 0/A from B/C. More importantly, the hypoxia score that can be obtained from the seven-gene set significantly correlates with survival and early recurrence.
Two other important microarray studies on HCC with prognostic value have also been reported (34, 41). These studies, however, concentrated on the surrounding liver tissue and showed a correlation with inflammatory pathways and tumor recurrence, whereas our seven-gene set is applicable on the tumor tissue itself. Indirect evidence suggests that our seven-gene set is liver-specific. First of all, the pattern of gene expression in the liver of these seven genes was different compared with other tissues (Appendix: Fig. 3). Second, the hypoxia score based on our seven-gene set had a poor performance in data sets from other malignancies. Third, if we repeat the same experimental sequence with gastrointestinal cancer cell lines, such as from colon and pancreas, we identified other sets of differently expressed genes (Supplementary Fig. S4).
All seven genes are, to a certain extent, related to either hypoxia, proliferation, or metabolism (see Supplementary Materials, Functions and interactions of the seven genes in relation to hypoxia). CCNG2 is a regulator of the cell cycle and is upregulated in response to DNA damage (42). Also WDR45L and RCL1 play a role in proliferation and RNA preprocessing, although the exact function of these two genes are unknown (43). EGLN3 and ERO1L are both upregulated in response to hypoxia. EGLN3 has been described as the main actor in the response to chronic hypoxia and the regulation of Hif1α (44). ERO1L is necessary in the proper protein folding under hypoxic conditions (45). FGF21 and MAT1A are downregulated in response to hypoxia and are involved in metabolic processes. FGF21 is involved in energy homeostasis, it has an endocrine function regulating glucose uptake in adipocytes (46). MAT1A is unbearable in normal liver function because it regulates most methylation processes (47).
As a clinical tool, the hypoxia score seems promising. It shows good correlation with survival as well as early recurrence in patients. It is a good alternative for the classification by Lee et al. (14, 15), with the advantage of requiring a much smaller gene set. In multivariate analyses, our hypoxia score had a better performance than, for example, tumor size and differentiation grade—two variables that are currently used to estimate the aggressiveness of tumors. The identification of these seven genes is the next step in predicting prognosis, and although preliminary, it might influence therapeutic decision-making in the future. Besides, it can be used in randomized trials for risk stratification.
Our findings might also have important implications for the use of antiangiogenic therapies. Nowadays, both chemoembolization and the anti–vascular endothelial growth factor receptor kinase sorafenib are standard treatments for patients with intermediate or advanced HCC. It has been previously noted that resistance emerges after antiangiogenic therapies, leading to untreatable aggressive tumors (48). Our findings on the role of hypoxia lead to the reconsideration of treatment strategies, but more importantly, our gene set could help in the search for new therapeutic targets.
In the last few years, the molecular classification of HCC has attracted a lot of attention. Based on gene expression, patients could be classified into several subgroups such as the β-catenin subgroup, the proliferation subgroup, and the inflammation subgroup. The exact prognostic and therapeutic implications of this categorization are still unclear. We analyzed our hypoxia signature in the studies of Chiang and Boyault, and there was a clear correlation with the proliferation cluster. This cluster consists of genes related to the mTOR pathway and several cell cycle genes, such as cyclins. In the proliferation cluster defined by Chiang, there was a significant increase in pRPS6 staining, indicating aberrant mTOR signaling. The mTOR signaling pathway regulates cell growth, cell proliferation, protein transcription, and survival by orchestrating several upstream signals. Recently, an important role for the mTOR pathway in HCC was shown (49, 50). Also, several studies have shown an interaction between mTOR and hypoxia (or HIF1A). There is evidence for an oxygen-independent induction of HIF1A by mTOR signaling and upregulation of targets such as vascular endothelial growth factor (51, 52). These facts warrant further examinations of the hypoxia score in combination with targeting the mTOR pathway in human studies.
In conclusion, we have identified a seven-gene set that correlates with prognosis in HCC patients. We were able to extract this very small gene set by combining four available microarray data and because we had identified in advance an in vitro gene set.
Our findings have potential implications in several areas:
Because our prognostic signature relies on a limited number of genes, its application should be broadly possible using generally available techniques such as real time PCR. However, it will be necessary to determine the criteria which a sample should fulfill and we recommend further cross-validation versus histology.
As a prognostic indicator, it can influence the therapeutic options that are available for a patient. Therefore, this signature needs to be further validated in new prospective studies to confirm its use.
The seven-gene set we have identified is based on an in vitro selection of genes that were altered during chronic hypoxia. These results suggest a role for hypoxia in the aggressive behavior of HCCs, although it still has to be proven that the gene expression of these seven genes in aggressive tumors is the direct consequence of chronic hypoxia.
The hypoxia score could possibly be further improved with genes that can be identified with the use of other cell lines.
The method we used to identify this limited gene set—the combination of a cell culture model and the global test method—could also be applied to other tumors. With this hypothesis-driven method, it is easier to extract the most important genes out of the large amount of information from the microarray technique, even when the number of patients that can be recruited into a single study is limited, as for HCC patients.
Although in vitro studies are never fully representative for the situation as it develops in an organ, the validation using four clinical data sets proves the value of our study beyond theoretical objections. These findings have prognostic implications for patients with HCC, and therefore, could be incorporated into the molecular classification of HCC.
Disclosure of Potential Conflicts of Interest
No potential conflicts of interest were disclosed.
The authors thank Dr. I. Bockx and P. Windmolders for their expert technical assistance, and Prof. J. Fevery and Prof. P. Van Hummelen for their support.
Grant Support: H. van Malenstein and A. Daemen are research assistants for the Fund for Scientific Research-Flanders (FWO-Vlaanderen). O. Gevaert is supported by the Institute for the Promotion of Innovation through Science and Technology in Flanders (IWT-Vlaanderen).
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.
- Received December 16, 2009.
- Revision received May 26, 2010.
- Accepted June 23, 2010.