Lactylation-driven molecular taxonomy of melanoma: linking epigenetic modifications to immune evasion and clinical outcomes

Abstract

Lactylation, a post-translational modification derived from elevated lactate levels, has gained attention as a potential regulator of melanoma's tumor metabolism and immune responses. Here, we combined single-cell RNA sequencing and bulk transcriptome profiling of cutaneous melanoma samples to establish a lactation-centric prognostic model. Our analyses revealed melanocytes as the most acetylation-enriched cell population and identified a six-gene lactylation signature that stratified patients into high- and low-risk groups with distinct survival outcomes. Mechanistically, high-risk tumors demonstrated significant immunosuppressive features characterized by M2 macrophage accumulation and depleted CD8+ T-cell activity, corresponding to reduced sensitivity to certain chemotherapeutic drugs. Pathway enrichment studies implicated DNA repair, Hedgehog, and JAK-STAT signaling in driving the aggressive phenotype of high-acetylation tumors. Additionally, pseudotime trajectory analyses highlighted developmental shifts in gene expression related to lactylation during melanocyte differentiation. The signature demonstrated robust predictive accuracy in training, testing, and external validation cohorts. Functional validation confirmed the critical role of RAN in promoting proliferation and migration in vitro. These findings unveil lactylation as a critical epigenetic factor influencing melanoma progression and immune evasion, offering a novel prognostic framework and potential therapeutic targets for precision medicine.

Keywords: Lactylation; immune microenvironment; melanoma; multi-omics; prognostic model; single-cell RNA sequencing.

Figure 1

Cell types and gene expression patterns of melanoma. A. UMAP visualization of single-cell data colored by Leiden cluster assignment (12 clusters). B. Annotations of seven major cell lineages (melanocytes, T cells, B cells, fibroblasts, endothelial cells, macrophages, and keratinocytes) based on marker genes. C. Bubble plot of classic marker genes illustrating their expression across these cell types. Red indicates high expression, and blue indicates low expression. D. Pie chart showing relative proportions of the seven cell types. E. Heatmap displaying log2 fold changes for the top six differentially expressed genes (DEGs) in each subpopulation.

Figure 2

Lactylation-associated cell subtypes and cell-cell communication. A. Bubble plot, Violin plot and Scatter plots showing comparative profiling in melanocyte populations to identify lactylation-enriched cellular subtypes. B. Venn diagram indicating the intersection of differentially expressed lactate-related genes. C. Intercellular communication network among diverse cell populations, visualized by line widths representing interaction strength. D. Bubble plots depicting incoming and outgoing signaling pathways between melanocytes with T cells, B cells and other clusters.

Figure 3

SCENIC analysis in melanocytes. A. Heatmap of regulon activity scores in individual melanocyte cells, indicating distinct transcriptional states. B. Scatter plot ranking transcription factors by regulon specificity score (RSS) in high-expression melanocytes (Hexp). C. Scatter plot ranking transcription factors by RSS in low-expression melanocytes (Lexp).

Figure 4

Prognostic model construction and validation. A. Forest Plot of seven prognostic genes (P < 0.05) from an initial list of 27 lactylation-related genes using univariate Cox regression. B-D. LASSO regression screening to refine these genes into a six-gene signature. E, F. Kaplan-Meier survival curves showing significantly lower overall survival (OS) in the high-risk group compared with the low-risk group for both training and testing cohorts. Statistical significance was determined by the Log-rank test (P < 0.001). G, H. ROC curves demonstrating strong predictive accuracy at 2, 4, and 6 years in the training and testing sets. The AUC values range between 0.5 and 1, where 0.5 indicates no discriminative ability and 1 represents perfect discriminative ability. I. External validation in the GSE53118 GEO dataset, confirming poorer OS for the high-risk group.

Figure 5

Correlation of the risk model with immune microenvironment. A. CIBERSORT-based immune cell composition in high- and low-risk patient groups. B. Heatmap illustrating correlations among immune cell subsets. C. Boxplots highlighting statistically significant differences in immune cell proportions between high- and low-risk tumors. D. Lollipop plot showing a positive correlation of the risk score with immunosuppressive macrophages (M0, M2) and a negative correlation with cytotoxic T cells.

Figure 6

Nomogram model for individualized prognostic prediction. A. Nomogram integrating the lactylation-driven risk score with clinicopathological variables to estimate survival probabilities. B. Calibration curve comparing predicted and observed survival outcomes in melanoma patients. C. Time-dependent ROC curves (2-, 4-, and 6-year) indicating the nomogram’s discriminative accuracy (AUCs: 0.6784, 0.7042, 0.7006). Shaded regions denote 95% confidence intervals. D. Decision curve analysis (DCA) demonstrating a superior net benefit of the nomogram versus conventional staging.

Figure 7

Drug sensitivity and associated signaling pathways. A. Predicted chemotherapy response (IC50) using the oncoPredict R package, anchored by GDSC database profiles. High-risk patients showed reduced sensitivity (elevated IC50) to vinblastine_1004, paclitaxel_1080, and vincristine_1818 (P ≤ 0.01). B, C. GSVA and GSEA highlighting distinct enrichment of pathways (DNA repair, Hedgehog, JAK-STAT) in high- versus low-risk groups. D. A molecular interaction network illustrating cross-talk among these enriched pathways.

Figure 8

Model gene expression in single cells. A, B. DotPlot and Violin Plot visualizing expression patterns of the six model genes (LAP3, RBM39, THRAP3, RAN, DDX3X, S100A11) across melanoma single-cell clusters.

Figure 9

Pseudotime trajectory analysis of model genes. A. Pseudotime ordering of melanocytes, illustrating transitional states during differentiation. B. Distribution of cells across distinct states. C. Cluster assignments along the pseudotime axis, indicating transcriptional shifts among subpopulations. D. Heatmap of gene expression patterns at various stages, capturing early (BCAN, LIMD1, HMCN1) and late markers (MITF, KCNJ13, AFF3). E. Expression trends of model genes, highlighting mid-stage peaks (LAP3, RAN, S100A11) and U-shaped behaviors (RBM39, THRAP3).

Figure 10

RAN as a critical driver of melanoma proliferation and migration. A-F. qRT-PCR validation of the expression of four genes (LAP3, RBM39, S100A11, RBM39, THRAP3, DDX3X and RAN) in A375, SK-MEL-28, and NHEK cells. G. Box plots showing RAN expression in 461 skin cutaneous melanoma tissues vs 558 normal tissues based on GEPIA (http://gepia.cancer-pku.cn). H. Survival Curve plots demonstrating the impact of RAN on overall survival (OS) and disease-free survival (DFS) in melanoma based on GEPIA (http://gepia.cancer-pku.cn). I. The knockdown efficiency of RAN in A375 and SKMEL28 cells transfected with siRNA1/2 was detected using qRT-PCR. J. The knockdown efficiency of RAN in A375 and SKMEL28 cells transfected with siRNA1/2 was detected using Western Blotting. K. Proliferative capacity was detected by MTS assay. L. Wound healing assay was used to detect the migration ability. Data are shown as mean ± SD from three independent experiments. Significance is indicated by *P < 0.05.