1 Construction of scRNA-seq atlases of melanoma.The results identified 28 distinct cell clusters, including immune, stromal, and melanoma cells (Fig. 1A). Further analysis revealed that (Fig. 1B) the immune cells could be categorized into T cells (IL7R, CD8B/3D, TRAC, and NKG7), B cells (CD79A, BANK1, and MS4A1), and monocytes (LYZ and CD14/68). Additionally, the stromal cells were classified into endothelial cells (PECAM1 and VWF) and fibroblasts (COL1A2/3A1). Notably, these cell cluster compositions, particularly the melanoma cells and T lymphocytes, were distinct in different melanoma samples (Fig. 1C). Subsequently, we examined differentially expressed genes (DEGs) among different cell types and marked some of the key genes in melanoma (Figs. 1D–E).
2 The communication networks among melanoma, immune and stromal cell subtypes.
Following the classification of melanoma samples at the cellular level, our objective is to establish a communication network between subtypes of melanoma, immune, and stromal cells. The comprehensive outcomes of the CellChat analysis are visually represented in Fig. 2 through Sankey diagrams, dotplot, and chordal graphs. Subsequent examination of incoming communication among these cells unveiled a shared pattern between stromal and melanoma cells (Figs. 2A-B). Furthermore, Figs. 2C-E provide detailed insights into potential molecular interactions. We can observe functional connections between these clusters of cells.
3 Copy number variation (CNV) analysis of melanoma cell subtypes
The melanoma cells were divided into ten clusters using cluster analysis and were visualized according to the tumor samples (Fig. 3A). Moreover, we visualized the top ten marker genes for each melanoma cell cluster (Fig. 3B). Subsequently, the CNV status from various cell types was determined utilizing T and B lymphocytes as reference controls using InferCNV analysis (Fig. 3C). Consistent with the findings of Zhang et al. in pre- and post-treatment samples of a single patient who received immunotherapy, variations in CNV were found on chromosome 411. The CNV observed in tumor specimens displayed significant heterogeneity, with varying degrees of CNV accumulation among various patients and tissue types (Fig. 3D). Subsequently, we categorized all melanoma cells into high- and low-CNV groups. Compared with cutaneous melanoma (CM1) and melanoma cells derived from cutaneous melanoma lymphatic metastasis tissues (CM1-lym) as well as pre- (AM3-pre) and post-immunotherapy acral melanoma tissues (AM3-post), CM1-lym and AM3-post tissues exhibited significantly elevated CNV levels, suggesting a more aggressive phenotype (Fig. 3E). The Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses for the DEGs between both CNV tumor samples showed the related pathways and functions (Fig. 3F). Figure 3G–H depicts that differences in gene expression patterns were significant between CM1 and CM1-lym as well as AM3-pre and AM3-post for melanoma cells.
4 CAF classification and pseudo-temporal trajectories analysis of CAFs
The CAFs were classified into six distinct cell subtypes (CAFs1–6) depending on biological functions, cellular interactions, marker genes, and spatial distribution in the TME. These subtypes included vascular CAF (vCAF), pericyte, matrix CAF (mCAF), inflammatory CAF (iCAF), tumor-like CAF (tCAF), dividing CAF (dCAF), antigen-presenting CAF (apCAF), and epithelial-like CAF (epi-CAF) (Figs. 4A–B)[16]. Figure 4C shows the respective marker genes of different CAF subtypes. Subsequently, we found that the differentiation trajectory of different CAF cell subtypes was different through pseudo-temporal trajectories analysis (Figs. 4D–E). Then, we applied gene set variation enrichment analysis (GSVA) to reveal CAF functions, demonstrating the related pathways of CAFs (Fig. 4F). Besides, the markers of the four reclassified CAF subtypes based on the pseudotime trajectories analysis were shown in heatmaps (Fig. 4G). The main transcription factors for the CAF subtypes were evaluated (Fig. 4H). A multivariate COX regression analysis for the CAFs1–6 identified the risk factors depending on the individual gene expression levels (Fig. 4I).
5 Functional analysis of T cell subtypes.
Herein, we also focused on the immune cell subtypes and classified them into five groups of natural killer (NK) cells, seven groups of CD4 T cells, four groups of CD8 T cells, and two groups of cycling T cells depending on the marker gene expression levels (Figs. 5A–C). Moreover, we used Z-score to show the expression of major genes in T cells in melanoma samples (Fig. 5D). The percentage of T cell clusters between AM3-pre and AM3-post tissues was represented on a proportion chart (Fig. 5E). The scRNA-seq data of gene expression levels in T cells were collected from the GEO database for 33 melanoma patients with immunotherapy resistance (GSE115978). Furthermore, we divided the immunotherapy-related T cells into two clusters: CD4/8 T cells (Fig. 5F). Meanwhile, we extracted 706 marker genes associated with ICP therapy from these T cells, 509 markers of T cells from 11 melanoma tissues, and eventually identified 33 IRTGs (Fig. 5G).
6 The IRTG prognostic signature construction and validation
The results indicated 33 prognostic genes mainly enriched in some immune-related pathways using a KEGG analysis (Figs. 6A–B). The LASSO and COX regression analyses were conducted for these 33 genes to identify six signature genes: IL27RA, PIM2, PRDM1, LTB, GBP5, and PDE4B (Figs. 6C–E). The risk score was determined depending on the signature gene expression levels alongside their corresponding risk coefficient value as follows: [IL27RA expression level × (–0.254915516014614)] + [PIM2 expression level × (–0.333557686339427)] + [PRDM1 expression level × (0.61377881693067)] + [LTB expression level × (0.195208657137769)] + [GBP5 expression level × (–0.443593637428344)] + [PDE4B expression level × (–0.326705622270952)]. Depending on the risk score, 469 SKCM patients acquired from the TCGA database were allocated into high- and low-risk groups. The low-risk groups had higher ICP gene expression levels (Fig. 6F) as well as overexpressed IL27RA, PIM2, PRDM1, LTB, GBP5, and PDE4B, indicated by a heatmap (Fig. 6G). Furthermore, the high-risk group exhibited a heightened death incidence (Fig. 6H). Multivariate COX regression analysis elucidated that age, sex, tumor stage, and risk score could serve as independent predictive factors (Fig. 6I). The Kaplan-Meier (KM) curves for the training cohort from TCGA manifested that low-risk patients had a significantly extended overall survival (OS; p < 0.001), with 1-, 3-, and 5-year AUC values of 0.654, 0.659, and 0.683, respectively (Fig. 6J). The outcomes from the four validation cohorts, GSE54467 (p = 0.036, 1-, 3-, and 5-year AUC = 0.451, 0.559, and 0.595, respectively), GSE65904 (p < 0.001, 1-, 3-, and 5-year AUC = 0.615, 0.656, and 0.634, respectively), GSE22153 (p = 0.012, 1-, 3-, and 5-year AUC = 0.630, 0.658, and 0.539, respectively), and GSE59455 (p = 0.113, 1-, 3-, and 5-year AUC = 0.694, 0.583, and 0.619, respectively), indicated that patients with a low-risk score experience longer disease-free survival (DFS).
7 TME and immunotherapy efficacy between both risk groups
The study revealed a direct correlation between risk scores and M2/M0 macrophages and both resting memory CD4 T cells and NK cells. Conversely, an adverse association was found between risk scores and γδ T cells, follicular helper T cells, CD8 T cells, activated memory CD4 T cells, plasma cells, and M1 macrophages (Fig. 7A). Furthermore, the six signature genes were related to several instances of immune cell infiltrations (Fig. 7B). Figure 7C shows that the low-risk group exhibited elevated stromal and immunological scores. Further, the low-risk group that underwent various types of ICP inhibition treatment had significantly increased immunophenoscores (IPSs; Fig. 7D), indicating a more favorable immunotherapeutic response. The low-risk group had increased ICP gene expression levels (Fig. 7E), suggesting a possibility for enhanced responsiveness to immunotherapy. Moreover, we further confirmed the risk score effectiveness in the prediction of ICI responses in the PRJEB25780, iMvigor210, PRJEB23709, and GSE35640 cohorts. Patients with stable disease (SD)/progressive disease (PD) had a lower risk score than those with high risk (Fig. 7F).
8 The scRNA-seq analysis of the six signature genes
Herein, we observed the six signature gene expression levels in many immune cell types in ten SKCM immunotherapy-related single-cell sequencing datasets. The six signature genes were mainly expressed in immune cells compared to stromal cells (Fig. 8A). Furthermore, Figs. 8B–J show the signature genes significantly expressed in the CD4/8 T cells.
9 PED4B expression positively correlates with tumor CD8 + T cell infiltration in SKCM
Differential expression levels were observed in the six signature genes when comparing SKCM patients to normal controls (Fig. 9A). Patients exhibiting elevated PED4B expression have a prolonged OS, as the KM analysis shows (Fig. 9B). Interestingly, pathway analysis based on PDE4B expression elucidated that PDE4B was mainly involved in pathways related to programmed cell death (Fig. 9C). Combined with the localization of PDE4B in single cells (Fig. 8), we speculated that PDE4B might be involved in the CD8 + T cells effector function. The TCGA database suggested that high PDE4B expression was accompanied by higher CD8 + T cell infiltration (Fig. 9D). To validate our conjecture further, we performed multiple immunofluorescence (IF) staining of PDE4B and CD8 using tissue microarrays. Excitingly, there was a co-localized expression of PDE4B and CD8 in tumor tissues, and PDE4B and CD8 expressions were positively correlated (Fig. 9F). Collectively, our findings indicate that PED4B was correlated with increased CD8 + T cell infiltration in SKCM.