Identifying a Potential Key Gene, TIMP1, Associated with Liver Metastases of Uveal Melanoma by Weight Gene Co-Expression Network Analysis

Introduction

Uveal melanoma (UM) is the most common primary intraocular malignancy in adults and can be fatal due to metastasis.¹ Up to 50% of UM patients develop distant metastasis through hematogenous spread, and the most common site of metastasis is the liver, with survival time less than 6 months after the diagnosis of metastasis.² Many effective local therapies have been developed for primary tumors, such as transpupillary thermotherapy, radiation therapy, and enucleation, while the poor prognosis for patients with metastasis has not improved.^3,4 Several clinicopathologic features are associated with an increased risk of metastasis, including increased patient age, increased tumor thickness and diameter, epithelioid cell type, and 3-chromosome monomer.^5–8 Several studies have identified a gene expression-based classification of primary UM that accurately predicts metastatic death.^9,10 Studies focused on tumor microenvironment have indicated that increased growth factors in the liver were associated with metastatic progression of UM, such as hepatocyte growth factor (HGF),^11,12 insulin-like growth factor 1,^13,14 and stromal cell-derived factor 1 (SDF-1).¹⁵ In addition, C-X-C chemokine receptor type 4 (CXCR4) is associated with the organ-specific homing of UM to the liver and is involved in the process of metastasis.^15,16 Currently, as metastatic UM patients are not responsive to existing radiotherapy and chemotherapy, it is important to understand the potential mechanism of liver metastasis to improve patient prognosis.

With the development of microarray and high‐throughput sequencing technology, bioinformatics has been used in recent studies to analyze genetic changes in tumors.^17,18 However, most studies compare differences in gene expression levels between different samples to identify potential biomarkers, ignoring the relationship between genes and clinical traits.^19,20 Weighted gene co-expression network analysis (WGCNA) is a bioinformatics analysis method, which allows cluster co-expressed genes with similar expression profiles in different samples and investigates the relationship between gene sets and clinical features.²¹ By developing the co-expression network, genes with similar expression profiles are classified into the same module. Recently, WGCNA has been widely applied to identify the important genes involved in disease processes, such as pancreatic cancer,²¹ clear cell renal cell carcinoma,²² endometrial carcinoma,²³ and systemic sclerosis-related pulmonary hypertension.²⁴

In this study, the WGCNA algorithm was used to analyze the microarray data downloaded from public datasets about UM metastasis to explore the metastasis-related genes. The genes may help to predict the malignant potential and prognosis of UM.

Materials and Methods

Dataset Collection and Study Design

The data of GSE22138 were obtained from Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo)²⁵ containing 35 metastatic UM tissues and 28 primary uveal melanoma samples.²⁶ The platform of GSE22138 was GPL570 [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array. In this study, the data of 63 tissues were analyzed; the data analysis process was presented (Figure 1A).

Figure 1 Description of the study, DEGs and clustering dendrogram of tumor samples. (A) Flow Diagram of this study. (B) Volcano plot of DEGs in GSE22138 (35 metastatic UM and 28 UM). The red dot refers to up-regulated genes, the green dot refers to down-regulated genes, and the black dot refers to unchanged genes. DEGs: differentially expressed genes. (C) GO and KEGG enrichment analysis of DEGs. (D) Sample dendrogram of 63 UM samples and their clinical characteristics. The clustering is based on DEGs data in metastatic UM samples compared to UM tissues. The color intensity was proportional to metastasis, age, months to endpoint, tumor diameter, tumor thickness, gender and extra-scleral extension.

Data Preprocessing and DEGs Identification

The “Affy” R package was applied to preprocess the raw expression data.²⁷ First, the raw data were background corrected and standardized by Robust Multiarray Average (RMA) analysis. Then, median polish probe sets annotated by the Affymetrix annotation files were summarized. The “Limma” R package was applied to screen the DEGs between uveal melanoma tissues with metastasis and without metastasis.²⁸ The selection standard for identifying DEGs was |log2(fold change) |>0.585 and P value < 0.05.

GO and KEGG Pathway Enrichment Analysis of DEGs

For further insight into the underlying mechanism of the DEGs of UM, the DEGs were uploaded to the Database for Annotation, Visualization, and Integrated Discovery (DAVID).^29,30 Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were performed using the DAVID database. P value < 0.05 was used to determine statistical significance.

Construction of Weighted Gene Co-Expression Network

The “WGCNA” package in R was used to construct a co-expression network for the 552 DEGs from 63 UM samples and the corresponding clinical data.³¹ First, a weighted adjacency matrix was constructed using Pearson’s correlations between gene pairs. Next, a soft threshold parameter β was chosen based on standard scale-free networks, and adjacencies of differential genes were analyzed applying a power function.^32,33 Finally, a topological overlap matrix (TOM) was built based on the adjacency.

An average linkage hierarchical clustering was performed according to a TOM dissimilarity measure. The genes dendrogram was built with the parameters of a minimum size of 30. The module dendrogram and merged modules were constructed with a cut-line of 0.25.

Identification of Interesting Modules of UM

The module–trait relationship analysis was applied to calculate the correlation between module eigengenes (MEs) and clinical features. The ME was the main component of a module and reflected all the characteristics of module genes. By calculating the correlation between MEs and clinical features, the interesting modules were determined. The assessed clinical phenotypes included metastasis, age, months to endpoint, tumor diameter, tumor thickness, gender, and extra-scleral extension. The interesting module mostly correlated with the clinical phenotypes of metastasis was selected to analyze subsequently.

Hub Module Functional Analysis

Using the ClueGO (version 2.3.5) plugin in Cytoscape (version 3.6.1),^34,35 which can visualize the biological terms for large clusters of genes, GO enrichment analysis of the genes in hub module was performed. ClueGO analysis was performed a threshold of P value < 0.05.

Module Visualize and Identification of Hub Genes

A gene with high interconnection with nodes in the interesting module was considered to be functionally significant. The network visualization of the interesting module was performed by Cytoscape software. Meanwhile, the Molecular Complex Detection (MCODE) plug was used to select the core network with degree cutoff = 2, node score cutoff = 0.2, k‐core = 2, and maximum depth = 100.³⁶

Identification of Candidate Hub Genes in Key Modules

Hub genes with a high interconnective degree in a module were considered to play an important role in the network. Genes from hub module were mapped in the Search Tool for the Retrieval of Interacting Genes (STRING) database to construct protein–protein interaction (PPI) analysis.³⁷ Further, cytoscape software was applied to establish a PPI network.³⁵ A confidence score >0.40 was chosen as the cut-off criterion. The common hub genes chosen from co-expression network and PPI network were selected for further analysis.

Gene Set Enrichment Analysis (GSEA)

Based on the expression level of the hub gene, 63 UM samples of GSE22138 were divided into two groups. To explore the potential function of TIMP1 in the progression of UM, gene set enrichment analysis (GSEA) was conducted with the cut-off criteria of P value < 0.05.³⁸

Kaplan–Meier Survival Analysis of TIMP1

To assess the association between the expression of TIMP1 and clinical outcomes of UM, the Gene Expression Profiling Interactive Analysis database (GEPIA, http://gepia.cancer-pku.cn/), an online database based on The Cancer Genome Atlas (TCGA), was used to analyze overall survival. UM patients were divided into low-expression and high-expression groups according to the expression level of TIMP1. The hazard ratio (HR) with 95% confidence intervals and log-rank P value were calculated.

Quantitative Real-Time Polymerase Chain Reaction (qRT-PCR)

RNA was extracted from cells with TRIzol reagent (Invitrogen, USA). The complementary DNA was synthesized by utilizing Easy Quick RT Master Mix (Cwbiotech, Beijing, China). RT-PCR was performed by applying the ABI7500 Real-time PCR System (Applied Biosystem, CA, USA) and AceQ^® qPCR SYBR Green Master Mix (Vazyme, Nanjing, China). The relative expression of the target gene was calculated based on the formula 2^−ΔCt with GAPDH setting as the internal control. Primer sequences of TIMP1 were 5′-CTGTTGTTGCTGTGGCTGAT −3′(Forward); 5′-AACTTGGCCCTGATGACG −3′(Reverse). Primer sequences of internal control GAPDH were 5′- TCATTGACCTCAACTACATGG −3′(Forward); 5′-TCGCTCCTGGAAGATGGTG −3′ (Reverse).

Cell Culture and Plasmid Transient Transfection

The human normal retinal pigment epithelium cell line ARPE19 and UM cell lines C918 and OCM-1 were purchased from Procell (PROCELL Life Science & Technology Co., Ltd., Wuhan, China). ARPE19 cells were cultured in F-12 supplemented with 10% fetal bovine serum (FBS; Gibco, Carlsbad, CA, USA). The UM cell lines were cultured in Roswell Park Memorial Institute 1640 (RPMI-1640; Gibco, Carlsbad, CA, USA) medium with 10% FBS. Cell lines were maintained at 37°C in a humidified 5% CO₂ atmosphere. The full-length TIMP1 sequence was inserted into the pcDNA3.1 plasmid (Gentec, Shanghai, China) to build the pcDNA3.1-TIMP1 overexpression plasmid. For transfection, Lipo8000TM transfection reagent (Beyotime, Shanghai, China) was used following manufacturer’s instructions. The control group included cells transfected with empty pcDNA3.1 plasmid. Cells were transfected with 2 μg of the plasmid for 24 hours before conducting further experiments.

Cell Proliferation Assay

Cell proliferation was analyzed by applying the Cell Counting Kit-8 kit (Dojindo, Kumamoto, Japan). Cells (2×10³/well) were seeded into 96-well plates with 0.2 mL cell suspension. Cells of each well were cultivated with 0.01 mL CCK-8 solution for 2 hours before cell harvest. Cell proliferation was assessed at four time points (24, 48, 72, and 96 hours) by analyzing the absorbance at 450 nm by a microplate reader (Bio-Rad, USA).

Scratch Healing Assay

Cells were seeded into 6-well plates and cultured to form a cell monolayer with 85–95% confluence. A straight scratch wound was made at the center of each well using a 200-μL pipette tip. Wells were washed with PBS 3 times to get rid of the unattached cells and cultured in the medium without serum. A microscope was used to observe the wound healing and take photos at three time points (0, 24, and 48 hours). Scratch areas were measured by Image-J software.

Invasion Assay

A transwell chamber with a 24-well plate was applied to perform the invasion assay, and the inner membranes were covered with 25 μL of Matrigel (Becton Dickinson Bioscience) diluted by a medium. Then, 4 × 10⁴ transfected cells in 200-μL medium without serum were added to the upper, while RPMI-1640 medium with 10% FBS was added to the lower surface for 24 hours of incubation. Cells on the lower surface were fixed with methanol for 10 minutes and then stained with 0.1% crystal violet. The Matrigel and cells on the upper surface were wiped away with cotton swabs. Cells that invaded the lower surface were observed and photographed by a microscope.

Statistical Analysis

Data were expressed as mean ± standard deviation (SD). SPSS 20.0 software was used for a student T test to analyze the data between groups. A P value <0.05 was used to demonstrate a significant difference.

Results

Identifying DEGs

The RNA sequencing data of GSE22138, including 35 metastatic uveal melanoma samples and 28 uveal melanoma tissues, were analyzed. The “Limma” package of R software was used to identify DEGs, based on the threshold of |log2(fold change) |>0.585 and a P value < 0.05. A total of 552 genes were screened as DEGs, with 358 upregulated and 194 downregulated (Figure 1B).

GO and KEGG Pathway Enrichment Analysis of DEGs

To determine the function of DEGs of UM, GO and KEGG pathway enrichment analyses were conducted by DAVID database. The DEGs were enriched in biological processes (such as Notch signaling pathway, response to drug, oocyte maturation, cellular senescence, and response to interferon-gamma), cellular components (such as endoplasmic reticulum membrane, extracellular exosome, early endosome, cytosol and extracellular matrix), and molecular functions (such as receptor binding, protein binding, protein homodimerization activity, gamma-tubulin binding and identical protein binding). The KEGG pathway analysis indicated that the DEGs were mainly enriched in three pathways including glutathione metabolism, focal adhesion, and platelet activation (Figure 1C).

Construction of Weighted Gene Co-Expression Network

A co-expression network was built applying dataset GSE22138 including 63 UM samples and the corresponding clinical data. The “WGCNA” package in R was used to construct co-expression network for the 552 DEGs (Figure 1D). In this study, the power of β=8 (scale‐free R2=0.8520) was chosen to ensure a free scale network (Figure 2A). The DEGs were put into modules according to similar expression patterns by average linkage clustering, and five modules were identified completely (Figure 2B).

Figure 2 Identification of soft-thresholding power in WGCNA and modules associated with the clinical traits of UM. (A) Determination of soft-thresholding parameter in the weighted gene co-expression network analysis (WGCNA). Analysis of the scale-free fit index for various soft-thresholding powers (β) and the mean connectivity for various soft-thresholding powers. (B) Dendrogram of all DEGs clustered based on a dissimilarity measure (1-TOM). (C) Heatmap of the correlation between module eigengene and clinical traits of UM. Each cell contains the corresponding correlation and P value. (D) GO functional enrichment analysis of the genes in hub module by using ClueGO. Small nodes labeled red represent genes, big nodes represent genes GO terms. Different colors of big nodes represent the functional annotation of ontologies.

Identification of Interesting Modules of UM

The module significantly correlated with metastasis is of great value in predicting the progression of UM. The turquoise module showed the highest correlation with metastasis (r = 0.56, P = 2e-06, Figure 2C), which was selected for further analysis.

Hub Module Functional Analysis

The GO functional enrichment analysis of the genes in hub module was performed by “ClueGO.” The results suggested that the genes in hub module are enriched in peptidyl-tyrosine autophosphorylation, circadian regulation of gene expression, positive regulation of interleukin-8 production, protein kinase C binding, response to retinoic acid, regulation of organ growth, lymphocyte costimulation, regulation of cardiac muscle tissue growth, and T cell costimulation (Figure 2D).

Module Visualize and Identification of Hub Genes

Module genes with high interconnection with nodes were considered to be functionally significant. The module was visualized by using the Cytoscape software (Figure 3A). And the core network was selected by applying plug-in MCODE (Figure 3B).

Figure 3 Identification of hub genes via WGCNA and PPI network. (A) The network visualization of the turquoise module by Cytoscape software. (B) The most significant module of the network by MCODE. Nodes represent genes, and node color intensity is proportional to connectivity of the gene (positive correlation in red and negative correlation in blue). (C) PPI network of genes in the turquoise module analyzed by STRING database. Nodes represent genes, and node color intensity is proportional to connectivity of the gene by degree (positive correlation in red and negative correlation in blue).

Identification of Hub Gene in Key Module

A PPI network of genes from hub module was built using STRING database (Figure 3C). The top three hub genes were identified based on the degree of connectivity, which met the standard of both co-expression network and PPI network. They were TIMP1 (tissue inhibitor of metalloproteinase1), GPAA1 (glycosylphosphatidylinositol anchor attachment 1), and biglycan (BGN).

TIMP1 has been shown to play an important role in the progression of many types of human cancers. However, the biological function of TIMP1 in UM remains unclear. Next, the role of TIMP1 in UM cell lines was further studied.

Gene Set Enrichment Analysis

Four gene sets were enriched, including “P53 pathway,” “Myogenesis,” “Coagulation,” and “Estrogen response early” (Figure 4).

Figure 4 Gene set enrichment analysis (GSEA). “P53 pathway,” “Myogenesis,” “Coagulation,” and “Estrogen response early” were enriched in UM samples with highly expressed TIMP1.

Identification of TIMP1 as a Potential Biomarker for the Worse Prognosis of UM Patients

The overall survival (OS) of UM patients with high TIMP1 expression was significantly shortened (Figure 5A). The qRT-PCR analysis indicated that the expression of TIMP1 in UM cells was significantly increased compared with ARPE19 cells (Figure 5B). To investigate the biological function of TIMP1 in UM, we established a model of TIMP1 sufficiency in UM cell lines. Overexpression of TIMP1 after transfection with plasmid was confirmed by qRT-PCR analysis in C918 and OCM-1 cells (Figure 5C). CCK-8 assay results showed that overexpression of TIMP1 significantly inhibited the proliferation of UM cells compared with the controls (Figure 5D). Scratch healing assay results revealed that overexpression of TIMP1 inhibited the migration of UM cells (Figure 5E and F). Transwell assay results indicated that there was no significant difference in the invasion of UM cells between overexpression of TIMP1 group and the empty vector group (Figure 5G and H).

Figure 5 Identification of TIMP1 as a potential biomarker for the worse prognosis of patients with UM. (A) Kaplan–Meier survival curve of overall survival time for TIMP1 based on TCGA data in the GEPIA database. (B) TIMP1 expression in different UM cell lines was measured by qRT-PCR analysis. TIMP1 was up-regulated in UM cells compared with ARPE19 cells. (C) Overexpression of TIMP1 in C918 and OCM-1 cells was confirmed by qRT-PCR analysis. (D) Overexpression of TIMP1 inhibited the proliferation of UM cells assessed by CCK-8 assay. (E and F) Overexpressed TIMP1 inhibited the migration of UM cells assessed by scratch healing assay. The wound area was measured by ImageJ software. (G and H) Transwell assays revealed that overexpressed TIMP1 had no effect on UM cell invasion. Data were shown as mean ± SD. *p < 0.05; **p < 0.01; ***p < 0.001; ****p < 0.0001; ns indicates p>0.05.

Discussion

In this study, TIMP1 is up-regulated in the metastasis progression of UM and can be a potential prognosis biomarker based on the comprehensive bioinformatics analysis. A total of 552 genes were screened as DEGs. Then, GO and KEGG pathway enrichment analyses were performed for further insight into the underlying mechanism of the DEGs. The turquoise module showed the highest correlation with metastasis of UM by applying WGCNA. Finally, the common hub gene, TIMP1, was selected from co-expression network and PPI network for further verification.

Tissue inhibitor of metalloproteinase 1 (TIMP1) is a dual-domain protein composed of a larger amino-terminal domain and a smaller carboxy terminal domain, named for the function of mediating extracellular matrix transformation by inhibiting matrix metalloproteinases (MMPs).^39,40 The two main functions of TIMP1 are to regulate wound healing and participate in tumor progression.^41,42

Previous studies have shown that MMPs play a central role in the process of tumor invasion and metastasis,⁴³ and TIMP1, as an effective natural MMP inhibitor, can inhibit tumor process.⁴² Prior studies have shown that the high expression of TIMP1 is associated with poor prognosis of several cancer types, such as breast,⁴⁴ lung,⁴⁵ endometrial,⁴⁶ and pancreatic cancer.⁴⁷

Previous experimental studies demonstrate that overexpression of TIMP1 inhibits the tumorigenesis and migration of various types of tumors. In transgenic mouse hepatoma models, overexpression of TIMP1 inhibited T-antigen-induced hepatocyte proliferation.⁴⁸ When the B16-F10 melanoma cell line is subcutaneously implanted into mice, overexpression of TIMP1 inhibits the tumor growth.⁴⁹ It has been reported that overexpression of TIMP1 can increase intercellular contact and inhibit the migration of human hepatoma cells.⁵⁰ Similarly, by analyzing the sequencing data of clinical specimens, our results also showed that high expression of TIMP1 is associated with poor prognosis of UM. In addition, consistent with previous experimental studies, our results showed that overexpression of TIMP1 inhibited tumorigenesis and migration of UM.

There are several possible explanations for the conflicting results observed between clinical and experimental studies. First, the degradation of extracellular matrix (ECM) is an important step in the acquisition of tumor invasion and metastasis ability.⁵¹ MMPs are involved in this important process.⁵² The dynamic balance between activated MMPs and free TIMP levels affects the degradation degree of ECM and thus influence tumor progression.⁴⁵ Previous studies have shown that elevated MMP9 levels are associated with a higher risk of UM metastasis.⁵³ Therefore, the increased expression of TIMPs may reflect the body’s response to tumor invasion and metastasis, in order to maintain the dynamic balance between MMP and TIMP, so as to maintain the ECM integrity.⁵⁴ Second, experiments manipulating of TIMP-1 in tumor cell lines have a limited ability to explain the entire pathological process of the body. Recent studies have shown that TIMP1 has been identified as a pro-metastasis factor that creates a pre-metastasis microenvironment in the liver.^55,56 Hepatic stellate cells were activated by TIMP1 via CD63 signaling and attracted neutrophils that promote metastasis,⁵⁷ which greatly increases the liver’s susceptibility to circulating tumor cells.⁵⁶ These findings may explain the association between increased TIMP1 levels and poor prognosis in patients with tumors.

To further explore the possible mechanism of TIMP1 in UM metastasis, we conducted GSEA analysis and found that the high expression of TIMP1 was mainly related to the p53 signaling pathway. Previous studies have shown that the p53 signaling pathway plays an important role in promoting the occurrence, apoptosis, and progression of many types of tumors.^58,59 These findings suggested that TIMP1 might affect UM metastasis by regulating the p53 signaling pathway.

Conclusion

In conclusion, by constructing a co-expression network, our results showed that TIMP1 might be associated with metastasis in UM, which might have important significance for identifying patients with high risk of metastasis and predicting the prognosis of UM.