Integrative Analysis Reveals the Landscape of Hypoxia-Inducible Factor (HIF) Family Genes in Pan-Cancer

Inside the cancer microenvironment, reduced O2 concentration, termed as hypoxia, is a common phenotype and leads to cancer progression. However, little is known about how and when those HIF members are dysregulated in distinct cancers. Here, by integrating a full range of data of thousands of patients, we comprehensively analyzed the genetics, epigenetics, and transcriptomic level of HIF genes and further defined pathways triggered by disrupted hypoxia-inducible factors. We reveal the expression landscape of HIF family genes and further demonstrate that copy number variations underlie such dysregulation. Further analysis indicates that HIF genes associate with cancer hallmarks such as cell cycle and DNA damage response. Drug resistance analysis showed that HIF globally impacts drug effectiveness such as docetaxel. In summary, the overall analysis reveals the landscape of HIF genes in pan-cancer and may assist mechanism research about hypoxia.


Introduction
For cells to adapt to hypoxia, they must be able to sense changes in oxygen and respond to them [1,2]. Under such a scenario, hypoxia-inducible factors (HIFs) mediate adaptive physiological responses to hypoxia [3][4][5].
Particularly in cancer cells, the HIF pathway is widely accepted to be enhanced or upregulated [15]. Inside the tumor, reduced O 2 availability is a common phenotype and leads to cancer hallmarks such as angiogenesis, metabolic reprogramming, and epithelial-mesenchymal transition [15]. Yet, little is known about how and when these HIF members are dysregulated. Copy number variations (CNVs), genetic mutation, and promoter methylation will all exert large effects on gene expression; however, how these factors contribute to HIF enhancement is unknown. We hence seek to address the origin of such upregulation and quantify the most enriched pathways linked with such upregulation. Defining enriched pathways of hypoxia might help improve our understanding of how HIFs shape cancer cells and offer the promise of developing targeting drugs.
Here, by integrating multiple data of thousands of patients, we aim to explore how cancer cells reprogram the HIF transcriptional program. We comprehensively analyzed the genetics, epigenetics, and transcriptomic level of HIF genes and further defined pathways triggered by disrupted hypoxia-inducible factors. is proposed work will elucidate key pathways in response to hypoxia that may underlie critical steps in carcinogenesis, genetic instability, tumor progression, and resistance to cancer therapies. Identification of strategies to prevent or reverse these pathways may provide the basis for new approaches to cancer prevention and therapy.

Data Source.
As for the gene expression analysis, methylation analysis, single nucleotide variation analysis, copy number variation analysis, and pathway enrichment analysis, we downloaded data from e Cancer Genome Atlas (TCGA, https://portal.gdc.cancer.gov/). As for the normal tissue gene expression analysis, we downloaded data from GTEx (https://www.gtexportal.org/home/datasets). As for drug response analysis, we downloaded data from GDSC (http://www.cancerrxgene.org/).

Gene Expression Analysis.
As for the differential gene expression analysis, we used mRNA count data of cancer samples and the matched normal samples. We used Deseq2 (http://bioconductor.org/packages/devel/bioc/vignettes/ DESeq2/inst/doc/DESeq2.html) to perform the differential gene expression analysis. We set the cutoff as 2 and 0.05 for the fold change and FDR. We used R to perform the Spearman correlation analysis.

Methylation Analysis and Single Nucleotide Variation
Analysis. As for the differential methylation analysis, we used methylation probe data of cancer samples and the normal samples. We used Student's t test to perform the differential methylation analysis. We set the cutoff as 0.05 for the FDR. e mRNA expression and methylation data were matched according to the barcode of each sample. We used R to perform the Spearman correlation analysis between methylation and mRNA expression. We used maftools (https://github.com/PoisonAlien/maftools) to summarize the single nucleotide variation events (driver mutation).

Copy Number Variation Analysis.
We divided all CNV events into heterozygous and homozygous events, which means the CNV occurs on only one or two chromosomes. We use GISTICS (ftp://ftp.broadinstitute.org/genepattern/ modules_public_server_doc/GISTIC2.pdf) to analyze the CPV events. Only genes with >5% CNV in cancers were shown as a corresponding point on the figure. We used R to perform the Spearman correlation analysis between CNV and mRNA expression.

Pathway Enrichment Analysis.
We calculated the pathway score according to the reverse-phase protein array data in TCGA. RPPA data are centered on the median, and each component is normalized by the STD of all patients to obtain relative protein levels. Gene expression was divided into two groups according to the median value (high and low), and Student's t test was used to determine the difference in the pathway activity score (PAS). e heat map shows genes that are functional (suppressed or activated) in at least 5 cancer types. e pathway_a (red) represents the percentage of cancers that are activated by specific genes and is inhibited similarly to the pathway_a (blue).

Drug Sensitivity Analysis.
To analyze the correlation between gene expression and drug sensitivity, we downloaded drug dose-response curve (AUC) values and regions under the gene expression profile for all cancer cell lines. Fisher's Z transform was used to normalize the transcription level and the Pearson correlation coefficient of the area under the curve. e two-tailed distribution corrected by Bonferroni had a family error rate of z-scores less than 0.025.
e Spearman correlation coefficient of the agent target was compared with the same number of drug-gene pairs computed by random sampling correlation.

Hypoxia-Inducible Factors Are Dysfunctional across
Cancers Which Leads to Worse Patient Outcomes. We used the differential gene analysis of hypoxia-inducible factor family members according to the raw count of these HIF genes from ten cancer types downloaded from e Cancer Genome Atlas. We selected cancers withed matched normal samples including THCA, KIRP, LIHC, STAD, BRCA, COAD, UCEC, BLCA, KIRC, KICH, and PRAD. We found that the majority of the detected HIF genes are significantly changed ( Figure 1). Among them, HIF3A and EPAS1 are significantly downregulated, whereas ATNTL and ARNT2 show minimal expression changes in most cancer types ( Figure 1). We further found that HIF family gene expression varied in cancer subtypes ( Figure S1). Especially in breast cancer, all HIF genes are differentially expressed in each subtype. In addition, we examined whether HIF gene expressions impact the overall survival of cancer patients ( Figure 1). Our result demonstrates that all HIFs had either favorable prognostic or unfavorable prognostic impacts on overall survival. Nevertheless, the prognostic power varied by HIF factor types. Some HIFs, such as HIF1A, ARNT, EPAS1, and HIF3A, were linked with unfavorable outcomes for pancancer, while two HIFs including ARNTL and ARNT2 were associated with favorable outcomes in pan-cancer ( Figure 1).

Heterozygous Amplification Contributes to Upregulated
Hypoxia-Inducible Factors. We also sought to explore the impact of copy number variations (CNVs) on HIF gene expression. We first separated all CNV events into heterozygous or homozygous copy number variations. Heterozygous amplification/deletion was mainly correlated with transcriptional alterations such as upregulated HIFs (e.g., HIF3A, ARNT, and ARNTL). Downregulated HIFs (e.g., HIF1A and ARNT2) were associated with decreased copy numbers (Figure 2 and Figure S2). However, we further examined the effect of homozygous amplification or deletion on mRNA level expression. To our surprise, only homozygous amplification of the ARNT gene associates with its upregulation in only 4 cancers (BLCA, BRCA, LIHC, and UCEC).
is indicates that heterozygous copy number variations underlie the transcriptional changes of HIFs. To confirm our result, we perform the correlation analysis and observe the positive correlation between CNVs and mRNA expression in most HIF genes across cancers.

Mutation Profile of Hypoxia-Inducible Factors.
Based on the somatic mutation data from e Cancer Genome Atlas, mutations of HIFs widely exist (Figure 3). Even for the lowestranked HIFs like ARNTL, the mutation frequencies were around 17%. EPAS2 and HIF3A were mutated in ∼26% of samples. Most HIF mutations distribute in cancers with high mutation loads such as UCEC and STAD (Figure 3(a)). For example, especially in UCEC, 31% of patients show EPAS1 mutation (Figure 3(a)). Among all mutations, missense mutations predominate, and most of mutations are C to G ( Figure S3). To sum up, frequent mutations of HIFs also contribute to the abnormal profile of hypoxia factors in cancers.

Demethylation Contributes to the Abnormal Hypoxia-Inducible Factors.
Since the methylation of the gene promoters can induce the downregulation of genes, we downloaded the promoter methylation of HIF genes in different cancer types. We found the abnormal methylation of the hypoxia genes in multiple cancer types, especially in BRCA, PRAD, and LIHC (Figure 4(a)). Further, we utilized the Spearman correlation analysis to compute the correlation between the gene expression and methylation level of HIF genes across cancer types and observed that the methylation level of HIF1A, ARNT, EPAS1, ARNT2, HIF3A, and ARNTL are significantly negatively correlated with their expression level, respectively (Figure 4(b)), which may partially explain the abnormal expression of these HIF genes in cancer samples.

Abnormal Hypoxia-Inducible Factors Associates with Key
Cancer Hallmarks. We next explore the HIF-related cancer hallmark pathways, including apoptosis, cell cycle, DNA damage response, EMT, hormone AR, hormone ER, PI3K/ AKT, RAS/MAPK, RTK, and TSC/mTOR. We found that upregulated HIF1A, ARNT, EPAS1, ARNT2, HIF3A, and ARNTL associates with the EMT signalling pathway (Figure 5 and Figure S4). Also, this observation is in consistency with previously reported results [16]. We also found previously unreported results. For example, HIF3A, which is reported to promote the cancer metastatic phenotype, was observed to be highly linked with PI3K pathway activation across multiple cancer types. In summary, our systematic analysis not only ensured the previously identified HIFrelated pathways but also shed light on novel potential HIFassociated signalling pathways in human cancers.

Abnormal Hypoxia-Inducible Factors Associates with
Drug Sensitivity. Many clinically actionable genes are targeted by anticancer drugs as identified in the GDSC project (Methods). Hypoxia in the tumor microenvironment is a common phenotype in cancers, and hence, targeting those hypoxia-related genes may help improve the survival of patients. To further assess how hypoxia impacts drug response, we computed the Spearman coefficient between drug sensitivity of 265 clinically used agents and mRNA expression of 6 HIFs across 1,080 cancer cell lines in GDSC. To our surprise, we observed no significant associations between HIF1A/ARNTL and drug response ( Figure 6). In contrast, responses to 62 drugs are related to EPAS1 expression. For example, docetaxel, a chemotherapy medication used to treat many types of cancer, negatively correlates with EPAS1 expression. Taken together, these data demonstrate hypoxia genes, especially EPAS1, may affect drug sensitivity.

Discussion
Hypoxia in the tumor microenvironment can impact the cancer cell phenotypes [17]. Although previous studies have examined the role of HIF genes in some cancers, a systemlevel analysis is still lacking. HIF genes in cancers are highly context-dependent, making the role of HIFs complicated and changeable in distinct cancers, thus impeding the effective clinical utility of HIFs. Here, by integrating multiple data across cancer types, we seek to explore the landscape of hypoxia-inducible factors (HIFs) family genes.
First, we report the transcriptional landscape of HIF family genes in cancers. Hypoxia in the tumor microenvironment has been linked with cancer development, and  Journal of Oncology hypoxia status is highly associated with effectiveness of anticancer drugs [18]. Our study presents the expression landscape of HIF genes, where most of them are dysregulated across cancer types. For example, hypoxia-inducible factor 3-alpha (HIF3A) is globally downregulated in cancers especially in breast cancer. is gene also shows a high rate (26%) of mutation in all samples. e pathway analysis of this gene shows its association with activation of EMT and PI3K, indicating this gene might serve as an important target in the clinical settings. Our system-level analysis represents a systematic computation of the widespread alterations of hypoxia-inducible factors in thousands of tumors. Further, our research shows that HIF genes are associated with cancer hallmarks. Hypoxia contributes to the suppressive tumor microenvironment in many cancers by activating multiple pathways that allow tumor cells to escape from the innate and adaptive immune defenses [1,19,20]. To better understand how hypoxia impacts the pathways inside the tumor, we selected ten hallmark pathways. Our results demonstrated that most HIF genes negatively correlates with DNA damage response and cell cycle. is indicates that reduced O 2 in the microenvironment may contribute to DNA damage. Our paper highlights the potentially linked signalling pathways of HIFs across a broad spectrum of cancers.
Finally, our study presents the interaction and correlation between the clinically actionable genes and HIFs, indicating that reprogramming the O 2 concentration in the tumor microenvironment should be considered in cancer therapy [21]. In particular, endothelial PAS domain protein 1 (EPAS1) gene expression associates with responses to 62 drugs.
is indicates that enhancing/ blocking EPAS1 may serve as a promising combinational strategy in the clinical setting. However, we observe no significant associations between drug response and HIF1A/ARNTL gene expression, although they are core hypoxia regulators. In a nutshell, our paper provides more evidence that further efforts should be made to combine targeting hypoxia with existed drugs.
Our paper does have some limitations. ese results should be validated in different patient cohorts. For example, future research should look at whether ARNTL can predict the outcome of patients. Second, e Cancer Genome Atlas does not document the synonymous mutation data. It is still unclear about the synonymous mutation of HIF genes and how they impact patients' prognosis.
In summary, our study reveals the expression landscape of HIF family genes and further addresses the origin, copy number variations, of such dysregulation. We also found that HIF genes associates with cancer hallmarks such as cell cycle and DNA damage response. Further drug resistance analysis showed that HIF globally impacts drug effectiveness such as docetaxel. ese findings provide new insights into  Percent -20 0 20 4 Figure 5: e HIF genes are widely associated with the hallmark cancer pathways. e heat map with function (inhibition or activation) in at least 5 cancer types. Pathway_a represents the potential activation of this pathway, and pathway_i represents the potential inhibition of this pathway.
cancer hypoxia and unravel new mechanisms of HIF genes that may be further explored in the future.
Data Availability e dataset supporting the conclusions of this article is included within the article.

Conflicts of Interest
e authors declare that they have no conflicts of interest.

Authors' Contributions
Qian-Kun Yang and Xue-Xin Wang contributed equally. YQK and NN designed and conducted this research. YQK, WXX, WY, and NN analyzed and interpreted the data. YQK wrote the paper. NN reviewed the manuscript and worked on the language.      Figure S3: a review of HIF gene variants in human cancers. (a) Variant type, (b) variant classification, (c) single nucleotide variation, (d) variation per sample, (e) mutation classification, and (f ) mutation of HIF gene between cancer types frequency. Figure S4: the HIF genes are widely associated with the hallmark cancer pathways in different cancer types. e pie plots represent the percentage of cancer samples that are correlated with HIF gene expression. (Supplementary Materials)