Identification of Keratinocyte Differentiation-Involved Genes for Metastatic Melanoma by Gene Expression Profiles

Background Melanoma is the deadliest type of skin cancer. Until now, its pathological mechanisms, particularly the mechanism of metastasis, remain largely unknown. Our study on the identification of genes in association with metastasis for melanoma provides a novel understanding of melanoma. Methods From the Gene Expression Omnibus (GEO) database, the gene expression microarray datasets GSE46517, GSE7553, and GSE8401 were downloaded. We made use of R aiming at analyzing the differentially expressed genes (DEGs) between metastatic and nonmetastatic melanoma. R was also used in differentially expressed miRNA (DEM) data mining from GSE18509, GSE19387, GSE24996, GSE34460, GSE35579, GSE36236, and GSE54492 datasets referring to Li's study. Based on the DEG and DEM data, we performed functional enrichment analysis through the application of the DAVID database. Furthermore, we constructed the protein-protein interaction (PPI) network and established functional modules by making use of the STRING database. Through making use of Cytoscape, the PPI results were visualized. We predicted the targets of the DEMs through applying TargetScan, miRanda, and PITA databases and identified the overlapping genes between DEGs and predicted targets, followed by the construction of DEM-DEG pair network. The expressions of these keratinocyte differentiation-involved genes in Module 1 were identified based on the data from TCGA. Results 239 DEGs were screened out in all 3 datasets, which were inclusive of 21 positively regulated genes and 218 negatively regulated genes. Based on these 239 DEGs, we finished constructing the PPI network which was formed from 225 nodes and 846 edges. We finished establishing 3 functional modules. And we analyzed 92 overlapping genes and 26 miRNA, including 11 upregulated genes targeted by 11 negatively regulated DEMs and 81 downregulated genes targeted by 15 positively regulated DEMs. As proof of the differential expression of metastasis-associated genes, eleven keratinocyte differentiation-involved genes, including LOR, EVPL, SPRR1A, FLG, SPRR1B, SPRR2B, TGM1, DSP, CSTA, CDSN, and IVL in Module 1, were obviously downregulated in metastatic melanoma tissue in comparison with primary melanoma tissue based on the data from TCGA. Conclusion 239 melanoma metastasis-associated genes and 26 differentially expressed miRNA were identified in our study. The keratinocyte differentiation-involved genes may take part in melanoma metastasis, providing a latent molecular mechanism for this disease.


Introduction
Melanoma is the cells' neoplasm, starting in skin cells, namely, melanocytes [1]. Environmental factors, such as ultraviolet light exposure, are regarded as the main reason causing melanoma [2]. This tumor is mainly in the skin or adjacent to the skin and spreads throughout the body [3], with a dramatic increase in global incidence in the last few decades [4]. Metastasis, the most important characteristic of malignant tumors, is the primary reason for deaths in related to melanoma [5]. What is more, for patients with advanced melanoma, there are still no satisfactory treatments, which involve complicated changes in multiple genes and signaling pathways. Thus, investigating the latent molecular mechanisms of metastasis in melanoma is of great significance.
Recently, microarray technology has a wide range of applications in studying gene alterations in cancer recurrence, metastasis, tumorigenesis, drug resistance and identification of biomarkers for tumor prognosis, and diagnosing and treating the tumor [6][7][8]. By making use of RNA-sequencing analysis on many genes, as reported, RNAs, which contain long noncoding RNAs (lncRNAs), messenger RNAs (mRNAs), and miRNAs, and proteins have a critical influence on progression, melanoma initiation, and recurrence. MicroRNA is a type of noncoding RNA with a length of 18-25 nucleotides, which serves as a posttranscriptional regulator [9]. miRNA plays a role in human diseases by binding to the target mRNA in the 3-untranslated region . Previous reports have shown that miR-182 inhibits the proliferation of malignant melanoma cells through RECK [10]. Androgen receptor (AR) promotes melanoma metastasis by altering miRNA-539-3p/USP13/MITF/AXL signal [11]. In our study, we gained 3 mRNA microarrays and 7 miRNA microarrays aiming at analyzing the DEGs and DEMs between primary melanoma and metastatic melanoma tissue samples. We applied functional enrichment and network analysis. The results demonstrated latent molecular mechanisms on metastatic melanoma.

Screening DEGs and DEMs.
We conducted the DEG and DEM analysis through the application of the limma software package in Bioconductor package (http://www.bioconductor .org/packages/release/bioc/html/limma.html) [16] in R software. In DEG analysis, our cutoff value was p < 0:01 and |fold change ðFCÞ | >1:5, and we made use of FDR < 0:05 in DEM analysis [12]. Unique DEGs and DEMs were selected.
2.3. GO and Pathway Analysis of DEGs. As a predominant bioinformatics initiative, Gene Ontology (GO: http://www .geneontology.org/) [17] contains the most annotations under three headings: molecular function (MF), biological processes (BP), and cellular component (CC), while we performed the Kyoto Encyclopedia of Genes and Genomes (KEGG: http://www.genome.ad.jp/KEGG) [18] pathway enrichment analysis for the investigation of the signaling pathways that were in association with the unique DEGs. We made GO and KEGG pathway analysis through the Database for Annotation Visualization and Integrated Discovery (DAVID: http://www.david.ncifcrf.gov/) for the identification of the DEGs' biological significance [19]. We regarded FDR < 0:01 and gene count > 2 as statistical significance.

Protein-Protein Interaction (PPI) Network Construction.
We first mapped the DEGs to the Search Tool for the Retrieval of Interacting Genes (STRING) (http://www .cytoscape.org/) [20] aiming at assessing functional associations among them. Then, through making use of the Molecular Complex Detection (MCODE), an app of Cytoscape software, we identified the functional modules of the PPI network.
2.6. Validation of Differential Expression of Metastasis-Associated Genes. UALCAN, a comprehensive web resource for analyzing cancer OMICS data (TCGA and MET500), was used to validate the expression of 11 metastasisassociated genes enriched in "keratinocyte differentiation" of Module 1. The result was performed in boxplots, and we regarded p value < 0.05 to have statistical significance.  Computational and Mathematical Methods in Medicine package, we conducted a differential expression analysis. Totally, we identified 1300, 731, and 1829 DEGs to be significantly differentially expressed from GSE46517, GSE7553, and GSE8401, respectively. Finally, we screened out 239 genes in all 3 datasets, which were inclusive of 21 positively regulated genes and 218 negatively regulated genes in metastatic melanoma tissues in comparison with primary melanoma tissues ( Figure 1). Meanwhile, we identified 63 DEMs from the seven miRNA expression profiles, which were made up of 35 positively negatively regulated miRNAs and 28 negatively regulated miRNAs in metastatic melanoma tissues in comparison with primary melanoma tissues (Table S1).

Functional Enrichment Analysis.
We performed the functional enrichment analysis by making use of downregulated DEGs because the number of the DEGs (218/239) was large. We performed three categories of GO functional     Figure 5: DEM-DEG pair network. Green nodes stand for negatively regulated genes or negatively regulated miRNA in melanoma metastasis tissue. Red nodes stand for positively regulated genes or positively regulated genes miRNA in melanoma metastasis tissue. 5 Computational and Mathematical Methods in Medicine exosome and desmosome (Figure 2(b), Table S2). The GO MF analysis results presented that negatively regulated DEGs were obviously abundant in structural molecule activity and structural constituent of cytoskeleton (Figure 2(c), Table S2). We only identified 2 KEGG pathways, considerably enhancing in the amoebiasis signaling pathway and arrhythmogenic right ventricular cardiomyopathy (ARVC) pathway (Table S2).

Construction and Analysis of PPI Network.
We performed the PPI network analysis of 239 DEGs through the STRING database, and 225 protein interactions with combined scores > 0:7 were identified (Figure 3). We identified three modules by making use of MCODE plugin, an app in Cytoscape. Module 1 was constructed with 21 nodes and 209 edges including DSC1, DSC2, and DSC3. Module 2 was made up of 9 nodes and 32 edges, which were inclusive of KRT2, KRT6B, and KRT8. Module 3 includes 7 nodes and 12 edges, which were made up of KRT1, KRT5, KRT6A, and so on. Furthermore, we performed DEGs' functional enrichment analysis in these 3 modules. The results demonstrated that Module 1 was enriched in the keratinocyte differentiation, extracellular exosome (Figure 4(a), Table S3), and so on in GO analysis. Module 2 was enriched in epidermis development, extracellular exosome (Figure 4(b), Table S3), and so on in GO analysis. Module 3 was enriched in epidermis development in GO analysis (Table S3). Unfortunately, we did not identify significant KEGG pathways in all three modules (FDR < 0:01).

Prediction of DEM Target
Gene and Construction of DEM-DEG Network. We predicted 816 genes in total would be regulated by the 63 miRNAs identified above through making use of the TargetScan, PITA, and miRanda database. Among the 816 genes, 92 genes were overlapped with DEGs, being made up of 11 upregulated genes targeted by 11 negatively regulated DEMs and 81 genes downregulated genes targeted by 15 positively regulated DEMs. We presented the results in detail in Table S4.
Subsequently, as presented in Figure 5, the DEM-DEG network was constructed by Cytoscape software. We found that hsa-miR-181b [24], hsa-miR-211 [25,26], and hsa-miR-24 [27,28], which had the most potential target genes, had a paramount influence on melanoma progression and metastasis. Table S3, Module 1 was obviously enriched in the keratinocyte differentiation in GO analysis. Until now, there were no detailed studies on keratinocyte differentiation in melanoma metastasis. As a result, we evaluated the expression of 11 enriched genes in keratinocyte differentiation through making use of the UALCAN database, including LOR, EVPL, SPRR1A, FLG, SPRR1B, SPRR2B, TGM1, DSP, CSTA, CDSN, and IVL ( Figure 6). The analytic results demonstrated that all of the 11 genes were significantly downregulated in metastatic melanoma tissue in comparison to primary melanoma tissue. Some of the genes, such as LOR, SPRR1B, and IVL, were even undetectable in metastatic melanoma tissue.

Discussion
Poor outcome of melanoma is primarily caused by metastasis. The reason for melanoma metastasis is still poorly understood to date. The more we know about the pathogenesis of metastasis, the better the drug can be developed and this disease can be treated. More and more evidence demonstrated that gene expression profiling analysis is of great use for humans to study the cancer progression and metastasis. We have also identified an abundance of genes involved in many crucial cellular pathways of cancer progression and found their expression in melanoma was aberrant.
The value of studies on metastasis-associated genes for melanoma by gene expression profiling is relatively limited to date. In this research, we collected gene expression profiling datasets of melanoma and made a systemic meta-analysis to retrieve genes in association with metastasis. We identified 239 DEGs in total. Many of them were known to take part in various cancer types. For example, DIO2 was underexpressed in nearly all papillary thyroid carcinomas and known as a latent target for thyroid tumor treatment [29,30]. FGFBP1 was identified as a prometastasis gene in HCC [31], whereas S100 family genes' expressions, particularly S100A7, were high in primary melanoma samples but low in metastatic melanoma. S100A7's expression was in relation to tumor invasion, and it may enhance melanoma's early diagnosis [32][33][34]. Additionally, it has been suggested that many predicted targets of DEGs were associated with cancer development. Interestingly, we found that most of the DEGs were downregulated (218/239). We identified DEGs' expression in TCGA database, and the results were the same as our analysis. But the reason is unknown.  Computational and Mathematical Methods in Medicine Furthermore, we finished constructing the PPI network based on the DEGs. What is more, GO and KEGG pathway enrichment analyses were to further interpret their biological functions. For example, we found that 11 DEGs were enriched in keratinocyte differentiation GO term and close relation to melanoma metastasis (Figure 7). We found LOR, EVPL, SPRR1A, FLG, DSP, and CSTA had been reported in previous studies on melanoma [14,27,[35][36][37][38][39]. But the report of SPRR1B, SPRR2B, TGM1, CDSN, and IVL in melanoma was a gap. Overall, this systematic metaanalysis of gene profiling was a step forward in the investigation of the mechanisms underlying melanoma's metastasis, and the role of SPRR1B, SPRR2B, TGM1, CDSN, and IVL in keratinocyte differentiation GO term needs to be further studied.
This study has some limitations. First, the expression levels of genes related to keratinocyte differentiation were not verified by qRT-PCR. Secondly, the role of miRNA-mRNA in melanoma metastasis needs further experimental verification. In future research, we will carry out these explorations.
To sum up, we finished the identification of many biological genes, which may take part in melanoma's metastasis, by collecting the gene expression datasets from the open database. This work identified keratinocyte differentiationinvolved genes' role in melanoma metastasis for the first time, offering additional insights into this disease's complicated process.

Data Availability
All data analyzed during this study are obtained from published article or are available from the corresponding author on reasonable request.

Consent
Consent is not applicable.

Conflicts of Interest
The authors declare that they have no competing interests.

Authors' Contributions
NX designed the study. KL, SG, ST, and QS performed the study and drafted the article. KL, SJ, BQ, and YS conducted data acquisition, data analysis, and interpretation. All authors discussed the results and agreed to be accountable for all aspects of the work. All authors read and approved the final manuscript.