Identication of Metastasis-Associated Gene And Its Correlation With Immune Inltrates For Skin Cutaneous Melanoma

Background: Skin cutaneous melanoma is a malignant and highly metastatic skin tumor. As the most common cause of death in skin cancer, its morbidity and mortality are still rising worldwide. However, the molecular mechanisms of melanoma metastasis are unclear. Methods: Three Gene Expression Omnibus (GEO) datasets (GSE15605, GSE7553 and GSE8401) were downloaded to identify the differentially expressed genes (DEGs) between primary and metastatic melanoma samples. Gene Ontology (GO) analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment were performed to explore the functional of DEGs by Metascape. The protein-protein interaction (PPI) network was constructed using STRING tool and Cytoscape software. We used the cytoHubba plugin of Cytoscape to identify the most signicant hub genes by four topological analyses (Degree, MCC, DMNC, and MNC). Hub genes expression was validated using UALCAN website. Finally, we explored the association between metastasis-associated genes and immune inltrates through Tumor Immune Estimation Resource (TIMER) database. Results: In total, we obtained 196 DEGs including 12 upregulated and 184 downregulated genes. GO and KEGG enrichment results indicated that DEGs were mainly concentrated in epidermis development, cornied envelope, structural molecule activity, and p53 signaling pathway. Eight hub genes were identied to be closely related to melanoma metastasis, including SPRR1B, DSC1, PKP1, TGM1, DSG1, IVL, SPRR1A and DSC3. On the ULCAN website, all hub genes expression levels are lower in metastatic tissues than in primary cancers. Results from TIMER database revealed that DSC1 and TGM1 were signicantly related with most of immune cell inltration. Conclusions: SPRR1B, DSC1, PKP1, TGM1, DSG1, IVL, SPRR1A and DSC3 may be hub genes involved in the progression of melanoma metastasis and thus may be regarded as therapeutic targets in the future. of metastatic melanoma by regulating the tumor inltration of immune cells.


Introduction
Skin cutaneous melanoma (SKCM) is an aggressive and highly metastatic skin tumor, with a 5-year survival rate less than 15% in patients with advanced melanoma, which seriously threatens human health [1,2] . There has been an increasing trend towards the SKCM incidence and death. It is estimated that there will be approximately 106,110 newly diagnosed melanoma patients and 7,180 new deaths worldwide in 2021 [3]. When melanoma is diagnosed at an early stage, surgical resection is the most effective treatment. However, patients suffering from metastatic melanoma have a poor prognosis. The 5-year survival rate for localized melanoma is 99%, but if distant metastasis occurs, it is only 20% [4].
Therefore, it is necessary to identify new diagnostic biomarker of melanoma metastasis and adopt effective treatment strategies.
In recent years, with the in-depth study of tumor microenvironment, immunotherapy has gradually become an active and important area of cancer treatment. Melanoma is one of the most sensitive tumors to immune regulation. The application of immune checkpoint inhibitors, such as monoclonal antibodies against cytotoxic T-lymphocyte-associated protein (CTLA-4) and programmed cell death protein 1 (PD-1), has been found to improve outcomes in patients with advanced melanoma [5]. A recent study showed that 5-year overall survival rate for patients with treatment-naïve, recurrent, or unresectable stage III/IV malignant melanoma treated with niluzumab was 26.1% [6]. Therefore, it is crucial to identify new possible immunotherapeutic biomarkers for this disease to improve prognosis.
In this study, we aimed to understand the differences in gene expression between primary and metastatic melanoma and determine the potential prognostic value of hub genes and their association with immune in ltration. We compared the gene expression pro les of the primary and metastatic samples of SKCM patients from the Gene Expression Omnibus (GEO) database. Then, the differentially expressed genes (DEGs) were underwent Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis to explore the biological functions. We also established a protein-protein interaction (PPI) network to get a better understanding of the mutual interactions of DEGs. We further validated the expression levels of hub genes and evaluated the relationship between metastasisassociated genes and immune in ltration using Tumor Immune Estimation Resource (TIMER). In summary, a series of bioinformatics analysis was used to seek novel predictive biomarkers in metastatic melanoma and immunotherapeutic targets.

Data collection and processing
The gene expression pro les of SKCM were retrieved from the GEO database (http://www.ncbi.nlm.nih.gov/geo) using the keywords "melanoma" OR "skin cutaneous melanoma" OR "SKCM". "Homo sapiens" and "Expression pro ling by array"were searched next round. The limitation criteria included: (1) human SKCM tissue samples; (2) inclusion of primary and metastatic melanoma tissue samples; (3) studies included at least 30 samples. After a systematic review, three gene expression datasets (GSE15605 [7], GSE7553 [8] and GSE8401 [9]) were downloaded. The mRNA expression pro le of GSE15605 and GSE7553 were detected by using GPL570 platform (Affymetrix Human Genome U133 Plus 2.0 Array). GSE8401 was based on GPL96 platform (Affymetrix Human Genome U133A Array).
Among them, GSE15605 consisted of 46 primary melanoma tissues and 12 metastatic melanoma tissues, GSE7553 contained 14 primary tissues and 40 metastatic tissues, and GSE8401 included 31 primary samples and 52 metastatic samples. GEO2R is an interactive web tool (http://www.ncbi.nlm.nih.gov/geo/geo2r) using 'limma' package of R to analyze the gene expression data of the above-mentioned microarrays and nd the differential genes between primary and metastatic melanoma samples. adj. P < 0.05 and |log2FC| > 1 were set as DEGs cutoff criterion. The overlapping portions between upregulated and downregulated genes of three datasets were subsequently examined with venn software (http://bioinformatics.psb.ugent.be/webtools/Venn/).

GO and KEGG pathway analysis
To reveal functions and interactions of DEGs, the Metascape platform (https://metascape.org/gp/index.html#/main/step1) was used to perform GO analysis and KEGG pathway enrichment analysis. P value < 0.01, a min overlap genes = 3 and min enrichment factor > 1.5 were set as the cut-off criteria for identifying terms and pathways.

PPI network construction and hub genes selection and analyses
The Search Tool for the Retrieval of Interacting Genes (STRING) [10] (http://string-db.org) (version 10.0) was used to construct the PPI network of 196 common DEGs. Those with a score ≥ 0.4 as the cut-off criterion [10]. Then, the string interactions were imported into Cytoscape [11] (http://www.cytoscape.org) (version 3.6.1) for visualization. Speci cally, we applied the Cytoscape plugin CytoHubba [12] to identify hub genes via four models: maximal clique centrality (MCC), maximum neighborhood component (MNC), Degree, and density of maximum neighborhood component (DMNC). GeneMANIA website (http://www.genemania.org/) [13] was used to analyze coexpression gene lists and predicting the function of hub genes. Finally, the cBioPortal website (https://www.cbioportal.org/) was applied to query the variation of hub gene in SKCM-Metastasis.

Validation of hub gene expression
We analyzed the mRNA expression of hub genes using TCGA data, which was obtained from the UALCAN website (http://ualcan.path.uab.edu/analysis.html) [14]. A total of 473 samples including 1 normal tissue, 104 primary melanomas and 368 metastatic melanomas were contained. P < 0.05 were statistically signi cant.

Immune in ltration analysis
The relationship between metastasis-associated genes expression and the abundance of immune in ltration in SKCM-Metastasis was evaluated using TIMER [15]. TIMER included six immune cells: B cells, CD4+ T cells, CD8+ T cells, neutrophils, macrophages, and dendritic cells. P < 0.05 was identi ed to be signi cant. Gene expression levels were visualized with log 2 TPM.

Screening for DEGs
Three gene expression pro les (GSE15605, GSE8401, and GSE7553) were downloaded from GEO database. DEGs between primary and metastatic samples of SKCM were identi ed using GEO2R. Based on the criteria of adj. P value < 0.05 and |log2FC > 1|, a total of 1420 DEGs (596 upregulated and 824 downregulated) were identi ed in GSE15605, 1001 DEGs including 447 upregulated and 554 downregulated were screen out in GSE8401, and 825 DEGs including 143 upregulated and 682 downregulated were identi ed in GSE7553. The volcano plots of pro le GSE8401 and GSE15605 data were shown in Fig. 1AB. Finally, a total of 196 overlapping DEGs were signi cantly differentially expressed among three datasets, which included 12 upregulated genes and 184 downregulated genes (Fig. 1CD). These genes were treated as candidate DEGs and used for further analysis.

Functional enrichment analysis for DEGs
To further explore biological functions of these 196 DEGs, GO and KEGG pathways analysis were performed by Metascape ( Fig. 2A). In BP, DEGs were mainly concentrated on epidermis development, peptide crossing-linking, multicellular organismal water homeostasis and intermediate lament cytoskeleton organization. In terms of CC, DEGs were primarily enriched in corni ed envelope, intermediate lament anchoring junction, extracellular matrix. For MF analysis, DEGs were mainly enriched in structural molecule activity, structural constituent of epidermis, calcium ion binding, serinetype peptidase activity. The results of KEGG pathway analysis revealed that DEGs were dominated in amoebiasis, estrogen signaling pathway, phenylalanine metabolism and p53 signaling pathway etc. (Fig.  2C). Additionally, GO and KEGG enriched terms were closely linked and clustered into intact networks (Fig.   2B).

PPI network analysis and hub Gene screening
The PPI network was constructed using STRING database. The interacting genes were then introduced into Cytoscape for network visualization analysis, which contained 161 nodes and 806 edges. We ranked the top 20 genes of the whole network based on the Cytoscape plugin CytoHubba models: Degree MCC, DMNC, and MNC( Fig. 3A-D). Venn analysis was performed to get the intersection of these genes. Notably, the top 20 genes from topological analysis algorithms included eight common genes: SPRR1B, DSC1, PKP1, TGM1, DSG1, IVL, SPRR1A, DSC3, which may play important roles in the development of metastatic melanoma.

Analysis of hub Genes
The network and function of hub genes and their co-expression genes were analyzed by GeneMANIA platform. Eight hub genes showed the complex PPI network with the co-expression of 35.29%, predicted of 24.15%, physical interactions of 16.64%, co-localization of 16.23%, shared protein domains of 7.34%, and Genetic Interactions of 0.35% (Fig. 4A). The network and function of hub genes were enriched in skin development, peptide cross-linking, keratinocyte differentiation, epidermis development and cell-cell junction. Genetic alterations of the eight genes in metastatic melanoma were showed in Fig. 4B, missense mutation was the most common type of mutation in the downregulated genes.

Validation of hub genes expression
To further verify previous de ned hub genes, we obtained expression data of 473 TCGA-SKCM samples which including 1 normal tissue, 104 primary melanomas and 368 metastatic tissue samples from the UALCAN database. The results showed that all the eight hub genes were signi cantly downregulated in metastatic melanoma compared with primary melanoma. All results had a P value < 0.001, which was statistically signi cant (Fig. 5).

We analyzed the correlations between hub genes expression and immune in ltration levels of six immune cells (CD8+ T cells, CD4+ T cells, B cells, neutrophils, macrophages and dendritic cells) in SKCM-
Metastasis microenvironment using TIMER. The purity-corrected partial spearman correlation and P value were shown in each plot. P 0.05 was considered as signi cance. We found that DSC1 expression was negatively correlated with in ltration degree of B cells, CD4+ T cells (P <0.05). The expression of TGM1 identi ed negatively correlated to the in ltration level of four immune cells, including CD8 + T cells, CD4+ T cells, macrophage cells, neutrophil cells and dendritic cells (P < 0.05) (Fig. 6). The expression of the remaining six genes was not associated with immune in ltration levels of six immune cells ( Supplementary Fig. S1).

Discussion
Metastasis is one of the most common causes of death in SKCM. Recently, several studies have focused on predicting the occurrence of distant seeding in melanoma [16,17], indicating that early detection of metastasis is an effective way to improve clinical outcomes. In the present study, we identi ed that 196 DEGs, including 12 upregulated genes and 184 downregulated genes, are closely related to melanoma metastasis. GO analysis showed that DEGs mainly enriched in epidermis development, peptide crossinglinking, multicellular organismal water homeostasis, regulation of epidermis development, and intermediate lament cytoskeleton organization. The above biological processes are mainly related to cell proliferation and differentiation, we speculate that DEGs may promote melanoma metastasis by regulating the proliferation and division of cancer cells. Moreover, KEGG enrichment analysis revealed that DEGs mainly enriched estrogen signaling pathway, phenylalanine metabolism and p53 signaling pathway. Studies have found wild type p53 is expressed at high levels in melanoma and contributes to survival of melanoma cells [18]. In melanoma cell line, through activation of p53 pathway signaling, Spetasin can induce apoptosis and inhibit cell migration [19]. Theses nding had con rmed our research results.
Fritsche and Knopf (2017) illustrated that the role of the tumor suppressor p53 and its pathway of p53 in mucosal melanoma.
Among eight identi ed hub genes, SPRR1A, SPRR1B and PKP1 have been reported to be associated with melanoma metastasis [8,20,21]. SPRR1(SPRR1A, SPRR1B), belongs to the SPRR multigene family, is regarded as a main precursor component of the keratinocyte corni ed envelope, and SPRR1 plays essential roles in the processes of epidermal development, keratinocyte differentiation, cell-to-cell signaling and cell adhesion [22,23]. SPRR1A and SPRR1B expression were both higher in primary melanoma than in metastatic melanoma [8,24]. PKP1 is a member of the armadillo protein family. Lee. P et al. reported that RIPK-PKP1 signaling pathway as novel axis involved in skin strati cation and tumorigenesis [20]. Another study validated that PKP1 was related to the metastasis of melanoma and could even distinguish low-and high-grade of metastatic melanomas [25]. In addition, through pairwise comparisons of normal skin, primary and metastatic cutaneous melanoma, PKP1 may also serve as novel markers for discriminating whether cutaneous melanoma metastasizes [26].
IVL is known as a main component of keratinized envelope, which has been proved to be an early and speci c marker of differentiation of epidermal keratinocytes [27][28][29]. In oral squamous cell carcinoma, the expression patterns of IVL were different between carcinoma in situ and invasive carcinoma, implying IVL expression was closely related to the degree of differentiation of the tumors [30]. IVL expression presented remarkable different in head and neck squamous cell carcinoma samples with or without clinical lymph node metastasis [31]. However, there are few articles about the relationship between IVL and melanoma metastasis. Transglutaminase 1 (TGM1) is a major subtype of transglutaminase (TGM) genes which is expressed in the epidermis and strati ed squamous epithelium [32]. It has been reported that TGM1 expression levels were correlated with gastric cancer patient survival, and TGM1 promotes the stemness and chemoresistance of gastric cancer cells by regulating Wnt/β-catenin signaling [33]. The role of TGM1 in melanoma metastasis has not been reported previously The desmosomal cadherins, DSG1, DSC1 and DSC3, are transmembrane molecules that mediate the adhesion between apposing cells through their extracellular domains [34,35]. DSG1 involves in regulating the role of keratinocyte and melanocyte paracrine signaling, and contributes to the occurrence of early melanoma [36]. DSC1 and DSC3 are members of the E-cadherin superfamily and involve in intercellular adhesion processes and cell-extracellular matrix interaction. The lower expression of DSC1 protein is signi cantly associated with poor tumor differentiation in lung cancer and low expression of DSC1 and DSC3 was signi cantly correlated with poor clinical outcome in human lung cancer [37]. Moreover, reduced expression of DSC1 and DSC3 were observed in colorectal cancer liver metastases and low expression of DSC3 was associated with an increased risk of developing tumor recurrence after resection of colorectal liver metastases [38], suggesting that the desmosomal cadherins protein may be involved in tumor cells invasion and distant metastasis to surrounding tissues.
Tumor microenvironment (TME) is composed of immune cells, mesenchymal cells, endothelial cells, in ammatory mediators and extracellular matrix (ECM) molecules [39,40]. The composition and abundance of immune cells in the tumor microenvironment strongly in uence the progress of tumors and the effect of immunotherapy. Many studies have demonstrated the important role of in ltration of macrophage, CD4 T cell, and CD8 T cell in tumor prognosis [41][42][43]. We found that DSC1 expression was negatively correlated with B cells and CD4+ T cells in ltration (P < 0.05). The expression of TGM1 identi ed negatively correlated to the in ltration level of four immune cells, including CD8+ T cells, CD4+ T cells, macrophage cells neutrophil cells, and dendritic cells (P < 0.05). However, the functions and pathways of these genes in tumor-in ltrating need to be further investigated.

Conclusions
In summary, 196 DEGs between the primary and metastatic skin cutaneous melanoma samples were screened, and 8 hub genes namely SPRR1B, DSC1, PKP1, TGM1, DSG1, IVL, SPRR1A and DSC3 were identi ed that may be associated with the metastasis of melanoma. Their functions and pathways involved in melanoma metastasis were explored. Moreover, two gene (DSC1 and TGM1) play an important role in the microenvironment of metastatic melanoma by regulating the tumor in ltration of immune cells. Further study is needed to investigate the special roles of these genes in the regulation metastatic melanoma progression, as well as their potential causal relationship with immune cells in ltration.

Consent for publication
Not applicable.

Availability of data
The data used to support the ndings of this study are available from the Gene Expression Omnibus (https://www.ncbi.nlm.nih.gov) and Cancer Genome Atlas database (https://cancergenome.nih.gov).

Competing interests
The authors declare that they have no con icts of interest.

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download. g.S1.tif