Identification of Metastasis-Associated MicroRNAs in Metastatic Melanoma by miRNA Expression Profile and Experimental Validation

It is reported that microRNAs (miRNA) have paramount functions in many cellular biological processes, development, metabolism, differentiation, survival, proliferation, and apoptosis included, some of which are involved in metastasis of tumors, such as melanoma. Here, three metastasis-associated miRNAs, miR-18a-5p (upregulated), miR-155-5p (downregulated), and miR-93-5p (upregulated), were identified from a total of 63 different expression miRNAs (DEMs) in metastatic melanoma compared with primary melanoma. We predicted 262 target genes of miR-18a-5p, 904 miR-155-5p target genes, and 1220 miR-93-5p target genes. They participated in pathways concerning melanoma, such as TNF signaling pathway, pathways in cancer, FoxO signaling pathway, cell cycle, Hippo signaling pathway, and TGF-beta signaling pathway. We identified the top 10 hub nodes whose degrees were higher for each survival-associated miRNA as hub genes through constructing the PPI network. Using the selected miRNA and the hub genes, we constructed the miRNA-hub gene network, and PTEN and CCND1 were found to be regulated by all three miRNAs. Of note, miR-155-5p was obviously downregulated in metastatic melanoma tissues, and miR-18a-5p and miR-93-5p were obviously regulated positively in metastatic melanoma tissues. In validating experiments, miR-155-5p's overexpression inhibited miR-18a-5p's and miR-93-5p's expression, which could all significantly reduce SK-MEL-28 cells' invasive ability. Finally, miR-93-5p and its potential target gene UBC were selected for further validation. We found that miR-93-5p's inhibition could reduce SK-MEL-28 cell's invasive ability through upregulated the expression of UBC, and the anti-invasive effect was reserved by downregulation of UBC. The results show that the selected three metastasis-associated miRNAs participate in the process of melanoma metastasis via regulating their target genes, providing a potential molecular mechanism for this disease.

It is reported that microRNAs (miRNA) have paramount functions in many cellular biological processes, development, metabolism, differentiation, survival, proliferation, and apoptosis included, some of which are involved in metastasis of tumors, such as melanoma. Here, three metastasis-associated miRNAs, miR-18a-5p (upregulated), miR-155-5p (downregulated), and miR-93-5p (upregulated), were identified from a total of 63 different expression miRNAs (DEMs) in metastatic melanoma compared with primary melanoma. We predicted 262 target genes of miR-18a-5p, 904 miR-155-5p target genes, and 1220 miR-93-5p target genes. They participated in pathways concerning melanoma, such as TNF signaling pathway, pathways in cancer, FoxO signaling pathway, cell cycle, Hippo signaling pathway, and TGF-beta signaling pathway. We identified the top 10 hub nodes whose degrees were higher for each survival-associated miRNA as hub genes through constructing the PPI network. Using the selected miRNA and the hub genes, we constructed the miRNA-hub gene network, and PTEN and CCND1 were found to be regulated by all three miRNAs. Of note, miR-155-5p was obviously downregulated in metastatic melanoma tissues, and miR-18a-5p and miR-93-5p were obviously regulated positively in metastatic melanoma tissues. In validating experiments, miR-155-5p's overexpression inhibited miR-18a-5p's and miR-93-5p's expression, which could all significantly reduce SK-MEL-28 cells' invasive ability. Finally, miR-93-5p and its potential target gene UBC were selected for further validation. We found that miR-93-5p's inhibition could reduce SK-MEL-28 cell's invasive ability through upregulated the expression of UBC, and the anti-invasive effect was reserved by downregulation of UBC. The results show that the selected three metastasis-associated miRNAs participate in the process of melanoma metastasis via regulating their target genes, providing a potential molecular mechanism for this disease. Keywords: miRNA,melanoma,metastasis

INTRODUCTION
As the neoplasm of the cells, melanoma starts in skin cells called melanocytes (McComiskey et al., 2015). Environmental factors, for example, ultraviolet light exposure, are considered as the main cause of melanoma (Kanavy and Gerstenblith, 2011). This tumor is predominant in the skin or adjacent to the skin and spread throughout the body (Bakkal et al., 2015), with a dramatically increased global incidence over the past few decades (Azoury and Lange, 2014). Moreover, there are still no satisfactory treatments for patients with advanced melanoma because of its complex pathogenesis. Like other types of malignant tumors, the main cause of melanoma-related deaths is metastasis (Shaikh et al., 2016). Thus, the pathogenesis and inhibition of metastasis have become a focal point of melanoma.
Recently, microarray technology has been widely used in investigating gene alterations in metastasis, tumorigenesis, drug resistance, and cancer recurrence, as well as to identify biomarkers for tumor diagnosis, prognosis, and therapy (Hu et al., 2010;Chien et al., 2015;Tabaries et al., 2015). With the nextgeneration sequencing technology's experimental application, Gene Expression Omnibus (GEO) collects data exponentially and gradually plays an important role in bioinformatics analysis (Barrett et al., 2013). Increasing reliable and functional miRNAs, as reported, have a vital influence on melanoma initiation, progression, and recurrence (Lou et al., 2018). DNA chipbased sequencing technology analyzes TP53 germline mutations in pediatric tumor patients, and simultaneously analyzes all coding exons of TP53 (Harris and McCormick, 2010). The entire sequencing process and data analysis are carried out within 24 h. Tissue microarray (TMA) has been successfully used in the immunohistochemical study of cervical adenocarcinoma (Tawfik El-Mansi and Williams, 2006). The expression of osteoblastspecific factor 2 in the prostate cancer-related stroma was analyzed by laser capture microdissection of clinical specimens from prostate cancer patients by whole genome expression microarray technology (Wu et al., 2005).
In our study, we selected three metastasis-associated miRNAs from miRNA expression microarray GSE24996. Data mining, network analysis, and experimental validation were applied, and the analytical and experimental results showed potential molecular mechanisms on metastatic melanoma.

Data Collection
From Gene Expression Omnibus (GEO, https://www.ncbi. nlm.nih.gov/geo/) (Li et al., 2016), we obtained the miRNA expression profile GSE24996 on the basis of the platform of GPL6955 (Agilent-016436 Human miRNA Microarray 1.0). Eight metastatic melanoma samples were screened out and analyzed. We retrieved 15 primary melanoma samples as control.

Screening for DEMs
For the analysis of DEMs, we first normalized the miRNA expression data through making use of the normalizeBetweenArray function from R package LIMMA (Smyth et al., 2005). Before and after normalizing the data, we showed them in Supplementary Figure 1. Then, the normalized data were used to analyze DEMs through making use of the Limma software package in R software (www.bioconductor.org/packages/release/bioc/html/limma.html) (Kerr, 2003). We set the cut-off value as P < 0.05 and |fold change (FC)|>2 for DEM analysis (Li et al., 2016).

Survival Data From OncoLnc Website
The prognostic values of identified DEMs in metastatic melanoma were obtained on OncoLnc website, which linked TCGA survival data to the expression levels of miRNA, mRNA, or lncRNA(http://www.oncolnc.org/). We thought Logrank P < 0.05 had statistical significance. We selected the survival-associated miRNAs for further study.

Predicting Target Genes
Through making use of miRTarBase (http://mirtarbase.mbc.nctu. edu.tw/php/index.php), we predicted latent targeted genes of the selected miRNAs. The tool is a microRNA-target interactions database that is experimentally validated (Chou et al., 2018).

Analysis of GO and Pathway
As a major initiative about bioinformatics, Gene ontology (GO: www.geneontology.org) (Ashburner et al., 2000) encompasses the largest variety of annotations under three headings: cellular component (CC), biological processes (BP), and molecular function (MF). The Kyoto Encyclopedia of Genes and Genomes (KEGG:www.genome.ad.jp/KEGG) (Kanehisa and Goto, 2000) pathway enrichment analysis was applied aiming at investigating the signaling pathways in relation to the unique DEGs. We performed the analysis of GO and KEGG pathway by making use of the Database for Annotation Visualization and Integrated Discovery (DAVID: www.david.ncifcrf.gov/) aiming at identifying the biological significance of DEGs (Huang et al., 2007). We considered P < 0.05 as statistical significance.

Construction of PPI Network and miRNA-Gene Network
We first mapped three groups of predicted genes to the Search Tool for the Retrieval of Interacting Genes (STRING) (www.cytoscape.org) (Shannon et al., 2003) to assess functional associations among them, respectively. Next, we analyzed the degree of connectivity in the PPI networks through making use of NetworkAnalyzer module in Cytoscape software and the top 10 higher degree nodes were used as hub genes to construct the miRNA-hub gene network.

Human Melanoma Samples and Cell Line
We gained 15 metastatic melanoma samples as well as 18 human primary melanoma samples from the Longhua Hospital Affiliated to Shanghai University of Traditional Chinese Medicine. After resection, we froze all specimens in nitrogen, which was liquid at once, and preserved them at −80 • C until use. We kept human melanoma cell line SK-MEL-28 (Cat. No. HTB-72, ATCC) in our lab and grew them in  The top 10 enriched biological processes of miR-155-5p, miR-18a-5p, and miR-93-5p; (B,E,H) The top 10 enriched cellular components of miR-155-5p, miR-18a-5p, and miR-93-5p; (C,F,I) The top 10 enriched molecular function of miR-155-5p, miR-18a-5p, and miR-93-5p. (J-L) KEGG pathway analysis for miR-155-5p, miR-18a-5p, and miR-93-5p. The black lines stand for gene count and the red lines stand for-log 10 (P-value).
5% CO 2 in DMEM (Gibco BRL) to which 10% (v/v) FBS (Hyclone) served as a supplement, and we keep the temperature at 37 • C.
2 − CT method, we aimed at determining fold change in each sample's RNA level in comparison with the sample for reference.

Transwell Assay
In 24-well transwell chambers (Corning, USA), we conducted the cell invasion experiment. We suspended 1 × 10 5 transfected cells in 0.2 ml serum-free medium. We added them to the inserts, which were previously coated with Matrigel (BD Bioscience, USA) on the upper surface. Then, we added 0.6 ml medium with 20% FBS to the lower compartment. Forty-eight hours later, after incubation at 37 • C, we carefully removed the cells on the membrane's upper surface through making use of a cotton bud. Then we successively fixed cells on the lower surface through making use of 100% methanol and stained them through making use of 0.1% crystal violet. We randomly selected 5 ocular fields of 200× magnification of each insert. Under a light microscope (Olympus, Japan), we enumerated them.

Dual Luciferase Assay
The dual luciferase assay used Liofectamine 2,000 to co-transfect UBC 3 ′ -UTR-wt/mut reporter plasmid and miR-93-5p mimic or NC into SK-MEL-28 cells. For 24 h, we used FLUOstar Omega to measure firefly fluorescence Reninase and renin luciferase activity, and renin luciferase activity as a reference. Dual-Lucy Assay Kit (Cat Nos. D0010) was purchased from Beijing Solabo Technology Co., Ltd.

Western Blot Analysis
We collected 2 groups of SK-MEL-28 cells (miR-93-5p mimic and control). The total protein was extracted by conventional methods, and the protein concentration was determined by the BCA method. After boiled denaturation, we took the same amount of protein-like behavior, SDS-PAGE, pass the membrane, room temperature 5%, milk block, add antibody UBC (Abcam), and incubated overnight. After 30 min of rewarming, TBST was rinsed four times (10 min/time), horseradish peroxidase-labeled two anti-incubated for 1.5 h, and TBST was washed with ECL chemiluminescence method and Quantity One 4.6.2 software for image analysis. β-actin is the reference for this experiment.

Statistical Analysis
We showed the results as mean ± SD. Through making use of the unpaired Student's t-test, we estimated differences between groups. We considered a two-tailed value of P < 0.05 as statistical significance.

Identification of Survival-Associated DEMs and Their Target Genes
Aiming at identifying DEMs between samples of metastatic melanoma and samples of primary melanoma, we analyzed the differential expression through making use of limma software package. In total, we identified 63 DEMs that were significantly differentially expressed in metastatic melanoma tissues in comparison with primary melanoma tissues, consisting of 50 upregulated and 13 downregulated miRNAs (Supplementary Table 1). To screen out the survivalassociated miRNA, all the 63 DEMs were analyzed on OncoLnc website. Finally, we picked out miR-18a-5p (upregulated), miR-155-5p (downregulated), and miR-93-5p (upregulated) as survival-associated miRNAs (Figure 1).

GO Function and KEGG Pathway Analysis
We analyzed three kinds of GO functional annotation on the three groups of predicted genes, cellular component (CC), molecular function (MF), and biological process (BP) included. The GO BP analysis results demonstrated that miR-155-5p target genes were obviously abundant in cell-cell adhesion, transcription regulated positively/negatively from RNA polymerase II promoter, an apoptotic process regulated negatively, and so on (Figure 2A). The GO CC analysis indicated that miR-155-5p target genes were abundant in membrane, nucleus, nucleoplasm, and cytosol ( Figure 2B). The GO MF analysis results indicated that miR-155-5p target genes were obviously abundant in protein binding, cadherin binding participating in protein kinase binding, and cell-cell adhesion ( Figure 2C). We showed the GO analysis results of miR-93-5p and miR-18a-5p in Figures 2D-I and the detailed results were presented in Supplementary Table 3.
Aiming at further analyzing the enriched pathways of the three groups of target genes, KEGG pathway enrichment analysis was subsequently conducted (Supplementary Table 4). For miR-155-5p (Figure 2J), the enriched KEGG pathways included colorectal cancer, Hepatitis B, pathways in cancer, and so on. Cell cycle, FoxO signaling pathway, and colorectal cancer, and so on were included in the enriched KEGG pathways for miR-18a-5p ( Figure 2K). Hepatitis B, HTLV-I infection, pathways in cancer, and so on were included in the enriched KEGG pathways for miR-93-5p ( Figure 2L). Interestingly, FoxO signaling pathway was highlighted in our analysis. As we know, FoxO signaling pathway widely participated in cell autophagy, apoptosis, and proliferation. Whether FoxO signaling pathway participates in melanoma metastasis needs further study.
Then, we finished constructing the miRNA-hub gene network by utilizing Cytoscape software. Through miR-18a-5p, MYC and JUN were regulated and miR-93-5p, UBC was regulated through miR-18a-5p and miR-93-5p, PTEN and CCND1 were regulated through the three selected miRNAs, which were all shown in Figure 3. As results showed, PTEN and CCND1 might serve as key targets in correlation with melanoma metastasis.

Inhibition of miR-93-5p Could Reduce SK-MEL-28 Cells Invasive Ability Through UBC
According to the results of previous studies, there is no related report of miR-93-5p in melanoma. To validate the results of comprehensive analysis, we picked out miR-93-5p and its potential target gene UBC for further study. Using the expression data of melanoma from ENCORI (http://starbase.sysu.edu.cn/), it was found that miR-93-5p expression was obviously in inverse correlation with UBC in 449 melanoma patients ( Figure 5A). The predicted site in UBC 3 ′ -UTR that can be bound by miR-93-5p is illustrated in Figure 5B. We consequently explored whether UBC served as miR-93-5p's direct target in SK-MEL-28 cells. The luciferase reporter experiment indicated that miR-93-5p mimic obviously inhibited the luciferase activity in SK-MEL-28 cells with wt-UBC-3 ′ UTR vector, but not in mutant UBC-3 ′ UTR ( Figure 5C). Moreover, transfected miR-93-5p mimic gave rise to an obvious reduction of UBC protein ( Figure 5D). To verify whether the reduction of SK-MEL-28 cells' invasive ability by inhibiting miR-93-5p was dependent on UBC, siRNA targeting UBC was transfected into miR-93-5p-upregulating SK-MEL-28 cells. The results showed that the anti-invasive effect by miR-93-5p's downregulation was reduced when the expression of UBC was decreased by its siRNA (Figure 5E).

DISCUSSION
Poor outcome of melanoma mainly results from metastasis. Until now, the mechanism of melanoma metastasis remained unclear. Through previous studies and reports, it was found that most primary melanomas can be cured by local resection (Damato, 2012), but metastatic melanomas have historically had a poor prognosis, with a median survival time of 9 months and a long-term survival rate of 10% (Hall et al., 2000). Metastatic melanoma accounts for approximately 80% of skin cancer-related deaths (Aladowicz et al., 2013). The more understanding of the pathogenesis of metastasis, the better we can develop drugs and treat this disease. Increasing evidence has shown that miRNA expression profiling analysis is a useful tool to study cancer progression and metastasis. Abundant miRNAs took part in numerous pivotal cellular pathways in relation to the progression of cancer, which we have recognized to be expressed aberrantly in melanoma.
Up to now, we have done the research on miRNA associated with metastasis for melanoma by miRNA expression profiling, whose values are relatively limited (Figure 6). As for the research, we gathered melanoma's miRNA expression profiling dataset and systemically analyzed them, aiming at retrieving miRNAs associated with metastasis. A total of 63 DEMs were identified, many of which were known to be involved in melanoma. For example, melanoma cell-secreted exosomal miR-155-5p, which was selected as a survival-associated miRNA, could induce a proangiogenic switch of fibroblasts associated with cancer through improving proangiogenic factors' expression in recipient fibroblasts via SOCS1/JAK2/STAT3 signaling pathway (Zhou et al., 2018). MiR-141-3p, which was regulated negatively in our study, showed the inhibitory effects on melanoma cells' anchorage-independent growth ability, their potential of invasion, and expression of an embryonic-like, multipotent, aggressive cancer phenotype. Further data suggested that vasculogenic mimicry was regulated by miR-141-3p through phosphatidyl inositol-3-kinase (PI3K) (Verrando et al., 2016) and extracellular signal-regulated kinase 1/2 (ERK1/2). MiR-330-3p inhibited melanoma proliferation by targeting TRX2 (Yao et al., 2018). miR-198, regulated by hsa_circ_0025039, inhibited cell growth, glucose metabolism, and invasion in malignant melanoma via targeting CDK4 (Bian et al., 2018).
Furthermore, the building of the PPI network was finished on the basis of the predicted genes. Additionally, we made use of analyses of GO and KEGG pathway enrichment aiming at further interpreting their functions about biology. Then, we identified the top 10 hub nodes whose degrees were higher for each miRNA as hub genes. AKT1, TP53, and UBC had the highest node degrees among these genes. Through miRNA-hub gene network construction, it was found that PTEN and CCND1 were regulated by all three miRNAs. The results indicated that AKT1, TP53, UBC, PTEN, and CCND1 may be pivotal targets in correlation with melanoma metastasis. Previous studies suggested that AKT1 promoted the development of melanoma metastases (Cho et al., 2015), and its melanoma-derived mutation, AKT1 E17K , activated focal adhesion kinase and promoted melanoma brain metastasis (Kircher et al., 2019). In addition, p53, encoding by TP53 gene, was reported as a therapeutic target of melanoma . The other three key genes, UBC, PTEN, and CCND1, were also well studied in melanoma (Roh et al., 2016;Donigan et al., 2017;Uguen et al., 2017;Mu and Sun, 2018;Zhang et al., 2018;Zhu et al., 2018;Chen et al., 2019;Giles et al., 2019). Overall, systematically analyzing miRNA profiling took a step in the investigation of the mechanisms that underlay melanoma's metastasis.
In future studies, we will collect more clinical samples to explore the prognostic value of miR-18a-5p, miR-155-5p, and miR-93-5p in melanoma. In addition, we will further study the mechanism of miR-93-5p/UBC in vivo.
To sum up, through analyzing the miRNA expression profile from an open database, many biological miRNAs were identified, which may participate in melanoma's metastasis. This work verified the function of miR-155-5p, miR-18a-5p, and miR-93-5p in melanoma metastasis, providing additional insights into the complex process of this disease. Meanwhile, we identified miR-93-5p as a pro-invasion miRNA by regulating the expression of UBC in melanoma for the first time. Our research found that miR-18a-5p, miR-155-5p, and miR-93-5p play a key role in the mechanism of melanoma metastasis, and proved that miR-93-5p/UBC is a potential effective target for melanoma.