Identification of Signatures of Prognosis Prediction for Melanoma Using a Hypoxia Score

Melanoma is one of the most aggressive cancers. Hypoxic microenvironment affects multiple cellular pathways and contributes to tumor progression. The purpose of the research was to investigate the association between hypoxia and melanoma, and identify the prognostic value of hypoxia-related genes. Based on the GSVA algorithm, gene expression profile collected from The Cancer Genome Atlas (TCGA) was used for calculating the hypoxia score. The Kaplan–Meier plot suggested that a high hypoxia score was correlated with the inferior survival of melanoma patients. Using differential gene expression analysis and WGCNA, a total of 337 overlapping genes associated with hypoxia were determined. Protein-protein interaction network and functional enrichment analysis were conducted, and Lasso Cox regression was performed to establish the prognostic gene signature. Lasso regression showed that seven genes displayed the best features. A novel seven-gene signature (including ABCA12, PTK6, FERMT1, GSDMC, KRT2, CSTA, and SPRR2F) was constructed for prognosis prediction. The ROC curve inferred good performance in both the TCGA cohort and validation cohorts. Therefore, our study determined the prognostic implication of the hypoxia score in melanoma and showed a novel seven-gene signature to predict prognosis, which may provide insights into the prognosis evaluation and clinical decision making.

Introduction

Melanoma is one of the highly malignant cutaneous neoplasms with a rising incidence around the world (Hallberg and Johansson, 2013; Domingues et al., 2018), characterized by its strong metastasis rate and poor prognosis (Nakamura and Fujisawa, 2018). Although surgery, chemotherapy, immunotherapy, and radiation have been performed to treat malignant melanoma, the efficacy of therapies remains limited (Domingues et al., 2018). Therefore, investigating the underlying biological mechanism and identifying new therapeutic targets are demanded.

Tumor microenvironment (TME) refers to the biological environment where tumors initiate, locate, and progress (Brandner and Haass, 2013; Roma-Rodrigues et al., 2019). The interaction between tumor and its TME influence the survival, migration, and invasion of tumor cells (Whiteside, 2008). Hypoxia is one of the essential features in the TME, which originates from the proliferation of tumor cells and increased oxygen consumption (Manoochehri Khoshinani et al., 2016). Tumor Hypoxia results in the activation of hypoxia-inducible factor (HIF), which mediates the expression of genes regulating metabolic pathways, pH regulation, DNA replication, and protein synthesis (Al Tameemi et al., 2019). Thus, tumor hypoxia contributes to heterogeneous changes, genetic instability, angiogenesis, and resistance to treatments, which has become an adverse prognostic factor of tumor assessment (Walsh et al., 2014; Jing et al., 2019). Many studies have suggested that hypoxia is related to poor prognosis in solid tumors (Winter et al., 2007; Ward et al., 2013). Likewise, hypoxia is a critical molecular program in melanoma, promoting tumor growth, invasion, treatment resistance, and relapse through the stabilization of HIF and the regulation of hypoxia-related responses (Widmer et al., 2013; Qin et al., 2016). In light of the essential role of hypoxia in melanoma, the detection and assessment of tumor hypoxia plays a critical role in clinical practice.

Assessment of the oxygen concentration, report of physiologic processes involving oxygen markers, and evaluation of endogenous molecules expression are considered as three major groups to detect tumor hypoxia status (Walsh et al., 2014). Deeply understanding the gene characteristics to estimate the degree of hypoxia would help the prognostic evaluation and treatment options. Immunohistochemistry (IHC) and plasma protein assays were developed for determining hypoxia (Russell et al., 2009; Khan et al., 2013). Recently, bioinformatics has been utilized to determine broader signatures. Based on the 26-gene hypoxia signature (Eustace et al., 2013), hypoxia status classifier was administrated in head and neck cancer (Brooks et al., 2019), and hypoxia score was implemented in lung adenocarcinoma (Liu Z. et al., 2020). Up till now, the hypoxia score in melanoma has not been investigated in detail.

Here, we calculated the hypoxia score for the analysis of gene expression profiles of melanoma which were collected from The Cancer Genome Atlas (TCGA, https://cancergenome.nih.gov). The correlation between hypoxia and prognosis was investigated, and hypoxia-associated molecules were determined. A seven-gene signature was further conducted using the profiles from TCGA and verified in the GSE54467, GSE53118, and GSE22153 dataset, providing novel insights for the assessment, treatment, and prognosis of melanoma. The workflow presenting the design of the present research was shown in Figure 1.

FIGURE 1

Figure 1. The workflow of the research.

Materials and Methods

Data Collection

The clinical information and RNA-sequencing data of skin cutaneous melanoma (SKCM) were downloaded from the TCGA database¹.

Calculation of the Hypoxia Score

Hypoxia score was calculated based on the 26-gene hypoxia signature (Eustace et al., 2013) and a gene set variation analysis (GSVA) (Eustace et al., 2013; Hänzelmann et al., 2013). GSVA is a GSE method which estimates variation of pathway activity over a sample population in an unsupervised manner (Hänzelmann et al., 2013). Hence, we used the 26-gene hypoxia signature and evaluated the GSVA score of each sample using the GSVA algorithm. The GSVA score was recognized as the hypoxia score, which represented the hypoxia status of each sample. The cut-off value was identified according to the method of best separation in R package survminer, and patients were divided into high- and low-hypoxia score groups. Such grouping aims to minimize the P value of the survival curve. Additionally, T-test was used to judge the differences of clinical indexes between groups.

Definition of Differentially Expressed Genes (DEGs)

EdgeR package was used to identify DEGs between high- and low-hypoxia score groups. The fold change (|fold change| ≥ 1.5) and adj.p < 0.05 were considered significant. Pheatmap package was used to generate the heatmap.

Identification of Hypoxia-Associated Genes by the Weighted Gene Co-expression Network Analysis (WGCNA)

The top 9829 genes, based on standard deviation, were used for further investigation. Co-expression networks were performed by using the R package WGCNA (Langfelder and Horvath, 2008). Among all the soft threshold values, we chose the β that showed the highest mean connectivity (β = 3). As the module Eigengenes (ME) was recognized to define the interpretation of gene expression profile, we associated the ME with the hypoxia feature, which showed high and low hypoxia score. Module with the highest correlation was selected, and genes of which were named hypoxia-related genes.

The Protein-Protein Interaction (PPI) Network and Functional Annotation

The overlapping genes between DEGs and hypoxia-related genes were depicted by the online Venn diagram analysis². We used the STRING (version 11.0, Search Tool for the Retrieval of Interacting Genes) and Cytoscape software (version 3.7.0) to construct the PPI network (Shannon et al., 2003; Szklarczyk et al., 2015). Molecular Complex Detection (MCODE) was utilized to determine the interaction clusters. The R package clusterprofile was used to perform functional enrichment analysis and KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway analysis (Yu et al., 2012).

Survival Analysis and Construction of the Hypoxia-Related Signature for Melanoma

For survival analysis, we utilized Kaplan–Meier survival. Survival-related genes in the multivariate Cox regression analysis were inferred using the least absolutes shrinkage and selection operator (LASSO) by the R package glmnet. Risk scores were obtained according to genes expression multiplied by a linear combination of regression coefficient acquired from the multivariate Cox regression, and patients were divided into a high-risk group and low-risk group based on the optimal cut-off point of risk score using the R package survminer. The Kaplan–Meier analysis and the receiver operating characteristic (ROC) curve were carried out using the R package ROCR.

Interaction Network Between the 7-Gene Signature and the 26-Gene Hypoxia Signature

To investigate the association between the 7-gene signature and the 26-gene list, genes from these two gene lists were input to the Gene-Cloud of Biotechnology Information (GCBI) analysis platform³ for data analysis.

External Validation of the Hypoxia-Related Signature Model

The signature model was validated using the GSE54467, GSE53118, and GSE22153 dataset derived from the Gene Expression Omnibus (GEO) database⁴. Risk scores were calculated using the same formula, and Kaplan–Meier and ROC curve analyses were implemented.

Results

Evaluation of the Degree of Hypoxia

Hypoxia scores were distributed between −0.699 to 0.659. A total of 368 patients were divided into high- and low-score groups based on the optional cut-off point of hypoxia score (0.43, Supplementary Figure S1). As shown in Figures 2A–E, no obvious differences in hypoxia scores were detected in patients with different clinical features. Additionally, mutations were common in melanoma, including BRAF (50%), NRAS (30%), MAP2K1 (6%), and KIT (2.6%). So, we also plotted the distribution of hypoxia scores to the status of driver mutations and found they were not significant (P = 0.375, P = 0.100, P = 0.765, P = 0.145, Figures 2F–I).

FIGURE 2

Figure 2. Distribution of hypoxia score in melanoma. (A) Distribution of hypoxia score of patients with different TNM staging. (B) Distribution of hypoxia score of patients with different T stage. (C) Distribution of hypoxia score of patients with or without lymph node metastasis. (D) Distribution of hypoxia score of patients with or without distant metastasis. (E) Distribution of hypoxia score of patients younger than 65 and those older than 65 years of age. (F–I) Distribution of hypoxia score of patients with BRAF mutant and BRAF wildtype, patients with NRAS mutant and NRAS wildtype, patients with MAP2K1 mutant and MAP2K1 wildtype, and patients with KIT mutant and KIT wildtype, respectively. (J) Patients were divided into high- and low-hypoxia score groups based on the cut-off value. Patients with a high hypoxia score showed a better prognosis compared to patients with a low score ((P) = 0.007). (K) Heatmap of the DEGs of high-hypoxia score group vs. low-hypoxia score group. p < 0.05, |fold change| ≥ 1.5. DEGs, differentially expressed genes.

The effects of hypoxia on prognosis were analyzed. The Kaplan–Meier plot suggested that patients with high hypoxia scores had a poor prognosis (P = 0.007, Figure 2J). To further determine the correlation of gene expression with hypoxia scores, we did differential gene expression analysis between high and low hypoxia scores. Of the 415 differential expression genes (DEGs), 365 genes were upregulated, while 50 genes were downregulated. Heatmaps in Figure 2K inferred distinct gene expression profiles of cases belong to high- vs. low-hypoxia scores.

Identification of the Most Relevant Module Genes for Hypoxia in Melanoma

We selected the top 9829 (of 19658) after sorting by the standard deviation (Figures 3A–C). The co-expression network was constructed, and 13 modules were determined. Correlation analysis between the module eigengenes and hypoxia scores showed that the yellow module (Figure 3D, Module–trait relationships = 0.43, P = 0.000) had the highest association with the degree of hypoxia. Then, 802 genes in the module were considered to be hub hypoxia-related genes for further investigation.

FIGURE 3

Figure 3. Determination of modules correlated with the hypoxia of melanoma in the WGCNA. (A) Analysis of the scale-free fit index and the mean connectivity for various soft-thresholding powers. (B) Checking the scale free topology when β = 3. Correlation coefficient = 0.9, which showed scale-free topology. (C) Dendrogram of genes clustered according to a dissimilarity measure (1-TOM). (D) Heatmap of the correlation between module Eigengenes and hypoxia. WGCNA, the weighted gene co-expression network analysis.

Protein-Protein Interactions and Functional Enrichment Analysis

A total of 337 genes were overlapped between DEGs and hypoxia-related genes (Figure 4A). To explore the interplay among 337 overlapping genes, the STRING tool with confidence > 0.7 was used to construct a PPI network. There were 10 modules in the network, including 195 nodes and 1173 edges. Modules with 10 or more nodes were selected for further analysis (Figure 4B). Based on the connection degree, we named these modules IVL, and FLG modules, respectively. In the IVL module, 528 edges involving 33 nodes were formed in the network. IVL, TGM1, LOR, SPRR1B, and PPL were the remarkable nodes, as they had the most connections with others. In the FLG module, FLG, DSG1, DSG3, PKP3, PKP1, KRT14, and DSC1 occupied the center of the module.

FIGURE 4

Figure 4. Analysis of DEGs. (A) Venn diagrams showing the number of commonly genes in DEGs and yellow module. (B) PPI networks of overlapping genes. A large node represented a higher degree. (C–E) Go enrichment analysis of biological process (BP), molecular function (MF), and cellular component (CC). (F) KEGG pathway enrichment analysis of the overlapping genes. DEGs, differentially expressed genes; PPI, the protein-protein interaction; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes.

To better understand the biological significance, we conducted enrichment analysis of the 337 overlapping genes. As shown in Figure 4, a total of 27 terms of biological process (BP), 8 terms of cellular component (CC), and 14 terms of molecular function (MF) were enriched (P < 0.05). Top GO terms comprised epidermis development, skin development and epidermal cell differentiation (Figure 4C), serine type endopeptidase activity, serine type peptidase activity, and serine hydrolase activity (Figure 4D), and cornified envelope and cell-cell junction (Figure 4E). Besides, KEGG analysis suggested overlapping genes were enriched in Staphylococcus aureus infection, estrogen signaling pathway, and IL-17 signaling pathway (Figure 4F).

Determination of Prognostic Molecules and a Prognostic Risk Model

We generated Kaplan–Meier survival curves to explore the independent prognostic impact of 337 overlapping genes and found that 29 genes were associated with prognosis in the log-rank test (P < 0.05). A total of 7 genes identified with the LASSO algorithms included ATP-binding cassette sub-family a member 12 (ABCA12), protein tyrosine kinase 6 (PTK6), fermitin family member 1 (FERMT1), gasdermin C (GSDMC), keratin 2 (KRT2), cystatin A (CSTA), and small proline rich protein 2F (SPRR2F), and constructed as a seven-gene signature model (Table 1). The risk score = 0.26084 * Expression (ABCA12) + 0.05797 * Expression (PTK6) + 0.14404 * Expression (FERMT1) + (−0.44473) * Expression (GSDMC) + (−0.09102) * Expression (KRT2) + (−0.02677) * Expression (CSTA) + 0.11245 * Expression (SPRR2F). The roles of these 7 genes in melanoma and hypoxia responses were described in Table 2. Also, we explored the relationships among genes from the 7-gene signature and 26-gene list. Although genes from the 7-gene signature were different from those of the 26-gene one, there were common regulators associated with hypoxic responses, including EGFR, ERBB2, and miR-125a (Figures 5A,B, Table 3).

TABLE 1

Table 1. The results of Univariate Cox regression analysis.

TABLE 2

Table 2. The roles of 7 genes in melanoma and hypoxia response.

FIGURE 5

Figure 5. Interaction network of molecules associated with the genes from the 7-gene signature and the 26-gene list. (A) Interaction network of genes associated with the 7-gene signature and the 26-gene list. Colors indicated types of genes: light blue, input genes; orange, activated genes; red, expressed genes; green, associated genes; dark blue, inhibited genes; yellow, the largest connection counts. Node size was adjusted according to the number of associated genes. (B) Interaction network of miRNAs correlated with the 7-gene signature and the 26-gene list. Colors indicated types of molecules: light blue, input genes; purple, targeted miRNAs; yellow, the largest connection counts. Node size was adjusted according to the number of associated miRNAs.

TABLE 3

Table 3. Common regulators and downstream effectors.

Kaplan–Meier curve and ROC were utilized to assess the prognostic capacity of the seven-gene signature model, and similar procedures were performed in the external data. The results showed that genes in the signature model performed well-predicting prognosis within the TCGA cohort (Figures 6A–G). Figures 6H–K suggested that patients with low-risk scores had significantly longer overall survival than those with high-risk scores in TCGA, GSE54467, GSE53118, and GSE22153 dataset (P < 0.001, P = 0.004, P = 0.017, P = 0.048). The AUCs were 0.716 (95% CI: 0.661–0.771), 0.667 (95% CI: 0.541–0.792), 0.648 (95% CI: 0.419–0.878), and 0.628 (95% CI: 0.406–0.849), respectively (Figures 6H–K).

FIGURE 6

Figure 6. Kaplan–Meier analysis, risk score analysis, and ROC analysis for the seven-gene signature. (A–G) Kaplan–Meier curves for overall survival of ABCA12, CSTA, FERMT1, GSDMC, KRT2, PTK6, and SPRR2F. (H–K) Kaplan–Meier curves for overall survival of risk score and ROC analysis for the seven-gene signature in TCGA cohort, GSE54467, GSE53118, and GSE22153. ROC, receiver operating characteristic; TCGA, The Cancer Genome Atlas.

Discussion

Hypoxia, one of the hallmarks of TME, is a biological condition present in most tumors (Jing et al., 2019). Tumor hypoxia exacerbates progression and metastasis through both physiological and genomic mechanisms (Akanji et al., 2019). Investigating crucial features of tumor hypoxia environment may facilitate clinical decision-making.

Previous studies identified several genes, long non-coding RNAs and miRNAs as promising therapeutic biomarkers in melanoma (Zhang et al., 2017; Wei et al., 2019; Xu et al., 2019). However, the differentially expressed signatures were explored between the normal and tumor samples, or between the primary and metastatic tissues, and molecules associated with the progression of cancer were not taken into consideration. Notably, the focus of our study was to estimate the degree of hypoxia according to the evidential basis for 26-gene hypoxia signature (Eustace et al., 2013), and high hypoxia score was demonstrated as a strong predictor of poor clinical outcome. Subsequently, we identified the promising hypoxia-related genes associated with prognosis.

Based on bioinformatics methods and databases, hypoxia score was calculated, and patients were divided into high- and low-score groups. DEGs were collected using differential gene expression analysis. At the same time, WGCNA analysis was performed to select the modules with the strongest relationship between genes in the modules and the module traits. The overlapping 337 genes of the above two clusters were determined as the hypoxia strongly associated genes related to melanoma. Functional analysis showed these 337 genes to be closely related to the development of melanoma, like via cell-cell junction. Cell junction was reported to be relevant for the metastatic process (Knights et al., 2012). Also, the process of epidermis development and epidermal cell differentiation were enriched. Previous studies showed the hyperplastic epidermal region was accompanied by aberrant expression of keratin 14, and melanoma cells were able to increase expression of keratins 8, 19 (Kodet et al., 2015). keratin 8, 14,19 were also observed in the FLG module from PPI network. These results showed that epidermis surrounding melanoma performed hyperplastic features, and indicated the possible interaction between melanoma cells and keratinocytes. KEGG analysis highlighted the estrogen signaling pathway and the IL-17 signaling pathway. Several studies pointed out that the estrogen signaling pathway relied on the balance between estrogen receptor (ER) α and ERβ expression, and the levels of ERβ regulated the capacity of melanoma invasion (Marzagalli et al., 2016; Rajabi et al., 2017). Additionally, IL-17/IL-17RA pathway stimulated cell proliferation of mouse B16F10 and human A375 and A2058 cell lines (Chen et al., 2019). IL-17 and IL-23 immunohistochemistry expression were increased in the melanoma tissues, possibly enhancing VEGF expression and angiogenesis (Ganzetti et al., 2015). Therefore, these analyses supported the hypothesis of the importance of hypoxia microenvironment in the regulation of the biological behavior of tumor cells and surrounding non-tumor cells.

Based on the log-rank test identifying the genes associated with prognosis, LASSO was performed, and seven characteristic variables were extracted. ABCA12 was upregulated in ovarian carcinoma and colorectal cancer, which was recognized as a promising candidate marker (Hlavata et al., 2012; Elsnerova et al., 2016). Mutations in ABCA12 were related to malignant melanoma (Natsuga et al., 2007). PTK6, a non-receptor type tyrosine kinase, was involved in breast, pancreatic cancer and metastatic skin cancer. It was recognized that PTK6 regulated proliferation and migration (Gotoh et al., 2014; Ito et al., 2017; Liu G. et al., 2020). FERMT1, encoding Kindlin-1, was correlated with metastasis and poor prognosis in several solid tumors (Liu et al., 2016; Sarvi et al., 2018). GSDMC functioned as an oncogene, enhancing cell proliferation and tumorigenesis in lung adenocarcinoma and colorectal carcinogenesis (Miguchi et al., 2016; Wei et al., 2020). It presented high in malignant melanoma but undetectable in normal epithelial cells, which might be associated with the metastasis of cells (Xia et al., 2019). CSTA, one of the tumor suppressors, had the anti-apoptotic effect and maintaining cell-cell adhesion. It was upregulated in several epithelial-derived malignancies, including squamous cell carcinoma (Gupta et al., 2015; Ma et al., 2018). KRT2 was found to form a mechanically resilient cytoskeleton and contribute to the skin homeostasis (Fischer et al., 2016). SPRR2F, a cross-linked envelope protein of keratinocytes, providing the protective barrier function (Cabral et al., 2001). Although there was no report of KRT2 and SPRR2F as a prognostic molecule of tumors, KRT2 and SPRR2F might function as promising biomarkers in melanoma. The consistency of our findings regarding ABCA12, PTK6, FERMT1, GSDMC and CSTA with previous studies suggested our method to be reliable, and thus supported the reliability of these potential prognostic and therapeutic targets to a certain extent.

Previous studies inferred that the expression of PTK6 were up-regulated, and FERMT1 were down-regulated in response to the hypoxia condition (Hlavata et al., 2012; Regan Anderson et al., 2013; Lin and Liu, 2019). PTK6 expression depended on both HIF-1α and HIF-2α, which were reported to have a direct regulation of PTK6 transcription. In the analytic process of investigating the effect of hypoxia on the vhl-deficient cells, HIF-regulated genes were obtained. FERMT1 was one of the 214 downregulated DEGs. Additionally, the increased expression of CSTA was detected in hypoxic A431 cells (Park et al., 2010). Although there was no common gene between the 7- and 26-gene signatures, a total of 3 genes, including epidermal growth factor receptor (EGFR), erb-b2 receptor tyrosine kinase 2 (ERBB2), and miR-125a, were identified as common regulators and effectors in these two gene lists in a context-dependent manner. PTK6 was reported to enhance EGFR signaling by direct phosphorylation of EGFR and inhibition of its degradation (Li et al., 2012), and EGFR might promote the cellular response to hypoxia by increasing HIF-1α expression (Swinson and O’Byrne, 2006). Through the split ubiquitin (Ub)-based membrane yeast two-hybrid assay, EGFR was reported to be physically associated with aldolase (ALDOA) and triosephosphate isomerase 1 (TPI1), respectively (Deribe et al., 2009). However, the potential functions of ALDOA and TPI1 need to be further explored. Furthermore, ERBB2, also known as HER2, was recognized as a regulator of HIF-2α and a driver of hypoxic responses (Jarman et al., 2019). PTK6 was coamplified with ERBB2 to promote cell proliferation (Xiang et al., 2008). Additionally, ERBB2 and keratin 17 (KRT17) were found to locate in the same chromosome region, which might have the following tumor associations (Zhang et al., 2013). Apart from regulating the expression of genes, hypoxia-regulated microRNAs (miRNAs) were identified. MiR-125a was a direct target of HIF-1α and drove the reduction of vascular endothelial growth factor A (VEGFA) (Dai et al., 2015; Pan et al., 2018). Based on the map of human miRNA interactome, enolase 1 (ENO1), FERMT1, and TPI1 were observed in the interaction sites of miR-125a and further examinations were demanded (Helwak et al., 2013).

Saxena and Jolly summarized different extents of hypoxia (Saxena and Jolly, 2019). Under acute hypoxia, HIF-1α levels stayed high to regulate acute response, while HIF-2α levels were stabilized later and played a crucial role during chronic hypoxia. Besides, cyclic hypoxia enhanced the expression of HIF-1α instead of HIF-2α. Several factors implicated in these hypoxia conditions were determined, including HSP-70, HAF, H3, H4, REST, and miR-429. Although genes identified in our study have been reported to function in hypoxic responses, there was no report of them to make a distinction of conditions of hypoxia, and further experimental verification is required.

Considering the accuracy of these prognostic genes, a seven-signature model was established based on the combination of genes. Cases in the low-risk group inferred obviously better survival than patients in the high-risk group. The prognosis predictive performance of the model was relatively good not only in the TCGA melanoma cohort but also in the GSE54467, GSE53118, and GSE22153 cohort. Additionally, we investigated whether the clinical features were correlated with the degree of hypoxia, and the results showed that no apparent differences in hypoxia score were observed. BRAF mutation was found to increase HIF-1α expression and influenced survival in previous studies (Kumar et al., 2007; Zerilli et al., 2010). KIT mutant was reported to require HIF-1α to transform melanocytes into melanoma cells (Monsel et al., 2010). In our cohort, the hypoxia score in the BRAF-mutant or KIT-mutant group was slightly higher than that of the wildtype group, but it was not statistically significant. It could be because of an inevitable limitation, the sample size. There were two other limitations to our study. Firstly, data were collected from TCGA, where the potential for selection bias could not be excluded, but we validated the results in the GEO database and demonstrated the reliability to some extent. Secondly, analysis in our study was descriptive, further research in vitro and in vivo could enhance our understanding of the critical genes.

In conclusion, we applied the hypoxia score to determine the degree of hypoxia in TME and identified the prognostic role of hypoxia score. Furthermore, using bioinformatics and machine learning methods, we determined the seven-gene prognostic signature as a potential prognostic predictor and therapeutic targets for melanoma.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi.nlm.nih.gov/geo/, GSE53118; https://www.ncbi.nlm.nih.gov/geo/, GSE54467; https://www.ncbi.nlm.nih.gov/geo/, GSE22153.

Author Contributions

XZ, FL, YY, and JX contributed to the design of this study. YS, LY, and YY contributed to the analysis of this study. YS contributed to drafting the text and preparing the tables and figures. All authors participated in the data collection, critical review, revision of this manuscript, contributed to the article, and approved the submitted version.

Funding

This research is supported by a grant from The National Natural Science Foundation of China (81673060).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2020.570530/full#supplementary-material

Supplementary Figure 1 | The process of finding optimal cut-off value to divide the patients into high- and low-hypoxia score groups.

Footnotes

1. ^ https://tcga-data.nci.nih.gov/
2. ^ http://bioinformatics.psb.ugent.be/webtools/Venn/
3. ^ https://www.gcbi.com.cn
4. ^ https://www.ncbi.nlm.nih.gov/geo

References

Akanji, M. A., Rotimi, D., and Adeyemi, O. S. (2019). Hypoxia-inducible factors as an alternative source of treatment strategy for cancer. Oxid. Med. Cell Longev. 2019:8547846. doi: 10.1155/2019/8547846

PubMed Abstract | CrossRef Full Text | Google Scholar

Al Tameemi, W., Dale, T. P., Al-Jumaily, R. M. K., and Forsyth, N. R. (2019). Hypoxia-modified cancer cell metabolism. Front. Cell Dev. Biol. 7:4. doi: 10.3389/fcell.2019.00004

PubMed Abstract | CrossRef Full Text | Google Scholar

Brandner, J. M., and Haass, N. K. (2013). Melanoma’s connections to the tumour microenvironment. Pathology 45, 443–452. doi: 10.1097/PAT.0b013e328363b3bd

PubMed Abstract | CrossRef Full Text | Google Scholar

Brooks, J. M., Menezes, A. N., Ibrahim, M., Archer, L., Lal, N., Bagnall, C. J., et al. (2019). Development and Validation of a Combined Hypoxia and Immune Prognostic Classifier for Head and Neck Cancer. Clin. Cancer Res. 25, 5315–5328. doi: 10.1158/1078-0432.Ccr-18-3314

PubMed Abstract | CrossRef Full Text | Google Scholar

Cabral, A., Voskamp, P., Cleton-Jansen, A.-M., South, A., Nizetic, D., and Backendorf, C. (2001). Structural Organization and Regulation of the Small Proline-rich Family of Cornified Envelope Precursors Suggest a Role in Adaptive Barrier Function. J. Biol. Chem. 276, 19231–19237. doi: 10.1074/jbc.M100336200

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, Y.-S., Huang, T.-H., Liu, C.-L., Chen, H.-S., Lee, M.-H., Chen, H.-W., et al. (2019). Locally Targeting the IL-17/IL-17RA Axis Reduced Tumor Growth in a Murine B16F10 Melanoma Model. Hum. Gene Ther. 30, 273–285. doi: 10.1089/hum.2018.104

PubMed Abstract | CrossRef Full Text | Google Scholar

Dai, J., Wang, J., Yang, L., Xiao, Y., and Ruan, Q. (2015). miR-125a regulates angiogenesis of gastric cancer by targeting vascular endothelial growth factor A. Int. J. Oncol. 47, 1801–1810. doi: 10.3892/ijo.2015.3171

PubMed Abstract | CrossRef Full Text | Google Scholar

Deribe, Y. L., Wild, P., Chandrashaker, A., Curak, J., Schmidt, M. H. H., Kalaidzidis, Y., et al. (2009). Regulation of epidermal growth factor receptor trafficking by lysine deacetylase HDAC6. Sci. Signal 2:ra84. doi: 10.1126/scisignal.2000576

PubMed Abstract | CrossRef Full Text | Google Scholar

Domingues, B., Lopes, J. M., Soares, P., and Populo, H. (2018). Melanoma treatment in review. Immunotargets Ther. 7, 35–49. doi: 10.2147/ITT.S134842

PubMed Abstract | CrossRef Full Text | Google Scholar

Elsnerova, K., Mohelnikova-Duchonova, B., Cerovska, E., Ehrlichova, M., Gut, I., Rob, L., et al. (2016). Gene expression of membrane transporters: importance for prognosis and progression of ovarian carcinoma. Oncol Rep. 35, 2159–2170. doi: 10.3892/or.2016.4599

PubMed Abstract | CrossRef Full Text | Google Scholar

Eustace, A., Mani, N., Span, P. N., Irlam, J. J., Taylor, J., Betts, G. N. J., et al. (2013). A 26-Gene Hypoxia Signature Predicts Benefit from Hypoxia-Modifying Therapy in Laryngeal Cancer but Not Bladder Cancer. Clin. Cancer Res. 19, 4879–4888. doi: 10.1158/1078-0432.Ccr-13-0542

PubMed Abstract | CrossRef Full Text | Google Scholar

Fischer, H., Langbein, L., Reichelt, J., Buchberger, M., Tschachler, E., and Eckhart, L. (2016). Keratins K2 and K10 are essential for the epidermal integrity of plantar skin. J. Dermatol. Sci. 81, 10–16. doi: 10.1016/j.jdermsci.2015.10.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Ganzetti, G., Rubini, C., Campanati, A., Zizzi, A., Molinelli, E., Rosa, L., et al. (2015). IL-17, IL-23, and p73 expression in cutaneous melanoma: a pilot study. Melanoma Res. 25, 232–238. doi: 10.1097/CMR.0000000000000151

PubMed Abstract | CrossRef Full Text | Google Scholar

Gotoh, N., Ono, H., Basson, M. D., and Ito, H. (2014). PTK6 promotes cancer migration and invasion in pancreatic cancer cells dependent on ERK signaling. PLoS One 9:e96060. doi: 10.1371/journal.pone.0096060

PubMed Abstract | CrossRef Full Text | Google Scholar

Gupta, A., Nitoiu, D., Brennan-Crispi, D., Addya, S., Riobo, N. A., Kelsell, D. P., et al. (2015). Cell cycle- and cancer-associated gene networks activated by Dsg2: evidence of cystatin A deregulation and a potential role in cell-cell adhesion. PLoS One 10:e0120091. doi: 10.1371/journal.pone.0120091

PubMed Abstract | CrossRef Full Text | Google Scholar

Hallberg, Ö, and Johansson, O. (2013). Increasing Melanoma—Too Many Skin Cell Damages or Too Few Repairs? Cancers 5, 184–204. doi: 10.3390/cancers5010184

PubMed Abstract | CrossRef Full Text | Google Scholar

Hänzelmann, S., Castelo, R., and Guinney, J. (2013). GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics 14:7. doi: 10.1186/1471-2105-14-17

PubMed Abstract | CrossRef Full Text | Google Scholar

Helwak, A., Kudla, G., Dudnakova, T., and Tollervey, D. (2013). Mapping the Human miRNA Interactome by CLASH reveals frequent noncanonical binding. Cell 153, 654–665. doi: 10.1016/j.cell.2013.03.043

PubMed Abstract | CrossRef Full Text | Google Scholar

Hlavata, I., Mohelnikova-Duchonova, B., Vaclavikova, R., Liska, V., Pitule, P., Novak, P., et al. (2012). The role of ABC transporters in progression and clinical outcome of colorectal cancer. Mutagenesis 27, 187–196. doi: 10.1093/mutage/ger075

PubMed Abstract | CrossRef Full Text | Google Scholar

Ito, K., Park, S. H., Katsyv, I., Zhang, W., De Angelis, C., Schiff, R., et al. (2017). PTK6 regulates growth and survival of endocrine therapy-resistant ER+ breast cancer cells. npj Breast Cancer 3, 1–7. doi: 10.1038/s41523-017-0047-41

CrossRef Full Text | Google Scholar

Jarman, E. J., Ward, C., Turnbull, A. K., Martinez-Perez, C., Meehan, J., Xintaropoulou, C., et al. (2019). HER2 regulates HIF-2α and drives an increased hypoxic response in breast cancer. Breast Cancer Res. 21:10. doi: 10.1186/s13058-019-1097-1090

CrossRef Full Text | Google Scholar

Jing, X., Yang, F., Shao, C., Wei, K., Xie, M., Shen, H., et al. (2019). Role of hypoxia in cancer therapy by regulating the tumor microenvironment. Mol. Cancer 18:157. doi: 10.1186/s12943-019-1089-1089

CrossRef Full Text | Google Scholar

Khan, R. H., Ahmad, Y., Sharma, N. K., Garg, I., Ahmad, M. F., Sharma, M., et al. (2013). An Insight into the Changes in Human Plasma Proteome on Adaptation to Hypobaric Hypoxia. PLoS One 8:e67548. doi: 10.1371/journal.pone.0067548

PubMed Abstract | CrossRef Full Text | Google Scholar

Knights, A. J., Funnell, A. P., Crossley, M., and Pearson, R. C. (2012). Holding tight: cell junctions and cancer spread. Trends Cancer Res 8, 61–69.

Google Scholar

Kodet, O., Lacina, L., Krejčí, E., Dvořánková, B., Grim, M., Štork, J., et al. (2015). Melanoma cells influence the differentiation pattern of human epidermal keratinocytes. Mol. Cancer 14:1. doi: 10.1186/1476-4598-14-11

CrossRef Full Text | Google Scholar

Kumar, S. M., Yu, H., Edwards, R., Chen, L., Kazianis, S., Brafford, P., et al. (2007). Mutant V600E BRAF Increases Hypoxia Inducible Factor-1α Expression in Melanoma. Cancer Res. 67, 3177–3184. doi: 10.1158/0008-5472.Can-06-3312

PubMed Abstract | CrossRef Full Text | Google Scholar

Langfelder, P., and Horvath, S. (2008). WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics 9:559. doi: 10.1186/1471-2105-9-559

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, X., Lu, Y., Liang, K., Hsu, J. M., Albarracin, C., Mills, G. B., et al. (2012). Brk/PTK6 sustains activated EGFR signaling through inhibiting EGFR degradation and transactivating EGFR. Oncogene 31, 4372–4383. doi: 10.1038/onc.2011.608

PubMed Abstract | CrossRef Full Text | Google Scholar

Lin, H., and Liu, H. (2019). Hypoxia on vhl-deficient cells to obtain hif related genes through bioinformatics analysis. bioRxiv [Preprint] doi: 10.1101/863662

CrossRef Full Text | Google Scholar

Liu, C. C., Cai, D. L., Sun, F., Wu, Z. H., Yue, B., Zhao, S. L., et al. (2016). FERMT1 mediates epithelial–mesenchymal transition to promote colon cancer metastasis via modulation of β-catenin transcriptional activity. Oncogene 36, 1779–1792. doi: 10.1038/onc.2016.339

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, G., Li, C., Zhen, H., Zhang, Z., and Sha, Y. (2020). Identification of prognostic gene biomarkers for metastatic skin cancer using data mining. Biomed. Rep. 13, 22–30. doi: 10.3892/br.2020.1307

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, Z., Tao, H., Tong, L., and Li, H. (2020). Genome-wide analysis of the hypoxia-related DNA methylation-driven genes in lung adenocarcinoma progression. Biosci. Rep. 40:BSR20194200. doi: 10.1042/bsr20194200

PubMed Abstract | CrossRef Full Text | Google Scholar

Ma, Y., Chen, Y., Li, Y., Grün, K., Berndt, A., Zhou, Z., et al. (2018). Cystatin A suppresses tumor cell growth through inhibiting epithelial to mesenchymal transition in human lung cancer. Oncotarget 9, 14084–14098. doi: 10.18632/oncotarget.23505

PubMed Abstract | CrossRef Full Text | Google Scholar

Manoochehri Khoshinani, H., Afshar, S., and Najafi, R. (2016). Hypoxia: a double-edged sword in cancer therapy. Cancer Invest. 34, 536–545. doi: 10.1080/07357907.2016.1245317

PubMed Abstract | CrossRef Full Text | Google Scholar

Marzagalli, M., Montagnani Marelli, M., Casati, L., Fontana, F., Moretti, R. M., and Limonta, P. (2016). Estrogen receptor beta in melanoma: from molecular insights to potential clinical utility. Front. Endocrinol. 7:140. doi: 10.3389/fendo.2016.00140

PubMed Abstract | CrossRef Full Text | Google Scholar

Miguchi, M., Hinoi, T., Shimomura, M., Adachi, T., Saito, Y., Niitsu, H., et al. (2016). Gasdermin C Is Upregulated by Inactivation of Transforming Growth Factor beta Receptor Type II in the Presence of Mutated Apc, Promoting Colorectal Cancer Proliferation. PLoS One 11:e0166422. doi: 10.1371/journal.pone.0166422

PubMed Abstract | CrossRef Full Text | Google Scholar

Monsel, G., Ortonne, N., Bagot, M., Bensussan, A., and Dumaz, N. (2010). c-Kit mutants require hypoxia-inducible factor 1alpha to transform melanocytes. Oncogene 29, 227–236. doi: 10.1038/onc.2009.320

PubMed Abstract | CrossRef Full Text | Google Scholar

Nakamura, Y., and Fujisawa, Y. (2018). Diagnosis and management of acral lentiginous melanoma. Curr. Treatment Options Oncol. 19:42. doi: 10.1007/s11864-018-0560-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Natsuga, K., Akiyama, M., Kato, N., Sakai, K., Sugiyama-Nakagiri, Y., Nishimura, M., et al. (2007). Novel ABCA12 Mutations Identified in Two Cases of Non-Bullous Congenital Ichthyosiform Erythroderma Associated with Multiple Skin Malignant Neoplasia. J. Invest. Dermatol. 127, 2669–2673. doi: 10.1038/sj.jid.5700885

PubMed Abstract | CrossRef Full Text | Google Scholar

Pan, L., Zhou, L., Yin, W., Bai, J., and Liu, R. (2018). miR-125a induces apoptosis, metabolism disorder and migrationimpairment in pancreatic cancer cells by targeting Mfn2-related mitochondrial fission. Int. J. Oncol. 53, 124–136. doi: 10.3892/ijo.2018.4380

PubMed Abstract | CrossRef Full Text | Google Scholar

Park, J. E., Tan, H. S., Datta, A., Lai, R. C., Zhang, H., Meng, W., et al. (2010). Hypoxic Tumor Cell Modulates Its Microenvironment to Enhance Angiogenic and Metastatic Potential by Secretion of Proteins and Exosomes. Mol. Cell. Proteom. 9, 1085–1099. doi: 10.1074/mcp.M900381-MCP200

PubMed Abstract | CrossRef Full Text | Google Scholar

Qin, Y., Roszik, J., Chattopadhyay, C., Hashimoto, Y., Liu, C., Cooper, Z. A., et al. (2016). Hypoxia-Driven Mechanism of Vemurafenib Resistance in Melanoma. Mol. Cancer Ther. 15, 2442–2454. doi: 10.1158/1535-7163.Mct-15-0963

PubMed Abstract | CrossRef Full Text | Google Scholar

Rajabi, P., Bagheri, M., and Hani, M. (2017). Expression of Estrogen Receptor Alpha in Malignant Melanoma. Adv. Biomed. Res. 6:14. doi: 10.4103/2277-9175.200789

PubMed Abstract | CrossRef Full Text | Google Scholar

Regan Anderson, T. M., Peacock, D. L., Daniel, A. R., Hubbard, G. K., Lofgren, K. A., Girard, B. J., et al. (2013). Breast tumor kinase (Brk/PTK6) is a mediator of hypoxia-associated breast cancer progression. Cancer Res. 73, 5810–5820. doi: 10.1158/0008-5472.CAN-13-0523

PubMed Abstract | CrossRef Full Text | Google Scholar

Roma-Rodrigues, C., Mendes, R., Baptista, P., and Fernandes, A. (2019). Targeting Tumor Microenvironment for Cancer Therapy. Int. J. Mol. Sci. 20:840. doi: 10.3390/ijms20040840

PubMed Abstract | CrossRef Full Text | Google Scholar

Russell, J., Carlin, S., Burke, S. A., Wen, B., Yang, K. M., and Ling, C. C. (2009). Immunohistochemical detection of changes in tumor hypoxia. Int. J. Radiat. Oncol. Biol. Phys. 73, 1177–1186. doi: 10.1016/j.ijrobp.2008.12.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Sarvi, S., Patel, H., Li, J., Dodd, G. L., Creedon, H., Muir, M., et al. (2018). Kindlin-1 Promotes Pulmonary Breast Cancer Metastasis. Cancer Res. 78, 1484–1496. doi: 10.1158/0008-5472.Can-17-1518

PubMed Abstract | CrossRef Full Text | Google Scholar

Saxena, K., and Jolly, M. K. (2019). Acute vs. chronic vs. cyclic hypoxia: their differential dynamics, molecular mechanisms, and effects on tumor progression. Biomolecules 9:339. doi: 10.3390/biom9080339

PubMed Abstract | CrossRef Full Text | Google Scholar

Shannon, P., Markiel, A., Ozier, O., Baliga, N. S., Wang, J. T., Ramage, D., et al. (2003). Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 13, 2498–2504. doi: 10.1101/gr.1239303

PubMed Abstract | CrossRef Full Text | Google Scholar

Swinson, D. E., and O’Byrne, K. J. (2006). Interactions between hypoxia and epidermal growth factor receptor in non-small-cell lung cancer. Clin. Lung Cancer 7, 250–256. doi: 10.3816/CLC.2006.n.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Szklarczyk, D., Franceschini, A., Wyder, S., Forslund, K., Heller, D., Huerta-Cepas, J., et al. (2015). STRING v10: protein–protein interaction networks, integrated over the tree of life. Nucleic Acids Res. 43, D447–D452. doi: 10.1093/nar/gku1003

PubMed Abstract | CrossRef Full Text | Google Scholar

Walsh, J. C., Lebedev, A., Aten, E., Madsen, K., Marciano, L., and Kolb, H. C. (2014). The clinical importance of assessing tumor hypoxia: relationship of tumor hypoxia to prognosis and therapeutic opportunities. Antioxid. Redox Signal. 21, 1516–1554. doi: 10.1089/ars.2013.5378

PubMed Abstract | CrossRef Full Text | Google Scholar

Ward, C., Langdon, S. P., Mullen, P., Harris, A. L., Harrison, D. J., Supuran, C. T., et al. (2013). New strategies for targeting the hypoxic tumour microenvironment in breast cancer. Cancer Treatment Rev. 39, 171–179. doi: 10.1016/j.ctrv.2012.08.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Wei, C. Y., Zhu, M. X., Lu, N. H., Peng, R., Yang, X., Zhang, P. F., et al. (2019). Bioinformatics-based analysis reveals elevated MFSD12 as a key promoter of cell proliferation and a potential therapeutic target in melanoma. Oncogene 38, 1876–1891. doi: 10.1038/s41388-018-0531-536

CrossRef Full Text | Google Scholar

Wei, J., Xu, Z., Chen, X., Wang, X., Zeng, S., Qian, L., et al. (2020). Overexpression of GSDMC is a prognostic factor for predicting a poor outcome in lung adenocarcinoma. Mol. Med. Rep. 21, 360–370. doi: 10.3892/mmr.2019.10837

PubMed Abstract | CrossRef Full Text | Google Scholar

Whiteside, T. L. (2008). The tumor microenvironment and its role in promoting tumor growth. Oncogene 27, 5904–5912. doi: 10.1038/onc.2008.271

PubMed Abstract | CrossRef Full Text | Google Scholar

Widmer, D. S., Hoek, K. S., Cheng, P. F., Eichhoff, O. M., Biedermann, T., Raaijmakers, M. I. G., et al. (2013). Hypoxia contributes to melanoma heterogeneity by triggering HIF1alpha-dependent phenotype switching. J. Invest. Dermatol. 133, 2436–2443. doi: 10.1038/jid.2013.115

PubMed Abstract | CrossRef Full Text | Google Scholar

Winter, S. C., Buffa, F. M., Silva, P., Miller, C., Valentine, H. R., Turley, H., et al. (2007). Relation of a Hypoxia Metagene Derived from Head and Neck Cancer to Prognosis of Multiple Cancers. Cancer Res. 67, 3441–3449. doi: 10.1158/0008-5472.Can-06-3322

PubMed Abstract | CrossRef Full Text | Google Scholar

Xia, X., Wang, X., Cheng, Z., Qin, W., Lei, L., Jiang, J., et al. (2019). The role of pyroptosis in cancer: pro-cancer or pro-“host”? Cell Death Dis. 10:650. doi: 10.1038/s41419-019-1883-1888

CrossRef Full Text | Google Scholar

Xiang, B., Chatti, K., Qiu, H., Lakshmi, B., Krasnitz, A., Hicks, J., et al. (2008). Brk is coamplified with ErbB2 to promote proliferation in breast cancer. Proc. Natl. Acad. Sci. U.S.A. 105, 12463–12468. doi: 10.1073/pnas.0805009105

PubMed Abstract | CrossRef Full Text | Google Scholar

Xu, Y., Han, W., Xu, W. H., Wang, Y., Yang, X. L., Nie, H. L., et al. (2019). Identification of differentially expressed genes and functional annotations associated with metastases of the uveal melanoma. J. Cell Biochem. 120, 19202–19214. doi: 10.1002/jcb.29250

PubMed Abstract | CrossRef Full Text | Google Scholar

Yu, G., Wang, L. G., Han, Y., and He, Q. Y. (2012). clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS 16, 284–287. doi: 10.1089/omi.2011.0118

PubMed Abstract | CrossRef Full Text | Google Scholar

Zerilli, M., Zito, G., Martorana, A., Pitrone, M., Cabibi, D., Cappello, F., et al. (2010). BRAF(V600E) mutation influences hypoxia-inducible factor-1alpha expression levels in papillary thyroid cancer. Mod. Pathol. 23, 1052–1060. doi: 10.1038/modpathol.2010.86

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, E. Y., Cristofanilli, M., Robertson, F., Reuben, J. M., Mu, Z., Beavis, R. C., et al. (2013). Genome Wide Proteomics of ERBB2 and EGFR and Other Oncogenic Pathways in Inflammatory Breast Cancer. J. Proteome Res. 12, 2805–2817. doi: 10.1021/pr4001527

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, Q., Wang, Y., Liang, J., Tian, Y., Zhang, Y., and Tao, K. (2017). Bioinformatics analysis to identify the critical genes, microRNAs and long noncoding RNAs in melanoma. Medicine 96:e7497. doi: 10.1097/MD.0000000000007497

PubMed Abstract | CrossRef Full Text | Google Scholar