Prognostic significance of deregulated microRNAs in uveal melanomas.

Uveal melanoma (UM) represents the most frequent primary tumor of the eye. Despite the development of new drugs and screening programs, the prognosis of patients with UM remains poor and no effective prognostic biomarkers are yet able to identify high‑risk patients. Therefore, in the present study, microRNA (miRNA or miR) expression data, contained in the TCGA UM (UVM) database, were analyzed in order to identify a set of miRNAs with prognostic significance to be used as biomarkers in clinical practice. Patients were stratified into 2 groups, including tumor stage (high‑grade vs. low‑grade) and status (deceased vs. alive); differential analyses of miRNA expression among these groups were performed. A total of 20 deregulated miRNAs for each group were identified. In total 7 miRNAs were common between the groups. The majority of common miRNAs belonged to the miR‑506‑514 cluster, known to be involved in UM development. The prognostic value of the 20 selected miRNAs related to tumor stage was assessed. The deregulation of 12 miRNAs (6 upregulated and 6 downregulated) was associated with a worse prognosis of patients with UM. Subsequently, miRCancerdb and microRNA Data Integration Portal bioinformatics tools were used to identify a set of genes associated with the 20 miRNAs and to establish their interaction levels. By this approach, 53 different negatively and positively associated genes were identified. Finally, DIANA‑mirPath prediction pathway and Gene Ontology enrichment analyses were performed on the lists of genes previously generated to establish their functional involvement in biological processes and molecular pathways. All the miRNAs and genes were involved in molecular pathways usually altered in cancer, including the mitogen‑activated protein kinase (MAPK) pathway. Overall, the findings of the presents study demonstrated that the miRNAs of the miR‑506‑514 cluster, hsa‑miR‑592 and hsa‑miR‑199a‑5p were the most deregulated miRNAs in patients with high‑grade disease compared to those with low‑grade disease and were strictly related to the overall survival (OS) of the patients. However, further in vitro and translational approaches are required to validate these preliminary findings.

Introduction

Uveal melanoma (UM) is an uncommon type of cancer. However, it remains the most frequent primary tumor of the eye in adults. In the majority of cases, it affects the choroid (1). Several factors have been associated with the development of UM; however, its etiology remains unclear. Among these factors, a light skin color, eye color and chronic ultraviolet exposure represent the main risk factors for the development of this malignancy (2,3). It is known that light exposure plays an important role in UM development; however, UV-associated mutational spectrum for UM is different from that observed for cutaneous melanoma (4,5).

The molecular mechanisms responsible for the carcinogenesis of UM have been linked with oncogenic mutations involved in cell cycle and apoptosis. Accordingly, UMs exhibit an overexpression of the cell-cycle regulatory protein cyclin D (CCND), implicating the involvement of the RAF/mitogen-activated protein kinase kinase (MEK)/extracellular signal-regulated kinase (ERK) pathway. The activation of this pathway causes the phosphorylation and inactivation of the retinoblastoma tumor suppressor gene (6).

Mutations of the genes GNAQ and GNA11, affecting the RAF/MEK/ERK pathway, have been detected in approximately 45 and 30% of UM cases, respectively. Moreover, GNAQ is involved in endothelin signaling, which plays a crucial role for melanocyte survival and expansion (7). The mutation of the BAP1 gene has been also identified in metastatic UM cases (8). Finally, the inactivation of the p53 pathway, activation of the pro-survival phosphoinositide 3-kinase (PI3K)/AKT pathway, and defects in the Bcl-2 pathway are other molecular alterations associated with the phenomenon of apoptosis in UM (9,10).

The diagnosis of UM is based primarily on a clinical examination by biomicroscopy and indirect ophthalmoscopy. Diagnostic procedures rarely include a histological evaluation. Its classification is based according to the American Joint Committee on Cancer (AJCC) (11). However, further diagnostic analyses are mandatory. In some cases, a tumor biopsy is useful for molecular analyses that can provide a prognostic value. Although improvements in the diagnosis and treatment of UM have been made, there is still a need for the identification of additional early markers that can be used to define the diagnosis and prognosis of patients, as well as therapeutic targets. Among such markers, an important role is played by microRNAs (miRNAs or miRs), a class of small non-coding RNAs with an average length of 19–20 nucleotides, able to alter the expression levels of several genes, thus modulating several molecular and signal transduction pathways (12).

In particular, several studies have demonstrated that epigenetic modifications, such as miRNA de-regulation, methylation or microbiota composition, play a key role in the regulatory mechanisms of different cellular processes, as well as in the regulation of carcinogenic processes (13–16).

Over the past years, the development of novel high-throughput technologies for genetics and epigenetics analyses has allowed researchers worldwide to collect a huge amount of molecular data for different types of tumors. To date, the majority of these data are collected within the public database ‘The Cancer Genome Atlas (TCGA)’, a program created by the National Cancer Institute (NCI) and the National Human Genome Research Institute (NHGRI) that has generated comprehensive, multi-dimensional maps of the key genomic changes in 33 types of cancer (17). Furthermore, to manage the enormous amounts of molecular data, each year, several new bioinformatics software are created to allow users to analyze the numerous data contained in TCGA and to perform different correlation studies (18).

Upon these bases, in the present study, different computational approaches were used to analyze the data contained in the TCGA UM (UVM) database in order to identify miRNAs which are able to be used as prognostic biomarkers for the management of uveal melanoma.

Materials and methods

Data source of UM miRNA expression profiles

The UM dataset of miRNA expression profiles was obtained from the UCSC Xena Browser (https://xenabrowser.net/) containing all UM molecular profiling data deposited on The Cancer Genome Atlas data portal (TCGA). As regards UM, the TCGA UVM cohort was selected for this study (https://xenabrowser.net/datapages/). The ‘Phenotype’ and ‘miRNA mature strand expression RNAseq by Illumina Hiseq’ datasets were downloaded in order to analyze the miRNA expression levels according to the clinicopathological data contained in the ‘Phenotype’ dataset and relative to 80 patients with UM. The datasets analyzed did not contain any normal samples.

From the ‘miRNA mature strand expression RNAseq by Illumina Hiseq’, the expression levels of 80 patients with UM were analyzed. In particular, the dataset contained the expression levels of 1,938 different mature miRNAs. However, for some of the 80 patients with UM some miRNA expression levels were missing (NA value); therefore, only miRNAs with the expression data of at least 50% of the patients (n=40) were considered for further analysis in order to avoid the analysis of miRNAs with non-representative expression levels. Therefore, the differential analyses were performed on 795 mature miRNAs.

Differential analyses of miRNA expression levels according to tumor stage and vital status

For the differential analyses, the 80 patients with UM were clustered according to tumor stage (T3-T4 vs. T1-T2 or high-grade vs. low-grade) and according to the vital status (deceased vs. alive). Following patient stratification, the downregulated and upregulated miRNAs were identified by calculating the fold change value obtained through the differential analysis between the high-grade and deceased groups compared to the low-grade and alive groups, respectively. For further analyses, only the miRNAs differentially expressed with a P-value <0.01 were selected. Subsequently, as reported in our previous studies, the top 10 lists of the most significantly downregulated and upregulated miRNAs were performed in order to select the 20 most representative deregulated miRNAs (top 20) (19–21).

Differential analyses were conducted only taking into account the tumor stage and vital status of patients as other clinical data (e.g., the type of treatment or the recurrence-free survival) were absent or fragmentary.

The annotation of the TCGA UM miRNAs was performed using miRBase (version 22, http://www.mirbase.org/) by converting the miRNA IDs ‘MIMAT00’ in ‘hsa-miR-’. The previously obtained top 20 miRNAs relative to the tumor stage and vital status were merged using a bioinformatics tool, Draw Venn Diagrams of the Bioinformatics & Evolutionary Genomics (BEG) (http://bioinformatics.psb.ugent.be/webtools/Venn/), for the comparison of the two top 20. With this approach, only the miRNAs strongly downregulated or upregulated and with a higher prognostic significance were identified.

Kaplan-Meier estimate of overall survival (OS) in patients with downregulated and upregulated tumor stage-related miRNAs

Overall survival was calculated from the date of diagnosis to patient death, or to the end of follow-up, whichever occurred first. The times of follow-up were different from patient to patient up to a maximum follow-up time of 2,499 days. Survival curves were estimated by using the MedCalc® statistical software (version 18.11.3; http://www.medcalc.org). For this analysis and for future ones only the deregulated miRNAs related to tumor stage were considered.

Correlation analysis between miRNAs and gene expression of the TCGA UVM datasets. In order to identify the putative genes and their relative expression levels altered by the tumor stage top 20 deregulated miRNAs, the publicly available miRCancerdb tool was used. miRCancerdb (https://mahshaaban.shinyapps.io/miRCancerdb/) is a free database of miRNA-gene/protein expression correlation in cancer with a web interface based on the data contained in TCGA and TargetScan databases (22). This R package allows the users to establish the correlation between the expression of miRNAs and genes/proteins. In particular, through miRCancerdb, the correlation value of each miRNA with different genes was obtained. The lists of correlated genes to each miRNA were merged using the Draw Venn Diagrams tool in order to identify the shared correlated genes among all miRNAs.

Since miRCancerdb uses miRNA interaction data derived only from TargetScan, in this study, the bioinformatics tool microRNA Data Integration Portal (mirDIP, version 4.1.11.1; Database version 4.1.0.3; September, 2018) (http://ophid.utoronto.ca/mirDIP) (23) was also used. In particular, mirDIP integrates 12 different miRNA prediction datasets from 6 different microRNA prediction databases providing the levels of interaction between a specific miRNA and all its target genes. With this second approach, the levels of interaction between the selected 20 deregulated miRNAs and the predicted correlated genes selected through miRCancerdb were established.

Pathways analysis and Gene Ontology (GO)

To better understand the functional role of the 20 selected miRNAs, a pathway prediction analysis was performed. For this purpose, the bioinformatics tool DIANA-mirPath (version 3) was used (24). All genes and cancer-related molecular pathways altered by the selected miRNAs were identified by using this computational approach. Furthermore, GO enrichment analysis was performed using the tool GO PANTHER version 14.0 (http://pantherdb.org/). GO PANTHER analysis was conducted for the lists of genes obtained from miRCancerdb and DIANA-mirPath analyses.

Statistical analyses

The fold change values of miRNA expression levels were calculated through differential analysis. Student's t-tests were performed to select the differentially expressed miRNAs with a statistical significance as reported in Tables I and II. For the Kaplan-Meier survival analysis, log-rank non-parametric test was used to compare the survival distributions of the patients with UM according to the downregulation or overexpression of selected miRNAs. Data with a P-value ≤0.01 (two-tailed) were considered to indicate statistically significant differences.

Table I.

Differentially expressed miRNAs in high-grade compared to low-grade uveal melanoma.

miRNA	miRNA ID	FC high-grade vs. low-grade	P-value	Number NA
Downregulated miRNAs
hsa-miR-514a-3p	MIMAT0002883	−4.54	5.14E-03	0
hsa-miR-508-3p	MIMAT0002880	−4.38	9.44E-03	0
hsa-miR-509-3-5p	MIMAT0004975	−4.33	2.40E-03	3
hsa-miR-513c-5p	MIMAT0005789	−3.93	3.33E-03	4
hsa-miR-509-5p	MIMAT0004779	−3.77	1.80E-04	30
hsa-miR-513a-5p	MIMAT0002877	−3.03	7.84E-03	11
hsa-miR-507	MIMAT0002879	−2.74	5.22E-03	31
hsa-miR-211-5p	MIMAT0000268	−1.76	1.82E-03	0
hsa-miR-374b-5p	MIMAT0004955	−1.75	1.86E-03	0
hsa-miR-29c-3p	MIMAT0000681	−1.75	4.95E-03	0
Upregulated miRNAs
hsa-let-7b-5p	MIMAT0000063	1.48	5.83E-03	0
hsa-miR-143-5p	MIMAT0004599	1.55	5.79E-03	1
hsa-miR-452-5p	MIMAT0001635	1.57	8.59E-03	0
hsa-miR-224-5p	MIMAT0000281	1.62	6.47E-03	0
hsa-miR-592	MIMAT0003260	1.64	7.09E-03	10
hsa-let-7b-3p	MIMAT0004482	1.65	1.01E-03	0
hsa-miR-210-3p	MIMAT0000267	1.83	2.18E-03	0
hsa-miR-767-5p	MIMAT0003882	1.94	8.82E-03	15
hsa-miR-199a-5p	MIMAT0000231	2.39	5.74E-03	0
hsa-miR-1247-5p	MIMAT0005899	2.55	2.38E-03	7

The miRNAs in common with those obtained in the differential analysis between the deceased and alive patients are presented in bold font. Number NA indicates the number of missing expression values for each miRNA.

Table II.

Differentially expressed miRNAs in the deceased group compared to the alive group.

miRNA	miRNA ID	FC deceased vs. alive	P-value	Number NA
Downregulated miRNAs
hsa-miR-508-3p	MIMAT0002880	−21.96	3.02E-08	0
hsa-miR-509-3p	MIMAT0002881	−19.83	2.13E-08	0
hsa-miR-508-5p	MIMAT0004778	−14.66	3.12E-08	0
hsa-miR-514a-3p	MIMAT0002883	−14.13	8.95E-08	0
hsa-miR-506-3p	MIMAT0002878	−9.12	1.78E-07	1
hsa-miR-509-3-5p	MIMAT0004975	−8.56	7.81E-08	3
hsa-miR-513c-5p	MIMAT0005789	−7.94	1.66E-07	4
hsa-miR-513a-5p	MIMAT0002877	−6.70	4.10E-08	11
hsa-miR-514a-5p	MIMAT0022702	−4.75	1.83E-05	25
hsa-miR-513b-5p	MIMAT0005788	−4.45	1.58E-05	25
Upregulated miRNAs
hsa-miR-592	MIMAT0003260	1.78	1.73E-03	10
hsa-miR-887-3p	MIMAT0004926	1.81	6.29E-03	4
hsa-miR-708-5p	MIMAT0000646	1.85	1.03E-03	0
hsa-miR-142-5p	MIMAT0004951	1.90	3.99E-04	2
hsa-miR-212-3p	MIMAT0000269	1.98	1.66E-05	0
hsa-miR-155-5p	MIMAT0000433	2.08	8.77E-03	1
hsa-miR-199b-3p	MIMAT0000434	2.19	9.08E-03	0
hsa-miR-199a-3p	MIMAT0004563	2.63	3.45E-05	0
hsa-miR-10b-5p	MIMAT0000232	2.63	3.46E-05	0
hsa-miR-199a-5p	MIMAT0000231	3.08	1.18E-05	0

The miRNAs in common with those obtained in the differential analysis between the high-grade and low-grade uveal melanoma groups are presented in bold font. Number NA indicates the number of missing expression values for each miRNA.

Results

Identification of deregulated miRNAs according to tumor stage and vital status

After downloading the data matrixes related to the datasets ‘Phenotype’ and ‘miRNA mature strand expression RNAseq by Illumina Hiseq’, the 80 patients were clustered into the high-grade group (n=44) and low-grade group (n=36) when considering the tumor stage and in the deceased group (n=23) and alive group (n=57) regarding the vital status. After excluding miRNAs with missing expression values in >50% of patients, the number of miRNAs decreased from 1,938 to 795 miRNAs.

The differential analyses performed between the obtained groups, namely the high-grade vs. low-grade and deceased vs. alive groups, allowed for the identification of 20 deregulated miRNAs for each stratification. As shown in Tables I and II, among the 20 deregulated miRNAs according to tumor stage and vital status, 7 were in common between the two lists. In particular, the 5 miRNAs, hsa-miR-514a-3p, hsa-miR-508-3p, hsa-miR-509-3-5p, hsa-miR-513c-5p, hsa-miR-513a-5p, and the 2 miRNAs, hsa-miR-592 and hsa-miR-199a-5p, were respectively downregulated and upregulated with similar expression levels in both stratification groups, suggesting a potential greater involvement of these miRNA in defining patient prognosis.

Involvement of the 20 tumor stage-related altered miRNAs in patient survival

The Kaplan-Meier estimate of the OS of the patients with UM revealed that the alterations of the expression levels of 12 miRNAs out of the 20 computationally selected miRNAs were associated with a worse prognosis. In particular, among the 10 downregulated miRNAs, hsa-miR-211-5p (log-rank test P=0.0301) and hsa-miR-514a-3p, hsa-miR-508-3p, hsa-miR-509-3-5p, hsa-miR-513c-5p, hsa-miR-513a-5p (log-rank test P<0.01) were significantly associated with a negative prognosis of the patients when downregulated (Fig. 1A). Similarly, the overexpressed miRNAs, hsa-let-7b-5p (log-rank test P=0.0110), and hsa-miR-452-5p, hsa-miR-224-5p, hsa-miR-592, hsa-let-7b-3p and hsa-miR-199a-5p (log-rank test P<0.01) were associated to a worse OS, when upregulated (Fig. 1B).

Figure 1.

The overall survival (OS) of patients with uveal melanoma according to miRNA expression. (A) OS of patients with uveal melanoma according to the downregulation and upregulation of the top 10 downregulated related miRNAs related to tumor stage; (B) OS of patients with uveal melanoma according to the downregulation and upregulation of the top 10 upregulated miRNAs related to tumor stage. Only Kaplan-Meier estimates of OS with a log-rank test value of P<0.05 were reported.

Identification of the miRNAs targeted genes and correlation analysis

After defining the prognostic role of the selected miRNAs, the miRNA targeted genes and modulated pathways were identified by different bioinformatics approaches. At first, through the miRCancerdb tool, all the genes related to each single selected miRNA were downloaded, reporting the correlation value relative to each gene. In particular, for each miRNA, a list of genes that varied from 6,209 to 13,906 items was downloaded.

Subsequently, all the lists of genes related to each miRNA were merged. In total, 53 genes in common between all the 20 selected miRNAs were identified using this approach. As shown in Fig. 2, by stratifying miRNAs for expression levels and stratifying the genes for correlation levels, a clear division between genes positively and negatively correlated to the identified miRNAs was obtained (Fig. 2). By considering the miRNAs, the highest correlation levels were observed for the upregulated miRNA, hsa-miR-767-5p; while the highest levels of negative correlation were observed for the upregulated miRNAs, hsa-let-7b-5p, hsa-miR-143-5p, hsa-miR-452-5p, hsa-miR-224-5p, hsa-miR-592 and hsa-let-7b-3p. Taking into account the genes instead, the AK3L1 and SDK2 genes were those with the highest positive correlation value with upregulated miRNAs. By contrast, the RPL35A, EIF4A2 and DTWD1 genes had lower correlation levels with upregulated miRNAs.

Figure 2.

Heatmap of the miRCancerdb correlation analysis. The upper side of the heatmap reports the 20 computationally selected miRNAs. In bold are reported the miRNAs with prognostic significance for the definition of the overall survival of patients with uveal melanoma. The miRNAs in common between tumor stage and vital status stratification are marked with an asterisk. On the left side of the heatmap all the genes shared and correlated with all the 20 miRNAs are listed. The green squares indicate a negative correlation, the red squares a positive correlation.

To further elucidate the interaction between the 20 selected miRNAs and their correlated genes, mirDIP analysis was performed to establish the interaction levels. mirDIP analysis revealed that 6 out of the 53 genes did not interact with the selected miRNAs. Notably, among these 6 genes, AK3L1 exhibited no level of interaction with the 20 miRNAs, although the miRCancerdb analysis revealed a marked positive correlation with upregulated miRNAs (Fig. 3). Overall, mirDIP analysis demonstrated a medium interaction level among genes and miRNAs. However, interesting data were obtained concerning the SPTBN1, RREB1, MARCKS, SHE and SLC6A15 genes that generally exhibited elevated levels of interaction with almost all the miRNAs examined. Finally, among the miRNAs, hsa-miR-143-5p exhibited lower interaction levels, while the miRNAs, hsa-miR-211-5p, hsa-miR-513a-5p and hsa-miR-224-5p, exhibited a high interaction with several correlated genes (Fig. 3).

Figure 3.

mirDIP gene target analysis - interaction between selected miRNAs and the 53 genes identified through miRCancerdb. For each miRNA is reported the level of interaction with the 53 genes positively and negatively correlated. The intensity miRNA-gene interaction is highlighted with a color scale ranging from dark red (very high interaction) to yellow (low interaction).

Pathway prediction analysis and GO enrichment of selected miRNAs and target genes

For the pathway prediction analysis, all the 20 miRNAs were entered into the search bar of DIANA-mirPath. The analysis revealed that for the miRNAs, hsa-miR-508-3p, hsa-miR-507 and hsa-miR-210-3p, there were not modulated pathways and targeted gene according to the reference database TarBase Version 7.0 of DIANA-mirPath. By performing a cumulative analysis of all the molecular pathways altered by the selected miRNAs, it was possible to establish that all the miRNAs, apart from the 3 not reported miRNAs, were able to modulate several cancer pathways. In detail, the miRNAs are able to alter 53 different pathways and 743 univocal genes (Table SI). Among the altered pathways, 25 were directly related to cancer development and progression. Among these 25 pathways, the most affected pathways were the Proteoglycans in cancer (hsa05205), Adherens junction (hsa04520), FoxO signaling pathway (hsa04068), Pathways in cancer (hsa05200) and the PI3K-Akt signaling pathway (hsa04151) modulated by 17 miRNAs. Furthermore, Pathways in cancer (hsa05200), the PI3K-Akt signaling pathway (hsa04151) and Proteoglycans in cancer (hsa05205) were the pathways with the highest number of modulated genes (Table III). Therefore, it is evident that all the selected miRNAs were strictly involved in cancer development, including UM. The analysis of the repetition count of genes altered by the selected miRNAs revealed that the genes most frequently affected in different or repeated molecular pathways were CCND1 (36 counts), MDM2 (30 counts), TP53 (28 counts), IGF1R and MYC (27 counts), CDKN1A, CTNNB1, MAPK1, NRAS (25 counts), CDK6 (21 counts) and BRAF (18 counts). Other key genes in neoplastic processes were PIK3CA, GSK3B, BAX, KRAS, etc., all involved in tumors as discussed below (data not shown).

Table III.

Molecular pathways involved in cancer development modulated by the 20 selected miRNAs.

No.	KEGG pathway	P-value	No. of genes	No. of miRNAs
1	Proteoglycans in cancer (hsa05205)	7.25E-09	90	17
2	FoxO signaling pathway (hsa04068)	4.91E-05	65	17
3	Pathways in cancer (hsa05200)	0.000203	150	17
4	Adherens junction (hsa04520)	1.84E-08	42	17
5	PI3K-Akt signaling pathway (hsa04151)	0.040504	120	17
6	Viral carcinogenesis (hsa05203)	8.36E-11	89	16
7	Chronic myeloid leukemia (hsa05220)	5.00E-05	39	16
8	Glioma (hsa05214)	8.05E-05	32	16
9	TGF-beta signaling pathway (hsa04350)	0.000203	39	16
10	Prostate cancer (hsa05215)	0.000738	44	16
11	Hippo signaling pathway (hsa04390)	3.24E-10	70	15
12	Cell cycle (hsa04110)	6.49E-10	68	15
13	Endometrial cancer (hsa05213)	0.003501	25	15
14	HIF-1 signaling pathway (hsa04066)	0.018482	44	15
15	Bladder cancer (hsa05219)	0.031409	20	15
16	Non-small cell lung cancer (hsa05223)	0.035274	24	15
17	Melanoma (hsa05218)	0.043902	29	15
18	Renal cell carcinoma (hsa05211)	1.46E-06	35	14
19	p53 signaling pathway (hsa04115)	5.19E-05	38	14
20	Central carbon metabolism in cancer (hsa05230)	0.003672	28	14
21	Small cell lung cancer (hsa05222)	7.16E-05	45	13
22	Pancreatic cancer (hsa05212)	0.000512	35	13
23	Thyroid cancer (hsa05216)	0.009723	14	13
24	Transcriptional misregulation in cancer (hsa05202)	0.031409	66	13
25	Colorectal cancer (hsa05210)	3.93E-05	34	12

By using the Draw Venn Diagram, the comparison of these last data with the list of TCGA UVM gene correlated with miRNA expression (miRCancerdb analysis) revealed that only 4 genes were in common between the list of genes obtained from DIANA-mirPath and that of miRCancerdb, that are LTBP1, ADCY3, EIF4A2 and AGAP1.

The final step of the analysis consisted of the GO enrichment analysis by PANTHER

By entering the list of 53 genes, the software was able to recognize only 50 genes. As shown in Fig. 4, the selected genes were grouped according to molecular function, biological process, cellular component, protein class and pathway (Fig. 4). As regards the molecular function, it was observed that the majority of genes were involved in binding (GO:0005488) functions (42.9% of genes) (Fig. 4A). When considering the biological processes, 21.1 and 17.5% of genes were involved in metabolic processes (GO:0008152) and biological regulation (GO:0065007) (Fig. 4B). As regards the cellular component, the majority of genes belonged to the cell (GO:0005623), while the most represented protein classes were hydrolase (PC00121), transferase (PC00220), cell adhesion molecule (PC00069) and transporter (PC00227), in accordance with the cellular component (Fig. 4C and D). Finally, the analysis of pathway did not reveal strong evidence. Indeed, only 14 genes out of the 50 recognized were assigned to a molecular pathway (Fig. 4E). Similar results concerning the molecular function, biological process and cellular component obtained for the GO PANTHER analysis of the 743 DIANA-mirPath-derived genes (Fig. 5A-C). However, as shown in Fig. 5D and E, respectively, the most represented protein classes were nucleic acid binding (PC00171), enzyme modulator (PC00095), transferase (PC00220) and transcription factor (PC00218); the most represented pathways were the Integrin signaling pathway (P00034), the Gonadotropin-releasing hormone receptor pathway (P06664), CCKR signaling map (P06959) and the Wnt signaling pathway (P00057) (Fig. 5D and E).

Figure 4.

Gene Ontology enrichment analysis by PANTHER for the 53 genes identified through miRCancerdb. (A) Distribution of genes according to molecular function; (B) Distribution of genes according to biological process; (C) Distribution of genes according to the type of cellular component; (D) Distribution of genes according to protein class; (E) Distribution of genes according to the analysis of pathway. Beside each category, the percentage of gene frequency was reported. The number of assigned genes may be greater than the number of recognized genes as the same gene can be included in different categories.

Figure 5.

Gene Ontology enrichment analysis by PANTHER for the 743 genes identified through DIANA-mir-Path. (A) Distribution of genes according to molecular function; (B) Distribution of genes according to biological process; (C) Distribution of genes according to the type of cellular component; (D) Distribution of genes according to protein class; (E) Distribution of genes according to the analysis of pathway. Beside each category, the percentage of gene frequency was reported. For each Ontology the 15 most represented categories are displayed. The number of assigned genes may be greater than the number of recognized genes as the same gene can be included in different categories.

Discussion

Over the past decades, the incidence and mortality of malignant melanoma have been continuously increasing despite the development of novel pharmacological treatments and the adoption of screening programs (25,26). One of the main critical issues in the management of malignant melanoma, including that of the uvea, is the lack of specific prognostic biomarkers able to define the tumor aggressiveness and, consequently, the design of personalized therapeutic strategies (27). Several genetic and epigenetic alterations have been proposed as good diagnostic and prognostic markers for both cutaneous melanoma and UM; however, these studies were not conclusive in identifying effective biomarkers for all patients and studies on this matter are still ongoing (28–30). Over the past decade, a growing body of evidence has indicated that miRNAs may represent good diagnostic and prognostic biomarkers for several pathologies, including cancer, as well as for other clinical applications (31–33).

In recent years, with the development of high-throughput technologies and the development of the omics sciences, there has been a marked increase in genomic and molecular data concerning all types of tumors (34). In particular, The Cancer Genome Atlas consortium collected and publicly shared all the genomics, epigenetics, proteomics, and mutational data of 33 different type of cancer useful for researchers all over the world to perform various bioinformatics analyses (17). For the present study, the TCGA UVM database, containing 24 different datasets, was analyzed. To the best of our knowledge, only two recent studies have analyzed the miRNAs and gene expression data contained in the TCGA UVM repository, identifying, respectively, a signature of 9-tumor-specific miRNAs and hub genes related to UM relapse (35,36). In this study, through several integrated bioinformatics approaches, it was possible to identify a set of miRNAs strictly related to UM patient tumor stage and vital status. In particular, the differential analysis of miRNA expression levels revealed that among the top 20 lists of deregulated miRNAs according to tumor stage and vital status, 7 were in common between the two lists, i.e., the miRNAs hsa-miR-514a-3p, hsa-miR-508-3p, hsa-miR-509-3-5p, hsa-miR-513c-5p, hsa-miR-513a-5p, hsa-miR-592 and hsa-miR-199a-5p. These first data are relevant to those reported in literature as other research groups have well characterized the so-called miRNA-506-514 cluster whose alteration is associated with melanocyte transformation and promoting melanoma growth, however contrasting data were generated on this matter (37,38). In this study, the miRNAs contained in the miRNA-506-514 cluster were down-regulated in the high-grade and deceased patients compared to the low-grade and alive patients with UM, showing a discrepancy with the data obtained in the study by Streicher et al (38). However, the data obtained from the TCGA analysis were concordant with that obtained in the study by Stark et al (37), where the downregulation of hsa-miR-514a was associated with a higher proliferation rate, and therefore with a worse prognosis (37). However, what is important is that this cluster is potentially involved in the mechanisms of progression of UM and therefore the miRNAs it contains (miR-506, −507, −508, −509, −510, −513 and miR-514) can play a strong role as biomarkers for the prognosis of UM.

Furthermore, by considering the deregulated tumor stage-related miRNAs, the downregulated miRNAs, hsa-miR-211-5p, hsa-miR-374b-5p and hsa-miR-29c-3p, are also associated with melanoma development and metastasis (37,39,40). As regards the upregulated miRNAs, these are all related to melanoma development, but not to the uveal type (41,42).

After this first analysis, the prognostic value of each miRNA was established. Among the 20 miRNAs, the Kaplan-Meier estimate revealed that 12 miRNAs had a strong prognostic value in defining the OS of patients with UM. Notably, among these 12 prognostic miRNAs that were present, these all belonged to the miR-506-514 cluster.

The third step of the study, performed by using miRCancerdb, allowed the identification of 53 genes correlated and in shared between the 20 tumor stage-related de-regulated miRNAs. The genes, AK3L1 and SDK2, and the genes, RPL35A, EIF4A2 and DTWD1, were revealed as the most positive and negative correlated genes with upregulated miRNAs, respectively. While the SLC46A1 and KCTD17, and SLC25A30 and ZNF836 genes were the most negatively and positively correlated with downregulated miRNAs, respectively. However, the mirDIP analysis revealed that no gene-miRNAs interactions were deposited for the AK3L1 gene. Furthermore, none of the 9 highly correlated genes exhibited high interaction levels with the selected miRNAs, but only medium levels of interaction.

Subsequently, prediction pathway analysis and GO enrichment were performed in order to define the functional role of the 20 selected miRNAs and their 53 correlated genes. As expected, all the miRNAs and genes were involved in molecular pathways commonly altered in cancer. In particular, it was widely demonstrated that UM is the result of different altered pathways, such as the MAPK pathway, Hippo signaling pathway, hypoxia-inducible factor (HIF)-1 signaling pathway, p53 signaling pathway, PI3K-Akt signaling pathway, etc. (43–45). The DIANA-mirPath analysis revealed that all these pathways were strongly altered by the 20 the selected miRNAs. It is clear that the alterations of these pathways were associated with the development of UM, as well as that of other tumors, as demonstrated in other similar studies (19–21,46). In this regard, the DIANA-mirPath analysis revealed that the genes most frequently altered within these pathways are genes already associated with the processes of tumor progression and metastasis, including those related to UM. Among these genes, noteworthy are the genes, TP53, IGF1R, MYC, NRAS, BRAF, PIK3CA and GSK3B, for which several studies have described their involvement in the progression of UM, cutaneous melanoma and other cancers (28,47–50).

Finally, through GO enrichment, it was demonstrated that the 53 genes correlated with the selected miRNAs and the 743 univocal genes derived from DIANA-mirPath analysis were mainly involved in binding and catalytic activity molecular function and in biological regulation and metabolic processes. All these data demonstrated that the data mining of the TCGA UVM dataset can provide useful information in order to identify specific miRNAs whose alterations may be predictive of a worse prognosis for patients with UM.

Overall, through several integrated computational approaches performed on TCGA datasets, this study identified a set of 20 miRNAs strictly associated with the prognosis of patients with UM. Some of these miRNAs, i.e., hsa-miR-514a-3p, hsa-miR-508-3p, hsa-miR-509-3-5p, hsa-miR-513c-5p, hsa-miR-513a-5p, hsa-miR-592 and hsa-miR-199a-5p, were found to more highly deregulated in patients with high-grade compared to those with low-grade disease. Therefore, further in vitro and in vivo evaluations are required to validate these preliminary bioinformatics data in order to propose these miRNAs as useful prognostic biomarkers to be used in clinical practice to address high-risk patients towards more aggressive therapeutic regimens.