Abstract

   Melanoma is an aggressive type of skin cancer that arises from
   melanocytes, the cells responsible for producing skin pigment. In
   contrast to non-melanoma skin cancers like basal cell carcinoma and
   squamous cell carcinoma, melanoma is more invasive. Melanoma was
   distinguished by its rapid progression, high metastatic potential, and
   significant resistance to conventional therapies. Although it accounted
   for a small proportion of skin cancer cases, melanoma accounts for the
   majority of deaths caused by skin cancer due to its ability to invade
   deep tissues, adapt to diverse microenvironments, and evade immune
   responses. These unique features highlighted the challenges of treating
   melanoma and underscored the importance of advanced tools, such as
   single-cell sequencing, to unravel its biology and develop personalized
   therapeutic strategies. Thus, we conducted a single-cell analysis of
   the cellular composition within melanoma tumor tissues and further
   subdivided melanoma cells into subpopulations. Through analyzing
   metabolic pathways, stemness genes, and transcription factors (TFs)
   among cells in different phases (G1, G2/M, and S) as well as between
   primary and metastatic foci cells, we investigated the specific
   mechanisms underlying melanoma metastasis. We also revisited the
   cellular stemness and temporal trajectories of melanoma cell
   subpopulations, identifying the core subpopulation as C0
   SOD3 + Melanoma cells. Our findings revealed a close relationship
   between the pivotal C0 SOD3 + Melanoma cells subpopulation and
   oxidative pathways in metastatic tumor tissues. Additionally, we
   analyzed prognostically relevant differentially expressed genes (DEGs)
   within the C0 SOD3 + Melanoma cells subpopulation and built a
   predictive model associated with melanoma outcomes. We selected the
   gene IGF1 with the highest coefficient (coef) value for further
   analysis, and experimentally validated its essential function in the
   proliferation and invasive metastasis of melanoma. In immune
   infiltration analysis, we discovered the critical roles played by M1/M2
   macrophages in melanoma progression and immune evasion. Furthermore,
   the development and progression of malignant melanoma were closely
   associated with various forms of programmed cell death (PCD), including
   apoptosis, autophagic cell death, ferroptosis, and pyroptosis. Melanoma
   cells often resisted cell death mechanisms, maintaining their growth by
   inhibiting apoptosis and evading autophagic cell death. Meanwhile, the
   induction of ferroptosis and pyroptosis was thought to trigger immune
   responses that helped suppress melanoma dissemination. A deeper
   understanding of the relationship between melanoma and PCD pathways
   provided a critical foundation for developing novel targeted therapies,
   with the potential to enhance melanoma treatment efficacy. These
   findings contributed to the development of novel prognostic models for
   melanoma and shed light on research directions concerning melanoma
   metastasis mechanisms and therapeutic targets.

   Keywords: Melanoma, Tumor metastasis, Oxidative phosphorylation,
   Macrophages, Single-cell RNA sequencing, Programmed cell death

Introduction

   Melanoma, a malignancy originating from melanocytes, stands as one of
   the most prevalent and severe form of cutaneous carcinoma
   [[36]1–[37]3]. Often triggered by exposure to ultraviolet radiation,
   its incidence has demonstrated a persistent increase in recent years,
   with an estimated global annual incidence of 300,000 cases and
   mortality rates reaching up to 60,000, imposing substantial burdens on
   both public health and healthcare systems [[38]3, [39]4]. Presently,
   surgical excision remains the cornerstone of melanoma management,
   exhibiting favorable prognostic outcomes in the absence of metastatic
   spread [[40]5]. However, unfortunately, melanoma frequently progressed
   to form aggressive metastatic tumor tissues, leading to the demise of a
   majority of afflicted individuals [[41]6]. Patients with metastatic
   disease typically faced a median survival duration of merely
   6–9 months, with a grim 5-year survival rate below 5% [[42]7],
   indicative of the limited efficacy of surgical interventions at this
   stage, which are also prone to relapse [[43]8]. Tumor invasion and
   metastasis served as critical determinants of disease progression,
   recurrence, and mortality in melanoma patients [[44]9]. Hence, there is
   a pressing need to delve deeper into the biological underpinnings of
   melanoma invasion and metastasis, as well as to identify pertinent
   pathways or immunotherapeutic targets, in order to ameliorate the
   prognosis of individuals afflicted with melanoma.

   Macrophages serve as effector cells in regulating the body's immune
   response, exhibiting high plasticity and categorizing into M1
   macrophages, which promote inflammation progression, and M2
   macrophages, which facilitate tissue repair. They participate in
   shaping the tumor microenvironment (TME) and play crucial roles in the
   initiation, proliferation, and metastasis of tumors [[45]10].
   Specifically, M1 macrophages are involved in inhibiting tumor growth,
   whereas M2 macrophages contribute to its progression [[46]11, [47]12].
   The redox balance is fundamental for cell survival and normal
   functioning, with oxidative processes in tumor cells often surpassing
   those in normal cells [[48]13, [49]14]. Oxidative phosphorylation, a
   pivotal pathway in this context, generates reactive oxygen species
   (ROS), during metabolism, which pose a risk of damaging nuclear and
   mitochondrial DNA, thereby increasing cancer susceptibility [[50]15,
   [51]16]. Furthermore, oxidative phosphorylation itself plays a certain
   promoting role in tumor metastasis [[52]17]. However, the specific
   mechanisms of macrophages and oxidative phosphorylation in the process
   of melanoma metastasis remain relatively understudied.

   Based on this, our study investigated the functional relevance and
   specific mechanisms of various melanoma cell subpopulations at the
   cellular level. Furthermore, we examined the correlation between the
   primary and metastatic tumor tissues of Melanoma cells, considering
   cell staging and metabolism. Through analysis of cellular stemness and
   pseudo-temporal trajectories, we identified the crucial subpopulation,
   and designed a prognostic model connected to them. In addition, we
   explored the interplay with immune function. These investigations aimed
   to provide a deeper insight into the mechanisms driving melanoma
   metastasis, with the ultimate goal of designing effective and specific
   immune therapies for melanoma and developing prognostic models.

Materials and methods

Data download and processing

   The single-cell sample data from 8 melanoma patients, including primary
   and metastatic tumor tissues, were sourced from the NCBI Gene
   Expression Omnibus (GEO) database (Accession No. [53]GSE189889).
   RNA-Seq data and clinical information of melanoma cells were retrieved
   from The Cancer Genome Atlas (TCGA) database. The data were normalized
   using R software (version 4.3.0). Since the data were sourced from
   public databases, no ethical review was required.

Quality control

   The R package “DoubletFinder” [[54]18, [55]19] was used to purify,
   filter, and remove doublets and other low-quality data from the
   single-cell data obtained from the GEO database, which included 9
   melanoma tumor tissue samples. The selection criteria were as follows:
   (1) 200 < number of expressed genes per cell (nFeature) < 6000; (2)
   500 < total number of transcriptional genes per cell (nCount) < 80,000;
   (3) proportion of red blood cell genes < 5%; (4) proportion of
   mitochondrial genes < 25%.

Dimensionality reduction, clustering, and cell annotation

   Single-cell data were normalized through the “NormalizeData” function
   in Seurat [[56]20–[57]22]. Identified highly variable genes with the
   Find Variable Features function, focusing on the top 2000 based on
   variance and standard deviation [[58]23, [59]24]. Conducted gene
   centering and standardization using the “ScaleData” function, and
   evaluated cell cycle effects with the Cell Cycle Features function.
   Performed principal component analysis (PCA) on the standardized highly
   variable genes with the PCA function. Then, addressed batch effects
   using the Harmony algorithm and executed clustering with
   “FindNeighbors” and “FindClusters” [[60]25]. Obtained cell annotations
   for different subpopulations from the Cell Marker database, refining
   based on relevant literature [[61]26, [62]27].

Cell annotation and tumor tissue preference analysis

   Then, we reclustered tumor cells and identified subpopulations based on
   marker genes. Differential expression analysis between subpopulations
   and between primary and metastatic samples was conducted using the
   “FindAllMarkers” function [[63]28, [64]29]. To assess the cancer
   preference of melanoma cell subpopulations, odds ratios were calculated
   using the calculation method.

Subpopulation cell activity and metabolism analysis

   The “scMetabolism” package in the R package was used to screen and
   analyze the metabolic pathways of different phases, primary and
   metastatic tumor tissues and subpopulations. The metabolic pathways of
   enrichment were visualized using violin and UMAP plots.

Subpopulation enrichment analysis

   We identified differentially expressed genes (DEGs) [[65]30] among
   subpopulation cells, with criteria including detection in at least 25%
   of cells, false discovery rate (FDR) < 0.05, and | logFCfilter |> 1.
   The DEGs before and after filtering were displayed for each
   subpopulation. Based on the DEGs among subpopulation cells, we
   performed Gene Ontology (GO) [[66]31–[67]33] enrichment analysis for
   each subpopulation by using “Cluserprofiler” package.

Cell trajectory analysis

   Utilized the R package CytoTRACE to predict the differentiation states
   of melanoma subpopulations. Combined this with slingshot (v2.6.0) and
   monocle2 to infer cellular lineages and pseudotemporal variations
   [[68]34–[69]36]. Integrated synchronous principal and branching curves
   using the “getCurves” function to extract smoothed trajectory curves.

Subpopulation cell stemness and transcription factors (TFs)

   The AUCell method was used to identify cells with active gene
   expression by analyzing single-cell RNA-seq data. This approach
   evaluated gene expression profiles by taking gene sets as input and
   assessing the activity of each gene set within individual cells. In
   this research, AUCell was specifically employed to assess the stemness
   characteristics of different tumor cell subpopulations.

   We utilized pySCENIC (v0.10.0) in Python 3.7 to infer regulatory
   networks at the single-cell level and to cluster tumor cell
   subpopulations [[70]37]. Target genes of TFs were identified using
   “GRNBoost” and further refined through DNA motif analysis. The AUCell
   method evaluated regulon activities, identifying the top 5 TFs with the
   significant changes in expression.

Clinical implications and prognostic evaluation of the C0 SOD3 + melanoma
cells subpopulation

   Developed predictive models for prognostication in patients with
   melanoma [[71]38]. Excluded incomplete clinical data and merged the
   intersection genes with the standardized clinical dataset. Conducted
   single-variable Cox regression analysis using the “coxph” function from
   the survival package in R, and validated the results with least
   absolute shrinkage and selection operator (LASSO)-penalized Cox
   regression [[72]39]. Identified prognostic differential genes through
   multivariable Cox regression analysis.

   Calculated the risk score for each sample (risk score = ΣXλ × coefλ,
   where X represents the relative expression levels of prognostic-related
   genes and coef represents the coefficients and stratified samples into
   high-risk and low-risk groups based on the median score. Employed PCA
   to examine sample distribution, and visualized the expression profiles
   of prognostic genes in each risk group using a heatmap. Highlighted
   coef associated with DEGs for further experimental validation.
   Illustrated the survival status of high-risk and low-risk groups over
   time using Kaplan–Meier curves [[73]40–[74]42], and assessed the
   specificity and sensitivity with time-dependent receiver operating
   characteristic (ROC) curves 's area under the curve (AUC)
   [[75]43–[76]45].

   Explored the link between risk scores and prognostic gene expression,
   and incorporated age and melanoma staging information to refine the
   prognostic model, modeling it with the “rms” package in R [[77]46].
   Developed a nomogram chart by combining the risk score with clinically
   relevant factors, serving as a predictive tool for melanoma patient
   prognosis. Validated the accuracy using ROC curves.

Immunological and enrichment analysis

   Utilized the “Xcell” and “CIBERSORT” algorithms to analyze immune cell
   infiltration in high-risk and low-risk groups. Investigated
   correlations between risk scores, prognostic differential genes, and
   immune cells. Employed the R package “ESTIMATE” to calculate immune
   scores, stromal scores, and overall scores for the TME, compared Tumor
   Immune Dysfunction and Exclusion (TIDE) scores between the groups, and
   visualized results with violin plots. Filtered DEGs using the “limma”
   package with criteria of |log2 Fold Change|> 1 and FDR < 0.05.
   Conducted Gene Set Enrichment Analysis (GSEA) on the filtered genes
   [[78]47].

Cell culture

   The SK-MEL2 and A375 cell lines were sourced from the American Type
   Culture Collection (ATCC) and cultured in F12K and PRMI1640 media
   (Gibco BRL, USA), supplemented with 10% fetal bovine serum (Gibco BRL,
   USA) and 1% penicillin/streptomycin, under standard conditions (37 °C,
   5% CO2, 95% humidity).

Cell transfection

   Small interfering RNA (siRNA) constructs targeting IGF1 knockdown were
   obtained from GenePharma (Suzhou, China). Transfection was performed
   according to the Lipofectamine 3000 RNAiMAX protocol (Invitrogen, USA).
   Cells were seeded at 50% confluence in 6-well plates and transfected
   with negative control (si-NC) or knockdown siRNAs (Si-IGF1-1 and
   Si-IGF1-2). Each transfection utilized Lipofectamine 3000 RNAiMAX
   reagent.

Cell viability assay

   Cell viability of transfected SK-MEL2 and A375 cells was evaluated
   using the CCK-8 assay. Cells were seeded at 5 × 10^3 per well in a
   96-well plate and incubated for 24 h. Afterward, 10 μL of CCK-8 reagent
   (A311-01, Vazyme) was added per well and incubated for 2 h at 37 °C in
   the dark [[79]48]. Absorbance was measured at 450 nm on days 1, 2, 3,
   and 4 using a microplate reader (A33978, Thermo), and the average OD
   values were plotted.

Western blot analysis

   Cells were cultured and transfected until reaching 70% confluence. Cell
   lysates were prepared using RIPA buffer, clarified by centrifugation at
   12,000 rpm for 15 min, and subjected to SDS-PAGE. Proteins were
   transferred onto PVDF membranes, blocked with 5% BSA at room
   temperature for 1.5 h, incubated with primary antibodies overnight at
   4 °C, followed by corresponding secondary antibodies for 1 h. Protein
   bands were visualized using ECL Western Blot [[80]49] Substrate.

Quantitative real-time polymerase chain reaction (qRt-PCR)

   The process of RNA extraction was performed using Trizol reagent,
   followed by reverse transcription with the PrimeScript™ Kit. qRt-PCR
   was then conducted with SYBR Green premix for detecting amplification
   [[81]50].

EdU proliferation assay

   Transfected SK-MEL2 and A375 cells were seeded at a density of 5 × 10^3
   cells per well in 6-well plates and cultured overnight. Cells were
   treated with 10 mM EdU solution in serum-free medium for 2 h at 37 °C.
   After fixation with 4% paraformaldehyde and permeabilization with
   glycine (2 mg/ml) and 0.5% Triton X-100, cells were incubated with 1X
   Apollo and 1X Hoechst 33,342 for 30 min at room temperature. Cell
   proliferation was quantitatively analyzed under a fluorescence
   microscope.

Wound healing assay

   Transfected SK-MEL2 and A375 cells were grown in 6-well plates until
   95% confluence. A straight scratch was made using a sterile pipette
   tip, followed by washing with PBS. The medium was replaced, and cells
   were incubated further. Images of the scratch were taken at 0 and 48 h
   to measure wound width.

Transwell assay

   Prior to the experiment, cells were serum-starved for 24 h in
   serum-free medium. After coating with Matrigel (BD Biosciences, USA),
   cell suspensions were added to the upper chambers of Costar Transwell
   plates, while serum-containing medium was added to the lower chambers.
   Cells were incubated in a humidified chamber for 48 h. Following
   incubation, cells were fixed with 4% paraformaldehyde, stained with
   crystal violet, and assessed for invasive capability.

Statistical analysis

   R and Python were employed for database data analysis, while
   experimental data were analyzed with GraphPad Prism version 8.0.1.
   Two-tailed p-values were applied in all analyses, with statistical
   significance determined as p-values less than 0.05.

Results

Single-cell analysis of melanoma tumor tissues

   Figure [82]1 illustrated the overall workflow of this study. At the
   beginning, we obtained single-cell data from 9 melanoma samples
   ([83]GSM5708993-[84]GSM5709001) representing tumors from 8 patients,
   including 5 primary and 4 metastatic tumor tissues. Employing rigorous
   data filtration and dimensionality reduction clustering techniques, we
   subdivided the melanoma tumor tissue into 34 subpopulations
   (Fig. [85]2A). Investigation into the origins and metastatic
   attributions of these cells revealed heterogeneous distribution
   patterns across the subpopulations, indicating diverse cellular
   compositions within the TME. The study mainly focused on primary acral
   melanoma tissues (AM2, AM3, AM4, AM6, AM8-Primary) and metastatic acral
   melanoma tissues (AM1, AM5, AM7, AM8-Mets, Mets indicated Metastatic).
   Notably, two types of samples were extracted from a single patient
   (Fig. [86]2B). Subsequent annotation and classification of the 34
   clusters, delineated 9 distinct cellular populations within the
   melanoma tissue (Fig. [87]2C, D). These populations comprised Melanoma
   cells (18,083), NK-T cells (10,548), Endothelial cells (ECs) (1418),
   B-Plasma cells (1295), Myeloid cells (1535), Fibroblasts (1651),
   Pericytes (552), Epithelial cells (EPCs) (225), and Mast cells (79).
   Notably, melanoma cells constituted the largest cluster, serving as the
   predominant cellular entity within the melanoma TME, followed by NK-T
   cells, indicative of their immunological relevance. To investigate
   cellular migratory capabilities, we analyzed the distribution and
   proportions of cellular origins in primary and metastatic tumor tissues
   across the major cell populations (Fig. [88]3A–C). Cell cycle
   distribution (Fig. [89]3D–F) revealed a prominent G2/M phase
   predominance in EPCs, indicative of heightened proliferative activity.
   These findings underscored the necessity for further subtyping of
   melanoma cells to elucidate mechanisms underlying tumor progression and
   metastasis.

Fig. 1.

   [90]Fig. 1
   [91]Open in a new tab

   Workflow of single cell sequencing analysis for melanoma. We used
   single-cell sequencing technology, combined with bulk sequencing, in
   vitro animal experiments, and a series of analyses to comprehensively
   investigate melanoma. The aim was to identify novel therapeutic targets
   and improve patient prognosis

Fig. 2.

   [92]Fig. 2
   [93]Open in a new tab

   Clustering of melanoma tumor tissues. A Following single-cell
   clustering, a total of 34 clusters were identified from tumors and
   surrounding tissues of 8 melanoma patients. B These samples originated
   from 9 specimens of 8 patients. It could be divided into 5 primary
   acral melanoma tissue and 4 metastatic acral melanoma tissues. C, D
   Annotation of different cell populations based on distinct cell surface
   markers, including Melanoma cells, NK-T cells, ECs, B-Plasma cells,
   Myeloid cells, Fibroblasts, Pericytes, EPCs, and Mast cells

Fig. 3.

   [94]Fig. 3
   [95]Open in a new tab

   The distribution of melanoma cell populations and the cell cycle
   phases. A–C The UMAP plots illustrated the spatial distribution of
   different cells within primary and metastatic melanoma tissues. Pie
   charts were used to depict the proportions of primary and metastatic
   melanoma tissues from different groups within each cell type. D–F The
   UMAP plots illustrated the spatial distribution of various cells within
   primary and metastatic melanoma tissues. Pie charts were used to
   describe the proportions of primary and metastatic melanoma tissues
   from different cell types across different cell cycle phases (G1, G2/M,
   S)

Subpopulation and primary-metastasis analysis of melanoma cells

   Further exploration was done to understand the relationship between
   heterogeneity and subpopulations in melanoma patients. Through further
   clustering of melanoma cells, we subdivided them into 8 subpopulations
   (Fig. [96]4A). Investigation into the cellular origins of these
   subpopulations revealed their derivation (Fig. [97]4B). According to
   the expression levels of cell-specific genes, we reclassified them into
   8 subpopulations, including C0 SOD3 + Melanoma cells, C1
   MTRNR2L12 + Melanoma cells, C2 PMEL + Melanoma cells, C3
   IGF2 + Melanoma cells, C4 COMP + Melanoma cells, C5 CXCL9 + Melanoma
   cells, C6 THBS2 + Melanoma cells, C7 PECAM1 + Melanoma cells
   (Fig. [98]4C). To delve deeper into distribution of subpopulations and
   metastatic relationships of each subpopulation, we evaluated the cell
   cycle phases and primary metastatic status of each subpopulation.
   Remarkably, the majority of cells in the C0 SOD3 + Melanoma cells
   subpopulation was found to be in the G1 phase (Fig. [99]4D–F). In
   addition, the proportion of metastatic tumor tissues in the C0
   SOD3 + Melanoma cells subpopulation was relative higher
   (Fig. [100]4G–I).

Fig. 4.

   [101]Fig. 4
   [102]Open in a new tab

   Comprehensive distribution landscape of melanoma cell subpopulations. A
   The melanoma cells underwent dimensionality reduction and clustering,
   resulting in the identification of 8 distinct subpopulations. B The 8
   subpopulations were derived from different melanoma tissues, primarily
   categorized into primary and metastatic tissues. C The melanoma cells
   were distinguished into 8 subpopulations based on the different
   characteristic genes. D–F The UMAP plots illustrated the distribution
   characteristics of the 8 subpopulations across different cell cycle
   phases. Pie charts depicted the proportions of each phase within the
   various subpopulations. G–I The UMAP plots showed the distribution
   characteristics of the 8 subpopulations based on their tissue origins.
   Pie charts displayed the proportions of primary and metastatic melanoma
   tissues within the different subpopulations

   To further investigate the interplay between cell cycle phases and
   primary metastasis among different subpopulations, we visualized the
   expression of key genes across all subpopulations (Fig. [103]5A).
   Strikingly, the key genes expression profile of the C0 SOD3 + Melanoma
   cells subpopulation closely mirrored that of metastatic tumor tissues,
   suggesting an association between the C0 SOD3 + Melanoma cells
   subpopulation and tumor metastasis. Analysis of the cell cycle in each
   subpopulations revealed that the C0 SOD3 + Melanoma cells subpopulation
   showed distinct differences compared to other subpopulations,
   exhibiting a considerably higher proportion of cells in the G1 phase,
   reaching 50.7%, while other subpopulations showed more uniform
   distributions across G2/M and S phases (Fig. [104]5B, C). Furthermore,
   our analysis of the cell proportions involved in primary and metastatic
   tissues across different subpopulations indicated that within the C0
   SOD3 + Melanoma cells subpopulation, metastatic tumor tissue cells
   represented the higher proportion, accounting for 78.1% (Fig. [105]5D,
   E). Further analysis of Ro/e preference analysis across subpopulations
   revealed that within all subpopulations, the C0 SOD3 + Melanoma cells
   subpopulation exhibited the highest Ro/e value of metastatic tumor
   tissue cells and G1 phase cells (Fig. [106]5F, G). Hence, we
   hypothesized that certain life activities or metabolic pathways of the
   C0 SOD3 + Melanoma cells subpopulation during the G1 phase could be
   associated with tumor metastasis, potentially influencing other
   subpopulations through specific metabolic pathways to promote tumor
   proliferation and metastasis.

Fig. 5.

   [107]Fig. 5
   [108]Open in a new tab

   Phase and distribution of subpopulations within primary and metastatic
   tumor tissues. A The bubble diagram showed the expression levels of
   different marker genes across the 8 subpopulations and various tissues.
   Pie charts illustrated the cell cycle proportions for different cell
   subpopulations and tissues, while violin plots presented the G2/M
   Score, S Score, and nCount-RNA levels for different cell subpopulations
   and tissues. B, C The bar graphs showed the proportion of each
   subpopulations in the different phases. D, E The bar graphs showed the
   proportion of each subpopulations in the primary and metastatic groups.
   The proportion of primary and metastatic melanoma tissues within
   subpopulations C0 to C7. F Ro/e values of 8 melanoma cell
   subpopulations in metastatic and primary tumor tissues G Ro/e values of
   8 melanoma cell subpopulations were assessed across different cell
   cycle phases (G1, G2/M, and S phases)

Metabolic pathways and activity analysis of subpopulations

   To better understand the mechanisms of melanoma cell subpopulations in
   the context of staging and metastasis, we assessed the metabolic
   pathways across distinct subpopulations, cell cycle phases, as well as
   primary and metastatic tumor tissues (Fig. [109]6A–C). We observed that
   the pathways of oxidative phosphorylation, glutathione metabolism, and
   glycolysis/gluconeogenesis consistently ranked high across different
   groups. To further explore and compare the scores of these pathways
   across different subpopulations, cell cycle phases, and primary and
   metastatic tumor tissue cells, we visualized them using UMAP and violin
   plots (Fig. [110]6D–O). Notably, there was a significant upregulation
   in the expression of these three metabolic pathways in the C0
   SOD3 + Melanoma cells subpopulation. In addition, we found that the
   above three metabolic pathways exhibited higher expression levels in
   the G1 phase and metastatic tumor tissues. This suggested a potential
   correlation between these pathways and the activity and metastasis of
   melanoma tumor cells. Interestingly, the expression levels nFeature-RNA
   and nCount-RNA were high in the C0 SOD3 + Melanoma cells subpopulation.
   Metastatic tumor tissue cells and those in the G1 phase also exhibited
   significantly higher expression compared to other cell types and cell
   cycle phases (Fig. [111]6P–W). We hypothesized that the C0
   SOD3 + Melanoma cells subpopulation might influence the activity and
   metastasis of tumor cells through oxidative-related pathways and
   potentially promote tumor proliferation by affecting other
   subpopulations' cells via differentiation or signaling pathways.

Fig. 6.

   [112]Fig. 6
   [113]Open in a new tab

   Cell metabolic pathway analysis landscape. A–C Top 20 metabolic
   pathways analyzed across different subpopulations, phases, as well as
   primary and metastatic tumor tissues. D–O UMAP and violin plots
   depicted the distribution and scores of oxidative phosphorylation,
   glutathione metabolism, glycolysis/gluconeogenesis pathways in
   different subpopulations, phases, as well as primary and metastatic
   tumor tissues. P–S UMAP and violin plots showed nCount-RNA distribution
   and scores across different subpopulations, phases, as well as primary
   and metastatic tumor tissues. T–W UMAP and Violin plots displayed
   nFeature-RNA distribution and scores across different subpopulations,
   phases, as well as primary and metastatic tumor tissues. (*P < 0.05,
   **P < 0.01, ***P < 0.001, ****P < 0.0001, "ns" was used to say that
   there was no significant difference)

Enrichment analysis of subpopulations

   Then, we computed the DEGs between each subpopulation of melanoma cells
   and showcased the top 5 upregulated and downregulated genes for each
   comparison (Fig. [114]7A–H). Subsequently, we proceeded with GO
   enrichment analysis on these DEGs (Fig. [115]7I–P). We discovered that
   the enriched pathways associated with the C0 SOD3 + Melanoma cells
   subpopulation included “extracellular matrix organization”,
   “extracellular structure organization”, “external encapsulating
   structure”, and “cell-substrate adhesion organization”. Existing
   studies had shown that the “extracellular matrix organization” pathway
   was closely related to cancer cell metastasis [[116]51]. Therefore, we
   inferred that the C0 SOD3 + Melanoma cells subpopulation was
   intricately linked to melanoma metastasis and might achieve this
   through pathways involving extracellular matrix and tissue adhesion.

Fig. 7.

   [117]Fig. 7
   [118]Open in a new tab

   Enrichment analysis of melanoma cell subpopulations. A–H DEGs within
   each melanoma cell subpopulations. I–P Functional enrichment analysis
   of GO-BP terms for differential genes within each melanoma cell
   subpopulation

Cell trajectory analysis of the subpopulations

   For a clearer insight into the developmental trajectories of melanoma
   cell subpopulations, we used pseudotime analysis to infer their
   progression. We utilized CytoTRACE to predict the cellular stemness of
   each subpopulation of melanoma cells (Fig. [119]8A, B). We observed
   that the C0 SOD3 + Melanoma cells subpopulation exhibited the strongest
   cellular stemness, followed by the C2 PMEL + Melanoma cells
   subpopulations. Additionally, we utilized slingshot to elucidate the
   temporal ordering of subpopulations (Fig. [120]8C, D), which revealed
   four distinct trajectories corresponding to various cellular
   differentiation pathways. Notably, all four trajectories originated
   from the C0 SOD3 + Melanoma cells subpopulation, suggesting that this
   subpopulation exhibited higher level of cellular stemness and functions
   as the primary source for the differentiation of the other
   subpopulations. Subsequently, we validated these results through
   monocle analysis and examined the temporal expression order of
   stemness-associated genes (Fig. [121]8E). Additionally, we analyzed the
   subpopulation trajectories calculated by monocle (Fig. [122]8F, G),
   revealing that C0 SOD3 + Melanoma cells was likely positioned at the
   initial phase of the temporal order, consistent with the results
   obtained through slingshot and CytoTRACE. Through our investigations,
   we identified the C0 SOD3 + Melanoma cells subpopulation as the key
   subpopulation, characterized by its strong cellular stemness, metabolic
   activity, and high correlation with metastasis.

Fig. 8.

   [123]Fig. 8
   [124]Open in a new tab

   Pseudotime analysis of subpopulation cells. A Using CytoTRACE,
   predicted ordering was performed on the 8 subpopulations, with higher
   scores indicating greater differentiation potential. B Distribution
   characteristics of predicted order and corresponding cell
   subpopulations of melanoma cells. C, D Prediction of time trajectories
   for subpopulations using slingshot, revealing 4 distinct trajectories,
   with UMAP plots illustrating each trajectory.
   (Lineage1:C0 → C1 → C2 → C3; Lineage2:C0 → C1 → C4;
   Lineage3:C0 → C1 → C5; Lineage4:C0 → C1 → C7. “ → ” represented the
   direction of pseudotemporal trajectory differentiation. E Dynamic
   expression trend of different stemness genes over time based on monocle
   analysis. F, G Violin and UMAP plots illustrated the pseudotime
   expression levels of different melanoma cell subpopulations

Stemness and TFs analysis

   To further investigate the crucial role played by tumor development,
   cell proliferation, metastatic behavior, and drug resistance, we
   investigated the stemness-associated genes of each cell subpopulation,
   different cell cycle phases, as well as primary and metastatic tumor
   tissues. (Fig. [125]9A–D). We observed significantly higher cellular
   stemness in the C0 SOD3 + Melanoma cells subpopulation, G1 phase, and
   metastatic tumor tissues. Moreover, we examined the expression patterns
   of stemness-associated genes across subpopulations, cell cycle phases,
   as well as primary and metastatic tumor tissues (Fig. [126]9E–H).
   Notably, TWIST1, MYC, HLFA, CD44, exhibited significant expression
   across C0 SOD3 + Melanoma cells subpopulation, and their expression
   levels were visualized using UMAP plots (Fig. [127]9I–L). Furthermore,
   stemness-associated genes showed pronounced expression in the C0
   SOD3 + Melanoma cells subpopulation, G1 phase, and metastatic tumor
   tissues, suggesting their involvement in the activity of the C0
   SOD3 + Melanoma cells subpopulation and potentially in the metastatic
   process.

Fig. 9.

   [128]Fig. 9
   [129]Open in a new tab

   Stemness analysis of cells. A–C The AUC value of cell stemness within
   each subpopulation, phase, as well as in the primary and metastatic
   tumor tissues. D The UMAP plot demonstrated the AUC density
   distribution characteristics of stemness across the cell
   subpopulations. E The expression levels of stemness genes in the 8
   subpopulations, as well as in primary and metastatic melanoma tissues,
   were analyzed. F–H Expression levels of stemness genes within each
   subpopulation of cells, across different phases, as well as within
   primary and metastatic tumor tissues. I–L Expression distribution
   characteristic of stemness genes TWIST1, MYC, CD44, and HIF1A within
   UMAP plots

   To delve deeper into the specific mechanisms underlying tumor
   metastasis, we investigated the TFs associated with primary and
   metastatic tumor tissues, computed and selected the top 5 TFs, and
   visualized them using UMAP plots (Fig. [130]10A–N). We found that TFs
   related to metastatic tumor tissues were predominantly expressed in the
   region associated with the C0 SOD3 + Melanoma cells subpopulation,
   while TFs related to primary tissues were distributed among various
   subpopulations. Additionally, we observed that TCF7L2 ranked highly in
   specificity scores within melanoma metastasis samples. Consequently, we
   conducted an in-depth analysis of this TF, focusing on its
   visualization across different melanoma subpopulations, cell cycle
   phases, and various classifications (Fig. [131]10O–R). The results
   indicated that TCF7L2 exhibited expression levels in the C0
   SOD3 + Melanoma cells subpopulation and metastasis samples. Notably,
   TCF7L2 was linked to the advancement of various cancers, particularly
   gastric cancer, where it played a key role in regulating upstream and
   downstream genes during transcription, enhancing tumor cell drug
   resistance, and significantly promoting apoptosis and metastasis.
   Moreover, the presence of TCF7L2 was often observed in cases where
   gastric cancer patients had poor prognoses [[132]52].Therefore, we
   might hypothesize that TCF7L2 played an integral part in the metastatic
   process of melanoma. This TF could serve as a potential therapeutic
   target.

Fig. 10.

   [133]Fig. 10
   [134]Open in a new tab

   TFs analysis. A, H TFs distribution characteristics of primary and
   metastatic melanoma tissues. B, I The scatter diagrams displayed the
   specificity scores of TFs within regulons in primary and metastatic
   melanoma tissues. C–G, J–N UMAP plots showed distribution of top 5 TFs
   in primary and metastatic tumor tissues. O UMAP plot displayed the
   distribution of TCF7L2 within subpopulations. P–R Violin plots depicted
   expression levels of TCF7L2 within each different subpopulations, cell
   cycle phases, as well as primary and metastatic tumor tissues.
   (****P < 0.0001. "ns" was used to say that there was no significant
   difference)

Construction and analysis of the prognostic model

   To further explore the patients' prognostic status, we performed an
   intersection analysis between the DEGs from the TCGA database and the
   signature genes of the C0 SOD3 + Melanoma cells subpopulation in
   melanoma. Using univariate Cox regression analysis, we identified 38
   risk genes associated with prognosis (Fig. [135]11A). The stability and
   reliability of these genes were confirmed through LASSO Cox regression
   analysis (Fig. [136]11B, C). Furthermore, through multivariable Cox
   regression analysis, we ultimately identified 15 prognostic genes and
   performed a coef analysis. The figure presented the coef values of
   these prognostic genes (Fig. [137]11D). Given that IGF1 had relatively
   higher coef values, we hypothesized its important role in melanoma
   metastasis and progression. Consequently, we selected this gene for
   further experimental validation. Among them, IGF1, CD81, FBLN1, GPC3,
   and APCDD1 were classified as high-risk genes, while APOD, PLSCR4, MGP,
   MRPS6, GADD45A, SERPING1, SPRY1, MT2A, RGS16, and IFITM1 were
   categorized as low risk genes. Risk scores were used to categorize
   patients into high-risk and low-risk groups (Fig. [138]11E), and their
   survival status over time was assessed (Fig. [139]11F). We observed a
   concentration of deceased patients in the high-risk group, with a
   significant decrease in the number of survivors over time.
   Additionally, we displayed the expression levels of prognostic genes in
   the high-risk and low-risk groups using a heatmap (Fig. [140]11G). The
   ROC curve's area under the curve values at 1-year, 3-years and 5-years
   were 0.69, 0.69, and 0.74, respectively (Fig. [141]11H), indicating the
   stable and reliable predictive capability of our model. Kaplan–Meier
   survival curves (Fig. [142]11I) revealed a substantial difference in
   survival rates between the high-risk and low-risk groups, with a
   statistically significant p-value of less than 0.0001. We also explored
   the correlation between genes and risk scores, highlighting genes with
   strong correlations (Fig. [143]11J, K). The expression levels of
   APCDD1, IGF1, CD81, and FBLN1 were positively correlated with risk,
   while MRPS6, PLSCR4, SERPING1, and IFITM1 were negatively correlated
   with risk. In other words, higher expression levels of APCDD1, IGF1,
   CD81, and FBLN1 were associated with an poor prognosis.

Fig. 11.

   [144]Fig. 11
   [145]Open in a new tab

   Independent prognostic analysis. A Univariate regression analysis
   uncovered 38 genes connected to prognosis. (HR < 1: protective factors,
   HR > 1: risk factors). B The coef spectrum distribution of 15
   prognostic genes as determined by LASSO regression analysis. C
   Parameter selection in optimal cross-validation LASSO regression. D The
   coef values of genes associated with prognosis. E Sorting patients into
   high-risk and low-risk groups according to risk scores. F Scatter
   diagram showed the temporal trends in the occurrence of death and alive
   events. G The heatmap visualized the expression levels of prognostic
   genes in high-risk and low-risk groups. H Time-dependent ROC curves
   with AUC values of 0.69, 0.69, and 0.74 for 1-year, 3-years, and
   5-years survival, respectively. I Kaplan–Meier survival analysis curves
   for high-risk and low-risk groups. J Correlation analysis between genes
   and risk scores. K Scatter diagrams illustrated the correlation of 8
   prognostic genes with risk scores

   Furthermore, we investigated the association between clinical variables
   like age, tumor stages, and prognosis (Fig. [146]12A–D). In addition,
   we performed risk assessments based on patients' risk groups, race,
   age, and tumor stages. Forest map revealed some independent risk
   factors (P < 0.05) (Fig. [147]12E). Based on the prognostic genes and
   clinical factors, we constructed a nomogram chart to predict the
   survival rates of melanoma patients at 1-year, 3-years, and 5-years
   (Fig. [148]12F). The accuracy of our nomogram chart model was validated
   using ROC curve, showing AUC values of 0.72, 0.79, and 0.79 at 1-year,
   3-years, and 5-years, respectively (Fig. [149]12G), indicating its high
   specificity and sensitivity.

Fig. 12.

   [150]Fig. 12
   [151]Open in a new tab

   Clinical relevance analysis. A–D Comparative analysis of risk scores
   between different prognostic factor groups. E The forest plot
   demonstrated the results of multivariate cox regression analysis
   integrating risk scores and clinical factors (age, race and tumor
   clinical stage T, M and N). F Nomogram showed the prediction of 1-year,
   3-years, and 5-year of OS based on race, tumor clinical stage (T, M,
   and N), age, and risk score. G Time-dependent ROC curves with AUC
   values of 0.72, 0.79, and 0.79 for 1-year, 3-years, and 5-years
   survival, respectively

Immuno-related and enriched functional analysis

   Additionally, we carried out an in-depth analysis to explore the role
   of immune infiltration in disease progression. We conducted an analysis
   of immune infiltration in both high-risk and low-risk groups
   (Fig. [152]13A, B) and identified immune cell populations with higher
   proportions (Fig. [153]13C), where macrophages constituted a
   significant portion. Further comparison of immune cell compositions
   between the two groups (Fig. [154]13D) revealed a higher proportion of
   M1 macrophages in the low-risk group and a higher proportion of M2
   macrophages in the high-risk group. Additionally, M1 macrophages
   exhibited a negative correlation with the risk score, while M2
   macrophages showed a positive correlation (Fig. [155]13E). In summary,
   we hypothesized that macrophages might be involved in the proliferation
   and metastasis of melanoma and thus affect the prognosis of patients,
   but the roles played by M1 and M2 macrophages were different.

Fig. 13.

   [156]Fig. 13
   [157]Open in a new tab

   Immune-related analysis. A The heatmap displayed the differences in
   immune infiltration expression of various factors between the high-risk
   and low-risk groups. B Differential proportions of immune cells between
   high-risk and low-risk groups. C Comparison of immune cell infiltration
   and estimated proportions. D Differential estimated proportions of
   immune cells between high-risk and low-risk groups. E Correlation
   analysis between risk scores and different immune cells. F Heatmap
   analysis of the correlation between gene expression, risk score, OS,
   and immune cell infiltration. G–I Differential immune scores, stromal
   scores, and ESTIMATE scores of TME in high-risk and low-risk groups. J
   Violin plot displayed TIDE expression levels in high-risk and low-risk
   groups. (*P < 0.05, **P < 0.01, ***P < 0.001, and ****P < 0.0001.)

   We also explored the interplay between prognostic genes, risk, OS and
   immune cells (Fig. [158]13F) and found an association between the OS
   and both M1 and M2 macrophages, displaying consistent with previous
   findings, M1 macrophages were positively correlated with OS, while M2
   macrophages were negatively correlated with OS. M2 macrophages might
   have been associated with tumor progression. Regarding TME-related
   scores, the low-risk group exhibited higher scores for stromal, immune,
   and estimate scores (Fig. [159]13G–I). While TIDE scores were also
   higher in the low-risk group, suggesting that patients in the high-risk
   group might have a poorer response to treatment with immune checkpoint
   inhibitors (Fig. [160]13J).

   Furthermore, GSEA revealed notable enrichment of pathways connected to
   "epidermis development", “skin development” and "keratinization" in the
   positive expression trend group, which were associated with skin
   metabolism. Conversely, pathways related to "positive regulation of
   leukocyte cell–cell adhesion", "regulation of lymphocyte mediated
   immunity," and "adaptive immune response" showed a negative expression
   trend (Fig. [161]14A–F).

Fig. 14.

   [162]Fig. 14
   [163]Open in a new tab

   GSEA depicted enriched functional pathways. A–C Enrichment analysis of
   pathways associated with positive trend expression. D–F Pathways
   enrichment analysis related to pathways with negative expression trends

In vitro experiment of IGF1

   Based on analysis of prognostic genes, we selected IGF1 as the target
   gene for our experiments. We conducted experiments using two melanoma
   cell lineages, employing methods to compare negative control and
   knockdown infection groups. In the cell viability assay (Fig. [164]15A,
   B), CCK-8 tests indicated significantly reduced cell viability after
   knockdown. In colony formation assays (Fig. [165]15C, D), the colony
   numbers formed by IGF1 knockdown cells was notably lower in contrast to
   the negative control group. Cell proliferation assays (Fig. [166]15E,
   F) showed a decreased cell count in both melanoma cell lines following
   IGF1 knockdown compared to the negative control group. In the wound
   healing assays (Fig. [167]15G, H), the width of scratches made on IGF1
   knockdown cells was significantly wider than those on the negative
   control group after 48 h. In the Transwell assays, the migration of
   both cell lines treated with IGF1 knockdown showed significantly
   smaller stained areas against the negative control group
   (Fig. [168]15I, J). Similarly, in the invasion assay (Fig. [169]15K,
   L), the two melanoma cell lines with IGF1 knockdown exhibited markedly
   reduced invasion compared to the negative control group. Therefore,
   through these experiments, we found that knockdown of IGF1 in melanoma
   cells led to decreased proliferation, migration, and invasion
   capabilities, indicating that IGF1 probably played a positive role in
   the proliferation, migration, and invasion of melanoma, thereby
   promoting its progression.

Fig. 15.

   [170]Fig. 15
   [171]Open in a new tab

   In vitro assays of prognostic gene IGF1 on proliferation and invasion
   metastasis. A, B CCK-8 assays revealed significantly reduced cell
   viability after IGF1 knockout. C, D Colony formation assays
   demonstrated a notable decrease in colony numbers in cells with IGF1
   knockdown against the negative control group. E, F EdU staining results
   indicated suppressed proliferation of SK-MEL2 and A375 cells upon IGF1
   knockdown. G, H Wound healing assays revealed markedly decelerated
   migration of SK-MEL2 and A375 cells following IGF1 knockdown. I, J
   Transwell assays showed significantly decreased migration of SK-MEL2
   and A375 cells upon IGF1 knockdown. K, L Transwell assays demonstrated
   a significant reduction in invasion of SK-MEL2 and A375 cells following
   IGF1 knockdown

Discussion

   Melanoma is a fast-growing and aggressive tumor, often linked to a poor
   prognosis. Its heterogeneity and the complexity of the TME posed
   significant challenges for treatment [[172]53]. In recent years,
   single-cell sequencing technology, as an innovative tool, provided
   precise insights into tumor heterogeneity and the dynamic changes of
   the microenvironment. By analyzing gene expression, immune
   characteristics, and metabolic states at the single-cell level, this
   technology enabled the identification of key pathological cell
   populations and their functional properties, offering a scientific
   basis for developing personalized therapeutic strategies.

   This study leveraged single-cell sequencing technology to
   comprehensively investigate the TME and cellular heterogeneity in
   melanoma [[173]54]. The aim was to elucidate the mechanisms of key
   pathological cell populations. The findings of this research not only
   contributed to optimizing treatment strategies and improving patient
   prognosis but also advanced the clinical application of single-cell
   technology in melanoma treatment, providing theoretical and data
   support for the development of personalized therapeutic approaches.

   In tumor types, metastasis and dissemination often signify poor
   prognosis for patients, with approximately 90% of cancer-related deaths
   attributed to metastasis [[174]55]. Melanoma, as one of the most
   malignant forms of skin neoplasm, lacked effective treatment options
   once metastasized, owing to the incomplete elucidation of its
   underlying mechanisms and molecular targets [[175]56]. Therefore,
   investigating and refining the mechanisms and treatment targets of
   melanoma metastasis, as well as identifying new therapeutic directions,
   was paramount. Hence, in this study, we selected key subpopulations
   associated with metastasis, explored their intrinsic metabolic
   pathways, constructed prognostic models relevant to them, and
   investigated their relationship with immunity.

   Our analysis of metabolic pathways revealed that oxidative-related
   pathways, particularly oxidative phosphorylation, ranked highest across
   all subpopulations. This pathway showed the strongest expression in the
   C0 SOD3 + Melanoma cells subpopulation and was notably elevated in
   metastatic tumor tissues compared to primary tumor tissues. Hence, we
   hypothesized that oxidative-related pathways, particularly oxidative
   phosphorylation, played a crucial role in melanoma metastasis.
   Literature review revealed the vital importance of oxidative
   phosphorylation in tumor biology, as inhibiting this pathway's activity
   could reduce the anti-tumor effect and cell proliferation of melanoma
   cells [[176]57–[177]59]. Intriguingly, certain studies suggested that
   tumor stem cells with high metastatic potential were more dependent on
   the oxidative phosphorylation process, and inhibiting this pathway
   could improve tumor recurrence [[178]60, [179]61], aligning with our
   findings of significantly elevated oxidative phosphorylation pathway
   expression in metastatic cells. Current studies [[180]62] had
   demonstrated that inhibiting oxidative phosphorylation could suppress
   ARNT and PDK1-deficient melanoma cell migration and invasion.
   Additionally, the oxidative phosphorylation process was accompanied by
   ROS production, and further studies had [[181]63] found that carbon
   dots-induced ROS promoted melanoma cell growth and invasiveness. What
   is more, other studies had indicated that ROS exert inhibitory effects
   on melanoma dissemination and metastasis [[182]64, [183]65].

   We believed that this discrepancy could be correlated with ROS, a
   hypothesis that necessitates further experimental investigation and
   confirmation. Subsequently, to further investigate the evolutionary
   trajectory of melanoma subpopulations, we employed pseudotime analysis.
   We observed that, regardless of lineage, the C0 SOD3 + Melanoma cells
   subpopulation consistently emerged as the initial subpopulation with
   the highest predicted ordering by CytoTRACE. Thus, compared to other
   subpopulations, the C0 SOD3 + Melanoma cells subpopulation was more
   primitive, with lower differentiation and higher proliferation
   capacity, indicative of higher malignancy and tumor-promoting
   potential. The unique characteristics of the C0 SOD3 + The melanoma
   cell subpopulation indicated that it might play a crucial role in
   melanoma treatment. Furthermore, when combined with the trajectory
   analysis provided by Slingshot, which traced the developmental pathways
   of melanoma subpopulations, our findings further emphasized the central
   role of the C0 SOD3 + Melanoma cells subpopulation in the early stages
   of melanoma progression. Slingshot’s reconstruction of the lineage
   relationships reinforced the notion that this subpopulation might have
   been pivotal in initiating malignant transformation, making it an ideal
   candidate for targeted therapeutic strategies.

   Stemness analysis revealed that the C0 SOD3 + Melanoma cells
   subpopulation exhibited a higher stemness score, indicating a more
   undifferentiated state. Higher stemness scores are typically associated
   with increased self-renewal capacity and the potential for
   differentiation into various cell types. This suggests that the C0
   SOD3 + Melanoma cells subpopulation not only maintains a primitive
   state but also contributes to tumor heterogeneity and progression.
   Cells with elevated stemness are often linked to poor prognosis, as
   they can drive metastasis and resist conventional therapies. Therefore,
   targeting the stemness characteristics of the C0 SOD3 + Melanoma cells
   subpopulation could provide new therapeutic strategies for melanoma
   treatment. In the analysis of stemness gene expression, key genes such
   as TWIST1, MYC, CD44, and HIF1A were predominantly expressed. TWIST1
   was a TF involved in the epithelial-mesenchymal transition (EMT)
   process [[184]66], playing a crucial role in the initiation and
   metastasis of melanomas [[185]67]. During EMT, TWIST1 suppressed EPCs
   characteristics while inducing the expression of mesenchymal markers,
   thereby endowing cancer cells with enhanced migratory and invasive
   capabilities. TWIST1 was implicated in the progression of various
   cancers, such as breast cancer, lung cancer, and gastric cancer, often
   correlating with aggressive tumor behavior, metastasis, drug
   resistance, and poor prognosis [[186]68–[187]70]. MYC was a well-known
   oncogene that played a pivotal role in cancer cell invasion and
   migration, thus promoting tumor progression [[188]71]. It was involved
   in the expression of multiple oncogenic non-coding RNAs, and when
   aberrantly expressed, MYC disrupted normal signaling pathways, leading
   to crosstalk between oncogenic pathways and possibly inducing
   replication stress, which could, in turn, lead to cancer development
   [[189]72, [190]73].CD44 was a transmembrane glycoprotein that acted as
   an adhesion molecule on the cell surface, promoting cell-extracellular
   matrix interactions. Its adhesive properties and ability to degrade
   hyaluronic acid made it a key driver of tumor invasion and metastasis
   [[191]74]. Additionally, CD44 was recognized as a cancer stem cell
   marker in various tumors [[192]75], closely associated with tumor
   recurrence, metastasis, and chemoresistance. While CD44 had both
   oncogenic and tumor-suppressive roles, in this study, we focused on its
   oncogenic effects. By binding to its ligands, such as hyaluronic acid,
   CD44 strengthened the connection between tumor cells and the
   extracellular matrix, thereby promoting tumor cell invasiveness and
   metastatic potential [[193]76]. HIF1A was a TF upregulated under
   hypoxic conditions, regulating cellular adaptation to hypoxia, which
   was especially critical in the TME. Hypoxia was a key element in tumor
   survival and growth, and the TME was frequently hypoxic. Activation of
   HIF1A helped tumor cells adapt to low oxygen conditions, promoting
   immune evasion by releasing exosome-mediated signaling, thus
   contributing to tumor progression [[194]77, [195]78].In summary, these
   genes and proteins regulated tumor cell proliferation, survival,
   migration, and metastasis through distinct mechanisms. They were not
   only key drivers of tumor progression but also potential therapeutic
   targets.

   To gain a comprehensive understanding of melanoma patient prognosis, we
   constructed a prognostic model to analyze genes that might contribute
   to poor prognosis, to intervene and improve patient prognosis. Based on
   current data, the IGF1 gene scored the highest coef values. Existing
   studies suggested that IGF1 was closely associated with drug
   resistance, and targeting the IGF1 gene could enhance drug sensitivity
   to some extent, leading to improved therapeutic outcomes for patients
   [[196]79]. IGF1 was often associated with poor prognosis in various
   cancers, such as breast cancer and Wilms tumor [[197]80, [198]81]. It
   contributed significantly to the initiation, progression, and prognosis
   of many cancers by activating different signaling pathways through the
   IGF1 receptor, promoting tumor cell proliferation, survival, invasion,
   and metastasis. For example, in colorectal cancer, IGF1 exerted its
   effects on downstream genes via the PI3K/AKT/HIF1α pathway, enhancing
   tumor cell proliferation and inhibiting apoptosis, thereby increasing
   cancer cell invasiveness and metastatic potential [[199]82].In
   conclusion, while research on IGF1 in melanoma was limited, based on
   these findings, we could infer that the IGF1 gene played a significant
   role in melanoma. The IGF1 gene was poised to be a critical target for
   improving patient prognosis.

   Therefore, to verify the role of IGF1 in tumor cells, we conducted
   experiments to support our conclusions. Consistent with the
   aforementioned findings, in assays such as cell viability, wound
   healing assays, and invasion tests, our experiments exhibited that IGF1
   was instrumental in promoting the proliferation, invasion, and
   migration of melanoma cells. These experimental results further
   reinforced that IGF1 was a promising therapeutic target for improving
   patient prognosis. Therefore, more attention should be given to the
   IGF1 gene, and by targeting this gene, inhibiting its activity could
   potentially lead to improved patient outcomes.

   In our investigation of immune-related phenomena, we observed a
   relatively high infiltration of macrophages, among which both M1 and M2
   macrophages exhibited significant correlations with risk scores, albeit
   in opposing directions. Typically, M1 and M2 macrophages could
   interconvert, with M1 primarily associated with inflammation and
   exerting anti-tumor activity, while M2 was associated with tissue
   repair and often promoted tumor progression and invasion
   [[200]83–[201]85]. Consistent with our findings, M1 was positively
   associated with patient prognosis, while M2 was indicative of a poorer
   prognosis. Further literature review revealed that M1 macrophages
   infiltrating the immune microenvironment could inhibit tumor metastasis
   through phagocytosis and immune regulation [[202]86–[203]88], while M2
   macrophages could promote tumor progression and metastasis through
   mechanisms such as secreting CCL5, upregulating CRYAB expression, and
   activating the ERK1/2/Fra-1/slug signaling pathway [[204]89–[205]91].
   Studies had shown that inhibiting M2 macrophage polarization could
   effectively suppress melanoma metastasis [[206]92]. Intriguingly,
   existing research had [[207]93] found that activation of
   oxidative-related pathways could lead to immune suppression and
   polarization of macrophages from M1 to M2. Mouton [[208]94] also
   discovered that M2 macrophage activity and conversion were highly
   dependent on the tricarboxylic acid cycle and oxidative
   phosphorylation. Additionally, ROS, besides participating in tumor
   formation, could induce macrophage polarization toward M2 and promote
   the production of IFN and IL-6 by macrophages, leading to a suppressive
   tumor immune microenvironment [[209]95, [210]96], ultimately resulting
   in M2 macrophage polarization and consistent immune-related functions,
   thereby increasing the likelihood of tumor escape and melanoma
   metastasis. Therefore, we proposed that the activation of oxidative
   phosphorylation and related oxidative pathways in the C0
   SOD3 + Melanoma cells subpopulation and metastatic tumor tissue cells
   might induce M2 macrophage polarization, thereby triggering a
   suppressive immune microenvironment and enhancing the likelihood of
   melanoma metastasis.

   We conducted an in-depth study on melanoma using data obtained through
   single-cell sequencing technology. We analyzed the heterogeneous
   subpopulations of cells within melanoma tumors, providing
   high-resolution insights into the distinct functions of diverse cell
   groups within the TME, particularly the key mechanisms involved in
   tumor metastasis. This cell-level information aid in identifying
   critical subpopulations related to malignant metastasis and prognosis,
   such as the C0 SOD3 + subpopulation identified in the study. Apart from
   examining tumor cells, we extended our analysis to the distinct roles
   of M1 and M2 macrophages in the TME, revealing their contributions to
   melanoma progression and investigating their associations with
   prognostic genes. This offered potential targets for the development of
   novel immunotherapies [[211]97]. In melanoma, PCD pathways had critical
   roles in disease progression and therapeutic resistance [[212]98].
   Melanoma cells frequently evaded apoptosis and autophagy, which allowed
   sustained tumor growth and resistance to conventional therapies.
   However, inducing ferroptosis and pyroptosis showed promise in
   stimulating anti-tumor immune responses, potentially limiting melanoma
   dissemination. By disrupting redox balance and releasing
   pro-inflammatory factors, these PCD pathways facilitated immune
   activation and reshaped the TME. Targeting PCD in melanoma thus
   provided innovative strategies to overcome therapeutic resistance and
   improve patient outcomes. However, there were certain limitations to
   our study. Firstly, the dataset was limited. Although the study
   analyzed single-cell samples from 8 melanoma patients, the sample size
   was relatively constrained, which might restrict the generalizability
   of the results, particularly for different types or stages of melanoma.
   Additionally, our experimental focus was primarily on in vitro cellular
   studies, with a lack of in vivo animal model experiments. Therefore, it
   is essential to conduct a more in-depth exploration.

Conclusion

   Our study delved into the detailed mechanisms of melanoma metastasis at
   the single-cell level, exploring the roles of oxidative stress pathways
   and the immune system therein. Firstly, we identified the critical
   subpopulation, C0 SOD3 + Melanoma cell subpopulation, which played a
   pivotal role in disease progression. These cells exhibited antioxidant
   properties and immune modulation capabilities, contributing to tumor
   plasticity and resistance to oxidative stress. Targeting this
   subpopulation could have disrupted key pathways that facilitated
   melanoma progression and metastasis. Therefore, understanding the
   functional characteristics of C0 subpopulation provided valuable
   insights for developing novel therapeutic strategies against melanoma.
   Additionally, we constructed a prognostic model associated with these
   findings. Thus, key prognostic genes were identified, and for this
   purpose, we analyzed them. What is more, our findings showed that the
   IGF1 gene was crucial for prognosis, with its expression associated
   with patient outcomes. Furthermore, experimental validation confirmed
   the critical role of the IGF1 in the proliferation and invasive
   metastasis of melanoma. These discoveries contribute to a deeper
   understanding of the underlying mechanisms of melanoma metastasis,
   offering new therapeutic avenues and targets for the treatment of
   melanoma. Besides, they also provide a basis for creating new
   prognostic models for melanoma diagnosis and prognosis. Furthermore,
   the outcomes of this study demonstrate crucial clinical implications,
   particularly for guiding therapy and improving prognostic predictions.
   Understanding melanoma cell heterogeneity and the role of immune
   infiltration in disease progression can lead to more personalized
   treatments and better outcome predictions.

Acknowledgements