Hsa-let-7c-5p, hsa-miR-130b-3p, and hsa-miR-142-3p as Novel miRNA Biomarkers for Melanoma Progression

This study aimed to screen miRNA biomarkers for melanoma progression. Raw melanoma data were downloaded from the Gene Expression Omnibus (GSE34460, GSE35579, GSE18509, and GSE24996) and the Cancer Genome Atlas (TCGA). Then, all differentially expressed miRNAs (DEmiRNAs) between benign vs. primary, metastatic vs. benign, and metastatic vs. primary groups were obtained in the GSE34460 and GSE35579 datasets, and the miRNAs related to disease progression were further screened. Then, the miRNA-gene network was constructed, followed by enrichment, survival, and cluster analyses. Differentially expressed genes (DEGs), tumor-infiltrating immune cells, and tumor mutation burden (TMB) between subtypes were analyzed. miRNAs were verified in the GSE18509 and GSE24996 datasets. A total of 132 and 209 DEmiRNAs were obtained in the GSE34460 and GSE35579 datasets, respectively, and 27 DEmiRNAs related to disease progression were screened. hsa-miR-106b-5p, hsa-miR-27b-3p, and hsa-miR-141-3p had a higher degree and were regulated by numerous genes in the miRNA-gene network. Moreover, four miRNAs were associated with prognosis: hsa-let-7c-5p, hsa-miR-130b-3p, hsa-miR-142-3p, and hsa-miR-509-3p. Furthermore, the bidirectional hierarchical clustering of 27 miRNAs was classified into three subtypes, and TMB and four types of immune cells, including activated dendritic cells, naïve CD4 T cells, M1 macrophages, and plasma cells, showed significant differences among the three subtypes. The expression levels of most miRNAs in the GSE18509 and GSE24996 datasets were consistent with those in the training dataset. These miRNAs, including hsa-let-7c-5p, hsa-miR-130b-3p, and hsa-miR-142-3p, and activated dendritic cells, naïve CD4 T cells, M1 macrophages, and plasma cells may play vital roles in the pathogenesis of melanoma.


Introduction
Melanoma is a highly aggressive malignant tumor originating from neural crest melanocytes [1,2]. Melanoma cells can infiltrate the tissues, lymph, and blood vessels. Melanoma accounts for 3% of all malignant tumors of the skin [3], with its incidence rate increasing rapidly worldwide [4,5]. In 2017, the incidence rate of melanoma in China was 0.9/10 million by age [6]. Melanoma is not sensitive to radiotherapy, chemotherapy, and biological immunotherapy, and its invasiveness and early metastasis make it the deadliest among skin malignant tumors, with only a 6% fiveyear survival rate of stage IV patients [7]. us, it is necessary to screen for reliable and effective biomarkers for the diagnosis and prognosis of melanoma.
MicroRNAs (miRNAs) are gene regulatory factors with excellent functions, ranging in size from 19 to 25 nucleotides [8]. miRNAs participate in the posttranscriptional control of gene expression by binding to the 3′-untranslated region (3′-UTR) of the target mRNA and then inhibits or degrades the target gene [9]. is molecular-level imbalance in biological processes is a potential mechanism for the development of many diseases. Numerous miRNAs have been reportedly involved in melanoma development. Lu et al. found that five miRNAs, including miR-25, miR-204, miR-211, miR-510, and miR-513c, could be used as prognostic biomarkers for patients with cutaneous melanoma cancer [10]. Zhu et al. revealed that miRNA-378a-3p, miRNA-23b-3p, and miRNA-200c-3p are associated with the occurrence of skin cutaneous melanoma [11]. Li et al. suggested that miR-300 and miR-629 are associated with melanoma [12]. However, these studies screened the miRNAs related to melanoma based on a single dataset with a small sample size, and few studies have identified the miRNAs involved in the progression of melanoma. erefore, miRNA biomarkers for melanoma progression should be identified based on multiple datasets.
is study aimed to screen miRNA biomarkers for melanoma progression. Using the data downloaded from the Gene Expression Omnibus (GEO) and the Cancer Genome Atlas (TCGA) databases, miRNAs related to melanoma progression were obtained, the miRNA-gene network was constructed, and prognosis-related miRNAs were screened. Moreover, cluster analysis was conducted based on miRNAs related to melanoma progression, and the differentially expressed genes (DEGs), tumor-infiltrating immune cells, and tumor mutation burden (TMB) between subtypes were analyzed. is study provided a powerful basis for identifying novel therapeutic targets for melanoma patients. A flowchart of this study is shown in Figure 1.

Data
Processing. For the GSE34460, GSE18509, and GSE24996 datasets, the raw data were preprocessed using limma [15] in the R package. e data preprocessing processes included data reading, RMA background adjustment, quantile normalization, and normalization. e expression matrix of probes was obtained, and the probes were annotated to miRNA ID according to the platform annotation information; then, the miRNA ID was converted using miRBase website [16], and the obtained miRNAs were used for the subsequent analysis. For the GSE35579 dataset, the processed probe expression matrix was downloaded, which was log-transformed, and the probes were annotated to miRNA ID according to the platform annotation information. en, the miRNA ID was converted using miRBase website, and the obtained miRNAs were used for the subsequent analysis. For the TCGA-SKCM dataset, samples with an OS. time > 0 were selected for subsequent cluster analysis.

Differential Expression
Analysis. Based on the GSE34460 and GSE35579 datasets, the typical Bayesian test provided by the R package limma package was used to screen the differentially expressed miRNAs (DEmiRNAs) between benign vs. primary, metastatic vs. benign, and metastatic vs. primary, with a cutoff value of P < 0.05 and |logFC| > 0.5. en, all overlapping DEmiRNAs screened between benign vs. primary, metastatic vs. benign, and metastatic vs. primary were obtained in each dataset, and the common DEmiRNAs were obtained from the two datasets. Based on the common DEmiRNAs, the expressed values of the common DEmi-RNAs in each sample of the GSE34460 and GSE35579 datasets were extracted. For the same group of samples, the average of each miRNA was calculated as the average level of the miRNA in the sample group. After comparing the expression levels of DEmiRNAs in the benign, primary, and metastatic groups, the DEmiRNAs were divided into four types: expression continuous upregulation (up-up), expression continuous downregulation (down-down), upregulated first and then downregulated (up-down), and downregulated first and then upregulated (down-up). e miRNAs that overlapped in the same types were considered miRNAs related to disease progression for subsequent analysis.

Prediction of Target Genes of miRNA and Enrichment
Analysis. Based on the six databases, miRWalk, Microt4, miRanda, miRDB, RNA22, and TargetScan, the miRNAgene interaction pairs were screened using the miRWalk 2.0 tool [17]. e miRNA-gene interaction pairs overlapped in the six databases were identified. Moreover, the "Skin Cutaneous Melanoma" was used as the keyword to screen the disease-related genes in the GeneCards database [18]. e disease-related genes were intersected with target genes of miRNA that overlapped in the six databases, and the miRNA-gene interaction pairs were further obtained, followed by visualization using Cytoscape [19]. Additionally, clusterProfiler [20] in the R package was used to perform enrichment analysis of the genes. e P value was corrected using Benjamini and Hochberg (BH), and the adjust P < 0.05 was regarded as a statistically significant difference.

Survival Analysis.
Based on the miRNAs related to disease progression, the expression values of the miRNAs in each sample were extracted from the miRNA expression data in 2 Genetics Research TCGA-SKCM samples. e patients were then allocated to high and low-miRNA expression groups according to the median expression value, and the Kaplan-Meier curve method was used to evaluate the association between miRNAs and prognosis.

Cluster Analysis of miRNAs.
Based on the miRNAs related to disease progression and the expression value of miRNAs in each sample in TCGA-SKCM, bidirectional hierarchical clustering was conducted based on the centered Pearson correlation algorithm using pheatmap [21] to identify different subtypes. e Kaplan-Meier curve method was used to evaluate the association between subtypes and prognosis.

Identification of DEGs between Subtypes.
e gene expression matrix corresponding to each subtype sample was extracted, and the DEGs were identified between the two subtypes using limma in the R package. e P value was corrected using BH, and the adjust P < 0.05 and |logFC| > 0.5 were set as the cutoff values. Moreover, the DAVID tool [22] was used to perform enrichment analysis on the common DEGs with a threshold value of P < 0.05 and count ≥ 2.
2.8. Tumor-Infiltrating Immune Cells. Based on the gene expression level of samples in each subtype, the CIBERSORT algorithm [23] was utilized to evaluate the proportion of immune cells in subtypes, and the differences in the proportion of immune cells in subtypes were analyzed using the Kruskal-Wallis test with a threshold value of P < 0.05.

TMB Analysis.
e Oncoplot function in Maftools [24] was used to visualize the top 20 genes with the highest mutation frequency in each subtype based on the genetic mutation information in each sample. e TMB value was calculated, and the differences in TMB values in subtypes were analyzed using the t-test with a threshold value of P < 0.05.

Validation Analysis.
To verify that the miRNAs were related to melanoma, the GSE18509 and GSE24996 datasets were used as the validation datasets. e expression values of miRNAs in each sample in the GSE18509 and GSE24996 datasets were extracted, and the differences in expression levels between benign and metastatic groups were analyzed using the t-test.

Identification of DEmiRNAs.
Based on the GSE34460 and GSE35579 datasets, DEmiRNAs between benign vs. primary, metastatic vs. benign, and metastatic vs. primary were identified. As given in Table 1, 80, 107, and 24 DEmiRNAs were obtained between benign vs. primary, metastatic vs. benign, and metastatic vs. primary in the GSE34460 dataset, respectively, and 138, 127, and 78 DEmiRNAs were screened between benign vs. primary, metastatic vs. benign, and metastatic vs. primary in the GSE35579 dataset, respectively. us, a total of 132 and 209 overlapping DEmiRNAs were obtained in the GSE34460 and GSE35579 datasets, respectively. en, totally 61 common DEmiRNAs were screened (Figure 2(a)). A total of 27 DEmiRNAs that overlapped in the same types were considered as miRNAs related to disease progression, including 18 and 9 miRNAs in down-down-and up-up types, respectively ( Figure 2

Prediction of Target Genes of miRNA and Enrichment
Analysis. A total of 2330 miRNA-gene interaction pairs containing 26 miRNAs and 1490 genes were screened using miRWalk 2.0. en, the disease-related genes screened in the GeneCards database were intersected with target genes of Genetics Research 3 miRNA, and a total of 604 miRNA-gene interaction pairs were further obtained, including 26 miRNAs and 377 genes. e miRNA-gene network is shown in Figure 3(a). hsa-miR-106b-5p, hsa-miR-27b-3p, and hsa-miR-141-3p had higher degree and were regulated by numerous genes (Table 2). e Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis showed that the target genes of 14 miRNAs were involved in the MAPK and PI3K-Akt signaling pathways, and the top5 enriched pathways are shown in Figure 3(b).

Cluster Analysis of miRNAs.
Bidirectional hierarchical clustering was conducted based on the expression value of the 27 miRNAs in each sample in TCGA-SKCM. As shown in Figure 5(a), the bidirectional hierarchical clustering of the 27 miRNAs was classified into three subtypes (clusters 1, 2, and 3), and most miRNAs were highly expressed in cluster 2; however, miRNAs showed low expression in clusters 1 and 3. e principal component analysis (PCA) of three subtypes is shown in Figure 5(b). Moreover, survival analysis showed that clusters 1 and 2 were associated with poor OS, and cluster 3 correlated with better OS (Figure 5(c)). Moreover, the Kaplan-Meier curve was drawn between the two clusters, and the results showed that clusters 3 and 1 and clusters 3 and 2 had significant differences (all P < 0.05, Figures 5(d) and 5(e)).
3.6. Tumor-Infiltrating Immune Cells. e CIBERSORT algorithm was used to evaluate the proportion of subtypes of immune cells, and the differences in these proportions were analyzed using the Kruskal-Wallis test.
e results showed significant differences in four subtypes of immune cells, including activated dendritic cells (P < 0.05), naïve CD4 T cells (P < 0.05), M1 macrophages (P � 0.02), and plasma cells (P � 0.02). A box diagram of the proportion of 22 subtypes of immune cells is shown in Figure 7.

TMB Analysis.
e waterfall plot of the top 20 genes with the highest mutation frequency in each subtype is shown in Figure 8(a), which revealed that the gene mutation frequency in cluster 2 was low, the gene mutation frequency in cluster 1 was high, and TTN, MUC16, BRAF, and DNAH5 in the three subtypes had high mutation frequencies. In addition, the differences in TMB values in subtypes were analyzed, and the TMB value in cluster 2 was significantly lower than that in the other two subtypes (P < 0.05) (Figure 8(b)).

Validation Analysis.
e GSE18509 and GSE24996 datasets were used as validation datasets to verify that the 27 miRNAs were related to melanoma. As none of the probes in Table 2    Genetics Research the GSE18509 dataset matched hsa-miR-424-3p and hsa-miR-23b-5p, 25 miRNAs were verified. e heatmap showed that 25 miRNAs were well distinguished between samples from the benign and metastatic groups (Figure 8(c)), suggesting that 25 miRNAs were related to melanoma. e differences in expression levels between benign and metastatic groups were evaluated using the t-test, and the box diagram showed that the expression levels of most miRNAs were significantly different, and the regulation trend was consistent with the results of the training dataset (Figure 8(d)). Moreover, the expression of 27 DEmiRNAs was analyzed using the GSE24996 dataset.
Because the two probes in the GSE24996 dataset were not matched, 25 miRNAs were verified. e results showed that the expression of 22 DEmiRNAs was significantly different among the benign, primary, and metastatic groups (Figures 8(e) and 8(f)), suggesting that the above results were reliable.

Genetics Research 7
public databases. hsa-miR-106b-5p, hsa-miR-27b-3p, and hsa-miR-141-3p had a higher degree and were regulated by numerous genes in the miRNA-gene network, and the target genes of miRNAs were involved in the MAPK and PI3K-Akt signaling pathways. Moreover, four miRNAs were associated with prognosis: hsa-let-7c-5p, hsa-miR-130b-3p, hsa-miR-142-3p, and hsa-miR-509-3p. Moreover, the bidirectional hierarchical clustering of the 27 miRNAs was classified into three subtypes, and TMB and four subtypes of immune cells, including activated dendritic cells, naïve CD4 T cells, M1 macrophages, and plasma cells, showed significant differences among the three subtypes. Sixty-one common DEmiRNAs were screened between benign vs. primary, metastatic vs. benign, and metastatic vs.       [33]. us, miRNAs may be involved in the progression of melanoma via the MAPK and PI3K-Akt signaling pathways. ree miRNAs, hsa-let-7c-5p, hsa-miR-130b-3p, and hsa-miR-142-3p, were related to both OS and DFF, and hsa-miR-509-3p was associated with OS. Liu et al. found that hsa-let-7c-5p was associated with melanoma metastasis [34]. Fu et al. suggested that hsa-let-7c-5p restrains breast cancer cell proliferation and facilitates apoptosis by targeting ERCC6 [35]. Murria et al. identified miRNAs associated with melanoma survival, and the results suggested poorer survival in melanomas with hsa-miR-130b-3p overexpression [36]. Tembe et al. indicated that hsa-miR-142-3p is related with the prognosis of patients with metastatic melanoma [37]. Li et al. found that GPC6 is an early biomarker for metastatic progression of melanoma, which can be regulated by hsa-miR-509-3p [38]. Patil et al. revealed that hsa-miR-509-3p inhibits the migration, invasion, and proliferation of osteosarcoma cells [39]. Our results were in line with those reported in previous studies, suggesting that hsa-let-7c-5p, hsa-miR-130b-3p, hsa-miR-142-3p, and hsa-miR-509-3p are associated with the prognosis of patients with melanoma. Tumor-infiltrating immune cells play a vital role in promoting or inhibiting tumor growth [40]. e bidirectional hierarchical clustering of the 27 miRNAs was classified into three subtypes, and four subtypes of immune cells, including activated dendritic cells, naïve CD4 T cells, M1 macrophages, and plasma cells, showed significant differences among the three subtypes. Among the immune cells involved in cancer immunity, properly activated plasmacytoid dendritic cells play a vital role in building a bridge between inherent and adaptive immune responses and directly eliminating cancer cells [41]. Dendritic cells with potent antigen-presenting abilities have long been considered a critical factor in antitumor immunity [42]. Su et al. found that the abundance of naïve CD4 + T cells is closely related to Tregs, suggesting that the prognosis of breast cancer patients is poor [43]. Tumor-associated macrophages could be used as treatment targets in oncology [44,45]. Plasma cells are the main effectors of adaptive immunity and are responsible for producing antibodies to protect the body against pathogens. Numerous studies have reported that plasma cells are related to the pathogenesis of cancer, such as triple-negative breast cancer [46], lung cancer [47], and colorectal cancer [48]. Moreover, TMB could be considered an immunotherapy biomarker for cancer [49][50][51]. In this study, the TMB value in cluster 2 was significantly lower than that in the other two subtypes, suggesting that these miRNAs might play important roles in the pathogenesis of melanoma by involving tumor-infiltrating immune cells and TMB.
However, this study has certain limitations. First, the data analyzed in this study were downloaded from public databases, and external validation is required. Moreover, the sample size should be expanded in future studies. Relevant experiments must be performed to verify the miRNAs identified in our bioinformatics analyses.
Data Availability e data supporting the findings of this study are available from the GEO database (GSE34460, GSE35579, GSE18509, and GSE24996 datasets).

Ethical Approval
e animal experiments in this study were approved by the Ethics Committee of the First Affiliated Hospital of Harbin Medical University.

Disclosure
Xuerui Wu and Yu Wang are the co-first authors.

Conflicts of Interest
e authors declare that there are no conflicts of interest.

Authors' Contributions
XW, LC, and SZ conceptualized and designed the study. XW, CC, YX, and YW acquired data. LC, XW, and SZ analyzed and interpreted data, XW, CC, YX, and YW performed statistical analysis. LC obtained fund and revised the manuscript for important intellectual content. XW and SZ drafted the manuscript. All authors have read and approved the final manuscript. Xuerui Wu and Yu Wang contributed equally to this work.