3.1 Identification of DEGs and construction of a prognostic neddylation-related genes signature in CM

The complete clinical information and corresponding transcriptome data from 455 CM patients were obtained from the TCGA database, including 455 CM tissues and 1 nontumor tissue. Transcriptome data of 812 normal skin specimens were obtained from the GTEx database. We conducted a comparison of the expression levels of 234 NRGs in a dataset comprising both the GTEx and TCGA databases, consisted of 813 nontumor tissues and 455 tumor tissues. Based on the predefined criteria of |logFC|> 1.0 and adj. P < 0.05, we identified a total of 182 DEGs, which have the potential to serve as prognostic genes. Out of these, 100 genes showed upregulation, while 82 genes exhibited downregulation in the tumor group. Subsequently, 455 CM cases were randomly partitioned into a training set (n = 273) and a testing set (n = 182) using an approximately equal split of 6:4. Table 2 presents the epidemiological and clinical characteristics of each CM individual diagnosed in the training, testing, and total sets. The prognostic potential of these NRGs was assessed by using univariate cox regression analysis on the overall survival (OS) data of CM patients in the training set, consequently, we identified 39 NRGs with significant prognostic value in CM patients (Fig. 2A). Out of these, 18 NRGs were identified and classified as “protective” genes, whereas the rest of the NRGs were categorized as “risk” genes. In order to further explore the association between the identified 39 NRGs and prognosis, we constructed a network diagram of prognostic genes (Fig. 2B). To alleviate multicollinearity and pinpoint the 39 prognostic NRGs exhibiting a substantial association with the survival outcome in CM patients, we employed a LASSO Cox regression analysis. According to the AIC value, eight distinct NRGs (ASB7, DTL, FBXO27, KLHL2, PSMA6, PSMB9, SPSB1, and WSB2) to construct a risk score model for forecasting the OS of CM patients, derived from the TCGA training set. The cvfit and lambda curves, providing a visual interpretation of the analysis outcomes, were shown in Fig. 2C and D. Then, the risk score of each CM individual was calculated using the following formula: Risk score = ASB7 * (−0.94848) + DTL * 0.46097 + FBXO27 * 0.19656 + KLHL2* (−0.45256) + PSMA6 * (−0.63426) + PSMB9 * (−0.25661) + SPSB1 * 0.27065 + WSB2 * 0.54281.

Table 2 The clinical characteristics of CM patients in the total, training and testing setsFig. 2

Identification of neddylation-related genes signature in CM. A The forest plot of prognostic NRGs (Protective genes: ASB6, BTBD1, COPS2, CUL2, DCAF6, FBXW7, KLHL2, KLHL9, NUB1, PSMA3, PSMA4, PSMA6, PSMB8, PSMB9, PSMC6, PSME1, PUM2, UBA3; Risk genes: ASB6, BIRC5, BTRC, CCNF, DCAF4, DTL, FBXL18, FBXO17, FBXO27, FBXW8, FBXW9, GAN, KLHL25, LRR1, PSMB7, PSMD2, PSMD9, SPSB1, UBA52, WDR5, WSB2). B The network diagram of prognostic genes. C LASSO coefficient profiles of NRGs. D The partial likelihood deviance with changing of log(λ). (* p < 0.05, ** p < 0.01, *** p < 0.001, **** p < 0.0001, ns, no significance). Abbreviation: CM, cutaneous melanoma; NRGs, neddylation-related genes

3.2 Evaluation of the prognostic model

The risk score of each CM individual in each dataset was calculated using the formula of the risk score and organized in a sequence from lowest to highest. Following this, based on whether their scores were either above or below the median values derived from the training set, CM individuals from every dataset were respectively classified into two distinct groups: high-risk and low-risk. The survival outcomes, risk classifications, and expression intensities of NRGs for each patient were respectively represented in Fig. 3A–I. This was done separately for the training set, testing set, and total set. The findings revealed a gradual increase in patient mortality corresponding to an escalating risk score. The Kaplan–Meier survival analysis uncovered a noticeably improved OS in the low-risk group in comparison to the high-risk group across all datasets (Fig. 3J–L). Our model demonstrated superior predictive precision for CM survival, as substantiated by the ROC curve analysis. This analysis revealed AUC values of 0.850, 0.722, and 0.746 for the timeframes of 1 year, 3 years, and 5 years, respectively (Fig. 3N). To assess its predictive potency with enhanced accuracy, congruent observations were noted in both the testing set and the total set (Fig. 3M and N).

Fig. 3

Validation of the NRGs signature in the total, training and test sets. A–I The distribution of the risk scores, the distributions of the overall survival status and heatmap of the expression of 8 NRGs in the total, training and test sets. J–L Kaplan–Meier curves for the overall survival of patients in the high- and low-risk groups in the total, training and test sets. M–O The area under the time-dependent ROC curve (AUC) of the time-dependent receiver operating characteristic curve (ROC) verified the prognostic accuracy of the risk score in predicting 1-, 3-, and 5-year OS in the total, training and test sets. Abbreviation: NRGs, neddylation-related genes; OS, overall survival

Next, in an endeavor to ascertain whether the risk score could serve as an independent prognostic factor for CM patients, we engaged in univariate and multivariate cox regression analysis, incorporating clinical factors. The univariate cox regression analysis revealed a significant correlation between OS and factors such as age (HR, 1.020; CI, 1.006–1.033; P = 0.005), stage (HR, 1.585; CI, 1.272–1.975; P < 0.001), T (HR, 1.373; CI, 1.170–1.611; P < 0.001), N (HR, 1.490; CI, 1.237–1.794; P < 0.001) and risk score (HR, 1.602; CI, 1.411–1.819; P < 0.001) in the training (Fig. 4B). Moreover, the multivariate cox regression analysis signified that N (HR, 1.677; CI, 1.244–2.261; P < 0.001) and risk score (HR, 1.623; CI, 1.404–1.877; P < 0.001) in the training, independently prognosticated OS in CM patients (Fig. 4E). Then, the predictive specificity and sensitivity of risk scores in CM patients were evaluated using the AUC in the train set. This analysis revealed higher values for risk scores compared to other clinicopathological factors, with an AUC of 0.850 (Fig. 4H). Similar statistical outcomes were observed within both the testing set and the total set (Fig. 4A, D, G, C, F and I). These discoveries underscored the importance of eight NRGs, highlighting their pivotal role as independent determinants influencing the prognosis of CM. Furthermore, PCA was performed to discern differences between the two subgroups utilizing all genes, NRGs, and the risk score of 8 NRGs, as illustrated in Fig. 5A–C. The PCA results highlighted that the analysis of the risk score manifested a more marked bifurcation in opposing directions between the low-risk and high-risk group in comparison to the others. These findings suggest that the risk score effectively classified CM patients into two distinct groups (low-risk and high-risk) possessing entirely disparate statuses.

Fig. 4

The prognosis value of the NRGs signature. A–C The result of univariate cox regression analysis based on the total, training and testing sets, respectively. D–F The result of multivariate cox regression analysis based on the total, training and testing sets, respectively. G–I Area under the time-dependent ROC curve (AUC) of receiver operating characteristic curve (ROC) curves comparing the prognostic accuracy of the risk score and other prognostic factors in the total, training and testing sets. Abbreviation: NRGs, neddylation-related genes

Fig. 5

The result of PCA. A PCA of all genes expression. B PCA of NRGs expression. C PCA of the prognostic 8 NRGs signature. Abbreviation: PCA, principal component analysis; NRGs, neddylation-related genes

3.3 Construction and verification of the nomogram in CM

Our nomogram integrated a 8-NRGs signature, along with age, gender, stage, T, N, and M as predictive variables. Following this, the prognostic nomogram was harnessed to assess the prognostic outcomes for patients with CM at intervals of 1, 3 and 5 years post-diagnosis (Fig. 6A). Furthermore, the calibration curves displayed satisfactory calibration (Fig. 6B) and the 3-year DCA curves also determined that the nomogram had good performance in clinical practice (Fig. 6C).

Fig. 6

Construction of nomogram models for CM. A A nomogram combining clinicopathological variables and risk score predicts the 1-, 3-, and 5- year overall survival. B The calibration curves for 1‐, 3‐, and 5‐year OS. C 3-year DCA curves. (* p < 0.05, ** p < 0.01, *** p < 0.001, **** p < 0.0001, ns, no significance). Abbreviation: Abbreviation: CM, cutaneous melanoma; OS, overall survival; DCA, Decision Curve Analysis

3.4 Gene set enrichment analysis (GSEA)

According to the median risk score, GSEA was executed to pinpoint significantly enriched pathways in both the high-risk and low-risk groups. The GSEA findings demonstrated that within the high-risk group, out of 178 signaling pathways, a total of 106 were upregulated, and 16 signaling pathways exhibited significant enrichment at a significance threshold of nominal (P < 0.01 and FDR < 25%). Conversely, in the low-risk group, out of the same 178 signaling pathways, 72 were upregulated, and 23 pathways indicated significant enrichment at a similar significance level of nominal (P < 0.01 and FDR < 25%). The enriched hallmark tumor‐associated pathways were lysine degradation, aminoacyl-tRNA biosynthesis, pyrimidine metabolism, base excision repair, purine metabolism, RNA polymerase, Alzheimer’s disease, glycerolipid metabolism, melanogenesis, glycolysis gluconeogenesis, Amino sugar and nucleotide sugar metabolism, mTOR signaling pathway and glycerophospholipid metabolism in the high-risk group (Fig. 7A), and the enriched hallmark pathways were predominantly immune-related pathway or tumor-related inhibitory pathways, such as antigen processing and presentation, cell adhesion molecules (CAMs), cytokine-cytokine receptor interaction, natural killer cell mediated cytotoxicity, JAK-STAT signaling pathway, chemokine signaling pathway, Toll-like receptor signaling pathway, intestinal immune network for IgA production, NOD-like receptor signaling pathway, T cell receptor signaling pathway, complement and coagulation cascades and RIG-I-like receptor signaling pathway in the low-risk group (Fig. 7B).

Fig. 7

Gene set enrichment analysis (GSEA) of the high- and low-risk groups based on the risk score

3.5 Analysis of the tumor microenvironment and immune infiltration

While scrutinizing the immune cell infiltration of all CM patients in the TCGA database using the CIBERSORT algorithm, there were discernible variances concerning 22 types of infiltrating immune cells between the two groups (Fig. 8A), and our findings unveiled a substantial elevation in the levels of memory B cells, CD8 T cells, CD4 memory T cells, Follicular helper T (Tfh) cells, and M1 macrophages within the low-risk group relative to the high-risk group. In contrast, the quantities of M0 macrophages, M2 and resting mast cells were significantly reduced in the low-risk group (Fig. 8B). Subsequently, the ESTIMATE algorithm was used to assess tumor microenvironment characteristics, specifically focusing on immune and stromal components. The ESTIMATE score combines the stromal score (reflecting the presence of stromal cells) and the immune score (indicating immune cell infiltration) to provide an overall assessment of the tumor’s microenvironment. We highlight that both the stromal and immune scores, as well as the overall ESTIMATE score, were significantly higher in the low-risk groups compared to the high-risk group (Fig. 9A). This finding suggests a more favorable immune landscape in the low-risk groups. The negative correlation between these scores and the risk score indicates that as risk increases, the infiltration of immune and stromal cells decreases (R = −0.22, P = 3.5e−6; R = −0.37, P < 2.2e−16; R = −0.33, P = 3.5e−13) (Fig. 9B–D), which could have implications for tumor progression and patient prognosis.

Fig. 8

Analysis of the immune cell infiltration in CM patients. A The proportions of 22 tumor infiltrating immune cells in individual CM patients by CIBERSORT analysis. B Violin diagram showing the immune cell composition between the two groups by CIBERSORT analysis. (* p < 0.05, ** p < 0.01, *** p < 0.001, **** p < 0.0001, ns, no significance). Abbreviation: CM, cutaneous melanoma

Fig. 9

Analysis of the immune cell infiltration in CM patients by ESTIMATE. A Stroma, immune, and ESTIMATE scores in the high-risk and low-risk groups in CM patients. B–D Correlation analysis between Stroma, immune, and ESTIMATE scores and the risk score. (* p < 0.05, ** p < 0.01, *** p < 0.001, **** p < 0.0001, ns, no significance). Abbreviation: CM, cutaneous melanoma

The ssGSEA result demonstrated that specimens from the low-risk group consistently manifested significantly elevated activity scores across multiple pathways, including Antigen-Presenting Cell (APC) co-inhibition, APC co-stimulation, Chemokine Receptor (CCR), check point, cytolytic activity, Human Leukocyte Antigen (HLA), inflammation promoting, Major Histocompatibility Complex (MHC) class I, parainflammation, T cell co-inhibition, T cell co-stimulation, Type I Interferon (IFN) response and Type II IFN response (Fig. 10A). We then evaluated the variations in expression of common checkpoint genes across both groups and observed that many of these genes displayed significant differences between the two groups and many tested immunotherapy targets, such as Programmed Cell Death 1 (PDCD1) (PD-1) and Cytotoxic T-Lymphocyte-Associated Protein 4 (CTLA4) were highly expressed in the low-risk group (Fig. 10B).

Fig. 10

Analysis of immune-related pathways in CM patients between the low-risk group and high-risk groups based on the ssGSEA scores. A Immune functional differences between the two groups. B The difference of common immune checkpoints expression in the two groups. (* p < 0.05, ** p < 0.01, *** p < 0.001, **** p < 0.0001, ns, no significance). Abbreviation: CM, cutaneous melanoma

3.6 The application of the NRGs signature in treatment with targeted therapeutics and chemotherapy drugs

It is widely acknowledged that targeted therapy against the gene BRAF and first-line chemotherapeutic agents in clinical treatment can enhance the survival outcomes of patients with CM, targeted therapy against the BRAF gene include dabrafenib, nilotinib, sorafenib and trametinib, and first-line chemotherapeutic agents, including docetaxel, oxaliplatin, paclitaxel and temozolomide according to the NCCN guidelines insights of cutaneous melanoma [2]. Thus, by evaluating the relationship between the risk score and the efficacy of these widely recognized anticancer drugs, it was observed that CM patients in the high-risk group demonstrated enhanced responsiveness to BRAF-targeted therapies, including nilotinib, sorafenib, and trametinib, relative to the low-risk group, and individuals in the low-risk group exhibited heightened responsiveness to conventional chemotherapy drugs, such as paclitaxel and temozolomide, but demonstrated reduced receptiveness to docetaxel compared to their counterparts in the high-risk group (Fig. 11). This finding has potential implications for the implementation of personalized therapy for patients with CM.

Fig. 11

The chemotherapeutic responses of the two groups to common anticancer drugs (targeted therapeutics and chemotherapy drugs). A Dabrafenib. B Nilotinib. C Sorafenib. D Trametinib. E Oxaliplatin. F Docetaxel. G Paclitaxel. H Temozolomide

3.7 Investigation of the expression alterations of the 8 NRGs

The NRGs signature was composed of eight genes, we investigated the expression levels of these genes in human melanoma cell lines (SK5 and MNT-1) and human normal melanocyte cell line (PIG1). Although the expression difference of some genes was not statistically significant, the result demonstrated that a similar trend in expression between human melanoma cell lines and human normal melanocyte cell line (Fig. 12A–H). Thus, these findings underscore the robustness of our model and indicate that these NRGs may have substantial implications in CM.

Fig. 12

The expression alterations of 8 NRGs in human melanoma cell lines (SK5 and MNT-1) and human normal melanocyte cell line (PIG1). (* p < 0.05, ** p < 0.01, *** p < 0.001, **** p < 0.0001, ns, no significance). Abbreviation: NRGs, neddylation-related genes