Comprehensive analysis and identification of key genes and signaling pathways in the occurrence and metastasis of cutaneous melanoma

Introduction

Cutaneous melanoma (CM) is the most dangerous type of skin cancer. It accounts for approximately 232,100 new cases of CM around the world each year, including 55,500 deaths (Schadendorf et al., 2018), and it ranks 15th among the most common cancers in the world (Leonardi et al., 2018). CM is one of the most aggressive and metastatic human cancers, and compared with other cancer types, it can spread from a small primary tumor to multiple sites throughout the body (Bostel et al., 2016). Although primary cutaneous melanoma can be removed and cured through the operation, when a few millimeters thick skin lesion is found, it represents an advanced stage, and there is a high chance of distant visceral metastasis (Wei, 2017). Once metastatic foci are established in distant organs, the 5-year overall survival rate of melanoma patients drops sharply to less than 10% (Braeuer et al., 2014). Therefore, it is urgent to identify the mechanisms that drive the occurrence and metastasis of CM, and to develop effective therapeutic strategies. Studies had classified the somatic mutations and expression profiles of metastatic melanoma (Cancer Genome Atlas Network, 2015), but the mechanisms of CM evolution and metastasis have not been fully understood (Shain et al., 2018). Understanding the gene expression changes during the development of CM will help to develop new biomarkers and therapeutic targets for the diagnosis and treatment of patients.

Gene expression microarray technology can be used to understand the biology associated with cancers, gene mutations and abnormal biological pathways, as well as to predict the diagnosis, treatment, prognosis or metastasis of patients (Tao et al., 2017). The results of microarray technology provide a wealth of information; thus, the data stored in public databases can be reintegrated and bioinformatics analyzed to search for new clues about the pathological mechanisms of cancers through computers rather than laboratories (Wei et al., 2019). In recent years, a large number of studies have predicted the key genes, signaling pathways and protein functions (Le et al., 2019) of many cancers by analyzing the patients’ genetic profiles. For example, the pathogenesis and metastasis mechanism of colorectal cancer (Huang et al., 2018), prostate cancer (Wang et al., 2018), breast cancer (Bertucci et al., 2019) and other cancers had been explored. A number of studies have done bioinformatics analysis of CM. Chen et al. (2020) compared the DEGs between normal skin and melanoma, then used bioinformatics methods to analyze and identify the pathogenesis of CM. Wang et al. (2020) compared the expression of CD38 in the tissues of healthy people and melanoma patients in the TCGA database, and analyzed the occurrence of subtypes and promoter methylation, so as to conclude that CD38 may be a potential biomarker for CM. Some studies obtained differentially expressed non-coding RNAs by analyzing the microRNA and lncRNA of melanoma, which proposed more possibilities for its occurrence and development mechanism (Li et al., 2019b). Meanwhile, some articles have reported the metastasis of CM. For example, Chen et al. (2019) analyzed the gene expression of primary and metastatic melanoma in a database and obtained some candidate genes for metastasis. Wang et al. (2019b) comprehensively analyzed the gene expression of PM and MM in TCGA and constructed a competitive endogenous RNA (ceRNA) network, then proposed a new idea that non-coding RNA and mRNA may act together on the metastasis of melanoma.

However, previous studies have mainly focused on the analysis of PM and N samples. Several studies have also explored the metastasis of CM, but the entire progression of melanoma has not been analyzed from the perspective of occurrence and development, nor has it been compared by combining multiple data sets. In this study, we downloaded three gene expression profiles (GSE7553, GSE15605, and GSE46517) from the GEO, which all included N, PM and MM samples. DEGs among N, PM and MM were determined by gene expression profiling. Subsequently, GO functions, KEGG pathways and PPI network analyses were performed on DEGs. Finally, verification and survival analysis were performed on identified hub genes, which may be potential biomarkers and therapeutic targets in the occurrence and transfer of CM. The flow chart is shown in Fig. 1.

Materials & Methods

Results

Identification of DEGs

R software was used to compare the gene expression of samples from GSE7553, GSE15605 and GSE46517 data sets, and the DEGs of N, PM and MM were obtained in each data sets (Table 1, table s1-9, available at https://doi.org/10.6084/m9.figshare.13019600.v1). Then, the overlaps of 134, 317, and 147 DEGs between PM and N, MM and N, and MM and PM are obtained from the three data sets, which were shown by Venn disgram (Fig. 2). Among them, there were 12 upregulated genes and 122 downregulated genes in the PM compared to N (table s10, available at https://doi.org/10.6084/m9.figshare.13019600.v1), 153 upregulated genes and 164 downregulated genes between MM with N (table s11, available at https://doi.org/10.6084/m9.figshare.13019600.v1), and MM had 79 upregulated genes and 68 downregulated genes compared with PM (table s12, available at https://doi.org/10.6084/m9.figshare.13019600.v1).

Table 1:

The information of GEO dataset and DEGs.

GEO datasets	Platform	Normal skin (N)	Primary melanoma (PM)	Metastatic melanoma (MM)	Differentially expressed genes (DEGs)
					PM and N	MM and N	MM and PM
GSE7553	GPL570	4	16	40	843	1129	929
GSE15605	GPL570	16	46	12	2605	2829	1560
GSE46517	GPL96	8	31	52	668	813	2466

DOI: 10.7717/peerj.10265/table-1

GO functions and KEGG pathways enrichment analyses of DEGs

DAVID was used for GO functions and KEGG pathways enrichment analysis (Figs. 3–5, table s13-15, available at https://doi.org/10.6084/m9.figshare.13019600.v1). GO functions analysis results showed that compared with N samples, upregulated genes of PM were enriched in the collagen catabolic process (BP), while downregulated genes were enriched in transcription from RNA polymerase II promoter (BP), plasma membrane (CC) and structural molecule activity (MF). The upregulated genes between MM and N samples were enriched in negative regulation of neuron apoptotic process (BP), spindle microtubule (CC) and protein kinase activity (MF), and the downregulated genes were enriched in epidermal development (BP), extracellular exosome (CC) and structural molecule activity (MF). In MM and PM samples, the upregulated genes mainly included cell division (BP), cytoplasm (CC), and protein binding (MF), and the downregulated genes mainly included epidermis development (BP), extracellular exosomes (CC), and structural molecule activity (MF).

The KEGG pathways of the overlapped DEGs were analyzed, the upregulated genes between PM and N were significantly enriched in the transcriptional misregulation in cancer, while the downregulated genes were enriched in the metabolic pathways. In MM and N, the upregulated DEGs enriched in the pathway in cancer and the transcriptional misregulation in cancer, and the downregulated DEGs enriched in the arachidonic acid metabolism and steroid biosynthesis. Complement and coagulation cascades was the top enriched term for upregulation genes of MM and PM, while the p53 signaling pathway was the top enriched term for downregulation genes.

PPI network construction and the most meaningful module analysis

The PPI network of DEGs was constructed using the STRING (Fig. 6) and the most important modules were obtained by Cytoscape (Fig. 7). The GO functions and KEGG pathways enrichment analysis showed that the important modules of PM and N were enriched in the cholesterol biosynthetic process, mitochondrion, oxidoreductase activity and metabolic pathway (Fig. 8A, table s16, available at https://doi.org/10.6084/m9.figshare.13019600.v1). The module genes between MM and N were mainly enriched in cell division, extracellular exosomes, protein binding and oocyte meiosis pathway (Fig. 8B, table s17, available at https://doi.org/10.6084/m9.figshare.13019600.v1). In the MM and PM, the module genes were enriched in keratinocyte differentiation, cytoplasm, structural molecule activity, protein binding and mismatch repair (Fig. 8C, table s18, available at https://doi.org/10.6084/m9.figshare.13019600.v1).

Discussion

Recently, many studies had carried out gene expression profiling and bioinformatics analysis on the molecular mechanism of CM occurrence, but the biological mechanism of its development and metastasis were still unclear. In this study, we downloaded three gene expression data sets from GEO and used a comprehensive bioinformatics method to directly compare the gene expression differences among N, PM and MM samples. 134, 317 and 147 DEGs, as well as 14, 18, and 18 hub genes were identified between PM and N, MM and N, MM and PM, respectively. Then, we used the online analysis website to verify the hub genes expression in TCGA samples and performed survival analysis on CM patients.

Through GO functions and KEGG pathways analyses of DEGs, we found that biological processes of upregulated genes mainly concentrated in cell division, spindle microtubule, protein kinase activity and the pathway of transcriptional misregulation in cancer. The downregulated gene mainly occurred in epidermis development, extracellular exosome, structural molecule activity, metabolic pathways and p53 signaling pathway. Studies have shown that the occurrence and metastasis of melanoma need to be realized through the promotion of cell mitosis and the growth of anti-aging and anti-apoptosis (Bennett, 2016). In addition, spindle microtubules can accelerate the proliferation and transfer of cells, and the regulation of metabolic pathways such as protein synthesis and transcriptional disorders can promote cell division (Soltani et al., 2005). The results showed that downregulated genes were associated with skin epidermal development, melanoma cells were produced in the basal layer of the epidermis and hair follicles, and epidermal keratinization could control the homeostasis of melanocytes (Orgaz & Sanz-Moreno, 2013). Studies have found that most exosomes mediate the tumor process in the progression of melanoma (Mannavola et al., 2019). Therefore, previous studies had confirmed our results.

In the 14 hub genes between PM and N, it had been found that matrix metalloproteinases (MMP) can participate in the skin matrix remodeling through degrading and rebuilding the matrix components, and affect the proliferation, survival, vascularization, protease expression and migration of melanoma cells (Napoli et al., 2020). MMP-9 knockdown can reduce the migration and invasion of melanoma cells and inhibit epithelial-mesenchymal transformation (EMT), thus being considered as a promising molecule for the CM treatment (Tian et al., 2019). Meanwhile, bone morphogenetic protein (BMPs) is involved in the regulation of MMPs and is an inevitable factor in the migration and invasion of melanoma cells (Rothhammer, Braig & Bosserhoff, 2008). Fibroblast growth factor receptor 3 (FGFR3) may promote melanoma growth, metastasis, and EMT behavior by influencing the phosphorylation levels of ERK, AKT, and EGFR (Li et al., 2019a). The loss of EphB6 may have deleterious immunological effects in cancer progression, while Hafner et al. (2003) found that its expression decreased gradually in N, PM and MM. These previous studies had suggested that these hub genes may have a potential role in the development of CM.

In MM and PM, we found 18 hub genes through PPI network, all of which were downregulated genes. Several were associated with keratinocyte differentiation and epidermal development, such as loricrin (LOR), involucrin (IVL), filaggrin (FLG), small proline-rich protein1 (SPRR1) keratin (KRT) and plakophilin (PKP1), which may result in loss of epidermal function. Among them, some hub genes had been found to be related to the production and metastasis of melanoma, such as SPRR1 (Koizumi et al., 1996) and PKP1 (Wang et al., 2019a). However, other genes had not been proven to be related to CM, such as LOR dysregulation was considered as an early indicator of potential malignant diseases, including oral submucosal fibrosis and leucoplasis (Nithya et al., 2019). IVL was a specific and sensitive marker of cell differentiation, the expression of IVL in head and neck squamous cell carcinoma patients with or without lymph node metastasis was significantly different (Jin & Qin, 2020). Loss-of-function mutations in FLG can lead to a decrease in epidermal filaggrin and its degradation products, and increase the sensitivity of CM (Thyssen et al., 2018). In MM and PM, the downregulated amplitude of FLG (log₂ FC = −5.404) was less than the amplitude in MM and N (log₂ FC = −8.586), while there was no significant difference between PM and N. It showed that FLG may be related to CM transfer, and this association needs to be verified by subsequent experiments. Studies had shown that KRT5 and KRT14 were involved in HNSCC differentiation and apoptosis as the epithelial proliferative markers (Wang et al., 2014). Meanwhile, immunohistochemical staining of KRT14 and KRT16 in PM and MM were mostly negative, and the positive distribution contributed to the diagnosis of poorly differentiated squamous cell carcinoma (Safadi et al., 2016). CDSN was expressed in hair follicles and keratinized epithelial cells, played an important role in intercellular adhesion, and was related to skin barrier function and epidermal defense pathway (Mondon et al., 2017). Studies had found that mutations in the CDSN gene could cause excessive keratosis of the skin, and lead to peeling skin disease and hypotrichosis simplex of the scalp (Van der Velden et al., 2020), therefore, down-regulation of CDSN may accelerate the development of CM by slowing down epidermal development.

In MM and PM, there were some hub genes associated with EMT. EMT was normally associated with embryogenesis and wound healing, but in tumor cells, it promoted tumor metastasis by enabling cells to leave the epithelium and acquire mesenchymal specificity (Hodorogea et al., 2019). This process increased the aggressiveness of the tumor by the loss of the epithelial phenotype (E-cadherin, desmosin, laminin-1) and the acquisition of the mesenchymal marker (N-cadherin) (Kalluri & Robert, 2009). Hub gene desmoglein1(DSG1) could control the role of keratinocytes, and contribute to the page-like behavior in the development of melanoma (Arnette et al., 2020). Furthermore, desmocollin1(DSC1) and desmocollin3(DSC3) are members of the E-cadherin superfamily, involved in cell–cell adhesion and cell-extracellular matrix interaction. Benign melanocytes expressed high levels of E-cadherin, and during the transition to melanoma cells, E-cadherin was down-regulated and N-cadherin was up-regulated (Jaeger et al., 2007). Studies had shown that desmoplakin (DSP) was a desmosomal protein involved in cell–cell adhesion. Desmosome formation was characteristic of cell differentiation and intercellular adhesion, and the loss of desmosome might accelerate the occurrence and early migration of tumor cells (Walia et al., 2014).

There are also some hub genes between MM and PM had been found to be related to CM, for example, the down-regulation of cystatin A (CSTA) expression has become an important feature to distinguish N, PM and MM (Wallin et al., 2017). Several S100 family genes had been found to be highly expressed in PM, but low level in MM (Xiong, Pan & Li, 2019). In particular, the loss of S100A7 is highly correlated with the metastasis progression score (Bhalla et al., 2019). Lentini et al. (2008) observed that the anti-invasion effect of transglutaminase(TGM) might lead to the post-translational modification of some components of the cell basal membrane, thereby interfering with the metastasis of melanoma cells. The activity of TGM2 had a protective effect on the progression of melanoma in vivo (Facchiano et al., 2013), but no studies had been conducted to prove the relationship between TGM1 and CM, so TGM1 may be a new potential marker.

We found that most of the above hub genes had been reported to be closely related to the generation and metastasis of CM. Moreover, through prognostic analysis, most hub gene expression differences in CM patients were connected with overall survival, which proved the reliability of our study. There are still a few genes that had not been reported or experimentally confirmed to be associated with CM, but some of them are related to the occurrence and development of other cancers. So they might be potential biomarkers of CM, and a large number of experiments are needed to confirm. Most hub genes, such as LOR, IVL, FLG, DSG3, TGM1, KRT16, SPRR1A, KRT14, DSP and CSTA, showed no difference in the expression of PM and N, but significantly decreased in MM and PM, suggesting that these genes might be potential predictors of CM metastasis. The expression of some genes, such as CDSN, DSG1, DSC3, DSC1 and DSP, was downregulated in all three groups, which might be relate with the occurrence and progression of CM. It is worth noting that three genes, SPRR1B, PKP1 and S100A7, were upregulated in PM and N, but downregulated in MM and PM, which were likely to be used as novel markers to distinguish whether CM was metastatic or not. These assumptions need to be tested experimentally.

Compared with previous studies, we used more samples, and compared N, PM and MM in pairs, then took the intersection, so as to make the experimental results more reliable. Besides using GEO data sets, we also used TCGA data for verification, which increased the sample size and accuracy. However, the limitation of this study is the lack of experimental verification of hub genes. Therefore, to understand whether the hub genes are really closely related to the generation and metastasis of CM, a large number of subsequent experiments are needed to explore.

Conclusions

In conclusion, through bioinformatics analysis, we obtained some potential key genes and biological pathways related to the occurrence and metastasis of CM, providing directions for future research. In addition to the genes that have been reported and demonstrated, other identified key genes may be potential prognostic markers and therapeutic targets for CM occurrence and metastasis.