Abstract

Introduction

Melanoma is a highly malignant and aggressive skin cancer that metastasizes early through lymphatic or hematogenous routes, leading to poor prognosis[1] . Despite accounting for only 2% of all skin cancers, melanoma has the highest mortality rate among them. Globally, Skin Cutaneous Melanoma (SKCM) accounts for more than 80% of deaths related to skin cancer annually, making it one of the most lethal cancers threatening human health. Both the incidence and mortality rates of SKCM have steadily increased over the past few decades, with a growing prevalence among younger populations[2] . Originating from melanocytes in the skin's basal layer, SKCM is highly invasive and ranks as the fourth leading cause of cancer-related deaths worldwide[2] . SKCM can be classified into benign tumors, primary malignant tumors, and metastatic melanoma. The five-year survival rate for patients with primary melanoma exceeds 95%, while those with metastatic melanoma rarely survive beyond one year[3] . Surgical resection is the primary treatment for early-stage SKCM, but advanced cases present significant challenges, often requiring radiotherapy and chemotherapy. Early detection through biomarkers is crucial for improving outcomes. Despite progress in diagnosing and prognosticating SKCM using markers like BRAF and NRAS mutations and immune checkpoints such as PD-L1, identifying additional molecular characteristics for advanced SKCM remains a challenge[4] .

The centrosome, a critical organelle in human cells, serves as the microtubule-organizing center in animal cells and is composed of two orthogonal centrioles embedded in the pericentriolar material (PCM). It plays an essential role in cellular processes, including signal transduction, cell polarity, division, and migration[5] . Defects in centrosome structure and function are associated with cancer. The most common defect is centrosome amplification (CA), which refers to the presence of more than two centrosomes or abnormally large centrosomes in a cell, including numerical and structural amplification[6] . Numerous studies have identified CA as a hallmark of cancer, often linked to abnormal tumor karyotypes and poor clinical outcomes. Mechanistically, CA impairs mitotic fidelity, leading to chromosomal instability, which underpins tumor initiation and progression [7] . Studies have indicated that CA is associated with the excessive replication of centrioles in human melanoma[8] . However, research into the molecular mechanisms of CA in melanoma remains limited. Identifying the genetic mechanisms underlying centrosome amplification-related genes (CA-RGs) in SKCM will contribute to understanding the mechanisms driving this cancer’s development.

This study employs SKCM transcriptome data and CA-RGs to conduct bioinformatics analyses, aiming to identify prognostic markers and elucidate their molecular mechanisms, thereby offering new perspectives for SKCM treatment..

Materials and Methods

Analysis of gene expression differences between SKCM and normal skin tissues

To identify differentially expressed genes (DEGs) between SKCM and non-tumor tissues, we downloaded gene transcriptome data and clinical information from the GEO dataset GSE15605 (https://www.ncbi.nlm.nih.gov/geo/database). The dataset comprises 58 SKCM samples and 16 normal skin tissue samples. Differential expression analysis was performed using the R package limma [9] . DEGs were filtered based on the criteria of P < 0.05 and |log2FC| > 0.05. A volcano plot was generated using the R package ggplot2 [10] , and a heatmap was created using the R package pheatmap [11] .

Identification of centrosome amplification-related genes in SKCM and functional analysis

To identify differentially expressed centrosome amplification-related genes (CA-RGs) in SKCM tissues, we obtained a set of 76 CA-RGs from the Gene Ontology database (http://geneontology.org/). The intersection of these genes with the differentially expressed genes identified in SKCM provided the differentially expressed CA-RGs. A Venn diagram was created using the R package ggvenn[12] . Next, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were conducted on these CA-RGs using the R package clusterProfiler[13] to identify potential functional pathways involved in melanoma development. A pentagon bubble plot was generated using the R package ggplot2[10] , and a chord diagram was created with the R package GOplot[14] . Additionally, to explore the protein-protein interaction (PPI) network among the differentially expressed CA-RGs, we used the STRING database (https://cn.string-db.org/) with a confidence score >0.4, and visualized the interaction network using Cytoscape software[15] .

Construction and validation of centrosome amplification-related risk models

To construct a risk model based on the differentially expressed CA-RGs, we downloaded the GDC TCGA-SKCM dataset from the UCSC Xena database (https://xena.ucsc.edu/), which includes gene expression data, clinical information, and survival data for 457 SKCM samples. First, we conducted univariate Cox regression analysis using the R package survival[16] , selecting genes with P < 0.05 as input features for LASSO Cox regression. The R package glmnet [17] was used to fit a generalized linear model with penalized maximum likelihood, applying L1 regularization to remove features with a penalty score of 0. The remaining non-zero features were identified as candidate genes for constructing the risk model. To ensure the validity and reliability of the Cox regression model, we tested the proportional hazards (PH) assumption for the candidate genes using the R package survival. Genes with P > 0.05 were considered to pass the PH test. Next, we performed multivariate Cox regression on the genes that passed the PH test, using the step function for stepwise Cox regression to identify prognostic genes. The regression coefficients for each gene were obtained from the model. A forest plot displaying the results of the univariate and multivariate Cox regression analyses was generated using the R package forestplot[18] .The risk score for each sample in the TCGA-SKCM cohort was calculated using the following formula: Risk score = ∑ Exp (mRNA) ∗ Coef (mRNA), where Exp represents the mRNA expression level of each prognostic gene, and Coef represents the corresponding Cox regression coefficient. The 457 samples in the TCGA-SKCM cohort were divided into high-risk and low-risk groups based on the optimal cut-off value of the risk score. The distribution of risk scores and survival status of the samples were visualized using the R package ggplot2. We then compared the survival probabilities between the high-risk and low-risk groups using the R package survival and plotted Kaplan-Meier survival curves with the R package survminer[19] . Receiver operating characteristic (ROC) curves were generated with ggplot2, and the area under the curve (AUC) for 1-year, 3-year, and 5-year survival was estimated using the R package survivalROC[20] to assess the sensitivity and specificity of the risk model. A heatmap showing the expression of prognostic genes across TCGA-SKCM samples was created with the R package pheatmap. To validate the robustness and effectiveness of the risk model, we downloaded an independent SKCM dataset, GSE100797, from the GEO database (https://www.ncbi.nlm.nih.gov/geo/database), which includes gene expression data and survival information for 25 SKCM samples.

Construction and evaluation of prognostic models

To evaluate the prognostic capability of combining risk scores with clinical factors, we constructed a prognostic model in the TCGA-SKCM cohort. First, we used the Wilcoxon test or Kruskal-Wallis test to assess differences in risk scores across six clinical subgroups: age, sex, stage, T stage, N stage, and M stage. A P value < 0.05 was considered statistically significant, and the results were visualized as boxplots using the R package ggplot2. We then analyzed the survival differences between high-risk and low-risk patients within subgroups showing significant differences in risk scores. Next, univariate Cox regression analysis was performed based on risk scores, age, sex, stage, T stage, N stage, and M stage in the TCGA-SKCM cohort. Factors with P < 0.05 were subjected to the proportional hazards (PH) assumption test, and those with P > 0.05 were included in the multivariate Cox regression model. Factors with P < 0.05 in the multivariate analysis were incorporated into the final prognostic model. The R package survival was used for survival analysis, univariate Cox regression, PH assumption testing, multivariate Cox regression, and model construction. Kaplan-Meier survival curves were plotted using survminer, while forest plots were generated using forestplot to display the results of both univariate and multivariate Cox analyses. A nomogram was drawn using the R package regplot[21] .Subsequently, we assessed the performance of the prognostic model using calibration curves, ROC curves, and decision curve analysis (DCA). Model calibration was performed using the R package rms[22] . The specificity and sensitivity of the ROC curves were computed with timeROC [23] , and decision curve data were calculated using ggDCA[24] .

Gene set enrichment analysis

To further explore the biological pathways involving prognostic genes, we first calculated the Spearman correlation coefficients between genes in the TCGA-SKCM cohort and ranked the genes based on these coefficients. Next, KEGG pathway enrichment analysis was conducted using the R package clusterProfiler, with significance thresholds set at p.adjust < 0.05 and |NES| > 1. The top five most significantly enriched pathways were visualized using the R package enrichplot.

Functional enrichment analysis

To investigate pathway alterations between the high- and low-risk groups in the TCGA-SKCM cohort, we used Gene Set Variation Analysis (GSVA) to estimate the biological functions and signaling pathways of the prognostic genes. Specifically, molecular hallmark scores were calculated for each patient using the Hallmarks database in the R package GSVA[25] . Differential pathway analysis between the high- and low-risk groups was performed using the R package limma. All pathways were visualized using ggplot2, with P < 0.05 considered statistically significant.

Immune infiltration analysis

To investigate immune cell infiltration in SKCM, we used the R package xCell[26] to calculate infiltration scores for 34 types of immune cells in the TCGA-SKCM cohort. Differences in cell scores between high- and low-risk groups were compared, with cells showing P < 0.05 identified as significantly different. The correlations between these different immune cells and their associations with risk scores were also analyzed. ggplot2 was used to visualize cell infiltration scores as bar plots and correlation lollipop charts, while ggpubr[27] was used to generate box plots for differential analysis. Additionally, corrplot[28] was employed to display correlation heatmaps. To further analyze the relationship between immune and stromal cell infiltration and prognostic genes across the high- and low-risk groups, we first calculated immune, stromal, and composite scores for each sample in the TCGA-SKCM cohort using the R package estimate[29] . Differences in scores between the two risk groups were then compared, and the Spearman correlation coefficients between these scores and risk scores were calculated. The ggplot2 package was used to generate violin plots for differential score analysis and scatter plots for correlation analysis.

Analysis of immune response capacity and immune status

To assess immune response capacity and immune status between high- and low-risk groups, we first downloaded melanoma clinical data from the TCIA database (https://tcia.at/home) and extracted patients’ Immune Phenotype Scores (IPS). The IPS differences between the high- and low-risk groups in the TCGA-SKCM cohort were then compared. Additionally, we extracted 42 immune checkpoints from the literature (PMID: 37215879) and compared the expression levels of these immune checkpoints between the high- and low-risk groups in the TCGA-SKCM cohort. The ggplot2 package was used to generate violin plots to visualize differences in IPS scores and box plots for immune checkpoint expression differences.

Drug sensitivity analysis

To evaluate the sensitivity of SKCM patients to common chemotherapeutic agents, we utilized the GDSC database (https://www.cancerrxgene.org/) to predict the half-maximal inhibitory concentration (IC50) values for 198 chemotherapeutic and targeted therapy drugs in the TCGA-SKCM cohort. Using the pRRophetic package[30] , we compared the IC50 values between high- and low-risk groups for each drug. The most significant differences in IC50 values were visualized using violin plots generated with the ggplot2 package.

Somatic variants analysis

To investigate the differences in gene mutations between high- and low-risk groups, we utilized the TCGAmutations package[31] to download somatic mutation data for melanoma. Subsequently, we employed the maftools package[32] to generate waterfall plots depicting the top 20 most frequently mutated genes in the TCGA-SKCM cohort for both risk groups.

Analysis of prognostic genes regulatory mechanisms

To further explore the interaction between prognostic genes and other genes or elements, we first accessed the miRNet database (https://www.mirnet.ca/) to identify transcription factors (TFs) targeting the prognostic genes. We selected the JASPAR database to construct a TF-gene regulatory network. Subsequently, we utilized the GeneMANIA database (https://genemania.org/) to identify genes associated with the prognostic genes, examining their interactions, including protein-protein, protein-DNA, and genetic interactions, as well as pathways, reactions, gene and protein expression data, and protein domain information. The interaction network and the TF-gene regulatory network were visualized using Cytoscape software.

Single-cell transcriptome sequencing analysis

To analyze the expression distribution of prognostic genes at the single-cell level, we used the single-cell dataset GSE215120 to explore SKCM pathogenesis. Raw data from three SKCM samples were downloaded from the GEO database (https://www.ncbi.nlm.nih.gov/geo/). Cells were filtered based on UMI counts (<500 or >20,000), mitochondrial content (>20%), and features (<200 or >6,000). Data normalization was performed using the NormalizeData function (scale factor=10,000), and variable genes were identified using the FindVariableGenes function (vst method; top 2,000 features). Features were normalized and centered using the ScaleData function. Principal component analysis (PCA) was conducted using the RunPCA function (50 principal components retained), with significance assessed via the JackStraw function. Cell clusters were identified using the FindClusters function (SNN modular optimization; resolution=0.8), and dimensionality reduction was performed with the RunTSNE function, visualized using DimPlot. Clusters were annotated using the CellMarker 2.0 database and manually curated marker genes. The top three cell types with the highest expression of prognostic genes were identified as key cells. All analyses were performed using the Seurat package[33] . To explore receptor-ligand interactions, we referenced the CellChatDB.human database and inferred interactions using the CellChat package[34] . The c Communication networks were visualized with network diagrams and heatmaps. Signaling pathways were inferred using the netAnalysis_computeCentrality function, and signal senders/receivers were visualized with the netAnalysis_signalingRole_scatter function. Pseudotime analysis was conducted using the monocle package[35] to assess dynamic changes in key cells. Cell trajectories were displayed using the plot_cell_trajectory function, while prognostic gene expression changes were visualized with the plot_genes_branched_heatmap and plot_genes_in_pseudotime functions.

Analysis of prognostic genes expression

To analyze the expression of prognostic genes, we first compared the transcriptomic expression differences between SKCM tissues and normal skin tissues in the dataset GSE15605, using the ggplot2 package to present the results as boxplots. To validate gene expression at the protein level, we downloaded immunofluorescence images of prognostic genes from the HPA database (https://www.proteinatlas.org/) for normal skin melanocytes and melanoma cells, and we compared the protein staining intensity between the groups.

Results

Analysis of gene expression differences between SKCM and normal skin tissues

We identified a total of 6,399 differentially expressed genes, comprising 3,527 upregulated genes and 2,872 downregulated genes and presents the expression levels of the top 10 upregulated and downregulated genes (Supplementary Figure 1).

Identification of centrosome amplification-related genes in SKCM and functional analysis

We intersected the previously identified 6,399 differentially expressed genes with 76 CA-RGs, yielding 17 differentially expressed CA-RGs for further analysis (Supplementary Figure 2A). The results of the GO annotation indicated that these genes were markedly augmented in the centrosomal functions and pathways, including centrosome duplication, centrosome cycle, centriole assembly, microtubule organizing center organization, and tubulin binding (Supplementary Figure 2B-D). As evidenced by the KEGG enrichment analysis, these genes were mostly correlated with pathways of cell cycle (Supplementary Figure 2E). To investigate the protein-protein interaction (PPI) relationships among the differentially expressed CA-RGs, we constructed a PPI network using the STRING database (Supplementary Figure 2F).

Construction and validation of centrosome amplification-related risk models

To construct a risk model for centrosome amplification-related genes (CA-RGs) in SKCM, we performed univariate Cox regression analysis on the TCGA-SKCM cohort using 17 differentially expressed CA-RGs. This identified 7 survival-associated genes (P < 0.05): CDK2, RTTN, KAT2B, PKD2, CCP110, NUBP1, and CEP120. CDK2 and RTTN were risk factors (HR > 1), while KAT2B, PKD2, CCP110, NUBP1, and CEP120 were protective factors (HR < 1) (Figure 1A).  LASSO Cox regression with 10-fold cross-validation refined the model to 6 genes: CDK2, RTTN, KAT2B, PKD2, NUBP1, and CEP120 (Figure 1B). All 6 genes met the proportional hazards assumption (P > 0.05) (Supplementary Table 1). Multivariate and stepwise Cox regression analyses further narrowed the model to 4 prognostic genes: CDK2, KAT2B, NUBP1, and CEP120. Although CEP120 was not significantly associated with survival, all 4 genes were retained as essential features based on stepwise regression. These genes were used to calculate the risk score, forming the final risk model (Figure 1C).

Based on the optimal cutoff value of the risk score (1.340736), patients were stratified into high-risk (n = 374) and low-risk (n = 83) groups. Survival analysis revealed significantly lower survival probabilities in the high-risk group (P < 0.0001) (Figure 2A). ROC curve analysis yielded AUC values of 0.6, 0.61, and 0.63 for 1-year, 3-year, and 5-year survival, respectively (Figure 2B). Risk score distribution and survival status (Figure 2C, D) showed more patients and deaths in the high-risk group. A heatmap of prognostic gene expression (Figure 2E) indicated higher CDK2 expression in the high-risk group, while KAT2B, NUBP1, and CEP120 were downregulated.

To validate the risk model's adaptability and stability, we recalculated risk scores for patients in the GSE100797 cohort. Using the optimal cutoff (-1.674684), patients were stratified into high-risk (n = 13) and low-risk (n = 12) groups. The Kaplan-Meier (KM) curve showed significantly lower survival probabilities in the high-risk group (P = 0.026) ( Figure 2F). ROC curve analysis yielded AUC values of 0.61, 0.67, and 0.67 for 1-year, 3-year, and 5-year survival, respectively (Figure 2G). The risk scores and survival status of patients were visualized in Figure 2H, I. The heatmap confirmed higher CDK2 expression in the high-risk group and lower expression of KAT2B, NUBP1, and CEP120 (Figure 2J).

Construction and evaluation of prognostic models

To evaluate the prognostic capability of the risk score alongside clinical factors, we constructed a prognostic model using the TCGA-SKCM cohort. Differences in risk scores across six clinical factors—age, sex, stage, T stage, N stage, and M stage—were assessed using Wilcoxon or Kruskal-Wallis tests. Significant differences were observed in T stage (P < 0.05) ( Supplementary Figure 3A-F). Survival analysis within T stage subgroups revealed significantly lower survival in the high-risk group of the T3/T4 subgroup compared to the low-risk group (Supplementary Figure 3 G, H).

Next, univariate Cox regression analysis was performed on the risk score, age, sex, stage, T stage, N stage, and M stage in the TCGA-SKCM cohort. Risk score, age, stage, T stage (excluding T2), N stage, and M stage were significantly associated with survival (HR > 1), indicating they are risk factors (Figure 3A). A proportional hazards (PH) assumption test confirmed that risk score, age, stage, N stage, and T stage met the assumption (P > 0.05). Multivariate Cox regression analysis revealed that risk score, N stage, and T stage (excluding T2) were significantly associated with survival and served as risk factors (Figure 3B). These three factors were integrated into a prognostic model to predict 1-year, 3-year, and 5-year survival probabilities. A nomogram was constructed to visualize the model, predicting survival probabilities for a patient with a total score of 129 based on their risk score, T stage, and N stage. The predicted probabilities were 0.968 (1-year), 0.795 (3-year), and 0.675 (5-year) (Figure 3C).

To evaluate the prognostic model’s performance, we generated a calibration curve, ROC curve, and decision curve analysis (DCA). The calibration curve showed close agreement between predicted and observed values (Figure 3D). The ROC curve achieved AUC values of 0.79, 0.79, and 0.75 for 1-year, 3-year, and 5-year predictions, respectively (Figure 3E-G). The DCA curve revealed that the prognostic model yielded the highest net benefit (Figure 3H).

Gene set enrichment analysis and functional enrichment analysis

We conducted independent Gene Set Enrichment Analysis (GSEA) for each gene, identifying the top five significantly enriched pathways. For CDK2, 193 pathways were enriched, with the top five being oxidative phosphorylation, lysosome, complement and coagulation cascades, cytokine-cytokine receptor interaction, and neuroactive ligand-receptor interaction (Supplementary Figure 4A). For KAT2B, 50 pathways were enriched, including autophagy (animal), phosphatidylinositol signaling system, T cell receptor signaling pathway, alcoholic liver disease, and NOD-like receptor signaling pathway (Supplementary Figure 4B). NUBP1 analysis revealed 143 pathways, with the top five being viral protein interaction with cytokine and cytokine receptor, graft-versus-host disease, allograft rejection, hematopoietic cell lineage, and systemic lupus erythematosus (Supplementary Figure 4C). For CEP120, 99 pathways were enriched, including nucleocytoplasmic transport, ubiquitin-mediated proteolysis, Polycomb repressive complex, biosynthesis of amino acids, and Vibrio cholerae infection (Supplementary Figure 4D).

We employed Gene Set Variation Analysis (GSVA) to estimate the biological functions and signaling pathways of prognostic genes, using Hallmarks as the background gene set. We identified 34 differential pathways, with 20 pathways upregulated and 14 pathways downregulated in the high-risk group (Supplementary Figure 4E).

Immune infiltration analysis

To investigate immune cell infiltration in SKCM, we calculated infiltration scores for 34 immune cell types in the TCGA-SKCM cohort using the xCell algorithm. Significant differences were observed in 25 cell types between high-risk and low-risk groups, with 23 showing lower infiltration scores in the high-risk group and 2 showing higher scores (Figure 4A). Spearman correlation analysis revealed predominantly positive correlations among immune cell infiltration scores (Figure 4B), while differential immune cell scores were negatively correlated with risk scores (Figure 4C). Using the ESTIMATE algorithm, we assessed immune cell, stromal cell, and combined scores, finding all three significantly lower in the high-risk group (Figure 4D). These scores also exhibited negative correlations with risk scores (Figure 4E).

Analysis of immune response capacity and immune status

We first analyzed the Immune Prognostic Score (IPS) in the TCGA-SKCM cohort. The results indicated that there was no significant difference in scores between the high-risk and low-risk groups when both CTLA4 and PD1 were suppressed. However, when either one of the immune checkpoints was suppressed or neither was suppressed, the low-risk group exhibited higher response scores (Supplementary Figure 5A). To further clarify the differences in immune checkpoint expression levels, we analyzed the transcriptional expression of 42 immune checkpoints between the high-risk and low-risk groups in the TCGA-SKCM cohort. The findings revealed that 40 immune checkpoints were expressed at lower levels in the high-risk group (Supplementary Figure 5B).

Drug sensitivity analysis

To assess the sensitivity of SKCM patients to commonly used chemotherapeutic agents, we predicted the IC50 values of 198 chemotherapeutic drugs within the TCGA-SKCM cohort and compared the differences between high-risk and low-risk groups. Among the top ten drugs presented, the high-risk group exhibited higher IC50 values (Supplementary Figure 6).

Somatic mutation analysis

We analyzed the mutation status of the top 20 genes in the TCGA-SKCM cohort. The results indicated that 94.86% of patients in the high-risk group (351 individuals) had somatic mutations, while 96.34% of patients in the low-risk group (79 individuals) exhibited somatic mutations. The most common type of mutation in both groups was multi-site mutations, followed by missense mutations. The gene with the highest mutation frequency was TTN, followed by MUC16. Notably, mutations in DNAH5 ranked third in the high-risk group but ranked tenth in the low-risk group (Supplementary Figure 7).

Analysis of prognostic genes regulatory mechanisms

We first utilized the miRNet database to predict the transcription factors targeting the prognostic genes and constructed a transcription factor-prognostic gene interaction network. CDK2 was targeted by 14 transcription factors. The transcription factor FOXC1 was found to target three prognostic genes (Supplementary Figure 8A). Subsequently, we employed the geneMANIA database to investigate the interactions of prognostic genes with other genes. The analysis revealed 20 interacting genes, with NUBP2 demonstrating the strongest interaction (Supplementary Figure 8B).

Single-cell transcriptome sequencing analysis

We utilized the single-cell dataset GSE215120 to explore SKCM pathogenesis. Supplementary Figure 9A displays the number of genes, counts, and mitochondrial proportions quality control, while Supplementary Figure 9B presents a volcano plot of the variable genes. Based on the JackStraw plot and elbow plot (Supplementary Figure 9C), we selected 15 clusters for analysis. Dimensionality reduction and clustering yielded 20 cell clusters, annotated into eight cell types: B cells, endothelial cells, fibroblasts, keratinocytes, malignant cells, monocytes, NK cells, and T cells (Figure 5A). Figure 5B illustrates marker gene expression for each cell type. Statistical analysis revealed malignant cells, T cells, and NK cells as the top three cell types (Figure 5C). Prognostic gene expression analysis identified malignant cells, keratinocytes, and monocytes as key players for pseudotime analysis (Figure 5D).

To investigate receptor-ligand interactions among the eight annotated cell types, we inferred intercellular communication. The network diagram illustrates the quantity and strength of communication networks (Figure 6A). Endothelial cells, fibroblasts, keratinocytes, NK cells, and malignant cells exhibited higher numbers of links, while B cells, T cells, and monocytes had fewer (Figure 6B). We categorized endothelial cells and fibroblasts as signal senders and malignant cells, keratinocytes, and NK cells as signal receivers. The bubble plot (Figure 6C) revealed that communication between these groups primarily involved genes such as COL1A1, COL1A2, COL4A1, COL4A2, COL6A1, COL6A2, and COL6A3, among others.

We performed pseudotime analysis based on cell expression pattern similarities, identifying a major path with two minor branches. This assigned pseudotime-dependent progression states to malignant cells, keratinocytes, and monocytes using Monocle2-reconstructed trajectories (Figure 6D). Prognostic gene expression patterns along the melanoma progression trajectory were analyzed, revealing two distinct clusters in the pseudotime heatmap (Figure 6E). The pseudo-time plot reveal that the expression levels of the prognostic genes CEP120 and NUBP1 remain relatively stable throughout the cell development process, with NUBP1 increasing only at the end. In contrast, CDK2 expression gradually rises during the coexistence of malignant cells, monocytes, and keratinocytes, accelerating when only malignant cells are present. KAT2B shows a slow increase in expression during the phase of coexistence among the three cell types, followed by a decline after monocytes completely transition. However, KAT2B expression rises sharply when only malignant cells remain. In summary, during the progression of melanoma, the changes in expression levels of CDK2 and KAT2B primarily regulate the transition of monocytes and keratinocytes into malignant cells. (Figure 6F).

Analysis of prognostic genes expression

We found that all four prognostic genes were upregulated in tumor tissues in the GSE15605 dataset (Supplementary Figure 10A). Subsequently, we utilized the Human Protein Atlas (HPA) database to download immunofluorescence images of prognostic genes in normal skin melanocytes and melanoma tissues. The results showed that CDK2 and NUBP1 exhibited high expression at the protein level in SKCM, while CEP120 did not display detectable fluorescence signals. KAT2B displayed high fluorescence signal; however, it was apparent that its expression was more intense in tumor tissues. In summary, the expression patterns of CDK2 and NUBP1 are consistent at both transcriptomic and protein levels (Supplementary Figure 10B).

Discussion

CA has been shown to be involved in the development and progression of human cancers and is also a valuable prognostic factor in various cancer types[7] . However, there are relatively few studies on centrosome amplification in SKCM. Therefore, we aimed to evaluate the effect of CA-RGs on cancer development, prognosis prediction, tumor immune microenvironment changes and treatment response in SKCM.

We identified 17 genes by intersecting normal-tumor differential genes from the GSE15605 database with 76 centrosome amplification-related genes (CA-RGs) from the Gene Ontology database, which were enriched in centrosomal functions and pathways. To refine CA-RGs in melanoma and investigate their association with prognosis in SKCM, we constructed a molecular signature based on 4 prognostic genes using Cox regression analysis. A CA-RGs risk score was calculated for each SKCM patient, and patients were stratified into high- and low-risk subgroups based on an optimal cutoff. The high-risk group showed a significantly worse overall survival (OS), consistent across stratified survival analyses in GEO cohorts. Risk scores were significantly associated with clinical T stages, with higher scores correlating with poorer survival in the T3/T4 stages. Univariate and multivariate Cox analyses confirmed the CA-related signature as an independent prognostic factor. To enhance its predictive performance and clinical applicability, a nomogram integrating the CA-related signature and two additional prognostic factors was developed. Calibration curves, multi-index ROC curves, and decision curve analysis (DCA) demonstrated its superior efficacy in predicting OS compared to traditional clinical features such as age and histological stage.

The most significant finding from this study is the identification of CDK2 as a key prognostic factor. Cyclin-dependent kinase 2 (CDK2) plays a critical role in the cell cycle by regulating the transition from the G1 phase to the S phase[36] , and is required for centrosome duplication[37] . Matsuura found that the phosphorylation of SMAD3 by CDK2–cyclin-E limits its transcriptional activity and eventually slows cell-cycle progression[38] . Our gene set enrichment analysis (GSEA) revealed that CDK2 is involved in oxidative phosphorylation and immune-related pathways as well, indicating its dual role in metabolic regulation and immune evasion[39] during melanoma progression. Lysine acetyltransferase 2B (KAT2B) ,also named as PCAF, plays a crucial role in the regulation of gene expression at the post-transcriptional level by acetylation[40] , and is associated with many types of cancer. It has been reported to be involved in regulating the acetylation and stability of HMGA2 to accelerate the growth of esophageal squamous cell carcinoma[41] . Moreover, induction of PGK1 enzymatic activity and cancer cell metabolism by KAT2B engagement-mediated acetylation plays an important role in liver cancer progression[42] . On the other side, one outcome demonstrated that PCAF blocked the growth of hepatocellular carcinoma via enhancing cell autophagy[43] . Li has found KAT2B degradation contributes to the development of pancreatic cancer via the Ras-ERK1/2 signaling pathway[44] . Our results show that KAT2B is a protective prognostic gene analyzed by Cox regression analysis and is engaged in autophagy, which is consistent with a previous study. In our analysis of the prognostic gene interaction network, we found that NUBP1 and NUBP2 exhibit the closest interaction. This is likely because both proteins are crucial for the synthesis and transport of intracellular iron-sulfur clusters, and they function in a coordinated manner to fulfill their roles in these processes[45] . NUBP1, NUBP2 and KIFC5A have been implicated in the regulation of centriole duplication with their depletion causing supernumerary centrioles[46] . This function is critical in maintaining genomic stability, as excessive CA is known to promote tumorigenesis by driving chromosomal instability[47] . The protective nature of NUBP1 identified in our study suggests that restoring its function could be a potential therapeutic strategy to counteract CA-induced tumor progression in SKCM. CEP120, although not significantly associated with survival, was included in our final risk model based on its role in centrosome biogenesis and cilia formation[48] . While its exact function in melanoma remains less clear, it has been implicated in centriole elongation and centrosomal microtubule organization[49] , both of which are essential for proper cell division . Further research is warranted to elucidate whether CEP120 plays a contributory role in SKCM progression through its involvement in CA or other mechanisms of tumor growth. Our results of GSVA analysis indicate that melanoma development is associated with multiple signaling pathways. For instance, the MYC_TARGETS_V1 and MYC_TARGETS_V2 pathways are linked to MYC, a transcription factor that regulates cell proliferation, growth, and metabolism, often overexpressed to promote tumor cell proliferation and survival[50] . Additionally, the INTERFERON_GAMMA_RESPONSE pathway involves interferon-gamma (IFN-γ), a crucial immune regulator that plays a key role in modulating immune responses, particularly in antiviral and antitumor reactions[51] . The downregulation of the IFN-γ-related pathway in melanoma indicates a reduced antitumor capacity, facilitating tumor progression.

Beyond the effects on tumor growth, CA-RGs also affect the composition of the tumor immune microenvironment. Tumor-associated macrophages (TAM) can be divided into two major subtypes: M1 TAMs counteract tumor progression while M2 TAMs promote tumor growth and inhibit tumor immunity[52] . In this study, we found that M1 macrophages infiltrated the low CA-related signature risk group more than the high CA-related signature risk group. The proliferation and cytotoxicity CD8+ T cells and their secretion of IFN-γ induce an antitumor response[53] . Tumor regression after ICB requires abundant infiltration of CD8+ T cells in proximity to tumor cells[54] . Consistent with this, our results show that the low-risk group harbors more CD8+ cytotoxic T lymphocytes. This suggests that the tumor microenvironment in the low-risk group may have stronger antitumor immune activity, while the high-risk group tends toward an immunosuppressive state. ESTIMATE analysis indicates that the low-risk group has significantly higher ESTIMATE, immune, and stromal scores compared to the high-risk group. This finding suggests that the low-risk group’s microenvironment is richer in immune and stromal components, potentially enhancing antitumor immune responses and contributing to better prognosis. Furthermore, the correlation analysis shows a negative correlation between the risk score and the ESTIMATE, immune, and stromal scores, indicating that as the risk score increases, immune and stromal component levels in the tumor microenvironment decrease, which indicates that a higher risk score may be associated with a microenvironment deficient in immune activity and stromal support, reflective of a more aggressive and immune-evasive tumor profile. To date, immunotherapeutic strategies have manifested curative therapeutic effects in unresectable or metastatic melanoma[55] . Unfortunately, the clinical benefits remain unsatisfactory owing to low objective remission rates and resistance to ICB in a substantial fraction of melanoma patients. Investigating a predictive marker for benefit from immunotherapy is thus urgently needed. Cytotoxic T-lymphocyte associated protein-4 (CTLA-4) and the Programmed Death Receptor 1 (PD-1) are immune checkpoint molecules that are well-established targets of antibody immunotherapies for the management of malignant melanoma[4] . In our study, across all checkpoint expression statuses, low-risk scores consistently associate with higher IPS scores in samples where at least one immune checkpoint is expressed. This indicates that high-risk score group may enhance immune evasion, while low-risk score group exhibit greater immune activity and higher potential responsiveness to immunotherapy. These findings support a stratified approach to immunotherapy, where patients with low-risk score may be more suitable candidates for PD1- or CTLA4-based checkpoint inhibition. We also evaluated the sensitivity of high-risk and low-risk patients to chemotherapeutics and targeted drugs to better guide clinical medication for SKCM patients and the high-risk group exhibited higher IC50 values, indicating that patients in the low-risk group have greater sensitivity to chemotherapy.

At the single-cell level, malignant cells, monocytes, and keratinocytes emerge as key cell types in melanoma progression. This process relies primarily on communication from fibroblasts to keratinocytes and malignant cells. Under the regulatory influence of prognostic genes CDK2 and KAT2B, this signaling promotes the transition of keratinocytes and monocytes toward malignant cell phenotypes.

Conclusion

Although our study has significant clinical implications for survival prediction and therapeutic decision-making, several limitations must be acknowledged. First, we identified four novel prognostic genes related to centrosome amplification and constructed high-performance risk and prognostic models. These models offer new, reliable tools for clinical survival prediction and patient risk assessment and provide potential therapeutic targets for melanoma, including drug and immunotherapy interventions. Furthermore, our findings reveal melanoma progression mechanisms and their relationships with prognostic genes at a higher-resolution level. However, given that this study is based on data analysis, there are inherent limitations. The molecular mechanisms through which these CA-related signature biomarkers contribute to disease progression and therapy resistance in patients with SKCM should be further investigated through in experimental studies.

Abbreviations

CA: Centrosome Amplification; CA-RGs: Centrosome Amplification-Related Genes; PCM: Pericentriolar Material; DEGs: Differentially Expressed Genes; GO: Gene Ontology; KEGG: Kyoto Encyclopedia of Genes and Genomes; PPI: Protein-Protein Interaction; TCGA: The Cancer Genome Atlas; ROC: Receiver Operating Characteristic; AUC: the Area Under the Curve; GSVA: Gene Set Variation Analysis; IPS: Immune Phenotype Scores; TME: Tumor Microenvironment; HPA: Human Protein Atlas

Supplementary Material

Supplementary methods, results, spectra, figures.

Acknowledgements

All authors would like to thank the TCGA, GEO datasets, which provided data for this study.

Author contributions

Wendong Chen contributed to the conceptualization, data curation ,methodology and software, Dawei Wang contributed to the formal analysis and investigation, Zhenyu Wu contributed to the validation and visualization, Shengrong Cheng contributed to supervision, Fei Zhu contributed to project administration and resources. Wendong Chen, Dawei Wang, Zhenyu Wu and Shengrong Cheng contributed to writing the original draft. Hailin Yao and Fei Zhu contributed to reviewing and editing.

Funding information

Not applicable.

Ethics approval and consent to participate

Not applicable.

Competing Interests

The authors declare that they have no existing or potential commercial or financial relationships that could create a conflict of interest at the time of conducting this study.

Data Availability

All data needed to evaluate the conclusions in the paper are present in the paper or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

References