Screening and identification of key genes and pathways in metastatic uveal melanoma based on gene expression using bioinformatic analysis.

The current study aimed to elucidate the molecular mechanisms and identify the potential key genes and pathways for metastatic uveal melanoma (UM) using bioinformatics analysis.Gene expression microarray data from GSE39717 included 39 primary UM tissue samples and 2 metastatic UM tissue samples. Differentially expressed genes (DEGs) were generated using Gene Expression Omnibus 2R. Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were performed using the online Database for Annotation, Visualization and Integrated Discovery (DAVID) tool. The web-based STRING tool was adopted to construct a protein--protein interaction (PPI) network. The MCODE tool in Cytoscape was used to generate significant modules of the PPI network.A total of 213 DEGs were identified. GO and KEGG analyses revealed that the upregulated genes were mainly enriched in extracellular matrix organization and blood coagulation cascades, while the downregulated DEGs were mainly related to protein binding, negative regulation of ERK cascade, nucleus and chromatin modification, and lung and renal cell carcinoma. The most significant module was extracted from the PPI network. GO and KEGG enrichment analyses of the module revealed that the genes were mainly enriched in the extracellular region and space organization, blood coagulation process, and PI3K-Akt signaling pathway. Hub genes, including FN1, APOB, F2, SERPINC1, SERPINA1, APOA1, FGG, PROC, ITIH2, VCAN, TFPI, CXCL8, CDH2, and HP, were identified from DEGs. Survival analysis and hierarchical clustering results revealed that most of the hub genes were associated with prognosis and clinical progression.Results of this bioinformatics analysis may provide predictive biomarkers and potential candidate therapeutic targets for individuals with metastatic UM.

Introduction

Melanoma is a life-threatening malignancy and the primary intraocular form is known as uveal melanoma (UM). Among primary intraocular tumors in the adult population, UM is the most common. UM may originate from the choroid, iris, or ciliary body, which are commonly known as the uvea. In approximately 90% of UM cases, the choroid is involved.[ The biological features and clinical behavior of UM are distinct from those of cutaneous melanoma. Currently, first-line treatment for UM includes resection, radiation, and eye enucleation. These therapy options are able to control the local disease but still did not reduce the risk of distant metastases, which is a key obstacle to improve the long-term survival of UM. Despite the emergence of novel treatment modalities, such as immune checkpoint blockade, gene-targeted therapy, and anti-angiogenic therapy, the survival rates of patients with UM have not changed in the past 40 years.[

Hematogenous metastases typically involve the liver in approximately 90% of metastatic cases, the lung(s) is involved in 24% of cases, and bone in 16%.[ The median time from initial diagnosis to metastasis is approximately 2 to 3 years; once metastases occur, prognosis is typically poor, with a median survival of 2 to 3 months.[ To improve the prognosis of metastatic UM, the mechanisms of how UM metastasizes and the prognostic factors that can predict the risk for metastasis have been extensively studied. Shields et al[ reported that the thickness of UM is positively associated with increasing risk for metastasis. Schmittel et al[ found that primary UMs with a largest diameter > 14 mm and ciliary body involvement have a poor prognosis. In addition, due to the advances in molecular biology, some researchers have found that noncoding RNAs,[ aberrant alterations in chromosomes 1, 3, 6, and 8,[ and loss-of-function mutations in the BAP1 gene[ are involved in metastasis. However, metastatic mechanisms in UM are particularly complicated, and there are no clinically applicable molecular biomarkers that can accurately predict metastatic risk.

In recent decades, advances in microarray technology and bioinformatics analysis have helped to identify key gene(s) and functional pathways involved in the progression and metastasis of cancers, which have offered new insights into the molecular mechanism of metastasis in UM. Thus, in the present study, messenger RNA (mRNA) microarray datasets from the Gene Expression Omnibus (GEO) database were obtained and analyzed to identify differentially expressed genes (DEGs) between patient-derived primary UM tissues and metastatic UM tissues, followed by Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis. Subsequently, a protein--protein interaction (PPI) network was constructed to interpret the biological interaction of DEGs. Module analysis of the DEGs was performed to identify key genes and pathways related to metastatic UM. Finally, a total of 213 DEGs and 14 hub genes were identified, which may be potential prognostic markers and therapeutic targets for metastatic UM.

Materials and methods

Microarray data

The gene expression dataset GSE397127[ was downloaded from the GEO database. The GEO (http://www.ncbi.nlm.nih.gov/geo)[ is a public database of high-throughput gene expression data, chips, and microarrays. GSE397127 was based on the GPL6098 platform (Affymetrix Illumina humanRef-8 version 1.0 expression beadchip), which contained 39 primary UM tissue samples and 2 liver-metastatic UM samples. The probes were converted into official gene symbols according to the annotation information of the platform. The ethical approval was not necessary for this study, as all datasets were retrieved from a public database.

Identification of DEGs

GEO2R (https://www.ncbi.nlm.nih.gov/geo/geo2r/)[ was used to screen DEGs between primary and metastatic UM tissue samples. GEO2R is an R-based application that enables users to identify DEGs in one or more datasets. LogFC (foldchange) > 3 and adj. A list of upregulated and downregulated DEGs were saved for subsequent analysis.

GO functional and KEGG pathway enrichment analyses

GO functional enrichment is a widely used approach for interpreting sets of genes.[ The KEGG database is a collection of pathway maps representing metabolism and various other biological functions.[ As a free online bioinformatics resource, the database for annotation, visualization, and integrated discovery (DAVID, https://david.ncifcrf.gov/) provides functional annotation and visualization of large-scale lists of genes.[ In this study, DAVID was used for the enrichment of GO functions and KEGG pathways for the systematic analysis of DEGs. Differences with P < .05 were considered to be statistically significant.

PPI network of DEGs

The web-based STRING tool (https://string-db.org/) was adopted to obtain PPI relationships for the DEGs followed by visualization using Cytoscape. PPIs with a combined score > 0.4 were selected. Cytoscape is an open source software for integrating biomolecular interaction networks with high-throughput expression data into a unified conceptual framework.[ The plug-in tool Molecular Complex Detection (MCODE) (version 1.5.1) of Cytoscape was adopted to detect strongly connected regions from the PPI network with the following parameters: degree cutoff = 2, k-core = 2, node score cutoff = 0.2, maximum depth = 100. MCODE is an application designed to find densely connected regions in a specific PPI network based on topology.[

Hub gene selection and analysis

Hub gene selection was performed using cytoHubba (version 0.1), a plug-in application of Cytoscape. CytoHubba computes 11 methods to identify important nodes in PPI networks.[ Genes appearing at least twice in the top 10 results of each computation method were considered Hub genes. The analyses of clinical prognosis including overall survival (OS) and disease-free survival (DFS) of hub genes were performed using Kaplan–Meier curve analysis and analyzed using GEPIA online platform[ (http://gepia.cancer-pku.cn/), and coexpression analysis of the hub genes were performed in both cBioPortal[ and Oncomine databases (https://www.oncomine.org/).[ Hierarchical clustering of hub genes was performed using the University of California Santa Cruz (UCSC) Xena platform (https://xenabrowser.net/).[

Statistical analysis

For identification of DEGs, the Student t test was adopted and Benjamini and Hochberg method was used to adjust the P value.[ For KEGG and GO analyses, the Fisher exact test was performed to determine whether differences were significant. For the above statistical methodologies, P value less than .05 was considered statistically significant.

Results

The microarray dataset GSE39717, deposited by Harbour et al,[ was downloaded from the GEO database. A total of 24,358 genes from 39 primary UM and 2 metastatic UM patient-derived tumor tissues were obtained. A total of 213 DEGs were identified between the primary and metastatic samples, including 70 (32.9%) upregulated and 143 (67.1%) downregulated genes.

GO and KEGG enrichment analyses of DEGs

On the basis of the enrichment analysis of DEGs using DAVID, a total of 115 GO terms of upregulated genes and 34 GO terms of downregulated genes were obtained. GO analysis revealed that changes in upregulated DEGs were significantly associated with extracellular region, matrix organization, space, and blood microparticles, while the downregulated DEGs were mainly related to protein binding, nucleus, negative regulation of ERK1 and ERK2 cascade, and covalent chromatin modification chromatin binding (Table 1). As shown in the KEGG pathway enrichment analysis (Table 2), upregulated DEGs were mainly involved in complement and coagulation cascades, extracellular matrix (ECM)-receptor interaction, amebiasis, focal adhesion, and protein digestion and absorption. The pathways enriched in the downregulated DEGs were mainly nonsmall cell lung cancer and renal cell carcinoma.

Table 1

Top 5 GO functional enrichment analyses of upregulated and downregulated DEGs.

Table 2

Top enriched KEGG pathways of upregulated and downregulated DEGs.

PPI network construction and module analysis

The PPI network of DEGs was constructed using Cytoscape (Fig. 1) and the most significant module was obtained using the MCODE application. As shown in Figure 2, the most significant module (MCODE Score = 12.5) contained 25 nodes and 120 edges. GO and KEGG enrichment analyses of genes involved in this module were conducted using DAVID. GO term enrichment analysis revealed that the genes in the above module were mainly involved in extracellular region, extracellular space, endoplasmic reticulum lumen, ECM organization, blood microparticle, ECM structural constituents, platelet degranulation, collagen catabolic processes, and extracellular exosome (Table 3). The results of KEGG pathway enrichment revealed that the genes were mainly related to ECM-receptor interaction, focal adhesion, protein digestion and absorption, amebiasis, PI3K-Akt signaling pathway, complement and coagulation cascades, platelet activation, small cell lung cancer, vitamin digestion and absorption, and proteoglycans in cancer (Table 4).

Figure 1

PPI network of DEGs was constructed using STRING and visualized in Cytoscape. Upregulated genes are marked in red; downregulated genes are marked in blue.

Figure 2

The most significant module of DEGs. The most significant module containing 25 nodes and 120 edges was extracted from PPI network using MCODE. Seed node was marked in yellow.

Table 3

Top 10 GO functional enrichment analyses in the most significant DEG module.

Table 4

Top 10 KEGG pathway enrichment analyses in the most significant DEG module.

A total of 14 genes were identified as hub genes from 33 candidate genes using the cytoHubba tool in Cytoscape (Supplemental Digital Content [Table S1]). The gene symbol, full name, and brief introduction of the functions for these hub genes are listed in Table 5. As illustrated in Figure 3, the survival analysis of the hub genes was performed using Kaplan–Meier curve analysis. In the OS analysis, the UM patients with high mRNA levels of FN1, VCAN, APOA1, and PROC genes demonstrated a worse prognosis (Fig. 3A). Meanwhile, UM patients with high mRNA levels of FN1, VCAN, SERPINC1, and ITIH2 demonstrated worse DFS (Fig. 3B). According to the results of cytoHubba analysis, FN1, APOB, F2, SERPINC1, and FGG were ranked highest, which suggested their potential role in UM metastasis. The prognosis analysis results indicated that the alteration of FN1 and VCAN mRNA levels was associated with worse OS and DFS. Nevertheless, APOPA1 and PROC worsen the OS while SERPINC1 and ITIH2 reduced DFS, although no statistical significance of the reduction was observed in OS affected by APOA1 and the reduction of DFS affected by FN1, SERPINC1, and ITIH2. Furthermore, hierarchical clustering analysis the UCSC Xena platform revealed that the mRNA levels of FN1, SERPINC1, SERPINA1, VCAN, PROC, and CDH2 were basically consistent with clinical grade (Fig. 4). Coexpression analysis using cBioPortal revealed that FN1 genes were highly coexpressed with VCAN in the UM tissue (Pearson correlation, 0.70; Spearman correlation, 0.85) (Fig. 5A). The Laurent Melanoma data in Oncomine revealed that the expression of VCAN was positively related to FN1 in 3 subtypes of UM (correlation index, 0.603) (Fig. 5B).

Table 5

Full name and functional roles of 14 hub genes.

Figure 3

(A) Overall survival and (B) disease-free survival analyses of hub genes were performed using GEPIA platform.

Figure 4

Hierarchical clustering of hub genes was performed using UCSC Xena platform. Upregulation of gene is marked in red, downregulation of genes is marked in blue.

Figure 5

FN1 and VCAN expression directly correlates in human uveal melanoma cancers. (A) Coexpression analysis of FN1 and VCAN via cBioPortal platform. (B) FN1 gene coexpressed with FN1 in Laurent Melanoma data in Oncomine database, which consists of three subtypes: 1. Epitheliod Cell Uveal Melanoma 2. Mixed Cell Uveal Melanoma 3. Uveal Melanoma. The reporters indicate probes used in the analysis. Upregulation of expression is marked in red and downregulation is blue, passing by white, with fluctuating color intensity.

Discussion

UM is one of the most common intraocular malignancies in adults; 62% of UM patients have confirmed melanoma metastasis at the time of death and 92% of metastatic sites are the liver.[ Currently, the management of liver metastasis from UM includes surgery, local chemotherapy, radiotherapy, and immune-embolization. Nevertheless, treatment of metastatic UM remains a daunting challenge in clinical practice due to the very poor prognosis of these patients.[ Benefitting from updated prognostication techniques, primary UM can be classified into distinct subgroups with various levels of metastatic risk based on gene expression profile.[ In 2004, Onken et al[ proposed that mRNA levels of PHLDA1, FZD6, and ENPP2 could be used as molecular signatures to predict prognosis. However, the oncogenic and metastatic mechanisms of UM remain controversial, and advances in the treatment of UM are not promising because survival of patients with UM has remained unchanged over the past 4 decades, from 1973 to 2013.[ Hence, identification of key genes and pathways of the metastatic mechanism of UM could contribute to the diagnosis and treatment of UM.

In the present study, gene expression profiles of 39 primary UM samples and 2 metastatic UM samples were obtained from the GEO39717 dataset. A total of 213 DEGs were identified, including 70 upregulated and 143 downregulated genes. To further understand the interactions of the DEGs, GO function and KEGG pathway analyses were performed using DAVID. The upregulated genes were mainly enriched in extracellular region, matrix organization, space and blood microparticle, complement and coagulation cascades, ECM-receptor interaction, amebiasis, focal adhesion and protein digestion and absorption, while the downregulated DEGs were mainly related to protein binding, nucleus, negative regulation of ERK1 and ERK2 cascade, covalent chromatin modification chromatin binding, nonsmall cell lung cancer, and renal cell carcinoma. In the most significant module generated by MCODE, DEGs were mainly enriched in extracellular region and space organization, blood coagulation process and the PI3K-Akt signaling pathway. According to previous studies, the extracellular environment is the key driver for both cancer development and progression.[ Blood coagulation pathways play a role in tumor progression and metastasis,[ phosphor-AKT protein levels are positively associated with a higher risk for metastasis in patients with UM,[ and ERK pathway promotes carcinogenesis and maintenance of UM.[ Thus, results of all of these studies support those of the current investigation.

Hub genes, namely FN1, APOB, F2, SERPINC1, SERPINA1, APOA1, FGG, PROC, ITIH2, VCAN, TFPI, CXCL8, CDH2, and HP, were identified from the PPI network of DEGs using the cytoHubba tool, indicating these genes may be vital in the metastatic process of UM. FN1 is involved in cell adhesion, cell motility, wound healing, and maintenance of cell shape. FN1 has been shown to promote metastasis in various types of tumors.[ Recently, Li et al[ reported that FN1 promotes cutaneous melanoma proliferation and metastasis by inhibiting apoptosis and regulating epithelial-to-mesenchymal transition, which is consistent with our results. APOA1 is the main protein constituent of high-density lipoprotein, which shuttles excess cholesterol from organs to the liver for excretion. APOA1 has been described to exert anti-apoptotic, anti-inflammatory, and antioxidant activities, which are involved in tumorigenesis.[ In a murine model of malignant melanoma, APOA1 also demonstrated anti-tumor effects.[ However, in the present study, we determined that APOA1 significantly increased in liver metastatic UM samples, indicating a stimulating role of APOA1 in UM metastasis. However, we cannot exclude the possibility that the increased APOA1 mRNA in liver metastatic UM was due to the fact that APOA1 mRNA levels are higher in the liver than any other tissue in the human body.[ APOB is a major protein constituent of chylomicrons, low-density lipoprotein and very-low density lipoprotein, lung cancer and colorectal cancer risk were increased with high APOB levels,[ whereas the role of APOB in UM remains unclear. SERPINA1 and SERPINC1 are members of the serpin family, SERPINA1 was found to improve nonsmall cell lung cancer cell migration, colony formation, and resistance to apoptosis,[ while knockdown of SERPINC1 was reported to inhibit neural progenitor cell proliferation via suppression of the PI3K/Akt/mTOR signaling pathway.[ ITIH2, also known as serum-derived HA-associated protein (SHAP), forms complexes with hyaluronan (HA) to regulate the localization, synthesis, and degradation of HA in serum. Elevated serum levels of the SHAP-HA complex indicate poor prognosis in endometrial and ovarian cancers.[ VCAN plays a role in intercellular signaling and in connecting cells with the ECM. It was reported that VCAN significantly increased in superficial spreading melanoma tissue and metastatic melanoma cell lines.[ Notably, another interesting finding of hub genes was that VCAN is highly relevant to FN1. Soikkeli et al[ reported that in melanoma lymph nodes, upregulation of POSTN, FN1, COL-I, and VCAN genes was confirmed in metastatic outgrowth, and all of these genes were inducible by transforming growth factor-beta, which indicated the activation of the transforming growth factor-beta signaling pathway. Although FN1 and VCAN have been found to be positively associated with UM metastasis, the downstream protein and concrete mechanism of VCAN, as well as the interaction of FN1 and VCAN, remains unclear. TFPI encodes a serine protease inhibitor that regulates the tissue factor dependent pathway of blood coagulation. TPFI contributes to the development of multiple drug resistance in breast cancer cells.[ CXCL8, also known as interleukin-8, is a member of the CXC chemokine family and is a mediator of the inflammatory response by regulating cancer stem cell proliferation and self-renewal, the CXCL8-CXCR1/2 axis may play an important role in tumor progression and metastasis.[ In melanoma, NFAT1 regulates CXCL-8/MMP3 and promotes tumor growth and lung metastasis.[ CDH2, also known as N-cadherin, belongs to the cadherin superfamily, and mediates calcium-dependent cell-cell adhesion. Elevated CDH2 is a well-known protein marker for the onset of epithelial-mesenchymal transition, which results in enhanced migratory capacity, invasiveness, and increased resistance to apoptosis in many types of cancers.[ In melanoma cells, increased N-cadherin expression contributes to proliferation and invasive potential by activating PI3/AKT, mTOR, and ERK kinase.[ The HP gene encodes haptoglobin, which combines with free plasma hemoglobin, thus enabling heme iron to be recycled in hepatocytes. Previous research found that cellular levels of HP are strongly associated with the recurrence rate of human head and neck cancers.[ A positive correlation between elevated serum haptoglobin level and the incidence of colorectal cancer was also observed.[

In addition, we performed hierarchical clustering and prognosis analysis for hub genes. The hierarchical clustering results illustrated that, as the clinical stage of UM increased, most of the hub gene mRNA levels also increased, indicating the consistency between hub gene expression and UM tumor progression. In addition, OS and DFS analysis of the hub genes demonstrated that high expression of FN1 and VCAN was related to worse OS and DFS, increased APOA1 and PROC reduced OS, while SERPINC1 and ITIH2 reduced DFS. Analysis of hub genes demonstrated that these genes may play an important role(s) in the progression, invasion, and metastasis of UM, and may be potential candidates for prognosis prediction and diagnostic biomarkers.

Finally, there were several limitations to the current study. First, all of the data were obtained from the GEO database rather than directly from UM patient tissues. Second, all conclusions were based on bioinformatics analysis; hence, caution must be exercised in interpreting the results, being aware that experimental verification is a better approach to confirm findings. Third, the GSE 39717 dataset consisted of 39 primary tumor samples and 2 metastatic samples, the imbalance between groups may have unintentionally introduced biases. In summary, larger-scale tissue samples derived from a primary and metastatic UM patient cohort with confirmatory experiments need to be performed to verify our conclusions.

Conclusion

The present bioinformatic analysis identified key genes and molecular pathways possibly involved in the metastatic process of UM. A total of 213 DEGs and 14 hub genes were identified to play crucial roles in the progression, invasion, and metastasis of UM, and could be potential candidates as diagnostic biomarkers.

Acknowledgments

We thank Wolters Kluwer for English language editing.

Author contributions

Conceptualization: Jiehui Xu, Jialu Xie.

Formal analysis: Zhenyu Wu.

Investigation: Jiehui Xu.

Project administration: Jiehui Xu.

Data curation: Jialu Xie, Xiaogang Xu.

Methodology: Jialu Xie, Zhenyu Wu.

Visualization: Jialu Xie.

Writing – original draft: Jialu Xie, Guanlu Liang.

Writing – review & editing: Jialu Xie, Xiaogang Xu.

Correction

When originally published, Dr. Xiagang Xu's degree appeared incorrectly as MD. It has been corrected to PhD.