Introduction:

Cutaneous melanoma (CM) is the most fatal and life-threatening form of malignant skin cancer. According to the World Health Organization, melanoma is responsible for approximately 57,000 deaths annually, making it one of the deadliest skin cancers (Arnold et al. 2022). Melanoma incidence has been increasing over the last 50 years. In 2020, more than 1.5 million new cases of skin cancer were reported. Malignant melanomas were responsible for one-fifth of these cancers. (Erdmann et al. 2013; Arnold et al. 2014). Its high metastatic capability is largely attributed to genetic and epigenetic alterations, including dysregulation of key oncogenic pathways (Rebecca et al. 2012). The incidence of melanoma is higher among the older generation, and the risk of the disease generally increases with age, with young adults being the most affected group (Fidler et al. 2017). In predominantly white populations, CM is responsible for the vast majority of melanoma cases, estimated at over 90%. Mucosal and uveal melanomas are significantly less common, accounting for less than 1–5% of cases, with variations depending on the country. In populations with darker skin tones, such as black, brown, and East Asian communities, cutaneous melanomas, particularly acral types, are the most common (Elder et al. 2020, 2018). Given these epidemiological findings, studying melanoma at a molecular level is essential for developing target therapies and markers which could improve patient treatment outcomes.

The mainstay of treatment for primary CM is excision with safety surgical margins, but once it metastasizes, treatment becomes more challenging and overall prognosis is poor. A favorable prognosis depends on early detection; however, inadequate treatment outcomes persist due to a lack of knowledge about the features of the TME in CM, as well as the expression and function of its regulators.

Furthermore, melanoma cells have a strong interaction with the tumor microenvironment and the immune system (Leonardi et al. 2018), indicating that the tumor microenvironment plays an important role in cancer initiation and metastasis (Visser and Joyce 2023). m6A is a key post-transcriptional modification that regulates RNA stability, splicing, and translation thereby influencing tumor progression. In melanoma, aberrant m6A modifications have been linked to altered gene expression patterns that promote tumor proliferation, immune evasion, and therapy resistance. This knowledge has led us to investigate post-transcriptional modifications, which have become a new trend in cancer research, with N6-methyladenosine (m6A) being the most common mRNA modification linked to the pathogenesis of many types of cancers (He et al. 2019; Wang et al. 2020). The three types of m6A regulators are writers, erasers, and readers, which are respectively responsible for m6A’s functions, methylation, and demethylation (Boulias and Greer 2023). Notably, m6A has been shown to play a critical role in melanoma oncogenesis and tumorigenicity. Specifically, m6A regulators such as METTL3 and FTO have been shown to modulate melanoma plasticity and drug sensitivity (Yu et al. 2021; Yang et al. 2019).

The function of tumor-infiltrating lymphocytes in the development and outcome of melanoma was recently thoroughly described by Dan Liu et al. (Liu et al. 2021). In the tumor microenvironment, Monica Marzagalli et al. uncovered the interplay between immune cells and melanoma (Marzagalli et al. 2019). Apart from tumor cells, TME cells comprise multiple cell types, including immune cells which have dual roles such as attacking tumor cells or, conversely, supporting tumor growth by creating an immunosuppressive environment. Furthermore, TME consists of stromal cells which signal molecules that promote tumor growth and facilitate metastasis and the non-cellular component which actively participate in cancer progression through angiogenesis and regulating cell behavior (Anderson and Simon 2020). Therefore, strategies that modulate the TME are being explored to enhance treatment outcomes. However, the cell–cell interaction between m6A mRNA modification-associated subtypes of TME cells and tumor cells has been the subject of a limited amount of research until recently. These findings prove that studying m6A modifications can provide novel insights into melanoma pathogenesis and potential therapeutic targets.

In this study, we examined the ramifications of m6A mRNA methylation on the primary TME cells, which include Melanocytes, T cells, B cells, EC, Fibroblasts, and Macrophages. Our analysis was based on 26,532 cells from three CM datasets that were sourced from public databases. 26 m6A RNA methylation regulators (Zaccara et al. 2019; Chong et al. 2021) were analyzed using NMF. It was observed that the m6A mRNA methylation sequences in each CM TME cell type subpopulation exhibited a substantial and distinctive link with tumor epithelial cells, and were associated with distinct immune and transcription characteristics, metabolic pathways, and prognosis. Essentially, our extensive single-cell analysis uncovers how m6A mRNA steers cell to cell interaction of TME cells with tumor cells and contributes to CM progression.

Materials and methods

Data collection

We obtained two skin cutaneous melanoma datasets from the GEO database: one single-cell dataset (GSE215121, N = 3) (Zhang et al. 2022) and one bulk-RNA dataset (GSE54467, N = 79) (Jayawardana et al. 2015). Additionally, we downloaded TCGA-SKCM dataset (N = 472) from UCSC (https://xenabrowser.net). The TCGA-SKCM dataset was used to establish a prognostic model of m6A-related clustering characteristic genes, while the GSE54467 dataset was used to further validate the expression profiles of characteristic genes. Since these public datasets were obtained from public databases, ethical approval was not required.

CM single-cell data analysis

First, the Read10X function in the Seurat package (v4.0.5) (Butler et al. 2018) in R software (v4.0.3) was used to read the data of each sample, and cells with mitochondrial expression greater than 20% were removed. This is because high mitochondrial expression usually indicates poor cell condition, possibly due to compromised cell membrane integrity, leading to a higher proportion of mitochondrial genes. Such cells may interfere with subsequent analyses, so it is necessary to filter them out based on the proportion of mitochondrial genes in the cells. We normalized the raw counts for each sample using Seurat function NormalizeData and identified 2000 variable genes using the FindVariableGenes function. Then, using the FindIntegrationAnchors function with default parameters and canonical correlation analysis (CCA) with 30 dimensions, we identified integration anchors between the three datasets. These anchors were integrated using the IntegrateData function with the first 30 CCA dimensions and default parameters. We scaled the data using the ScaleData function and performed principal component analysis (PCA) with the RunPCA function. Next, we performed unsupervised clustering on the data using the FindNeighbors and FindClusters functions. For the FindNeighbors function, we used the first 30 PCA dimensions to construct a shared nearest neighbor (SNN) graph for our dataset. We then clustered the cells using an SNN modularity optimization-based clustering algorithm with a resolution of 0.5. Finally, for visualization, we used the RunTSNE function with default parameters and 30 PCA dimensions. To identify cell groups, we used the FindAllMarkers function on the integrated dataset, which uses the Wilcoxon rank-sum test to identify representative genes for each cluster. These representative genes were used to determine the identity of each cluster and to identify cell types typically found in human skin based on literature markers. We used the R package BayesPrism (v2.0) (Chu et al. 2022), leveraging scRNA-seq as prior information to predict cell type composition from the TCGA melanoma dataset.

Cell–Cell communication analysis

We conducted an analysis of cell–cell communication analysis using the R CellChat package (v1.6.1) (Jin et al. 2021). This package contains a ligand-receptor interactions database for humans and mice, which allows for the analysis of cell–cell communication networks using scRNA-seq data annotated with various cell clusters. First, we used CellChat to assess the main signal inputs and outputs between cell clusters through CellChatDB.human database. Then, we used the netVisual_circle function to display the strength of cell–cell communication networks from the target cell cluster to all NMF cell clusters among different clusters. Finally, the netVisual_bubble function was employed to display bubble plots of significant ligand-receptor interactions between the target cell cluster and other clusters.

m6A gene set NMF clustering analysis

The 26 m6A gene sets were clustered across all TME cell types using the NMF R package, v0.25 (Gaujoux and Seoighe 2010). The different clusters of these cell types were determined based on the scRNA expression matrix. The execution of each of these steps was in line with previous studies (Chen et al. 2020; Puram et al. 2017).

Pseudotime analysis

We used the Monocle package (v2.22.0) in R to analyze the single-cell RNA data of fibroblast cell types (Qiu et al. 2017). Based on the following standards, highly variable genes were chosen: mean expression ≥ 0.1. The dynamic expression heatmap of the m6A genes was visualized using the plot_pseudotime_heatmap function, and dimensionality reduction was accomplished using the DDRTree method.

Transcription factor analysis of m6A gene set NMF clusters

Python-based SCENIC analysis software called pySCENIC (v0.9.0) was used to investigate the gene regulatory networks of transcription factors (TFs) in CM (Kumar et al. 2021). Two gene motif rankings from the RcisTarget database (hg19-tss-centered-10 kb and hg19-500 bp-upstream) were used to detect transcription start sites (TSS) and gene regulatory networks in the scRNA-seq data of CM. The threshold for transcription factors was set to Benjamini–Hochberg-False Discovery Rate (BH-FDR) < 0.05.

Construction of prognostic model of m6A cluster characteristic genes

Based on the FindAllMarkers function in the Seurat package, we generated characteristic genes between m6A gene set NMF clusters. Then, using multivariate Cox analysis based on overall survival, we further identified characteristic genes associated with patient prognosis and constructed a prognostic model with these genes. We used the survminer package (v0.4.9) in R to determine the model’s threshold, dividing all patients into low-risk and high-risk groups according to the threshold. Kaplan–Meier (K–M) survival curves were generated using the log-rank test to confirm prognostic differences between the two groups. The timeROC package (v0.4) in R was used to generate 1-year, 2-year, and 3-year ROC curves, and the corresponding area under the curve (AUC) was calculated.

Verification of prognostic characteristic genes

The GSE54467 dataset was used to externally validate the characteristic genes used to construct the prognostic model. The expression profiles of the characteristic genes included in the GSE54467 dataset were extracted. All patients in the validation cohort were divided into low-risk and high-risk groups using the method used to construct the prognostic model. The accuracy of characteristic genes in the validation cohort was evaluated using K-M survival analysis and ROC curves.

Utilizing the TIGER database for melanoma ICB datasets to investigate gene set behaviour in treatment outcomes

The TIGER database (The Tumor Immunotherapy Gene Expression Resource portal, http://tiger.canceromics.org/#/) (Chen et al. 2023a) contains an immunotherapy response module that includes melanoma data with immunotherapy information. This module was used to analyze the immunotherapy efficacy of different gene sets.

Results

Overall situation of TME in SKCM

A total of 26,532 cells from 3 CM samples were divided into 6 subgroups (Fig. 1A–C), namely melanocytes accounting for 55.75% (Melanocytes, MLANA, PMEL), T cells accounting for 31.48% (T cells, CD3D), B cells accounting for 4.82% (B cells, MS4A1, CD79A), fibroblasts accounting for 4.56% (Fibroblasts, DCN, COL1A1), endothelial cells accounting for 2.61% (EC, PECAM1, SELE), and macrophages accounting for 0.77% (Macrophages, LYZ, CD68, CD14), with melanocytes being the largest proportion. In the CellChat cell–cell communication network diagram, the thickness of the lines was determined by numerical values, which represented the number of ligand-receptor connections. A more significant number indicated more connections and stronger intercellular interactions. As shown in Fig. 1D, melanoma cells exhibited strong interactions with Fibroblasts, Macrophages, Endothelial Cells (EC), and T cells, suggesting that these four cell types played a crucial role in the tumor microenvironment of skin melanoma. They may directly or indirectly influence the occurrence and progression of melanoma within the tumor. Cell proportions estimated using BayesPrism based on the TCGA melanoma dataset showed that melanoma cells had the highest proportion (Fig. 1E), consistent with the single cell data. The proportions of other cell types also aligned well with the cell composition observed in single-cell dataset. We evaluated the average expression levels of 26 m6A regulatory factors in different cell types (Fig. 1F) and found that some regulatory factors are highly expressed in all cells and could serve as potential therapeutic targets in the TME. Specifically, HNRNPA1 and HNRNPA2B1 were highly expressed in melanocytes compared to other cell types.

Fig. 1
figure 1

Single-cell annotation of skin melanoma and the expression of the m6A gene set in single-cell. A tSNE plot showing the composition of skin melanoma cells; B bar chart of the proportions of cell composition in each sample; C violin plot of cell annotation markers; D interaction network diagram between cells; E TCGA cutaneous melanoma data cell proportions bar chat; F heatmap of the expression of the m6A gene set in various cell types

Full size image

Contribution of m6A-mediated fibroblasts to TME

Next, we analyzed fibroblasts, known as CAFs (cancer-associated fibroblasts) in tumors. First, we re-clustered them into three subgroups (Fig. 2A, B): pan-myCAFs, myofibroblast-like CAFs (TAGLN, ACTA2, MYLK); pan-dCAFs, desmoplastic CAFs (COL1A1, COL3A1); and pan-mCAFs, melanocytes CAFs (PMEL, MLANA). Then, using the NMF unsupervised clustering method, we clustered 26 m6A regulatory factors, resulting in four m6A-associated fibroblast clusters, namely YTHDF3 + CAF-C1 (n = 331), CBLL1 + CAF-C2 (n = 45), IGF2BP3 + CAF-C3 (n = 774), and METTL4-CAF-C4 (n = 60) (Fig. 2C). The selection criteria for clustering rank values and the heatmap can be found in the attachment (Fig. S1A, Fig. S1B).

Fig. 2
figure 2

Display of fibroblast subpopulations and m6A-related fibroblast clusters. A tSNE plot of fibroblast subpopulations; B heatmap of fibroblast subpopulation markers; C NMF clustering heatmap of m6A regulatory factors; D interaction results between m6A clusters and melanoma cells; E, F pseudotime analysis of m6A clusters; G heatmap of differential genes at various differentiation nodes in pseudotime analysis; H GOBP enrichment analysis of differential genes in m6A clusters; I correlation analysis between fibroblast subpopulations and m6A clusters; J expression of each gene set in m6A clusters

Full size image

Cell–cell communication analysis using the CellChat method revealed strong interactions between the YTHDF3 + CAF-C1 and IGF2BP3 + CAF-C3 clusters with melanoma cells (Fig. 2D).

Using monocle pseudotime analysis, we found that the cell differentiation direction was from IGF2BP3 + CAF-C3 to YTHDF3 + CAF-C1 (Fig. 2E, F), with m6A regulatory factors playing a key role in the trajectory of CAFs (Fig. 2G). GOBP functional enrichment analysis showed that the IGF2BP3 + CAF-C3 cluster was enriched in various biosynthetic processes, including secondary metabolite biosynthetic process, pigment biosynthetic process, and phenol-containing compound biosynthetic process, while the YTHDF3 + CAF-C1 cluster exhibited cell chemotaxis functions such as eosinophil chemotaxis and T cell chemotaxis (Fig. 2H). We also found that IGF2BP3 + CAF-C3 was closely related to pan-myCAF myofibroblasts, and YTHDF3 + CAF-C1 was closely related to pan-dCAF desmoplastic fibroblasts (Fig. 2I). Common signaling pathway genes such as collagen, extracellular matrix (ECM), matrix metalloproteases (MMPs), transforming growth factor beta(TGFB), Neo-Angio, Contractile, RAS, and pro-inflammatory genes were expressed differently in each m6A-associated fibroblast cluster (Fig. 2J).

Using the single-cell transcription factor analysis method pySCENIC, we found that the transcriptional activities of three transcription factors—ATF4, MXD3, and POLE4—were higher in IGF2BP3 + CAF-C3 compared to the other three m6A cell clusters, while RORA, KDM5B, ATF6, and ELK3 were higher in YTHDF3 + CAF-C1 compared to the other three m6A cell clusters (Fig. S2.A).

Contribution of m6A-mediated T/MAC cells to TME

First, T cells were isolated and subdivided into four subgroups: CD4 + T, CD8 + T, NK, and Treg (Fig. 3A, B). Then, using the NMF method combined with 26 m6A regulatory factors, T cells were subjected to unsupervised clustering, resulting in four clusters (Fig. 3C, D). The selection criteria for clustering rank values and the heatmap can be found in the attachment (Fig. S1C, Fig. S1D). The four clusters are METTL16 + T-C1 (n = 588), YTHDF1 + T-C2 (n = 143), VIRMA + T-C3 (n = 7052), and METTL4 + T-C4 (n = 444). It can be observed that the number of cells in C3 is the largest, while the proportion of each T cell subgroup within each m6A-related cell cluster does not differ significantly.

Fig. 3
figure 3

Display of T/macrophage subpopulations and m6A-related clusters. A tSNE plot of T cell subpopulations; B heatmap of T cell subpopulation markers; C NMF clustering heatmap of T cell m6A regulatory factors; D proportions of T cell subpopulations in each m6A cluster; E expression of each gene set in T cell subpopulations and m6A clusters; F heatmap of transcription factor differences in m6A clusters within T cell subpopulations; G interaction between m6A clusters and melanoma cells; H NMF clustering heatmap of macrophage m6A regulatory factors; I interaction between m6A clusters and melanoma cells; J heatmap of transcription factor differences in macrophage m6A clusters

Full size image

Through the cellchat cell communication analysis method, it was found that the number of ligand-receptor connections between melanoma cells and m6A-related T cell clusters varies (Fig. 3G).

Using the single-cell transcription factor analysis method pySCENIC, it was discovered that there are significant differences in transcription factor activity (average AUC value) among m6A-related cell clusters C1, C2, C3, and C4 within each T cell subpopulation (Fig. 3F). Additionally, a heatmap showed the average expression levels of these genes in CD4+, CD8+, Treg, and NK cells among the m6A cell clusters, highlighting significant differential characteristics. These include four T cell functional features (T exhaustion score, T effector score, and T escape score), as well as several immune stimulants, inhibitors, and T cell function marker genes (Fig. 3E).

We also isolated macrophages and used NMF clustering analysis to divide them into three classes (Fig. 3H, I): FTO + MAC − C1 (n = 150), METTL16 + MAC − C2 (n = 27), and METTL3 + MAC − C3 (n = 27). The selection criteria for clustering rank values and the heatmap can be found in the attachment (Fig. S1E, Fig. S1F). The pySCENIC analysis revealed that there was a significant difference in transcription factor activity between C1 and the other two m6A-related cell clusters (Kruskal–Wallis test, p < 0.001) (Fig. 3J).

Contribution of m6A-mediated EC Cells to TME

Since strong interactions between EC and melanoma cells were observed in Fig. 1D, we isolated endothelial cells (EC) for separate analysis. Subpopulation analysis was performed, dividing them into three subpopulations: Vascular EC, Melanocytes-like EC, and Fibroblast-like EC (Fig. 4A, B).

Fig. 4
figure 4

Display of endothelial cell subpopulations and m6A-related clusters. A tSNE plot of endothelial cell subpopulations; B heatmap of endothelial cell subpopulation markers; C NMF clustering heatmap of endothelial cell m6A regulatory factors; D interaction between m6A clusters and melanoma cells; E heatmap of transcription factor differences in endothelial cell m6A clusters

Full size image

Using NMF, EC were clustered into two m6A-related cell clusters (Fig. 4C). The selection criteria for clustering rank values and the heatmap can be found in the attachment (Fig. S1G, Fig. SH). The pySCENIC transcription factor analysis revealed significant differences between C1 and C2 (Wilcox-Wallis, p < 0.001) (Fig. 4E).

Construction of m6A cluster prognostic model

Differential analyses were performed on m6A clusters of fibroblasts, T cells, macrophages, and endothelial cells to identify genes for constructing a prognostic model. This model assesses whether the m6A clusters are meaningful and can positively impact tumor prognosis.

Using the Kruskal–Wallis test for differential analysis of the four fibroblast clusters, differential genes were identified as candidate gene sets. These genes were filtered using multivariate Cox regression to construct a prognostic model. The risk score value was calculated, and combined with survival data from TCGA-SKCM, the model’s prognostic capability was validated. Survival analysis revealed highly significant differences between high-risk and low-risk groups (p < 0.0001) (Fig. 5A), with ROC curve values at multiple time points exceeding 0.7 (Fig. 5D).

Fig. 5
figure 5

Survival analysis of various cell types and multi-timepoint ROC curve graphs. A TCGA data survival analysis for fibroblasts; B TCGA data survival analysis for T cells; C TCGA data survival analysis for macrophage cells; D TCGA data multi-timepoint ROC curve for fibroblasts; E TCGA data multi-timepoint ROC curve for T cells; F TCGA data multi-timepoint ROC curve for macrophage cells; G TCGA data survival analysis for endothelial cells (EC); H GSE54467 data survival analysis for fibroblasts; I GSE54467 data survival analysis for T cells; J TCGA data multi-timepoint ROC curve for endothelial cells (EC); K GSE54467 data multi-timepoint ROC curve for fibroblasts; L GSE54467 data multi-timepoint ROC curve for T cells; M GSE54467 data survival analysis for macrophage cells; N GSE54467 data survival analysis for endothelial cells (EC); O GSE54467 data multi-timepoint ROC curve for macrophage cells; P GSE54467 data multi-timepoint ROC curve for endothelial cells (EC)

Full size image

Using the Kruskal–Wallis test, differential analyses were conducted separately on the four T cell clusters and three macrophage clusters to identify differential genes as candidate gene sets. These genes were filtered using multivariate Cox regression to construct prognostic models. The risk score values were calculated, and combined with survival data from TCGA-SKCM, the prognostic capabilities of the models were validated.

Analysis revealed that in the prognostic model constructed using the differential genes of the T cell clusters as candidate genes, the differences between the high-risk and low-risk groups were extremely significant (p < 0.0001) (Fig. 5B), with ROC curve values at multiple time points all exceeding 0.65 (Fig. 5E). Similarly, in the prognostic model constructed using the differential genes of the macrophage clusters, the differences between the high-risk and low-risk groups were extremely significant (p < 0.0001) (Fig. 5C), with ROC curve values at multiple time points all exceeding 0.7 (Fig. 5F).

Using the Wilcox-Wallis test, differential analysis was conducted on the two EC cell clusters to identify differential genes as candidate gene sets. These genes were filtered using multivariate Cox regression to construct a prognostic model. The risk score values were calculated, and combined with survival data from TCGA-SKCM, the prognostic capability of the model was validated.

Analysis revealed that the differences between the high-risk and low-risk groups were extremely significant (p < 0.0001) (Fig. 5G), with ROC curve values at multiple time points all exceeding 0.7 (Fig. 5J).

At the same time, the prognostic models constructed for each cell m6A cluster were revalidated using the GEO dataset GSE54467. The results showed that the m6A cluster survival analysis for fibroblasts was significant (p = 0.018) (Fig. 5H), with ROC curve values at multiple time points all above 0.6 (Fig. 5K); the m6A cluster survival analysis for T cells was significant (p = 0.0028) (Fig. 5I), with ROC curve values at multiple time points all above 0.6 (Fig. 5L); the m6A cluster survival analysis for macrophages was not significant (p = 0.065) (Fig. 5M), with ROC curve values at multiple time points all below 0.6 (Fig. 5O); the m6A cluster survival analysis for endothelial cells was significant (p = 0.015) (Fig. 5N) but the 1 and 2 year ROC curve values at multiple time points were below 0.6 (Fig. 5P).

This indicates that utilizing m6A-mediated fibroblast and T cell subpopulations has great prognostic significance for the tumor microenvironment.

Interaction analysis between melanoma cells and fibroblast m6A clusters

Interaction analysis was conducted between melanoma cells and fibroblast m6A clusters (Fig. 6A). Ligand-receptor gene pairs were divided into high and low expression groups based on their median expression levels in the TCGA-SKCM dataset, followed by survival analysis. The results showed that the survival analysis for COL3A1 was notably different (p < 0.05) (Fig. 6B), and COL3A1 expression was higher in fibroblasts compared to melanoma cells (Fig. 6C). This indicates that the COL3A1 gene is likely an important marker linking melanoma cells and fibroblasts and may be a critical target for advancing melanoma treatment.

Fig. 6
figure 6

Interaction between melanocytes and fibroblast m6A clusters. A Scatter plot of ligand-receptor pairs interactions between melanocytes and fibroblast m6A clusters; B Kaplan–Meier curve for the COL3A1 gene; C violin plot showing COL3A1 gene expression in fibroblast m6A clusters and melanoma cells

Full size image

7. Construction of prognostic models between gene sets of TME cells and TIGER melanoma datasets to probe response group (R) and non-response group (NR) disparity

Prognostic models were constructed for gene sets associated with Fibroblasts, Macrophages, Endothelial Cells (EC), and T cells, and the melanoma dataset from TIGER was used to examine differences between the response group (R) and non-response group (NR). As shown in Fig. S3, there were significant differences in the gene sets of these cell types between the immunotherapy response groups.

Discussion

Thus far, a number of investigations have demonstrated a connection between the pathophysiology of cutaneous melanoma and RNA m6A modification (Liu et al. 2023; Wang et al. 2023a; Ran et al. 2023). Nevertheless, the prospective neoplastic role of m6A-modified single cells has been the subject of only a small number of studies. To our knowledge, this is the first study to systematically investigate m6A modification regulators of primary cell types in the tumor microenvironment (TME) of cutaneous melanoma. Moreover, the composition and characteristics of the TME can profoundly influence treatment responses. For example, a TME rich in immunosuppressive cells like certain macrophages can diminish the effectiveness of immunotherapies. Similarly, dense ECM can hinder the penetration of therapeutic agents, reducing their efficacy (Wang et al. 2023b). We have also discovered, at the single-cell sequence level, the diversity of cell–cell interactions between tumor cells and TME cell subtypes associated with m6A. Through this novel and distinctive viewpoint, we were able to comprehend the ways in which the RNA m6A modification of these various TME cellular constituents impacts the prognosis of individual cutaneous melanoma patients.

The majority of tumor tissue in CM is composed of melanoma cells, which are responsible for tumor expansion. Concurrently, the prognosis of patients and the responsiveness of melanoma cells to treatment are dictated by their heterogeneity. Besides melanoma cells, TME cells, including vascular endothelial cells, infiltrating immune cells, and various types of stromal cells, all contribute to the growth and immune evasion of cutaneous melanoma. Based on the single-cell analysis, we discovered that the TME cells, which included fibroblasts, T cells, macrophages, and endothelial cells, exhibited a wide range of m6A regulatory patterns and a significant amount of communication with melanoma cells. Additionally, the analysis of cellphones revealed that ligand-receptor pairs facilitated communication between melanoma cells and TME cell subtypes associated with m6A.

One of the essential parts of stromal cells, cancer-associated fibroblasts (CAFs), were categorized as pan-myCAFs and pan-dCAFs according to particular molecular characteristics (Galbo et al. 2021). We also discovered melanocytes CAFs, or pan-mCAFs, for the first time. The possible significance of RNA m6A modification in CAFs has not yet been fully explored in many studies. Our study revealed that fibroblasts mediated by m6A had more interactions with melanoma cells than fibroblasts not mediated by m6A. Moreover, fibroblasts exhibited higher expression of desmoplastic CAF Collagen Type III Alpha 1 (COL3A1) in comparison to melanoma cells, and there was a significant difference observed in its survival analysis. Cell interaction analysis revealed that fibroblasts and melanoma cells communicate through the COL3A1-ADGRG1 receptor-ligand pair. Additionally, COL3A1 has a significant prognostic impact on survival in TCGA-SKCM skin melanoma, suggesting that COL3A1 plays a crucial role in tumor progression and metastasis by affecting melanoma cells. It holds potential therapeutic and prognostic value. Future studies could further explore the molecular mechanisms and clinical applications of COL3A1, providing a potential target for melanoma metastasis treatment. Pathogenic variants resulting from COL3A1 mutations are associated with poor prognosis and metastasis of cutaneous melanoma (Zhao and Pritchard 2016). Excessive COL3A1 expression was affiliated with melanoma progression and outcome (Cao et al. 2023; Yuan et al. 2017). Furthermore, COL3A1 is upregulated in Triple-negative breast cancer (TNBC), a rare form of invasive breast cancer, and inhibiting it significantly aided in the metastatic regression of TNBC (Yang et al. 2022). Therefore, it was postulated that RNA m6A modification CAF COL3A1 could lead to an immunosuppressive interaction with melanoma cells, thereby helping the cancer to progress.

In recent years, studies have uncovered the substantial role of RNA m6A methylation in the modulation and development of immune cells (Wu et al. 2024). Obesity-associated protein FTO, an m6A eraser, increased melanoma cell proliferation, migration, and invasion in vitro as well as in vivo in a study by Seungwon Yang et al. and colleagues (Yang et al. 2019). They proposed that the tumorigenicity of melanoma was suppressed by FTO inhibition. Upregulation of FTO also promotes CAF-mediated angiogenesis, a key component in tumor development, according to research by Liao et al. (2024). Molecular clustering of m6A regulators was employed using Non-negative matrix factorization (NMF) which offers advantages over consensus clustering and k-means clustering, including stronger interpretability, sparsity, and greater flexibility. NMF can be combined with other methods incorporating sparse constraints or regularization to enhance clustering performance further. Moreover, since single-cell transcriptomic expression matrices are typically sparse, they align well with the sparsity characteristic of NMF, making it a suitable clustering method. Through NMF analysis, we discovered that all m6A-mediated subtypes of macrophages showed extensive interactions with melanoma cells; however, pySCENIC analysis showed a significant difference in transcription factor activity between FTO + MAC − C3 and the other two m6A-related cell clusters, indicating that m6A-mediated macrophages, particularly for the FTO + MAC − C3 subtype, clearly displayed activation of melanoma tumorigenesis. Furthermore, we discovered that melanoma cells also exhibited extensive interaction with m6A—mediated T cells (CD8+, CD4+, NK, and Treg), B cells, and EC. Moreover, distinct T cell active and inactive features were displayed by m6A-mediated subtypes of the four primary T cell populations. These findings all point to the importance of RNA m6A methylation in immune evasion and the tumor-promoting consequences of macrophages and T cells.

Transcriptional factors (TFs) were analysed at the single-cell level in order to identify gene regulatory networks that are specific to individual cells.In general, each subtype of CAF, T cells, macrophages, B cells, and EC exhibited distinct TF characteristics. For CAFs, IGF2BP3 + CAF-C3 and YTHDF3 + CAF-C1 both exhibited unique TF gene signature, i.e.; ATF4, MXD3, POLE4 and RORA, KDM5B, ATF6, ELK3 respectively. Previous research has demonstrated the connection between m6A modification and the expression of ATF4, RORA, KDM5B, ATF6, and ELK3 (Mo et al. 2023; Han et al. 2021; Huang et al. 2022; Guo et al. 2021; Hua et al. 2021; Zhang et al. 2021), indicating m6A’s function in CAF regulation (Chen et al. 2023b). Moreover, we discovered unique TF traits of m6A-mediated cell subtypes for macrophages, T cells, and EC. To sum up, the TME can be reshaped and reprogrammed by m6A-mediated cell subtypes through their modulation of different TF regulatory networks. Ultimately, melanoma cells were in close communication and connection with these m6A-mediated TME cells, as demonstrated by cell network analysis. Given that either m6A-mediated CAFs or immune cell subtypes demonstrated enhanced interaction with melanoma cells, it is clear that the development of an immunosuppressive tumor microenvironment may be partially impacted by RNA m6A methylation.

With consideration for the complex intrinsic patterns of RNA m6A methylation in TME cells, we compiled in detail the relationships between the scores of these sub-clusters and the immune response and prognosis from the public bulk RNA-seq cohorts. It was evident that patients with varying degrees of TME cell m6A regulator domination had significantly different melanoma prognoses and significantly different immune responses for patients receiving ICB therapy, particularly for macrophages and CAFs. These findings highlight the crucial role that TME m6A plays in the immune response of melanoma patients and warrant further investigation.

To validate these observations, potential wet lab experiments such as MeRip-Seq (m6A RNA immunoprecipitation sequencing) on melanoma cells with perturbed m6A regulators, map changes in m6A methylation sites followed by comparing results to the study’s bioinformatics predictions to confirm overlap in target transcripts. Additionally, utilizing CRISPR-Cas9 to target silencing specific m6A regulators identified as key players in melanoma clusters then validate knockdown efficiency via qRT-PCR (mRNA) and Western blot (protein) can aid in validation of m6A regulator expression and function (Du et al. 2021). In the future, we aim to include wet lab experiments in an effort to comprehensively validate observations.

The low depth of scRNA-seq and the insufficient sample size were the main limitations of our preliminary study, and additional patient verification is required to confirm our findings. A potential source of bias in our study’s clustering method is the fact that, when compared to bulk RNAseq, the scRNA-seq of some m6A regulators in melanoma would mostly be minor and have more zero observations. However, since our study focuses on cutaneous melanoma samples, we utilized three samples from CM single-cell dataset (GSE215121). Single-cell data differs from conventional bulk RNA sequencing, as it allows for gene expression analysis at the single-cell level. Even with a small number of samples, it provides rich information, revealing cellular heterogeneity and overcoming the limitations of averaging cell populations in traditional methods. Additionally, its technical advantages enable efficient single-cell isolation and capture, reducing sample requirements. Through mRNA amplification, sufficient data can be obtained even from minimal RNA amounts. In the future, to broaden our investigation, we aim to include more cutaneous melanoma samples (e.g., GSE81383) in our research. We also plan to expand our sample types to include acral melanoma samples (GSE189889) or melanoma brain metastases and leptomeningeal metastases (GSE174401), enabling a more comprehensive study of melanoma. Moreover, the scRNA-seq analysis still offers us a fresh perspective on the properties of m6A methylation regulators in different TME single cells to lessen melanoma assortment, which is crucial in advanced clinical practice.

Conclusion

Taking into account the intricate patterns of RNA m6A methylation within the cells of the tumor microenvironment (TME), it became evident that variations in m6A regulatory dominance in TME cells resulted in significant differences in the prognosis of cutaneous melanoma (CM) patients and showcased distinct immune responses, particularly among cancer-associated fibroblasts (CAFs) highlighting COL3A1 gene as a critical link and potential therapeutic target in CM. These findings provide new perception on the crucial role of m6A modifications in the TME offering valuable insights into potential biomarkers and targeted treatment strategies for CM.