Integrative Analysis of Hypoxia-Associated Signature in Pan-Cancer.

Hypoxia is serving crucial roles in cancers. This study aims to comprehensively analyze the molecular features and clinical relevance of a well-defined hypoxia-associated signature in pan-cancer using multi-omics data. Data were acquired from TCGA, CCLE, GDSC, and GEO. RNA expression pattern, copy number variation (CNV), methylation, and mutation of the signature were analyzed. The majority of the 15 genes were upregulated in cancer tissues compared with normal tissue, and RNA expression was negatively associated with methylation level. CNV occurred in almost all the cancers, whereas mutation frequency was low across different cancer types. The signature was also closely related to cancer hallmarks and cancer-related metabolism pathways. NDRG1 was upregulated in kidney cancer tissues as indicated by immunohistochemistry. Besides, most of the 15 genes were risk factors for patients' overall survival. Our results provide a valuable resource that will guide both mechanistic and therapeutic analyses of the hypoxia signature in cancers.

Introduction

Cancer is a major burden worldwide and is also a leading cause of death (Siegel et al., 2019). It is a multi-step disease and characterized by complex biological features. Hypoxia, which is defined as reduced oxygen levels, is one of the hallmarks of malignant tumors (Giaccia and Schipani, 2010). Most of the tumor cells have insufficient blood supply due to rapidly increasing cell proliferation and increasing tumor mass, resulting in a hypoxic tumor microenvironment (Foster et al., 2014), and hypoxia has a crucial role in shaping the behavior of cancers. It has been reported that hypoxia could cause tumor cells to acquire more aggressive phenotypes. Hypoxia can also induce epithelial-mesenchymal transition (EMT) (Muz et al., 2015), which is closely related to tumor migration and invasion (Mittal, 2018). Moreover, hypoxia promotes self-renewal capability and suppresses the differentiation of stem cells in certain tumors (Li et al., 2009). Besides, hypoxia is associated with abnormal vasculature (Carmeliet and Jain, 2011), metastasis, and radio-chemotherapy resistance (Muz et al., 2015).

Some approaches, including both invasive and noninvasive ones, have been developed to measure tumor oxygen levels in tumors (Vaupel et al., 2007). However, these procedures have limited applications. Therefore a robust and clinically applicable hypoxia gene signature has been introduced (Buffa et al., 2010). The signature is composed of 15 genes with known function: VEGF (Vascular endothelial growth factor A), PGAM1 (Phosphoglycerate mutase 1), ENO1 (Enolase 1), LDHA (Lactate dehydrogenase A), TPI1 (Triosephosphate isomerase 1), P4HA1 (Prolyl 4-hydroxylase, a-polypeptide I), MRPS17 (Mitochondrial ribosomal protein S17), ADM (Adrenomedullin), NDRG1 (N-myc downstream regulated 1), TUBB6 (Tubulin, b6), ALDOA (Aldolase A, fructose-bisphosphate), MIF (Macrophage migration inhibitory factor), SLC2A1 (Solute carrier family 2, member 1), CDKN3 (Cyclin-dependent kinase inhibitor 3), and ACOT7 (Acyl-CoA thioesterase 7). The signature is the best performer, proved by recent extensive researches that have evaluated and validated the power of different hypoxia signatures (Ye et al., 2019).

Despite the presence of studies about the hypoxia-associated signature, there is still a lack of researches about the expression, copy number variation (CNV), methylation, or mutations of these hypoxia-related genes in tumors. In this study, we aimed to comprehensively identify the molecular features and clinical relevance of these genes across a wide variety of cancer types using multi-omics data.

Results

Genomic Alterations of Hypoxia-Associated Signature in Different Cancer Tissues

A common and compact hypoxia metagene has been identified and selected (Buffa et al., 2010). We analyzed the molecular landscape of the 15-gene hypoxia signature in various cancers to provide basic information for future studies. We analyzed CNV, methylation, and mutation data from TCGA, CCLE, GDSC, and GEO (Database: TCGA, CCLE, GDSC. and GEO). As shown in Figure 1A, CNV was observed in all the hypoxia-related genes, among which NDRG1 and MRPS17 had the highest gain frequency across the 33 cancer types. The average mutation frequency of the 15 genes was low, ranging from 0 to 0.25 with an average value of 0.01 (Figure 1B).

Figure 1

Genomic Alterations of Hypoxia-Associated Signature in Different Cancer Tissues

(A) CNV of the hypoxia-associated signature in cancers. The upper half of each rectangle exhibits the loss frequency, and the bottom half exhibits the gain frequency.

(B) Mutation frequency of the signature among different cancer types.

(C) Heatmap of the methylation level of the 15 genes in 33 cancer types.

(D) Log2 fold change of the genes in cancer tissues versus normal tissues. Both TCGA (left) and GEO (right) data were used.

(E) Correlation plot of RNA expression and methylation level of each gene in each cancer type. The white grid means p > 0.05

(F and G) (F) Violin plot of each gene in normal tissues and cancer tissues in LIHC. (G) IHC result of NDRG1 expression in kidney cancer tissues.

Data were represented as mean ± SEM. ∗∗p < 0.01; ∗∗∗p < 0.001.

The methylation data from TCGA was also analyzed. The beta value was used to quantify methylation status. Methylation was observed in most of the 15 genes with the beta value range being 0–0.7 (Figure 1C). ACOT7, TUBB6, and NDRG1 were the most highly methylated genes as shown in Figure 1C. As gene methylation is closely correlated with mRNA expression, gene expression data from TCGA and GEO were then analyzed. Gene expression of normal tissues and cancer tissues was compared. Interestingly, both TCGA and GEO data showed that lung, kidney, and uterus cancer tissues had general higher hypoxia-associated gene expression compared with normal tissues (Figure 1D). Linear models were also used to adjust for tumor purity, ploidy, age, sex, and race (Table S1). To assess the relationship between methylation and gene expression, the correlation was analyzed, and it indicated that there was a reverse relationship between RNA expression and methylation level in most of the 33 cancer types (Figure 1E). In particular, RNA expression data from LUSC were further analyzed, and it showed that 14 of the 15 genes were significantly upregulated in cancer tissues compared with normal tissue with p < 0.001 (Figure 1F). As NDRG1 was significantly upregulated in kidney cancer tissues in both GEO and TCGA data, we validated the results in tissues by immunohistochemistry (IHC). It showed that NDRG1 was significantly higher in kidney cancer tissues (Figure 1G).

Genomic Alterations of Hypoxia-Associated Signature in Different Cancer Cell Lines

To further understand the molecular landscape of this 15-gene signature (Buffa et al., 2010), we also analyzed the data of cancer cell lines from CCLE and GDSC. Data of 715 cell lines of 24 cancer types in CCLE and 988 cell lines of 30 cancer types in GDSC were identified. Unexpectedly, CNV frequency of most cell lines in CCLE was quite low ranging from 0 to 0.65, whereas CNV frequency in GDSC was high on average ranging from 0 to 1 (Figures 2A and 2B).

Figure 2

Genomic Alterations of Hypoxia-Associated Signature in Different Cancer Cell Lines

(A–G) CNV frequency of the 15 genes in cell lines of CCLE (A) and GDSC (B). Mutation frequency of the 15 genes in cell lines of CCLE (C) and GDSC (D). Log2 fold change of genes' expression in metastatic cell lines versus primary cell lines in CCLE (E) and GDSC (F). Heatmap of the correlation between CNV of the signature and patient's OS. Red shows that a higher copy number was related to a worse OS, and blue shows a link with a better OS (G).

UCEC exhibited the highest mutation frequency among all the cancer types, which is consistent with the fact that UCEC has a higher global mutation burden. The mutation frequency of the genes in CCLE was similar to that of GDSC, both ranging from 0 to less than 0.5 (Figures 2C and 2D).

Hypoxia is a potent microenvironmental factor promoting metastatic progression (Rankin and Giaccia, 2016). To evaluate and validate the potential role of the 15-gene signature in metastasis, we further assessed the expression pattern of the signature in metastatic and primary cancer cell lines. The results revealed that the expression of these genes was on an average equal to or higher in metastatic cell lines in most cancer types (Figures 2E and 2F), indicating a potential role of these genes in cancer metastasis.

Unexpectedly, there was some inconsistency between CCLE and GDSC databases, which could be due to the use of different array platforms, experimental artifacts, measurement error, or divergence of the cell lines (Haibe-Kains et al., 2013). Extra attention should be paid for future researchers when using these databases.

Hypoxia-Associated Signature and Cancer Pathways

We hypothesized that a variety of pathways interacted with the hypoxia-related signature. Moreover, genes do not function in isolation (Li et al., 2019); we thus first investigated the expression correlation among the hypoxia gene signature. To our expectations, the majority of these genes were positively related to each other (Figure 3A). In particular, TPI1 and ENO1 were most closely correlated. GSVA analysis was adopted to assess the involvement of the signature in cancer-related hallmarks. Cancer hallmark gene sets from MSigDB collections, which summarize and represent specific well-defined biological states or processes (Subramanian et al., 2005), were applied in GSVA analysis. We then examined the correlation between the expression of individual genes of the signature and the activity of cancer hallmark-related pathways. We found that the expression of the signature was correlated with the activation (correlation coefficient > 0, p < 0.01) or inhibition (correlation coefficient < 0, p < 0.01) of multiple oncogenic pathways (Figure 3B). E2F target pathway, MYC target pathway, glycolysis pathway, hypoxia pathway, and unfolded protein response pathway were the most frequently and highly activated pathways, whereas myogenesis pathway, hedgehog signaling pathway, heme metabolism pathway, bile acid (BA) metabolism pathway, and Kras signaling pathway were the most highly inhibited pathways (Figure 3B). Besides, the activated pathways greatly outnumber the inhibited pathways (Figure 3C), indicating a cancer promotion role of these genes.

Figure 3

Hypoxia-Associated Signature and Cancer Pathways

(A) RNA correlation plots of the 15 hypoxia-associated genes.

(B) Network of cancer hallmarks and the hypoxia-associated signature.

(C) Number of activated and inhibited cancer hallmarks of each gene. The red sector represents activated hallmarks, and the blue sector represents inhibited hallmarks.

(D) Network of metabolites and the hypoxia-associated signature.

(E) Log2 fold change of cellular fitness after knocking out each gene in cells.

(F) Correlation between clinically actionable genes and hypoxia-associated genes.

(G) Average correlations between the gene of interest and functional states in different cancers.

To further study the role of hypoxia-related genes in the regulation of metabolites, metabolomics data were retrieved from CCLE. The correlativity between metabolomics and transcriptomics data of the 15-gene signature were then analyzed. RNA-metabolite pairs with |R|≥0.6 and p < 0.05 were shown in Figure 3D. We further investigated the metabolic pathways that were engaged in the network. It showed that about 30 metabolism pathways were included in the network. Interestingly, most of these pathways were identified as common cancer-related pathways like BA biosynthesis, biopterin metabolism, and lysine metabolism (Table S2).

Hypoxia promotes cancer malignancy, so the 15-gene hypoxia signature is supposed to contribute to cancer cells' progression. We adopted cellular fitness analysis, which could predict the functions of genes in various biological processes (Kurata and Lin, 2018; Perez et al., 2017), to validate the roles of the 15 genes in normoxic cancer cells. Cell fitness data from Sanger genome-wide CRISPR essentialities were used in our study to evaluate the function of these genes in cells. These essentiality profiles provide a measure of the loss or gain of cellular fitness resulting from knocking out a gene from targeted disruption via a single guide RNA. There was a dramatic decrease in cellular fitness after knocking out ALDOA, LDHA, TPI1, and PGAM1 in most cancer types, whereas only a slight increase in cellular fitness after knocking out TUBB6, ACOT7, NDRG1, CDKN3, VEGFA, and MIF (Figure 3E). These data suggested that the 15-gene signature mainly acts as cancer promoters. In addition, clinically actionable genes were collected and downloaded from a previously published article (Van Allen et al., 2014) to study the potential clinical and medical use of the 15-gene hypoxia signature. Interactions between these genes with the hypoxia-associated signature were evaluated, and 80 of 121 genes were found to be associated with the signature (Figure 3F), which indicated that the hypoxia signature might interact with these genes and be used as drug targets. Moreover, correlations between the gene list and functional states in different cancers were analyzed. As shown in Figure 3G, most of the genes were positively related to malignant behaviors, especially hypoxia, invasion, and EMT.

Prognostic Role of the Hypoxia-Associated Signature in Cancers

The hypoxia-associated signature was playing important roles in cancers based on the aforementioned findings. Therefore we speculated that the signature would also affect clinical outcomes. Cox analysis was adopted to evaluate the effect of the 15 genes on overall survival (OS) of patients. The results demonstrated that the majority of these genes were risk factors for OS in various cancer types, especially in liver hepatocellular carcinoma (LIHC) and lung adenocarcinoma (LUAD) (Figure 4A). Therefore, LIHC cases, as an example, were divided into two groups according to the median risk score (Risk Score = ). Based on the risk score pattern of the signature, two subgroups of patients with LIHC cancer were identified (Figure 4B). As shown in Figure 4C, cases of a low risk score had a much better OS than those of a high risk score. The data demonstrated that the signature is powerful to predict OS at least for LIHC. Time-dependent receiver operating characteristic analysis was then used to evaluate the accuracy of the signature-based prognostic prediction model in LIHC. Area under the curve = 0.69 showed that the model had excellent predictive value (Figure 4D). The role of P4HA1 in the prognosis of cancers was also analyzed, and it indicated that P4HA1 was a risk factor in more than one-third of all cancers (Figure S1).

Figure 4

Prognostic Role of the Hypoxia-Associated Signature in Cancers

(A) Heatmap of the relationship between the expression of the signature and OS. Red shows that a higher expression was related to a worse OS, and blue shows a link with a better OS.

(B) Risk score value distribution in patients with LIHC.

(C) Survival plot of patients with LIHC according to the risk score.

(D) Receiver operating chanracteristic curve of risk score prediction model.

Discussion

In our study, we comprehensively analyzed multi-omics data of hypoxia-associated signature. Metabolomics, transcriptomics, epigenetics, and proteomics data are integrated to investigate the role of the hypoxia-related signature in cancers in an extensive manner. Multi-omics data enable us to predict novel functional interactions between molecular mediators at multiple levels. Also, these data have the potential to uncover crucial biological observations into hallmarks and pathways that would otherwise not be obvious through single-omics studies. Hypoxia has been paid wide attention recently. Ye et al. classified tumor samples into hypoxia score-high and score-low groups based on the same hypoxia signature and identified molecular alterations that correlated with drug responses to anticancer drugs (Ye et al., 2019). Also, Bhandari et al. quantified hypoxia in various cancer types and found elevated hypoxia associated with increased mutational load (Bhandari et al., 2020), as well as observed widespread hypoxia-associated dysregulation of microRNAs (Bhandari et al., 2019). However, molecular description of the 15-gene hypoxia signature in cancer tissues and cell lines is still lacking. Our study performed a comprehensive analysis of the hypoxia signature, aiming to provide a useful resource for future related researches.

CNVs are nearly ubiquitous in cancer (Heitzer et al., 2016; Zack et al., 2013) and have a great effect on cancer genome as a key type of genomic variations (Heitzer et al., 2016). Genes' CNV landscape is varied across different cancer types, and specific CNVs are related to cancer outcomes (Nibourel et al., 2017). In our study, we found that NDRG1 and MRPS17 are characterized by high gain frequency across the 33 cancer types. Further analysis showed that the copy number gains of NDRG1 and MRPS17 are prone to be risk factors in cancers (Figure 2G). Therefore, CNVs of the two genes might be used to predict prognosis in some cancers.

Not up to expectations, the general mutation frequency of these genes is quite low with a maximum frequency of 0.25, indicating that mutation is not the main alteration of these genes. Methylation is one of the common epigenetic mechanisms that involved in cell proliferation, autophagy, differentiation, and cell cycle in cancer (Pan et al., 2018). We also looked at the methylation level of the 15 genes across all the 33 cancer types. The methylation level of the signature is similar in 30 of the 33 cancer types, demonstrating the homogeneity of the signature in cancers. DNA methylation in the gene body is supposed to increase gene expression, whereas promoters' methylation is negatively related to gene expression (Yang et al., 2014). The negative correlation between RNA expression and methylation status for the majority of the genes reveals that promoter methylation is common for these genes, and methylation inhibitors might be used based on the presumed therapeutic target, particularly in the case of NDRG1, TUBB6, and ACOT7.

Data from TCGA and GEO showed alike RNA profiling in normal tissue and cancer tissue, proving the reliability of our analysis. Also, RNA data from cell lines invalidated the metastasis promoter roles for most of these genes. These data are consistent with previous studies (Bhandari et al., 2019, 2020; Ye et al., 2019) and further proved the roles of hypoxia in carcinogenesis and cancer development process. Interestingly, the RNA expression patterns of COAD are much the same as those of READ, indicating the general histological homogeneity of colon and rectum. On the other hand, the patterns are almost opposite between KIRC and KIRP, both of which are from renal tissues. It illustrates the discrepancy of clear cell cancer and papillary cancer of the kidney.

The majority of these 15 genes are upregulated in cancer tissues compared with normal tissue, implying that these genes could be positively correlated with each other. Data from RNA correlation analysis consistently validated the hypothesis. In particular, TPI1 and ENO1 are the most closely related. Both TPI1 and ENO1 are glycolytic enzymes, and studies have reported that TPI1 and ENO1 are overexpressed in various types of cancers, along with other glycolytic enzyme members like ALDOA and LDHA (Fedorova et al., 2019; Hamaguchi et al., 2008; Migneco et al., 2010). TPI1 and ENO1 would also be functionally correlated. Therefore, future studies about hypoxia or glycolysis could incorporate TPI1 and ENO1 into existing knowledge. Moreover, GSVA analysis also found that these 15 genes are more likely to be oncogenes as the number of activated pathways is much larger than that of inhibited pathways. Furthermore, a lot of important cancer hallmarks, such as apoptosis, p53 pathway, and transforming growth factor-beta signaling, interact with the signature. It is not surprising that the signature promotes cancer, but we listed the activated pathways in our results aiming to provide some useful clues for future studies. To our expectation, hypoxia is significantly related to the majority of these 15 genes. Of interest, there are some overlaps between the signature-related cancer hallmarks and metabolomics pathways among which are glycolysis and BA biosynthesis. Aerobic glycolysis has been closely linked with tumorigenesis for a long time ever since the “Warburg effect” was proposed (Warburg, 1956). Yet, the underlying mechanistic details regarding the causes and effect of the metabolic phenotype remained unclear (Ganapathy-Kanniappan and Geschwind, 2013). Our study might provide some clues to subsequent studies. BA biosynthesis has been serving crucial roles in multiple cancers and a large number of pathways linking BAs to cancer have been identified including oxidative stress with DNA damage and genomic instability, apoptosis, epigenetic factors, activation of nuclear receptors and metabolic and cellular homeostasis, and interactions with and changes of gut microbiota (Di Ciaula et al., 2017). Studies have showed that BA could be linked to glycolysis. BA interruption by BA sequestrants is able to affect Carbohydrate-response element-binding protein (ChREBP) transactivity, subsequently influencing aerobic glycolysis through p53 inhibition (Iizuka, 2017). Also, BAs could potently modulate microbiome (Tian et al., 2020), and microbial communities provide molecules and metabolites, which are essential for glycolysis (Belizario et al., 2018). It has been suggested that hypoxia could repress BA synthesis through inhibition of CYP7A1 expression (Moon et al., 2016). Moreover, hypoxia could negatively regulate BA homeostasis by regulating the biosynthesis, conjugation, secretion, and absorption of BAs via Farnesoid X receptor (FXR) (Fujino et al., 2009). Our results were consistent with those of the aforementioned studies that most of the genes from the 15-gene signature were negatively associated with BA metabolism pathway. Tyrosine, a crucial component of tyrosine metabolism and biopterin metabolism pathways, is widely involved in multiple cancer pathophysiologies. We also found that some genes of the signature were correlated with tyrosine metabolism. Previous studies showed that hypoxia could decrease tyrosine by directly or indirectly upregulating tyrosine hydroxylase (Kumar et al., 2003; Lim et al., 2015). Therefore, tyrosine might be used to prevent hypoxia-induced effects. In fact, it has been proved that tyrosine administration counters hypoxia-induced decrements in learning and memory (Shukitt-Hale et al., 1996). Future researches are warranted to develop drugs targeted on tyrosine to treat cancers.

Cellular fitness has been applied in many fields to uncover new oncogenes (Kurata and Lin, 2018) and study molecular structures (Perez et al., 2017). In our study, cellular fitness was used to measure the potential function of these genes in cells. The results also indicate that most of these genes are essential to maintain normal cancer cellular fitness. It is worth noting that the cellular fitness here was performed in cancer cells, which were cultured under normoxia. Multiple hypoxia-mediated gene regulatory networks will be initiated and produced under lasting or temporary hypoxic situations and finally lead to significant change of cellular functions and behaviors (Chen et al., 2020). Subsequently, irreversible processes develop, which could result in pathophysiological disorder consequences. Therefore, the cells that are cultured in normoxia cannot completely reflect the the functions of the signature in hypoxia. However, these data could at least provide and validate some information about roles of the 15 genes. Also, cells from CCLE and GDSC were cultured under normoxia, therefore data from these cells are not closely related to hypoxia in primary cancers. However, we aimed to provide a molecular landscape of the signature in cancer cells, which could aid future researchers to select suitable hypoxia cell models.

A previous study has found that hypoxia status could affect most clinically actionable genes and experimentally validated the predicted effects of hypoxia on the response to several drugs in cultured cells (Ye et al., 2019). In our study, we further evaluated the correlation of the 15-gene signature and 121 clinically actionable genes in a direct way and also showed that the hypoxia signature is closely related to the clinically actionable genes. Both the studies warrant that hypoxia-targeted therapy is a promising cancer treatment, probably as a component of combination therapy targeting clinically actionable genes.

Prognosis prediction of cancers is crucial, as it can aid the subsequent clinical management of patients. Despite the fact that remarkable improvement has been achieved in cancer research globally over the past several decades, it is still a big challenge to predict patients' prognosis. We also attempted to evaluate the role of the signature in cancer prognosis. The findings demonstrate that the genes are risk factors in about 65% gene-cancer pairs, which accords with the cancer promoter role of these genes in cancers. The risk score based on expression and hazard ratio could greatly separate LIHC cases, and the risk score is powerful to predict the prognosis of LIHC. This evidence warrants further study of hypoxia-associated signature in cancers. It could be also used as an effective prognosis prediction model.

There may be some possible limitations in this study. Data from multiple sources are difficult to be processed and analyzed under a unifying framework given the diversity of methodological work in this area (Wu et al., 2019). It may lead to minor inconsistencies between our study and others'. Besides, some data are collected from retrospective studies in which some critical parameters are not recorded, and remarkable biases could influence selection of controls. Therefore, further prospective clinical studies are needed to validate our findings.

In conclusion, we systematically analyzed the landscape of genetic, epigenetic, and metabolomics alterations and biological and clinical relevance of hypoxia-associated signature in cancers. It helps to better understand the dysregulation of hypoxia in cancer. The findings from our study can be readily applied to further studies and are highly promising for personalized treatment regimens and clinical management.

Limitations of the Study

Cells that were used to analyze the profile of the hypoxia signature were cultured in normoxia. This would not perfectly reflect the situation of cells that are cultured in hypoxia.

Resource Availability

Lead Contact

Qian Zhang, zhangqian1171@126.com.

Materials Availability

The study did not generate new unique reagents.

Data and Code Availability

The datasets and code used and/or analyzed during the current study are available from the corresponding author on reasonable request.

Methods

All methods can be found in the accompanying Transparent Methods supplemental file.