Comprehensive Analysis of The Signicance of Hypoxia-Related Genes in Ovarian Cancer

Background: Hypoxia-related genes have been reported to play important roles in a variety of cancers. However, their roles in ovarian cancer (OC) have remained unknown. The aim of our research was to explore the signicance of hypoxia-related genes in OC patients. Methods: In this study, 15 hypoxia-related genes were screened from The Cancer Genome Atlas (TCGA) database to group the ovarian cancer patients using the consensus clustering method. Principal component analysis (PCA) was performed to calculate the hypoxia score for each patient to quantify the hypoxic status. Results: The OC patients from TCGA-OV dataset were divided into two distinct hypoxia statuses (cluster.A and cluster.B) based on the expression level of the 15 hypoxia-related genes. Most hypoxia-related genes were expressed more highly in the cluster.A group than in the cluster.B group. We also found that patients in the cluster.A group exhibited higher expression of immune checkpoint-related genes, epithelial-mesenchymal transition-related genes, and immune activation-related genes, as well as elevated immune inltrates. PCA algorithm indicated that patients in the cluster.A group had higher hypoxia scores than that in in the cluster.B group. Conclusions: In summary, our research elucidated the vital role of hypoxia-related genes in immune inltrates of OC. Our investigation of hypoxic status may be able to improve the ecacy of immunotherapy for OC.


Introduction
Ovarian cancer (OC) is one of the most serious threats to female health. Currently, surgery combined with chemotherapy is the main treatment strategy. However, 70% of patients relapse within 3 years and median progression-free survival (PFS) is only 12 to 18 months [1,2]. Therefore, it is important to explore new therapeutic strategies to improve clinical e cacy and delay recurrence. Immune checkpoint inhibitors (ICI) targeting programmed death protein 1 (PD1), programmed death ligand (PD-L1), and cytotoxic T lymphocyte-associated protein 4 (CTLA-4) have been approved for the treatment of a variety of malignancies. Although PD-L1 is highly expressed in epithelial ovarian cancer (EOC), ICI monotherapy for PD-1/PD-L1 or CTLA-4 has only shown a response rate of 5.9-15% [3]. Therefore, the search for predictive biomarkers is crucial to improve the effectiveness of ICI treatment.
The progression and metastasis of tumors mainly depend on the bidirectional interaction between tumor cells and their microenvironment. The tumor microenvironment (TME) includes immune cells, broblasts, endothelial cells, mesenchymal stem cells, extracellular matrix, blood vessels, and cytokines. Hypoxia is one of the essential features of the TME [4]. The limitation of oxygen diffusion and the disordered vascular system prevent rapidly developing tumors from obtaining su cient oxygen, resulting in a hypoxic microenvironment [5,6]. Numerous studies have reported that hypoxia plays an important role in the occurrence and progression of a variety of cancers [7][8][9]. For example, Klemba et al. summarized that hypoxia is strongly associated with angiogenesis, chemoresistance, epithelial-to-mesenchymal transition (EMT), prognosis, acquisition of stem-like characters, and immune suppression of ovarian cancer [10].
However, there is still a lack of studies about the role of hypoxia in the progression and immune in ltrates of OC.
In this study, we performed a comprehensive analysis of the functions of 15 hypoxia-related genes in the progression and immune in ltrates of OC based on TCGA-OV dataset from The Cancer Genome Atlas (TCGA) database. The patients in TCGA-OV dataset were classi ed into two distinct hypoxia statuses (cluster.A and cluster.B) using the consensus clustering method based on the expression levels of the 15 hypoxia-related genes. We found that patients in the cluster.A group exhibited higher expression of immune checkpoint-related genes and immune activation-related genes, as well as elevated immune in ltrates, suggesting that they may be more sensitive to ICI treatment.

Consensus clustering of distinct hypoxic statuses based on 15 hypoxia-related genes
Consensus clustering is an algorithm used to identify different subgroups based on resampling. The "ConsensusClusterPlus" package in R was used to divide the OC patients into distinct hypoxic statuses based on the 15 hypoxia-related genes using consensus clustering. To ensure stability of the classi cation, the consensus clustering was repeated 1000 times [15].

Identi cation of differentially expressed genes (DEGs) between distinct hypoxic statuses
The "limma" package in R was used to screen for DEGs between the distinct hypoxic statuses based on the RNA transcriptome pro les from TCGA-OV dataset. P-value < 0.01 and logFC > 1 were selected as the screening criteria.

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) functional enrichment analyses
GO and KEGG enrichment analyses were applied to identify the possible mechanisms and pathways by which the DEGs are involved in the hypoxic status of OC using the "clusterPro ler" package in R software. The results were visualized with an enrichment circle diagram or bubble diagram. . The PCA algorithm maps high-dimensional data to low-dimensional space and retains some features of the original data. To quantify the hypoxic status, we used PCA to calculate the hypoxia score for each patient from TCGA-OV dataset. First, PCA was conducted to distinguish the hypoxic status based on the expression levels of the 15 hypoxia-related genes. Then, the hypoxia score was calculated according to the following formula: hypoxia score = ∑ (PC1 i + PC2 i ), where PC1 represents principal component 1, PC2 represents principal component 2, and i represents the expression levels of the 15 hypoxia-related genes in samples from TCGA-OV dataset [12,13].

Single sample gene set enrichment analysis
Single sample gene set enrichment analysis (ssGSEA) was performed to calculate the abundance of 28 immune cell types in samples from TCGA-OV dataset. First, we obtained the genes representative of each immune cell type. These genes were obtained from the study by Charoentong et al., in which various human immune cell subtypes were stored [17]. Next, we searched for these genes in the RNA transcriptome pro les of TCGA-OV dataset. The enrichment scores calculated by ssGSEA analysis based on the expression levels of the representative genes were used to represent the relative abundance of the 28 immune cell types in each OC sample.

Statistical analysis
Kruskal-Wallis tests were used to analyze differences between groups. All parametric analyses were based on two-tailed tests, the statistical signi cance for which was set at P < 0.05. The Kaplan-Meier curve and log-rank method were used to compare the survival outcome between different hypoxic statuses. The "maftools" package in R software was used to screen genes that were differentially mutated between different hypoxic statuses. The results were visualized with forest plots and cooncoplots. All statistical analyses were performed using R version 4.0.0.

Two distinct hypoxic statuses identi ed by consensus clustering
We performed clustering of the consensus matrix to obtain two matrices, 0 and 1, where 1 represents clustering together and 0 represents complete separation. The consensus clustering method was performed to identify distinct hypoxic statuses based on the 15 hypoxia-related genes using the "ConsensusClusterPlus" package in R software. As shown in Fig. 1A, we found that similar samples were all clustered together in the blue module. However, there was substantial noise in the blue module shown in Fig. 1B-D. Therefore, we believed that the clustering effect was optimal at k = 2, so we divided the ovarian cancer samples into two groups (cluster.A and cluster.B). The cluster.A group contained 192 patients, while the cluster.B group contained 187. Then, we illustrated the levels of differential expression of the 15 hypoxia-related genes between these two clusters using a histogram. With the exception of MRPS17, the other 14 hypoxia-related genes showed higher expression levels in the cluster.A group than in cluster.B (Fig. 1E). PCA indicated that the 15 hypoxia-related genes could completely distinguish the two hypoxic statuses (Fig. 1F).
We then applied ssGSEA to calculate the abundance of 28 immune cell types in OC samples from TCGA-OV dataset. We explored the differential immune cell in ltration between the cluster.A group and the cluster.B group, and the results indicated that patients in the cluster.A group had increased immune cell in ltration compared with those in the cluster.B group (Fig. 1G). The Kaplan-Meier curve showed that there was no signi cant difference in survival outcome between the cluster.A group and the cluster.B group (Fig. 1H). We also compared the expression levels of immune checkpoint-related genes and immune activation-related genes between the cluster.A group and the cluster.B group. Figure 2A shows that the immune checkpoint-related genes CD274, CD80, CD86, HAVCR2, IDO1, PDCD1LG2, and VTCN1 were more highly expressed in the cluster.A group ( Fig. 2A). Moreover, the immune activation-related genes CXCL10, CXCL9, GZMB, PRF1, and TNF had higher expression levels in the cluster.A group (Fig. 2B). Hypoxia can also induce epithelial-mesenchymal transition (EMT) of cancers [18]. In our research, we con rmed that the EMT-related genes COL4A1, PDGFRA, TGFBR2, TWIST1, VIM, and ZEB1 were more highly expressed in the cluster.A group (Fig. 2C).

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) functional enrichment analyses
A total of 227 hypoxia-related DEGs between the two hypoxic statuses were selected. To understand the possible enriched functions and pathways of these DEGs in OC, GO and KEGG functional enrichment analyses were performed. The GO functional enrichment analysis indicated that the genes were mainly enriched in GO:0006954 ~ in ammatory response, GO:0071347 ~ cellular response to interleukin-1, GO:0006955 ~ immune response, GO:0005125 ~ cytokine activity, and GO:0071356 ~ cellular response to tumor necrosis factor, all of which are related to in ammatory and immune responses (Fig. 3A). Moreover, the KEGG functional enrichment analysis indicated that the genes were mainly enriched in cytokine-cytokine receptor interaction, ECM-receptor interaction, NOD-like receptor signaling pathway, Jak-STAT signaling pathway, MAPK signaling pathway, and Rap1 signaling pathway (Fig. 3B).

Identi cation of two distinct hypoxia gene statuses
Consensus clustering was performed to divide the OC patients into different genomic subtypes using the "ConsensusClusterPlus" package in R software based on the expression levels of 227 hypoxia-related DEGs through differential analysis between cluster.A and cluster.B, to further con rm the validity of the initial analysis. We also obtained two distinct gene clusters (gene cluster.A and gene cluster.B) (Fig. 4A-D). Figure 3E indicates that the expression levels of most of the 15 hypoxia-related genes were higher in gene cluster.A than in gene cluster.B (Fig. 4E). PCA indicated that the 227 hypoxia-related DEGs could completely distinguish the two gene clusters (Fig. 4F). Figure 4G shows that all of the other 27 of 28 immune cell types were more highly in ltrated in the gene cluster.A group, with the exception of CD56bright natural killer cells. Moreover, the Kaplan-Meier curve indicated that the survival outcome did not differ signi cantly between the gene cluster.A group and gene cluster.B group (Fig. 4H). The expression levels of immune checkpoint-related genes CD274, CD80, CD86, HAVCR2, IDO1, PDCD1, PDCD1LG2, TNFRSF8, and VTCN1 were higher in the gene cluster.A group (Fig. 5A). With the exception of TBX2, all immune activation-related genes CD8A, CXCL10, CXCL9, GZMA, GZMB, IFNG, PRF1, and TNF were overexpressed in the cluster.A group (Fig. 5B). The EMT-related genes COL4A1, PDGFRA, TGFBR2, TWIST1, VIM, and ZEB1 were more highly expressed in the gene cluster.A group (Fig. 5C). All of the above results indicated that the gene cluster.A group corresponded to the cluster.A group, while the gene cluster.B group corresponded to the cluster.B group, suggesting that the consensus clustering method was stable in its grouping. We then considered the gene cluster.A group and gene cluster.B group as two distinct hypoxia gene statuses.

Principal component analysis (PCA) algorithm to quantify the hypoxic status
To quantify the hypoxic status, each patient from TCGA-OV dataset was assigned a hypoxia score using a PCA algorithm. The expression levels of most of the 15 hypoxia-related genes were higher in the highscore group (Fig. 6A). PCA indicated that the 227 hypoxia-related DEGs could completely distinguish the high-and low-score groups (Fig. 6B). We then compared the hypoxia score between the clusters identi ed using either the 15 hypoxia-related genes or the 227 hypoxia-related DEGs. The results showed that the hypoxia score in cluster.A or gene cluster.A was higher than that in cluster.B or gene cluster.B (Fig. 6C and  D). However, the hypoxia score showed no statistically signi cant difference among the different stages and grades (Fig. 6E and F). The relationship among hypoxic status, hypoxia gene status, and hypoxia score was visualized in a Sankey diagram (Fig. 6G). All of the 28 immune cell types were more highly in ltrated in the high-score group (Fig. 6H). The Kaplan-Meier curve indicated that the survival outcome did not differ signi cantly between the high-and low-score groups, which was in line with our expectations (Fig. 6I). Roughly similar to previous results, the expression levels of immune checkpoint-related genes (Fig. 7A), immune activation-related genes (Fig. 7B), and EMT-related genes (Fig. 7C) were higher in the high-score group than in the low-score group. Finally, considering that heterologous proteins produced by gene mutations can lead to an increase of immune cell in ltration, we analyzed the different gene mutations between the high-and low-score groups with different immune in ltrate statuses using the "maftools" package in R software, and visualized the results with forest plots (Fig. 8A) and cooncoplots (Fig. 8B). The results indicated that the mutation frequencies of BRCA1, PRUNE2, APOB, RELN, MDN1, ZSWIM8, MGAT4C, FAT2, and VWF were higher in the high-score group, while those of HECW1 and RB1 were higher in the low-score group.

Discussion
Increasing evidence has indicated that hypoxia-related genes are associated with the malignant progression of ovarian cancer [19][20][21]. However, most of the research focused on exploring the function of single hypoxia-related genes in OC. There have been fewer reports investigating the signi cance of hypoxia-related signature sets in OC. In this study, we extracted the RNA transcriptome pro les of 15 hypoxia-related genes from TCGA-OV dataset. Then, we classi ed the 379 OC patients into two groups (cluster.A and cluster.B) based on the expression levels of the 15 hypoxia-related genes using the consensus clustering method. We explored signi cant differences in immune in ltration level between the cluster.A and cluster.B groups. Notably, previous studies also identi ed the relationship between hypoxia and immune in ltrates in cancer. For example, Wang et al. revealed that hypoxia can aggravate immune cell in ltration to promote the development of hepatocellular carcinoma (HCC) by inducing myeloid derived growth factor (MYDGF) [22]. Moreover, Xiao et al. established a four-gene tumor microenvironment-related model based on TCGA-GBM dataset. These four genes were reported to potentially be associated with the hypoxic phenotype of glioblastoma [23]. Furthermore, Gong et al. divided the breast cancer patients from TCGA database into two groups (cluster1 and cluster2) based on the expression levels of 13 hypoxia-related genes. The results indicated that the immune in ltrating state differed between cluster1 and cluster2, suggesting the hypoxia can regulate immune in ltrates in breast cancer [24]. Our study further supports the relationship between hypoxia and immune in ltrates in OC, but the potential mechanism behind this should be investigated. We hypothesized that patients in the cluster.A group were more sensitive to immune checkpoint inhibitors (ICIs) therapy, given their higher level of immunoin ltration. We then compared the expression levels of immune checkpoint molecules between the cluster.A and cluster.B groups. Consistent with our expectations, the immune checkpoint molecules had higher expression levels in the cluster.A group than in the cluster.B group. We also found that the immune activation-related genes were upregulated in the cluster.A group, which again met our expectations. EMT is a biological process in which epithelial cells transform into mesenchymal ones through speci c processes. Previous studies reported that hypoxia was involved in the speci c processes of EMT [25][26][27]. Our research also supported this, in that patients in the cluster.A group with an activated hypoxic status showed higher expression levels of EMT-related genes.
To further con rm the hypoxic status, consensus clustering was performed to obtain two distinct hypoxia gene statuses (gene cluster.A and gene cluster.B) based on the 227 hypoxia-related DEGs. We found that the gene cluster.A group corresponded to the cluster.A group. The patients in the gene cluster.A group exhibited higher expression of immune checkpoint-related genes, EMT-related genes, and immune activation-related genes, as well as elevated immune in ltrates. Moreover, GO and KEGG functional enrichment analyses based on the 227 hypoxia-related DEGs further con rmed that hypoxia was associated with in ammatory response, immune response, and EMT of OC. To quantify the hypoxic state of patients and facilitate clinical interpretation, a PCA algorithm was used to calculate a hypoxia score for each patient. The patients with a high score corresponded to the cluster.A group and the gene cluster.A group. A high hypoxia score was associated with activated immune in ltrates, and the upregulation of immune checkpoint-related genes, EMT-related genes, and immune activation-related genes, suggesting that the hypoxia score calculated by the PCA algorithm can well quantify the hypoxic state of OC patients. All of the above results repeatedly indicate the important role of hypoxia in the immune in ltration and malignant progression of OC. The mutual validation of the analysis results led us to believe that the hypoxia score constructed based on the hypoxia-related genes would be bene cial for the clinical selection of patients sensitive to immune checkpoint inhibitors.
Cell mutations can be divided into germline mutations and somatic mutations. The former are found in all cells in the body, due to originally occurring in the sperm or egg of the individual's parent. In contrast, somatic mutations are acquired mutations limited to a certain part of the body, which can occur due to environmental factors [28]. Non-synonymous somatic mutations of tumors may produce immunestimulating neoantigens [29,30]. Therefore, we hypothesize that hypoxia can elicit an immune response by inducing mutations in speci c genes that produce neoantigens. To explore the genes that were differentially mutated between the groups with high and low hypoxia scores, the "maftools" package in R software was used. The results revealed that the mutation frequencies of BRCA1, PRUNE2, APOB, RELN, MDN1, ZSWIM8, MGAT4C, FAT2, and VWF were higher in the group with a high hypoxia score. From a search of the literature, no relevant research on hypoxia-regulating gene mutations of OC has been reported. Our results may provide a reference for future research in this eld. Nevertheless, the potential immunogenicity of BRCA1 mutation has been demonstrated, so it may act as a biomarker for ICI therapies of patients with high grade serous ovarian cancer [31]. This provides some support for the reliability of our study.

Conclusions
In conclusion, the current study revealed that 15 hypoxia-related genes were related to the immune in ltrates and malignant progression of OC. Although the speci c regulatory mechanism is still unknown, it is clear that these 15 hypoxia-related genes can play a vital role in guiding ICI therapy of OC.

Declarations
Ethical Approval and Consent to participate Not applicable.

Consent for publication
Not applicable.

Availability of supporting data
All data generated or analyzed during this study are included in this published article and its supplementary information les.

Competing interests
The authors declare that they have no competing interests.

Funding
Not applicable.

Author Contributions
Wancheng Zhao and Lili Yin conceived and designed the study, developed the methodology, analyzed and interpreted the data, and wrote, reviewed, and/or revised the manuscript.     (H) The Kaplan-Meier curve showed that there was no signi cant difference in survival outcome between gene cluster.A and gene cluster.B (*P < 0.05, **P < 0.01, ***P < 0.001, "ns" means no signi cant difference).  -mesenchymal transition (EMT)-related genes between gene cluster.A and gene cluster.B (*P < 0.05, **P < 0.01, ***P < 0.001, "ns" means no signi cant difference).  Phenotypes associated with hypoxia score in ovarian cancer. (A) Differential expression of immune checkpoint molecules between high-and low-hypoxia-score groups. (B) Differential expression of immune activation-related genes between high-and low-hypoxia-score groups. (C) Differential expression of epithelial-mesenchymal transition (EMT)-related genes between high-and low-hypoxia-score groups (*P < 0.05, **P < 0.01, ***P < 0.001, "ns" means no signi cance difference).