Gene expression data for 46 metastatic SKCM samples, 77 primary SKCM samples and 21 normal skin samples were obtained from the GEO database (https://www.ncbi.nlm.nih.gov/geo/), specifically from the GSE8401 (21), GSE15605 (22), GSE46517 (23) and GSE65904 (24) datasets. Raw data were log-transformed and normalized prior to downstream analyses (Fig. S1). A summary of patient characteristics is presented in Tables I and SI. Moreover, RNA sequencing (RNA-seq) data from the TCGA-SKCM dataset were retrieved from the University of California, Santa Cruz Xena platform (https://xena.ucsc.edu/) and converted to transcripts per million (TPM) using the ‘count2tpm’ R package (version 4.3.0; http://www.r-project.org/). Following quality control and filtering for complete survival data, 461/472 patients in the TCGA-SKCM dataset with matched RNA-seq and survival information were included in the subsequent analyses.

In addition, seven pairs of SKCM and adjacent normal tissue samples were obtained from patients who underwent surgical treatment at the Institute of Dermatology, Peking Union Medical College and Chinese Academy of Medical Sciences (Nanjing, China) between 2018 and 2021. A total of 7 patients, aged 26–92 years (4 males and 3 females) were included. All paraffin-embedded tissue sections were confirmed as cutaneous melanoma by two or more dermatopathologists, with cases declining to participate excluded. The fresh tumor and adjacent tissues were promptly frozen and stored at −80°C for subsequent analyses. The present study was approved by the Ethics Committee of the Institute of Dermatology, Peking Union Medical College and Chinese Academy of Medical Sciences (approval no. 2017-KY-044), and all procedures were performed in accordance with approved institutional guidelines. Written informed consent was obtained from all participants.

TCGA-SKCM cases were classified into primary and metastatic groups based on patient clinical information. Least absolute shrinkage and selection operator (LASSO) Cox regression analysis was performed using the ‘glmnet’ package in R (version 4.3.0; http://www.r-project.org/), with the optimal λ value determined using 10-fold cross-validation (25). This approach identified the most predictive gene set, which, combined with primary-metastatic status, was used to construct a risk score model. Individual risk scores were calculated for each patient, who were then stratified into high- and low-risk groups using the median risk score as the cutoff. Survival differences between the two groups were assessed using Kaplan-Meier analysis and compared using the log-rank test, with significance defined as P<0.05. The predictive accuracy of the risk model was evaluated by receiver operating characteristic (ROC) curve analysis, and model performance was quantified using the area under the curve. Internal robustness was assessed using 10-fold cross-validation of the LASSO-Cox model. The resulting risk scores for each fold and the coefficients selected at the λ_min and λ_1se penalty parameters are summarized in Table SII. To address potential violations of the proportional hazard rate assumption due to survival curve crossover for certain genes, the analyzed time range was restricted to periods preceding the crossover points. Survival analyses were re-performed within these truncated time intervals to ensure robustness of the log-rank test results (26).

Based on the risk scores generated by the machine learning model, the TCGA-SKCM cohort was stratified into high- and low-risk groups. Differential expression analyses were performed to compare gene expression profiles among primary tumors, metastatic tumors and normal tissues from the GEO datasets, as well as between the high- and low-risk groups in the TCGA cohort. For the TCGA-SKCM RNA-seq data, raw read-count matrices were analyzed using the ‘DESeq2’ package (R version 4.3.0; http://www.r-project.org/) and the median-of-ratios procedure was used to normalize library size and estimate dispersion, followed by Wald tests to identify DEGs. Genes with adjusted P<0.05 and |log₂FoldChange (FC)|>1 were deemed significant. TPM values, calculated using the count2tpm tool (R version 4.3.0; http://www.r-project.org/), were used only for downstream visualization (27). For the GEO datasets, raw or processed expression matrices were analyzed using the ‘limma’ package (R version 4.3.0; http://www.r-project.org/) with empirical Bayes moderation to identify DEGs, applying the same adjusted P-value and fold-change thresholds. DEGs identified from multiple GEO datasets were then integrated using the ‘RobustRankAggreg’ (RRA) package (R version 4.3.0; http://www.r–project.org/) (28) to obtain a robust DEG list.

For single-gene differential expression analyses, thresholds were set at |log2FC|>1 and adjusted P<0.05 to identify differentially expressed lncRNAs and mRNAs associated with each of the four candidate genes. Overlapping lncRNAs and mRNAs were determined using the ‘VennDiagram’ package in R (version 4.3.0; http://www.r-project.org/). Visualization of DEGs was achieved using volcano plots and heatmaps generated using the ‘ggplot2’ package (R version 4.3.0; http://www.r-project.org/).

To assess the biological processes and signaling pathways underlying SKCM metastasis, enrichment analyses were performed on DEGs identified from GEO datasets comparing primary and metastatic melanoma, and co-expressed DEGs associated with the risk-related genes [keratin (KRT)17, plakophilin 1 (PKP1), cadherin 3 (CDH3) and cellular retinoic acid binding protein 2 (CRABP2)]. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were performed using the Metascape database (http://metascape.org/) and the ‘clusterProfiler’ R package (version 4.3.0; http://www.r–project.org/) (29). The top 20 enriched GO and KEGG terms were visualized using network analysis in Cytoscape 3.6.1 (30). Additionally, gene set enrichment analysis (GSEA) was performed using clusterProfiler, with significance thresholds set at false discovery rate (FDR) <0.25 and adjusted P<0.05 (31).

Based on previous studies (32), 21 key m6A regulatory genes [FTO, ALKBH5, RBM15, RBM15B, METTL3, METTL14, METTL16, WTAP, VIRMA, zinc finger CCCH-type containing 13 (ZC3H13), HNRNPC, HNRNPA2B1, insulin like growth factor 2 mRNA binding protein (IGF2BP)1, IGF2BP2, IGF2BP3, RBMX, YTHDC1, YTHDC2, YTHDF1, YTHDF2 and YTHDF3] were selected to assess the potential association between m6A methylation regulation and metastatic SKCM. Pearson correlation analysis was performed to assess the relationships between the expression levels of these m6A regulators and key genes, with |R|>0.2 and P<0.05 considered statistically significant.

Potential miRNAs targeting key mRNAs were predicted using the miRWalk database (http://mirwalk.umm.uni-heidelberg.de/). Concurrently, lncRNA-miRNA interactions were identified using the TarBase v8 database (http://carolina.imis.athena-innovation.gr/diana_tools/web/index.php?r=lncbasev2%2Findex). Common miRNAs predicted by both approaches were determined using Venn diagram analysis using the VennDiagram R package. A lncRNA-miRNA-mRNA ceRNA network was subsequently constructed and visualized using Cytoscape (v3.6.1; http://cytoscape.org/release_notes_3_6_1.html). Key regulatory nodes within the network were identified using the ‘cytoHubba’ plugin. To further characterize the lncRNAs, sequences were retrieved from LNCipedia (https://lncipedia.org/), and their subcellular localizations were predicted using lncLocator (http://www.csbio.sjtu.edu.cn/bioinf/lncLocator/). The subcellular localization of mRNAs was determined using the DM3Loc web tool (http://dm3loc.lin-group.cn/). Potential miRNA binding sites within genomic sequences were identified using miRanda (https://www.bioinformatics.com.cn/local_miranda_miRNA_target_prediction_120) (33). Only lncRNAs and mRNAs predicted to localize mainly in the cytoplasm (lncLocator or DM3Loc probability ≥0.6) were retained, ensuring structural plausibility for ceRNA crosstalk.

Relative enrichment scores for 24 immune subsets were calculated from log₂(TPM+1) matrices with the GSVA package (v1.34.0; method=‘ssgsea’, kcdf=‘Gaussian’) (34,35). Absolute cell fractions were inferred using CIBERSORTx (https://cibersortx.stanford.edu/index.php) in ‘absolute mode’ with the LM22 signature and 1,000 permutations; samples with P≥0.05 were discarded. Deconvolution was repeated using the EPIC algorithm implemented in TIMER2.0 (http://timer.cistrome.org/TIMER2.0) to provide an additional independent estimate. To evaluate the immune and stromal components of the tumor microenvironment, the Estimation of STromal and Immune cells in MAlignant Tumor tissues using Expression data (ESTIMATE) algorithm was applied using the ‘estimate’ R package (version 4.3.0; http://www.r-project.org/). Immune scores, stromal scores and tumor purity were calculated for each sample according to the original algorithm. These packages were used to assess the correlation between the expression of m6A-associated genes (KRT17, PKP1, CDH3 and CRABP2) and the infiltration levels of 24 immune cell types. Pearson correlation analysis was applied, with |R|>0.2 and P<0.05 considered statistically significant.

Tissues were fixed in 10% neutral buffered formalin at 4°C overnight, followed by paraffin embedding and sectioning at a 4–5 µm thickness. Sections were sequentially deparaffinized in xylene, dehydrated and rehydrated through a graded ethanol series, and finally immersed in distilled water. Subsequently, 10% goat serum (Gibco; Thermo Fisher Scientific, Inc.) was applied to the sections and incubated at room temperature for 1 h to complete blocking. Sections were then incubated overnight at 4°C with primary antibodies against PKP1 (1:500; cat. no. ab154622; Abcam), KRT17 (1:200; cat. no. ab109725; Abcam), CRABP2 (1:1,000; cat. no. ab211927; Abcam) and CDH3 (1:250; cat. no. ab242060; Abcam). After washing, the slides were incubated with HRP-conjugated secondary antibody (1:2,000; cat. no. ab6759; Abcam) in a humidified chamber for 30 min, followed by color development using a DAB detection kit (Beyotime Institute of Biotechnology). A total of two independent pathologists evaluated the stained sections under an optical microscope (Olympus Corporation). The proportion of positively stained cells was scored as follows: 0, 0–5%; 1, 6–25%; 2, 26–50%; 3, 51–75%; and 4, >75%. Staining intensity was graded as follows: 0, negative; 1, light yellow; 2, yellow-to-brown; and 3, brown. The final immunohistochemical score was calculated as the product of the proportion and intensity scores, and classified into four categories: Negative (−; 0), weakly positive (+; 1–4), moderately positive (++; 5–8) and strongly positive (+++; 9–12).

Total RNA was extracted from SKCM and matched adjacent normal tissues using TRIzol™ reagent (Invitrogen™; Thermo Fisher Scientific, Inc.). cDNA synthesis was performed using the EVO M-MLV RT Mix Kit [cat. no. AG11728; Accurate Biotechnology (Hunan) Co., Ltd.] according to the manufacturer's instructions. qPCR was then performed using the SYBR Green Premix Pro Taq HS qPCR Kit [cat. no. AG11701; Accurate Biotechnology (Hunan) Co., Ltd.] on a LightCycler^® 480 Real-Time PCR System (Roche Diagnostics). Each 10 µl reaction mixture contained 1 µl cDNA, 5 µl 2X SYBR Green Pro Taq HS Premix, 0.5 µM of each primer and nuclease-free water. The qPCR cycling conditions were as follows: Initial denaturation at 94°C for 30 sec, followed by 40 cycles at 94°C for 5 sec, 60°C for 10 sec and 50°C for 30 sec. Relative gene expression levels were calculated using the 2^−ΔΔCq method (36). The primer sequences used in the present study are listed in Table II.

All statistical analyses, except those involving transcriptomic sequencing data, were performed using GraphPad Prism (version 8.0; Dotmatics) and SPSS software (version 23.0; IBM Corp.). For comparisons of immunohistochemistry scores between tumor and adjacent normal tissues, non-parametric tests were applied: The Mann-Whitney U test for two-group comparisons or the Kruskal-Wallis test for multiple groups, followed by Dunn's post hoc test with Benjamini-Hochberg correction for multiple comparisons. For comparisons among more than two groups with normally distributed data, one-way ANOVA was used, followed by Tukey's post hoc test for multiple pairwise comparisons. RT-qPCR data were analyzed using the Wilcoxon signed-rank test to compare expression levels between tumor and matched normal tissues. Unless otherwise specified, P-values were adjusted for multiple comparisons using the Benjamini-Hochberg false-discovery-rate (FDR) procedure. P<0.05 was considered to indicate a statistically significant difference.

To identify key genes associated with SKCM metastasis, a comprehensive analysis of multiple GEO datasets was performed. Prior to analysis, gene expression data were normalized (Fig. S1), and the ‘limma’ package in R was used to identify DEGs within each microarray dataset, using adjusted P<0.05 and |log₂FC|>1 as cutoff criteria (Fig. S2). Subsequently, DEGs from all datasets were integrated using the RRA package, resulting in the identification of 94 consistently differentially expressed mRNAs, 84 upregulated and 10 downregulated, in metastatic SKCM samples compared with in primary SKCM samples across the four datasets (Fig. 1A).

Subsequently, GO biological process enrichment analysis was performed using the ‘clusterProfiler’ package in R. The results indicated that pathways related to cornification, epidermis development, and keratinization were significantly enriched, suggesting their critical involvement in SKCM metastasis. Notably, genes such as cathepsin G (CTSG), serine peptidase inhibitor kazal type 1 and S100 calcium binding protein A8 (S100A8) were identified as key contributors within these pathways (Fig. 1B and C).

In addition, the correlations between the 94 metastasis-related DEGs and 21 m6A regulatory genes were evaluated using Pearson correlation analysis, with thresholds set at |R|>0.5 and P<0.05. DEGs whose expression levels were significantly correlated with ≥1 m6A regulator were defined as m6A-related mRNAs. Based on this criterion, a total of 45 mRNAs were identified that showed significant correlations with m6A (Fig. 1D).

A LASSO-Cox analysis of the upregulated pathways and related genes involved in SKCM metastasis was subsequently performed. The partial likelihood deviance in relation to log(λ) was derived from the LASSO-Cox regression model (Fig. 2A and B). Using the TCGA transcriptomic data and survival information, 12 SKCM metastasis regulators associated with overall survival (OS) were identified: KRT5, KRT17, PKP1, KRT15, filaggrin (FLG), cystatin A (CSTA), CRABP2, cystatin E/M (CST6), calmodulin like (CALML)5, CDH3, S100A8 and CTSG. The risk score and survival status distribution indicated that patient survival time decreased with an elevated risk score. Moreover, combined with the gene expression distribution, higher expression of these genes correlated with a higher risk score (Fig. 2C). Furthermore, survival analysis revealed that patients in the high-risk group exhibited significantly worse long-term survival than those in the low-risk group [hazard ratio (HR), 2.62; P<0.001; Fig. 2D]. Ultimately, the present study focused on four mRNAs strongly related to m6A regulation (CDH3, KRT17, PKP1 and CRABP2). These genes appear to serve crucial roles in SKCM metastasis.

According to previous studies, mRNA expression is modulated by miRNAs, whilst lncRNAs can competitively bind miRNAs and thus indirectly regulate mRNA expression. This lncRNA-miRNA-mRNA triple regulatory network is known as a ceRNA network (37). To clarify the regulation of metastasis-related mRNAs in SKCM by lncRNAs via miRNAs, the present study constructed a ceRNA network centered on CDH3, KRT17, PKP1 and CRABP2. First, single-gene differential expression analysis was performed for each of the four genes in the TCGA-SKCM cohort, identifying significantly co-expressed or differentially up or downregulated transcripts, and for common differentially expressed mRNAs and lncRNAs (Fig. S3). Ultimately, 494 differentially expressed mRNAs and 91 differentially expressed lncRNAs were identified (Fig. 3A).

Subsequently, miRNAs predicted to target the 91 lncRNAs (TarBase) and the 4 mRNAs (miRWalk) were screened, identifying the shared miRNA interactions between lncRNAs and mRNAs. Consequently, 13 lncRNAs, 20 miRNAs and 4 mRNAs were integrated into the ceRNA network (Fig. 3B). Kaplan-Meier survival analysis revealed that KRT17, PKP1, CDH3, CRABP2, SOX21-AS1, hsa-miR-3200-3p, hsa-miR-1976, hsa-miR-1276 and hsa-miR-151a-3p were all significantly associated with OS in SKCM. HR analysis also demonstrated that KRT17, PKP1, CDH3, CRABP2, SOX21-AS1, hsa-miR-1276 and hsa-miR-151a-3p are risk factors for SKCM, whereas hsa-miR-3200-3p and hsa-miR-1976 serve as protective factors (Fig. 3C).

Additionally, co-expressed differentially expressed mRNAs associated with KRT17, PKP1, CDH3, and CRABP2 were subjected to functional enrichment analysis using the Metascape online tool. The results revealed that these genes are closely associated with the estrogen signaling pathway, Ras signaling pathway, cytokine-cytokine receptor interaction, melanogenesis and chemical carcinogenesis, with the CALML3/5 and phospholipase A2 group IVA families serving important roles in these pathways (Fig. 3D).

To further elucidate the abnormal m6A regulatory mechanisms in SKCM, the expression of 21 m6A regulators in SKCM and normal tissues from TCGA were analyzed. A total of four downregulated regulators (METTL3, WTAP, YTHDC1 and YTHDC2) and 13 upregulated regulators (ALKBH5, VIRMA, RBM15, ZC3H13, IGF2BP1, IGF2BP2, IGF2BP3, YTHDF1, YTHDF2, YTHDF3, HNRNPC, HNRNPA2B1 and RBMX) were identified in SKCM (Fig. 4A). The present study also evaluated how these 21 m6A regulators vary across different clinical stages (I–IV), and the ANOVA results revealed that ALKBH5, METTL3, WTAP, RBM15, YTHDF2, HNRNPC, HNRNPA2B1, RBMX, YTHDC1 and IGF2BP3 displayed continuous patterns of increased or decreased mRNA expression with advancing SKCM stage (Fig. 4B). This suggests that certain m6A regulators may be closely associated with SKCM tumorigenesis and progression.

The associations between the four key metastatic prognostic mRNAs (KRT17, PKP1, CDH3 and CRABP2) and the aforementioned 21 m6A regulators were subsequently assessed. The results showed that the expression level of KRT17 was significantly associated with METTL3, RBM15, ALKBH5, IGF2BP3 and HNRNPC. The expression level of PKP1 was significantly associated with METTL16, ALKBH5, IGF2BP3, YTHDF1 and YTHDF2. The expression level of CDH3 was significantly associated with VIRMA, RBM15B, YTHDC2, IGF2BP1 and RBMX, whereas the expression level of CRABP2 was significantly associated with METTL3, RBM15B, FTO, ALKBH5, IGF2BP1 and YTHDF1 (Fig. 5).

To thoroughly assess the influence of KRT17, PKP1, CDH3 and CRABP2 expression on immune cell infiltration, we used single-sample (ss)GSEA to assess correlations between these genes and 24 tumor-infiltrating immune cell types in 472 SKCM samples. The results indicated that PKP1 is positively correlated with mast cells and Th17 cells (Pearson-value >0.2; P<0.05; (Fig. 6A and B); KRT17 is negatively correlated with Tγδ (Tgd) cells, T helper cells (Th cells), and Th1 cells (Pearson-value <-0.2; P<0.05); CDH3 is positively correlated with mast cells and Th17 cells but negatively correlated with B cells and T helper cells; and CRABP2 is positively correlated with mast cells, natural killer (NK) cells, immature dendritic cells and neutrophils, and negatively correlated with B cells (Fig. 6A and B). To corroborate the ssGSEA-derived landscape, the immune infiltration data was re-evaluated using ESTIMATE and CIBERSORTx. The results revealed that tumors with a high expression of PKP1, KRT17, CRABP2 or CDH3 displayed significantly lower StromalScore, ImmuneScore and composite ESTIMATEScore than their low-expression counterparts (all FDR <0.01; Fig. S4), mirroring the global immune-suppressed phenotype observed with ssGSEA. Consistent with the ssGSEA results, CIBERSORTx deconvolution analysis further characterized the tumor immune microenvironment. Specifically, high CRABP2 expression was associated with significantly reduced levels of CD8⁺ T cells (Figs. 6B and S4). By contrast, both high-PKP1 and high-CRABP2 tumors showed increased levels of NK cells (Fig. 6B). Furthermore, both tumor subtypes shared a common feature of an increased M0-to-M2 macrophage fraction (all FDR <0.05; Fig. S4). The concordant results from the three independent algorithms strengthen the conclusion that elevated expression of these keratin-related genes is associated with an immunologically ‘cold’ microenvironment in melanoma.

ROC curves revealed that KRT17, PKP1, CDH3 and CRABP2 demonstrated predictive power for OS in the TCGA cohort of patients with SKCM (Fig. 7A). The associations of these genes with clinical factors was then further assessed. KRT17, PKP1, CDH3, and CRABP2 demonstrated a positive association with pathological stage (Fig. 7B). For tumor (T)-node (N)-metastasis (M) staging, their expression levels were significantly associated with the T stage, but had no discernable association with N or M stages (Fig. 7C). Similarly, no significant association was revealed between their expression and patient age (Fig. 7D).

Representative immunohistochemistry images demonstrated markedly higher cytoplasmic staining for all four genes in tumor tissues compared with in adjacent skin (Fig. 8). Consistently, RT-qPCR analysis on the same sample set showed 2.3- to 4.5-fold upregulation at the mRNA level (Fig. S5). These concordant protein and transcript results corroborate the bioinformatic prediction of the present study that the four hub genes are overexpressed in metastatic SKCM.