Bioinformatics Analysis Identifies MicroRNAs and Target Genes Associated with Prognosis in Patients with Melanoma.

BACKGROUND Melanoma of the skin can be associated with early metastases and poor prognosis. This study aimed to identify microRNAs (miRNAs) and target genes associated with prognosis in melanoma using bioinformatics analysis. MATERIAL AND METHODS The Gene Expression Omnibus (GEO) database identified the microarray dataset GSE20994. Differentially expressed miRNAs (DE-miRNAs) were first identified using R language software and validated by GEO2R. Potential target genes of DE-miRNAs were screened, and their targets and prognostic role were evaluated in the miRTarBase database. Pathway enrichment and functional annotation analysis for target genes were established using the DAVID database. miRNA-hub gene networks and protein-protein interaction (PPI) networks were constructed and visualized using the STRING database and Cytoscape. Kaplan-Meier survival curves were constructed using transcriptome and survival data from the UALCAN web tool. RESULTS There were 132 upregulated and 134 down-regulated DE-miRNAs identified from human melanoma samples. From the top three upregulated miRNAs, there were 580 potential predicted target genes, and from the top three down-regulated miRNAs, there 543 potential predicted target genes. miR-300 was upregulated, and miR-629 was down-regulated in melanoma. Two pivotal bub genes, TP53 and GAPDH, were identified in the PPI network. Five out of ten hub genes were modulated by upregulated miR-580, and five by miR-629. Increased mRNA expression of DAPK2 was associated with increased OS, and increased mRNA expression of SKCM, TECPR2, and ZNF781 were associated with reduced OS. CONCLUSIONS Bioinformatics analysis identified miRNAs and target genes associated with melanoma that may represent potential prognostic and therapeutic biomarkers.

Background

Melanoma is a malignant tumor that most commonly arises in the skin, but can develop from melanocytes that originate in the embryonic neural crest at other sites. Worldwide, melanoma accounts for 90% of deaths from skin tumors [1]. In 2018 in the US, there were more than 90,000 reported cases of invasive melanoma, and the number of cases is predicted to increase to 116,000 by the year 2026 [2,3]. Melanoma is an aggressive form of malignancy that is associated with rapid growth, early metastasis, local recurrence, and poor prognosis. Melanoma is currently ranked as the sixth most common malignancy in the US, and has a 5-year survival rate of less than 70%, even with early radical surgery [4]. Therefore, early detection and recognition of melanoma are keys to improve survival. The standard treatment for primary melanoma of the skin is surgical excision [5]. Other treatments, such as immunotherapy, targeted therapy, chemotherapy, and radiation therapy, are applied to treat advanced melanoma. However, the results are not satisfactory, especially for melanomas that metastasize to lymph nodes or major organs. Therefore, it is important to investigate the underlying molecular mechanisms of melanoma and to identify more beneficial early diagnostic techniques and more reliable molecular markers for monitoring recurrence and evaluating prognosis.

MicroRNAs (miRNAs) are small noncoding RNAs between 21–25 nucleotides in length. They are found in diverse organisms and play important roles in regulating the translation and degradation of messenger RNAs through base pairing and binding to the 3′-untranslated region [6]. Expression profiling studies have identified several miRNAs that are involved in the proliferation, apoptosis, and migration of melanoma cells in vitro and the immune response to cutaneous melanoma in vivo [7]. Studies that have profiled the expression of miRNA in melanoma have identified several miRNAs that have a role in melanoma cell proliferation, migration, and invasion in vitro, and miRNAs involved in melanoma cell apoptosis and the immune response to cutaneous melanoma in vivo [8-10]. MiR-211, miR-200c, miR-205, miR-196a, and miR-218 participate in the initiation and progression of melanoma. MiRNAs, including miR-30b/d, miR-34a/b/c, miRNA-155, and miR-182, are activated in the regulation of immune response in melanoma. Invasion and metastasis associated with miRNAs have been widely studied, including Let-7a/b, miR-200a/c, miR-145, and miR-203. In terms of therapeutic options for melanoma, upregulation of miR-21 has been shown to augment chemosensitivity and radiosensitivity. Although previous studies have identified the role of individual miRNAs in melanoma, there have been few previous gene array studies to investigate the role of miRNAs and target genes. Gene expression microarrays provide a snapshot of genome-wide expression in health and disease and have become a method to study the biological behavior of human tumors.

Therefore, this study aimed to identify miRNAs and target genes associated with prognosis in melanoma using bioinformatics analysis. The study included screening differentially expressed miRNAs (DE-miRNAs), analysis of functional and pathway enrichment, and integration of protein-protein interaction (PPI) networks, to identify miRNAs associated with prognosis in melanoma.

Material and Methods

Microarray data in cases of melanoma

Microarray data were obtained from the National Center for Biotechnology Information (NCBI) Gene Expression Omnibus (GEO) database (). The microarray dataset GSE20994 was identified based on the GPL9040 platform (febit Homo sapiens miRBase 13.0), which included samples from 35 patients with melanoma and 22 normal controls, was downloaded and used for this study.

Screening for differentially expressed miRNAs (DE-miRNAs)

The downloaded platform and series of matrix files were converted for use with the R language software and annotation package. Data were normalized through the array function of the limma R package () [11]. The differences between the melanoma patient group and the controls were compared using Student’s t-test. The screening thresholds of DE-miRNAs were set as P <0.05, and a fold-change (FC), defined as the ratio of the difference between the final value and the initial value when divided by the initial value, was identified as >1. Further validated procedures for the DE-miRNAs were performed using the online GEO2R analysis tool from the GEO database.

Prognostic value of target genes of DE-miRNAs in melanoma

The miRTarBase online database was used to identify the prognostic value of target genes. The miRTarBase database contained thousands of miRNA-target interactions (MTIs) collected from filtered research articles and experimentally validated () [12]. In this study, the top ten upregulated and down-regulated DE-miRNAs were identified by miRTarBase, a database of curated experimentally validated miRNA-target interactions, and were analyzed to identify the target genes.

Gene Ontology (GO) annotation and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis of DE-miRNAs

The DAVID database () is an online tool for high-throughput functional analysis of genes. The functional and pathway enrichment of the candidate miRNAs were analyzed and annotated using the DAVID database [13]. GO annotations were performed using the DAVID online tool on the screened DE-miRNAs. The KEGG pathway analysis of DE-miRNAs was also performed by using the DAVID database. DE-miRNAs that were significantly upregulated and down-regulated were determined from microarray melanoma data, with P<0.05 as the threshold for statistical significance.

Integration of protein-protein interaction (PPI) networks

The STRING database () was used to analyze the PPI networks and to predict the genes that targeted the top ten most upregulated and down-regulated DE-miRNAs. The results were visualized using Cytoscape version 3.7.1, a software platform used to visualize complex networks [14]. After miRNA and hub gene networks were established, the top 25 hub genes were acquired and visualized using Cytoscape version 3.7.1.

Kaplan-Meier survival analysis associated with target genes for patients with melanoma

UALCAN is an interactive web resource for analyzing cancer transcriptome data. The target genes of DE-miRNAs were uploaded to UALCAN (). Kaplan-Meier survival curves were calculated for the distinct effect of target gene expression on patient survival. Using two-sided statistical analysis, significance was assumed to be P<0.05.

Results

Identification of differentially expressed microRNAs (DE-miRNAs) and the prognostic value of target genes in melanoma

The downloaded microarray dataset GSE20994 from the Gene Expression Omnibus (GEO) online database included 35 cases of melanoma (GSM524779-GSM524813) and 22 normal controls (GSM524757-GSM524778). Before further analysis, the data were normalized (Figure 1). Data were then analyzed using the unpaired t-test (P<0.05; log2 FC >1). There were 266 DE-miRNAs that were identified from the dataset (Figure 2), and approximately half the miRNAs were upregulated (132 miRNAs). Table 1 lists the top ten most upregulated miRNAs and the most down-regulated miRNAs. Following ranking of the miRNAs according to the degree of fold-change (FC), hsa-miR-300, hsa-miR-519e, and hsa-miR-580 were the three leading upregulated miRNAs, and hsa-miR-629, hsa-miR-601, and hsa-miR-29c were the three leading down-regulated mi RNAs. Following identification by the miRNA and target interaction database, 580 candidate target genes were identified for the three upregulated miRNAs that had prognostic value in melanoma. There were 543 possible targets that matched the three down-regulated miRNAs. Figure 3 shows the volcano plot and the scatter plot of all the DE-miRNAs.

Figure 1

Normalization of the GSE20994 microarray dataset obtained from the Gene Expression Omnibus (GEO) database. (A) Data before standardization. (B) Data after standardization.

Figure 2

Differentially expressed microRNAs (DE-miRNAs) in 35 cases of melanoma and 22 normal controls. (A) Heatmap of DE-miRNAs. (B) Details of the 57 cases.

Table 1

Top ten up-regulated and down-regulated differentially expressed miRNAs between melanoma and normal control.

miRNAs	Log2FC	Average expression	t	P value	adj. P value	B	Regulated
hsa-miR-300	2.7961	3.848519	5.562796	7.13E-07	1.03E-05	5.615772	Up-regulated
hsa-miR-519e	2.6663	4.096406	4.716727	1.56E-05	0.000137	2.636554	Up-regulated
hsa-miR-580	2.548238	3.961057	5.432352	1.16E-06	1.56E-05	5.145503	Up-regulated
hsa-miR-592	2.229587	4.232154	3.828597	0.000319	0.001585	−0.24548	Up-regulated
hsa-miR-520g	2.176182	4.115652	4.291374	6.85E-05	0.000445	1.217952	Up-regulated
hsa-miR-513b	2.13154	3.861382	4.286952	6.95E-05	0.000448	1.203549	Up-regulated
hsa-miR-1243	2.0813	4.89348	4.161777	0.000106	0.000641	0.799017	Up-regulated
hsa-miR-519a	2.07649	3.689535	4.416181	4.46E-05	0.000326	1.627542	Up-regulated
hsa-miR-873	2.05318	4.245304	3.383928	0.001288	0.004893	−1.55677	Up-regulated
hsa-miR-1245	1.974701	4.37269	4.068774	0.000145	0.000814	0.502595	Up-regulated
hsa-miR-629	−3.90545	6.058667	−6.55068	1.68E-08	6.05E-07	9.261758	Down-regulated
hsa-miR-601	−3.76249	3.126457	−6.49409	2.09E-08	6.68E-07	9.049936	Down-regulated
hsa-miR-29c	−3.73207	3.475553	−6.36988	3.36E-08	9.37E-07	8.585987	Down-regulated
hsa-miR-142-3p	−3.53267	4.662733	−5.60801	6.02E-07	8.82E-06	5.779536	Down-regulated
hsa-miR-328	−3.42793	5.918518	−6.60093	1.38E-08	5.20E-07	9.449997	Down-regulated
hsa-miR-942	−3.35367	3.757525	−5.76897	3.29E-07	5.62E-06	6.365488	Down-regulated
hsa-let-7b	−3.32995	3.329909	−5.62132	5.73E-07	8.69E-06	5.827828	Down-regulated
hsa-miR-186	−3.3159	6.216168	−8.04408	5.28E-11	5.29E-09	14.88909	Down-regulated
hsa-miR-483-3p	−3.22229	4.785116	−5.76707	3.32E-07	5.62E-06	6.358526	Down-regulated
hsa-miR-1262	−3.08486	2.046869	−5.65304	5.09E-07	8.00E-06	5.943025	Down-regulated

Figure 3

Volcano plot of the differentially expressed microRNAs (DE-miRNAs). The black points represent genes with no significant difference. The red points represent upregulated genes screened based on the fold-change >1.0 and a corrected P-value of <0.05. The green points represent down-regulated genes screened based on the fold-change >1.0 and a corrected P-value of <0.05.

Enrichment analysis of the KEGG pathway

GO functional annotation, and the KEGG pathway enrichment analysis were performed to analyze the pathways of candidate target genes further. Three GO categories, the biological pathway (BP), the cellular component (CC), and the molecular function (MF) were chosen for functional annotation. For the target genes associated with the top ten upregulated DE-miRNAs, the top ten GO terms included the regulation of gene expression, metabolism, and energy pathway in the BP category, lysosomes, plasma membranes, and exosomes in the CC category, and DNA binding, cytokine activity, and the regulation of translational activity in the MF category (Figures 4A–4C). KEGG pathway enrichment analysis showed that the regulation of DE-miRNAs was localized to three pathways, the Notch, Wnt, and N-cadherin signaling pathways (Figure 4D). Figure 5A–5D show the results of the analysis of the target genes of the leading ten most down-regulated DE-miRNAs.

Figure 4

Gene Ontology (GO) annotation and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis for the target genes of the top ten most upregulated microRNAs (miRNAs). (A) Enriched functional biological process (BP) of the target genes. (B) Enriched cellular component (CC) of the target genes. (C) Enriched molecular function (MF) of the target genes. (D) Enriched KEGG pathways of the target genes.

Figure 5

Gene Ontology (GO) annotation and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis for the target genes of the top ten most down-regulated microRNAs (miRNAs). (A) Enriched functional biological process (BP) of the target genes. (B) Enriched cellular component (CC) of the target genes. (C) Enriched molecular function (MF) of the target genes. (D) Enriched KEGG pathways of the target genes.

Mapping of the protein-protein interaction (PPI) network and the miRNA target network

Physiologically, proteins rarely function alone but function in networks. In this study, PPI networks were identified for potential target genes to identify the top ten most down-regulated and upregulated miRNAs using the STRING database of known and predicted PPIs. Cytoscape software was used to visualize the data, and 25 hub nodes were identified and evaluated by degree, as shown in Table 2. For upregulated miRNAs, TP53, HSPA8, MDM2, NOTCH1, CDKN1A, CREB1, HNRNPA2B1, SMAD2, CDN1B, and RAC1 were the top ten hub genes, and TP53 showed the highest node degree (node degree=83). For down-regulated miRNAs, the leading ten hub genes were GAPDH, VEGFA, PTEN, JUN, SITT1, CDC42, MDM2, CREBBP, MMP2, and FOS and GAPDH showed the highest node degree (node degree=64). Figure 6 shows the networks based on these screened hub genes, constructed using Cytoscape visualization software.

Table 2

Hub genes identified in the PPI networks.

Up-regulated miRNAs		Down-regulated miRNAs
Gene symbol	Degree	Gene symbol	Degree
TP53	83	GAPDH	64
HSPA8	42	VEGFA	63
MDM2	38	PTEN	52
NOTCH1	35	JUN	45
CDKN1A	32	SIRT1	39
CREB1	27	CDC42	39
HNRNPA2B1	25	MDM2	38
SMAD2	24	CREBBP	37
CDKN1B	22	MMP2	36
RAC1	21	FOS	35
RPS6KB1	19	MAPK14	34
SOX9	18	CASP8	29
PSMC1	18	CCNA2	29
PARP1	18	IGF1R	27
NCOA3	18	COL1A2	26
RPS16	18	DICER1	25
PSME3	17	HIST1H2AC	25
NIFK	17	CDK6	25
RPL12	17	COL1A1	24
NOLC1	16	HIST1H2BB	23
CTPS1	16	MCL1	23
PCF11	16	HIST1H2BH	23
WEE1	15	AKT2	23
THY1	15	MYCN	22
GNG12	15	PDGFRB	22

Figure 6

The mapped networks of the top 25 hub genes constructed using Cytoscape software. (A) The mapped networks for upregulated miRNAs. (B) The mapped networks for down-regulated miRNAs.

The miRNA hub gene network was mapped further, as shown in Figure 7A. Five out of ten hub genes (DMRT2, MTCH1, MTERF4, NHS, and MMGT1) could be potentially modulated by hsa-miR-580. Three hub genes and two hub genes could be potentially modulated by hsa-miR-300 and hsa-miR-520g, respectively. Also, hsa-miR-629 could potentially target eight out of ten hub genes (ASCAN22, ZNF781, ZFP1, ABHD2, IRF2BP2, TECPR2, CCDC43, and BNC2). Three hub genes were chosen as the most likely targets modulated by hsa-miR-29c (Figure 7B). Therefore, hsa-miR-580 and hsa-miR-629 were identified as being most likely to be associated with the pathogenesis and development of melanoma.

Figure 7

The mapped microRNA (miRNA) and hub gene networks. (A) The mapped networks for upregulated miRNAs. (B) The mapped networks for down-regulated miRNAs.

Survival analysis

Survival analysis of potential target genes showed that patients with melanoma who had high mRNA expression of DAPK2 had improved overall survival (OS) (Figure 8). However, patients with melanoma who had high mRNA expression for SKCM, TECPR2, and ZNF781 had reduced OS (Figure 9).

Figure 8

Kaplan-Meier survival curve analysis of the prognostic value of hub genes for upregulated microRNAs (miRNAs) in melanoma.

Figure 9

Kaplan-Meier survival curve analysis of the prognostic value of hub genes for down-regulated microRNAs (miRNAs) in melanoma.

Discussion

Melanoma is an aggressive malignancy that is heterogeneous in terms of behavior and prognosis. Melanoma of the skin is often diagnosed at a late stage and has a low response to chemotherapy. Surgery is the primary treatment for early-stage melanoma. Radiotherapy may not improve overall survival (OS) in patients with advanced melanoma. However, there is a need for diagnostic and prognostic biomarkers for melanoma. In the present study, differential expression analysis of micro RNA (miRNA) arrays downloaded from Gene Expression Omnibus (GEO) database identified 266 differentially expressed miRNAs (DE-miRNAs) in melanoma samples compared with normal controls. In melanoma, hsa-miR-300 was mainly upregulated, and hsa-miR-629 was mainly down-regulated.

Previous studies have shown that miR-300 expression was associated with several human malignancies. In gastric cancer, miR-300 expression was found to be significantly upregulated when compared with adjacent normal gastric tissue, and reduced cell proliferation in gastric tumorigenesis [15]. Also, serum levels of miR-300 were significantly increased in patients with osteosarcoma and hepatocellular carcinoma (HCC) when compared with healthy controls, and increased expression of miR-300 promoted epithelial-mesenchymal transition (EMT), cell proliferation, and cell migration in vitro through activation of the bromodomain-containing protein 7 (BRD7) gene, but in HCC it acted through the FAK/PI3K/AKT signaling pathway [16-18]. Similar findings were shown for colorectal and breast cancer, where miR-300 promoted cancer cell proliferation and cell migration by targeting p53 [19,20].

In the present study, hub gene and miRNA networks were constructed to identify the interactive relationships in melanoma. The results showed that six DE-miRNAs could regulate the top 25 hub genes from the top ten upregulated miRNAs in melanoma. Among these miRNAs, miR-629 was the leading targeted hub gene in eight hub genes identified. The results showed that miR-629 might be a negative modulator for the development of melanoma. Previous studies have shown that miR-629 functions as a tumorigenic factor in renal cell carcinoma, pancreatic cancer, and ovarian cancer [21-23]. By regulating TGF-β/Smad signaling, miR-629 has been shown to increase cell migration of renal cell carcinoma cells, [23]. Overexpression of miR-629 increased cell proliferation and migration of pancreatic carcinoma in both in vitro and in vivo by targeting FOXO3 [21]. Also, miR-629 has been shown to promote proliferation and migration of ovarian cancer cells by suppressing the tumor suppressor gene, testis-specific Y-like protein 5 (TSPYL5) [22]. In this study, although miR-629 was predicted to be a tumor suppressor in melanoma, the molecular mechanisms require further study, particularly given the heterogeneity of melanoma in terms of behavior and prognosis.

In this study, the DAVID database, Gene Ontology (GO) annotation, and the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis identified ten miRNAs. Kaplan–Meier survival curve analysis was performed for the hub genes for five down-regulated miRNAs and five upregulated miRNAs associated with prognosis in patients with melanoma. The most significant enrichment of upregulated miRNAs was found in the platelet-derived growth factor (PDGF) receptor signaling network. PDGF receptor signaling promoted tumor growth by several mechanisms, including promoting tumor angiogenesis by vascular endothelial growth factor A (VEGFA), and mitosis in cancer-associated fibroblasts [15,16]. The PDGF signaling pathway has previously been shown to promote angiogenesis in colorectal cancer (CRC) by targeting pericytes and vascular smooth muscle cells, and high levels of PDGF receptors were associated with tumor invasion and metastasis [24]. Clinical studies have shown that increased expression of PDGF was associated with more aggressive biological behavior in breast cancer and prostate cancer, and resulted in poor prognosis [25,26]. Our findings were with the findings from previously published studies that showed that PDGF also promoted more aggressive clinical behavior in melanoma [15].

Following the construction of the protein-protein interaction (PPI) networks in this study, TP53 and GAPDH were mapped as the hub genes with the highest degree of connectively. TP53, also known as tumor protein p53, is a tumor suppressor gene that encodes a DNA-binding transcription factor governing multiple cellular processes, including DNA repair, cell growth, cell metabolism, cell apoptosis, cell senescence, and cell death [27]. TP53 has a role in the pathogenesis of several human cancers. Wang et al. showed that miR-300 regulated the expression of p53 protein by binding to the 3′-UTR of P53 in human colorectal cancer cells to promote cell proliferation and cell migration of colorectal cancer cells [19]. In the present study, TP53 was identified in melanoma by PPI network analysis. However, this finding requires validated in further experimental studies. Glyceraldehyde-3-phosphate dehydrogenase (GAPDH) is a key housekeeping gene and reference gene for quantification of gene expression in human tumors [28]. A previous study showed that GAPDH was significantly expressed in melanoma and was associated with cell invasion, but its mechanism in melanoma remains unclear [29].

Kaplan-Meier survival curves were constructed using transcriptome and survival data from the UALCAN web tool. Although DAPK2 was not in the top 25 hub genes, this was the only gene associated with improved patient survival in melanoma (Figure 8). DAPK2 is regulated by both miR-520g and miR-519e, and is a tumor suppressor gene that has a role in programmed cell death, the regulation of autophagy, and diverse developmental processes [30,31]. Dysregulation of miR-520g expression was shown to be involved in hepatocellular carcinoma (HCC) and colorectal cancer (CRC) [30,31]. Increased expression of miR-520g reduced survival in patients with ovarian cancer and was associated with tumor progression and chemoresistance to platinum-based chemotherapy by down-regulating DAPK2 [32].

In the present study, TECPR2 and ZNF781 were associated with poor prognosis and reduced overall survival (OS) in patients with melanoma, regulated by miR-629. To highlight the heterogeneity of melanoma, in 2017, Zhang et al. used microarray data from a small number of patients with melanoma (N=21) and benign nevus (N=11) and showed that COL17A1 was enriched in melanoma [33]. The miRNA and gene regulatory network showed that miR-375 targeted CCL27, and miR-375 targeted IGF1R [33]. To our knowledge, there have been no previous studies on the association between the two genes, TECPR2 and ZNF781, and miR-629 in melanoma. However, as this study has shown, bioinformatics analysis and the identification of miRNAs and target genes can be identified that are potentially associated with melanoma.

Conclusions

This study aimed to identify microRNAs (miRNAs) and target genes associated with prognosis in melanoma using bioinformatics analysis. In this study, miR-300 was upregulated, and miR-629 was down-regulated in melanoma, two pivotal bub genes, TP53 and GAPDH, were identified in the protein-protein interaction (PPI) network, increased mRNA expression of DAPK2 was associated with increased OS, and increased mRNA expression of SKCM, TECPR2, and ZNF781 were associated with reduced OS. These preliminary findings may stimulate further research on the prognostic role of these miRNAs and target genes in patients with melanoma.