3.1 Establishment of the apoptosis-related tumor microenvironment signature (ATM) and its independent prognostic impact in melanoma patients

To investigate the characteristics of the tumor microenvironment in patients with distinct overall apoptosis statuses, we calculated the Apoptosis score by a single-sample gene set enrichment analysis (ssGSEA) [30] and validated our apoptosis score using experimental data from four distinct datasets with well-documented apoptotic states (GSE196610, GSE147256, GSE137574, GSE242151) (Fig. S1a). Then, weighted gene co-expression network analysis (WGCNA) was used to obtain the apoptosis score associated modules. We identified a module that exhibited a collective biological function significantly associated with apoptosis score, independent of other clinical characteristics (e.g., age, gender, stage, Clark level, Breslow depth, mitotic count rate, melanoma ulceration indicator, etc.). This selection was based on rigorous analysis of the data, as illustrated in Figs. 1b and S1b. The genes comprising this module were found to be significantly enriched in signaling pathways related to immunoglobulin complexes, antigen binding, Fcgr activation, and the initial triggering of complement, as determined through Geneset Enrichment Analysis (GSEA) (Fig. 1c). Moreover, we observed that the high module score group with a higher apoptosis score exhibited the highest hallmark scores in seven immune-related hallmarks, demonstrating that these module genes reflect immune-related changes in the tumor immune microenvironment (Figs. 1c and S2). Next, we constructed an apoptosis-related model comprising 19 genes using lasso regression, utilizing 216 genes from the apoptosis-related ME_black module. The coefficients of the model genes were optimized using multi-cox regression (Fig. 1d). Subsequently, the apoptosis-related tumor microenvironment signature (ATM) was calculated using the predict function of the survival R package [31]. The ATM score, which exhibited a negative correlation with the apoptosis score (R = -0.37, P = 1.1e-15), was utilized for stratifying melanoma patients (Fig. 1f). The results consistently indicated a negative association between ATM score and apoptosis score across multiple datasets (TCGA-SKCM, GSE19234, GSE22153, GSE54467, GSE65904) (Figs. 1e and S3).

To assess the impact of ATM on the prognosis of melanoma patients, we conducted survival analysis and performed multivariate Cox regression model analysis using the TCGA-SKCM cohort. The results revealed that subgroups with low ATM scores exhibited improved overall survival (OS) and ATM serves as an independent prognostic factor even when considering other potential confounding factors (e.g., age, gender, stage, Clark level, Breslow depth, mutation status) (Figs. 1f and S4a).

To ascertain the prognostic capability of the ATM in melanoma patients, we gathered four independent cohorts (GSE19234, GSE22153, GSE54467, GSE65904) for subsequent survival analysis. The results consistently demonstrated that melanoma patients in the low ATM group exhibited improved overall survival (OS), thus confirming the predictive value of ATM (Fig. 1g). Furthermore, we performed multi-cox analysis on these four datasets, encompassing a total of 389 patients. Remarkably, the results consistently revealed that ATM can serve as an independent prognostic factor for predicting the prognosis of melanoma patients, surpassing other factors such as age, gender, tissue location, mutation status, and stage (Figs. 1h and S4b). The consistent outcomes across multiple independent cohorts strengthen the validity and generalizability of ATM as a reliable predictor of melanoma patient outcomes.

3.2 Plasma cells as pivotal contributors in model interpretation: implications for enhanced prognosis and immune response

To uncover the molecular characteristics of the TME associated with apoptosis in melanoma patients, we conducted differential expression analysis comparing the high and low ATM groups within the bulk samples of the TCGA-SKCM cohort. This analysis revealed distinct patterns of gene expression, particularly in relation to chemokines such as CXCL13, CXCL11, CXCL10, as well as other factors including PTPRC, JCHAIN, IGKJ5, and IFNG. Notably, these genes exhibited significantly higher expression levels in the low ATM group, indicating a strong potential for antibody secretion within this subgroup (Fig. 2a). Next, we conducted a GSEA analysis to investigate the mRNA-level differences between the high and low ATM groups. Remarkably, the low ATM group exhibited significant enrichment in immune-related signaling pathways. Notably, these pathways included the B cell receptor signaling pathway, B cell activation, immunoglobulin complex, positive regulation for B cell activation, etc. (Figs. 2b and S1d). These findings highlight the pronounced involvement of B cell-related immune pathways in the low ATM group, suggesting a heightened B cell-mediated immune response in this subgroup. Conversely, the high ATM group primarily displayed enrichment in pathways associated with epidermal development, keratinization, keratinocyte differentiation, and skin development (Figs. 2b and S1c). We also analyzed the differential hallmark score between the two ATM score groups in the TCGA-SKCM cohort. We observed that the low ATM score group exhibited significantly higher scores in six immune-related hallmarks, demonstrating that ATM score reflects immune-related changes in the tumor immune microenvironment (Fig. S6). Collectively, these results provide valuable insights into the distinct molecular signatures and biological processes associated with different ATM groups. Besides, to investigate clinical features associated with ATM in melanoma patients, we evaluated ATM across different clinical features, showing that patients with Age < 60, Clark levels I–III and IV had lower ATM relative to those with Clark level V (Fig. S5).

Fig. 2

Plasma cells as pivotal contributors in model interpretation: implications for enhanced prognosis and immune response a Differential expression genes between the high ATM group and the low ATM group. b The top 10 enriched Gene Ontology (GO) signaling pathways associated with differential expression genes between the high ATM group and the low ATM group. c, d UMAP plot colored by various cell types or ATM score in GSE120575 and GSE123139 cohorts. e Kaplan–Meier curves for patients with high and low plasma ratio in the GSE12575 cohort show that patients with high plasma ratio(red) exhibited better overall survival. f The box plot illustrates that immunotherapy responders exhibited a higher plasma ratio in the GSE120575 cohort. g, h Violin diagram shows the ATM score across various cell types in GSE120575 and GSE123139 cohorts. i Spatial UMAP plot colored by Malignant signature score, plasma cell signature score, and ATM score in 10x Genomics Human Melanoma Spatial Gene Expression Data. j The correlation between ATM score and plasma cell signature score in 10x Genomics Human Melanoma Spatial Gene Expression Data

In addition to the insights gained from bulk sample analysis, we extended our investigation to explore the ATM at a single-cell resolution to better understand intra-tumor heterogeneity. Employing graph-based principal component clustering combined with marker-based annotations, we classified cells from the GSE120575 dataset into 14 clusters based on the RNA_snn_res.0.55 and cells from the GSE123139 dataset into 10 clusters based on the RNA_snn_res.0.15, revealing distinct cellular subpopulations (Figs. 2c, d and S7a, b). Notably, plasma cells exhibited the lowest ATM score, further highlighting their potential role as pivotal contributors in our model interpretation (Fig. 2c, d, g, h). Moreover, analyzing the spatial transcriptomic data of Human Skin Melanoma obtained from the official 10X platform consistently demonstrated that plasma cells had the lowest ATM score (Fig. 2i). Interestingly, plasma cells adjacent to tumor cells demonstrates its beneficial spatial location to play its further anti-tumor immunity.

Furthermore, within the GSE120575 cohort, we found that a higher plasma cell ratio was associated with improved survival and enhanced immune response (Fig. 2e, f). Through further analysis, it was determined that this particular cell population predominantly exhibited IgG expression, with a subset displaying IgA (Fig. S7e, f). Moreover, previous investigations corroborated the correlation between IgG+ plasma cells and enhanced survival [38, 39].

3.3 Plasma cells orchestrate a wide spectrum of immune activation state by interacting with cellular constituents of the immune microenvironment

Further, cell communication analysis provides the opportunity to study cell-cell interactions based on ligand-receptor binding. Cellchat is thus utilized to conduct an in-depth analysis of the interplay between plasma cells and other cells, elucidating their pivotal role in orchestrating an effective anti-tumor immune response. 10 individual samples exhibiting superior survival rates and the highest plasma ratio from the GSE120575 cohort were taken for an in-depth examination of cellular interactions. The findings unequivocally demonstrated the significant involvement of plasma cells as primary senders and influencers within the MIF signaling pathway network. Remarkably, plasma cells also assume a crucial role as senders and influencers within the SEMA4 signaling pathway network (Fig. 3a, b). It is widely acknowledged that the cytokine macrophage migration inhibitory factor (MIF), a proinflammatory cytokine, plays a key role in inflammatory diseases with chemokine receptors CXCR2, CXCR4, and CD74/CXCR4 complexes as functional receptors. By activating CXCR2, CXCR4, or CD74/CXCR4 complexes, MIF displays chemokine-like functions and acts as a major regulator of inflammatory cell recruitment such as monocytes, T cells, and B cells [40,41,42]. According to our findings, plasma cells within the MIF signaling pathway exhibit robust chemotactic capabilities towards B cells, pDC and Tfh (Fig. 3a, c), leading to their migration towards the effector site and fostering opportunities for cellular interactions. Moreover, plasma cells in the SEMA4 signaling pathway engage in interactions with B cells, promoting their aggregation through SEMA4D and CD72 ligand receptors [43] (Fig. 3b, c). Notably, our results highlight that B cells play a central role as the primary senders in the MHC-II signaling pathway, facilitating antigen presentation to M1, Monocyte, pDC, Tfh, and Treg (Fig. 3d). We utilize the ssGSEA [30] to estimate the score of MHC.II gene signatures (Fig. 3e) from Thorsson et al. [34]. It is noteworthy that germinal center B cells (GCB) demonstrate the most pronounced MHC-II signature within the B cell subpopulation (Fig. 3e). Several studies have indicated that follicular helper T (Tfh) cells are essential for germinal center formation, affinity maturation, and the development of most high-affinity antibodies and memory B cells [44]. Recent study has also highlighted that pDCs were identified as a component of the T cell zone of TLS, which are major regulators of adaptive antitumor immunity [45]. Collectively, these findings imply that plasma cell chemotaxis facilitates the aggregation of B cells and Tfh which contribute to germinal center formation, allowing GCB to present antigens to Tfh, thereby contributing to subsequent B cell differentiation and the formation of long-lived plasma cells.

Fig. 3

Plasma cells orchestrate a wide immune activation state by interacting with cellular constituents of the immune microenvironment a, b Heatmaps showing the relative importance of each cell population based on the computed network centrality measures of MIF signaling pathway network (a) and SEMA4 signaling pathway network (b) in GSE120575 cohort. c The significant ligand-receptor pairs that contribute to the signalings (MIF, SEMA4) sending from plasma to other cell populations. The dot color and size represent the calculated communication probability and p-values. p-values are computed from one-sided permutation test. d Heatmaps showing the relative importance of each cell population based on the computed network centrality measures of MHC-II signaling pathway network in GSE120575 cohort. e Violin diagram shows the MHC.II signature expression level across various cell types in GSE120575 B cell subpopulation. f The significant ligand-receptor pairs that contribute to the signalings (APRIL, BAFF) sending from other cell populations to plasma cell. g–j Heatmaps showing the relative importance of each cell population based on the computed network centrality measures of APRIL signaling pathway network (g), BAFF signaling pathway network (h), CADM signaling pathway network (i), and LAMININ signaling pathway network (j) in GSE120575 cohort. k Bar plot shows the cell types with significantly different cell proportions between two ATM groups which were revealed by the Wilcox test in TCGA-SKCM cohort and 3 validation cohorts (GSE65904, GSE19234, GSE54467). Asterisks denoted p-value. l Correlation between ATM score and immune checkpoints in TCGA-SKCM cohort and 4 validation cohorts (GSE65904, GSE19234, GSE54467, GSE22153). m The correlation between ATM score and TMB, TCR, and BCR in TCGA-SKCM cohort. Asterisks denoted p-value. (“*”p < 0.05; “**”p < 0.01; “***”p < 0.001; “****”p < 0.0001; ns was the abbreviation of no significance)

Additionally, plasma cell chemotaxis also influences CD8+Tem migration (Fig. 3a, c), thereby enabling CD8+Tem cells to exert their anti-tumor immune function at the effector site. Moreover, plasma cells exhibit chemotactic effects on monocytes (Fig. 3a, c). Our results also indicate that within the APRIL and BAFF signaling pathways, M1, M0, and pDC cells reciprocally interact with plasma cells as senders through the TNFSF13B-TNFRSF17 ligand-receptor interaction (Fig. 3f–h), thereby promoting the longevity of plasma cells [46]. Furthermore, plasma cells, as sender and receiver, played an important role in a series of cell adhesion-related signaling pathways (e.g. CADM, LAMININ) (Fig. 3i, j). Numerous investigations have provided evidence of the involvement of CADM1 (Cell adhesion molecule 1), an immunoglobulin superfamily member, in cell-cell interactions [47]. Furthermore, molecules such as LAMININ, which are part of the ECM-Receptor signaling pathway, offer specific opportunities for interactions between cells and the ECM. These intricate interactions exert direct or indirect control over cellular activities encompassing adhesion, migration, differentiation, and proliferation [48]. Notably, research has highlighted the significance of the adhesion phenotype of plasma cells within the B-cell compartment, contributing to their differentiation and homing [49]. Recent studies have also suggested that chemoattractant cytokines (chemokines), in conjunction with tissue-specific adhesion molecules, coordinate the migration of antibody-secreting cells (ASCs) from lymphoid tissues, where they undergo antigen-driven differentiation, to effector tissues [50]. Taken together, these suggest that this adhesion phenotype of plasma cells is in a “being ready” state for further differentiation or homing to target effector tissues to play anti-tumour immune effects.

These findings highlight the antitumor immune properties of plasma cells themselves, as well as their involvement in coordinating a broad spectrum of immune activation states through interactions with other cells within the immune microenvironment. Consistent results were observed in the immune infiltration analyses of the TCGA-SKCM, GSE54467, GSE65904, and GSE19234 datasets through cibersort [51]. Specifically, the low ATM group exhibited higher proportions of plasma cells, Memory B cells, M1 macrophages, activated CD4+ memory T cells, and CD8+ T cells, indicative of an immunoactivated state. In contrast, the high ATM group displayed higher proportions of M0 macrophages and M2 macrophages, suggesting an immunosuppressive immune microenvironment (Fig. 3k).

In order to investigate the association between ATM and immune checkpoints, we conducted a comprehensive analysis of the correlation between stimulatory immune checkpoint expression and ATM in TCGA-SKCM cohort, as well as four additional GEO cohorts (GSE65904, GSE22153, GSE54467, GSE19234). Remarkably, our findings consistently revealed a negative correlation between ATM and the majority of immune checkpoints (Fig. 3l, Table S5). Notably, the stimulatory immune checkpoints (CD40, IL2RA, CD96, ICOS, CD28, TNFRSF4, BTN3A1) exhibited a consistent negative correlation with ATM across all five datasets. These findings suggest a potential regulatory relationship between ATM and the expression of key immune checkpoints, indicating their coordinated involvement in shaping the tumor immune microenvironment. Subsequently, in order to elucidate the potential mechanism underlying the role of ATM in immunotherapy, we conducted an examination of the association between established immunotherapy biomarkers including tumor mutation burden (TMB), T-cell receptor (TCR), B-cell receptor (BCR), and ATM. We observed the negative correlation between ATM and these biomarkers, indicating a potential interplay between ATM and the immunotherapy response (Fig. 3m).

3.4 ATM serves as a predictor of immunotherapy efficacy

The aforementioned analysis reveals a robust anti-tumor immune response orchestrated by plasma cells in collaboration with other immune cell populations. To investigate the association between ATM and immunotherapy efficacy and prognosis, multiple datasets with anti-PD-L1/PD1/CTLA4 cohorts were collected in the study. The in-house anti-PD1 treatment melanoma patients cohort, two public melanoma patients cohorts with immune checkpoint therapy (the Riaz N cohort [21]/GSE91061: Anti-PD1-treated advanced melanoma (Nivolumab)), and the Van Allen, E. M cohort [22]: Anti–CTLA4-treated metastatic melanoma) and other cohorts with immune checkpoint therapy (the Balar AV cohort [23]/IMvigor210: Anti-PD-L1-treated locally advanced and metastatic urothelial carcinoma and the Braun DA cohort [25]: anti-PD-1-treated advanced clear cell renal cell carcinoma) were taken for analysis.

In Braun DA cohort, Balar AV cohort, and in-house cohort (579 patients), immunotherapy responders were shown to have lower ATM (p-value = 0.016 in the Braun DA cohort; p-value = 0.0043 in the Balar AV cohort; p-value = 0.0081 in the in-house cohort) and the results of the chi-square analysis showed that the proportion of responder would be higher in the low ATM group (p-value = 0.033 in the Braun DA cohort; p-value = 0.058 in the Balar AV cohort; p-value = 0.012 in the in-house cohort) (Fig. 4a, b). Three immunotherapy datasets showed better OS in the group with low ATM, and the Braun DA cohort showed better OS and PFS in the group with low ATM (Fig. 4c, d). All these results exhibited that ATM has a robust independent prognostic ability and predictive power of immunotherapy efficacy.

Fig. 4

ATM serve as a promising predictor of immunotherapy efficacy a Differentiation of ATM in response and non-response groups in two independent immunotherapy cohorts (Braun DA and Balar AV cohorts) and In-house cohort. b The proportion of response/non-response patients in high- and low-ATM groups in the Braun DA and the Balar AV cohorts and In-house cohort. c, d Kaplan–Meier curves for patients with high ATM and low ATM in four independent immunotherapy cohorts (the Van Allen, E. M; the Balar AV; Braun DA; and the Riaz N cohorts) and In-house cohort show patients with lower ATM (blue) exhibited better overall survival and/or progression-free survival

3.5 Influence of apoptosis-related tumor microenvironment signature (ATM) on anti-cancer drug response

Massive studies have demonstrated that apoptosis can affect drug response mostly in chemotherapy. In order to explore the potential efficacy of influencing other drugs’ sensitivity by apoptosis and develop novel therapeutic hypotheses, we comprehensively depicted the associations between apoptosis-related tumor microenvironment and drug response. We calculated the correlation between the ATM and imputed drug data of drugs in CTRP [28, 29].

A total of 51 genes are targeted by 46 drugs including 6 Food and Drug Administration (FDA)-approved drugs, 10 clinically used drugs, and 30 probes that are associated with ATM, most of which are targeted therapy (Fig. 5a). Among these drugs, the area under the curve (AUC value) of 45 drugs is negatively correlated with ATM which means that these drugs are more sensitive in patients with low ATM. CHIR-99021 which is the GSK3B inhibitor is negatively correlated with ATM indicating that it is more sensitive in patients with low ATM with better prognosis (Figs. 4b and 5a). A previous study has also demonstrated that it may function as an antagonist of MYC degradation pathways to transiently elevate MYC levels and then confer chemosensing within a narrow window [52].

Fig. 5

Influence of apoptosis-related tumor microenvironment signature (ATM) on anti-cancer drug response a The drug names, targeted gene symbol, and drug type of the three classes of drugs whose drug response correlates with the ATM. b Signaling paths targeted by drugs whose drug sensitivity correlates with the ATM. Orange (negative correlation) or blue (positive correlation). c Differentiation of AUC value in high- and low-ATM groups in Apoptosis inducer agents (SZ4TA2; gossypol) and agents targeting VEGF signaling (linifanib; tivozanib; quizartinib; vandetanib) was revealed by the Wilcox test. Asterisks denoted p-value. (“*”p < 0.05; “**”p < 0.01; “***”p < 0.001; “****”p < 0.0001; ns was the abbreviation of no significance)

Notably, six of these drugs are majoring in the VEGF signally and the AUC of targeting VEGF signaling drugs such as linifanib, tivozanib, quizartinib, and vandetanib are significantly higher in the group with low ATM (Fig. 5b, c). In addition, two apoptosis inducer agents (SZ4TA2, gossypol) are also significantly higher in the group with low ATM (Fig. 5c).

3.6 Apoptosis-related tumor microenvironment signature (ATM) is significantly associated with tertiary lymphoid structures (TLS), exhibiting stronger patient stratification ability compared to classical “hot tumors”

To explore the association between ATM score and tumor categorization, we compiled genes related to cold and hot tumors from the study by Dong Wang et al., encompassing 12 hot tumor-related genes (CXCL9, CXCL10, CXCL11, CXCR3, CD3, CD4, CD8a, CD8b, CD274, PDCD1, CXCR4, and CCL5), and 3 cold tumor-related genes (CXCL1, CXCL2, and CCL20) [53]. In the TCGA-SKCM cohort, our analysis revealed that the ‘Hot’ group exhibited significantly improved overall survival, aligning with prior research findings. Notably, we further stratified patients based on ATM score and Hot tumor signature. Intriguingly, we observed that irrespective of tumor type (hot or cold), the low ATM group consistently displayed superior survival outcomes compared to the high ATM group. Most importantly, the high ATM score has the capability to identify high-risk patient groups even within ‘hot’ tumors, thus offering crucial insights for patient stratification compared to the conventional ‘hot’ and ‘cold’ tumor classifications. (Fig. 6b, c).

Fig. 6

Apoptosis-related tumor microenvironment signature (ATM) is associated with tertiary lymphoid structures (TLS), exhibiting stronger patient stratification ability compared to classical “hot tumors” a Kaplan–Meier curves for group hot and cold patients based on the average expression of 12 hot tumor–related genes (CXCL9, CXCL10, CXCL11, CXCR3, CD3, CD4, CD8a, CD8b, CD274, PDCD1, CXCR4, and CCL5) in the TCGA-SKCM cohort show that group hot patients (red) exhibited better overall survival. b Kaplan–Meier curves for 4 groups of patients (Cold-High group, Cold-Low group, Hot-High group, and Hot-Low group) based on the average expression of 12 hot tumor and ATM score in the TCGA-SKCM cohort. c sankey diagram for four groups of patients stratified by hot tumor score and ATM score in the TCGA-SKCM cohort. d, e The top 10 enriched Gene Ontology (GO) signaling pathways associated with differential expression genes between the Hot-ATM low group and Hot-ATM high group (d) or between the Cold-ATM low group and Cold-ATM high group (e). f–h The correlation between ATM score and TLS score in TCGA-SKCM cohort and 4 validation cohorts (GSE65904, GSE19234, GSE54467, GSE22153). (ssGSEA was applied to calculated the TLS_9 score (CD79B, CD1D, CCR6, LAT, SKAP1, CETP, EIF1AY, RBP5 and PTGDS); TLS_12 score (CCL2, CCL3, CCL4, CCL5, CCL8, CCL18, CCL19, CCL21, CXCL9, CXCL10, CXCL11 and CXCL13); TLS_29 score (IGHA1, IGHG1, IGHG2, IGHG3, IGHG4, IGHGP, IGHM, IGKC, IGLC1, IGLC2, IGLC3, JCHAIN, CD79A, FCRL5, MZB1, SSR4, XBP1, TRBC2, IL7R, CXCL12, LUM, C1QA, C7, CD52, APOE, PTLP, PTGDS, PIM2, and DERL3 genes)

Furthermore, we observed that among patients with either hot or cold tumors, the enriched signaling pathways of the differential gene functional analysis between patients with low ATM scores and those with high ATM scores contained a greater representation of B cell, cell-cell adhesion, tissue remodeling and humoral immune-related pathways (e.g., B cell receptor signaling pathway, humoral immune response mediated circulating immunoglobulin, adaptive immune response based on somatic recombination of immune receptors, regulation of B cell differentiation, lymphocyte costimulation, regulation of tissue remodeling, positive regulation of leukocyte cell-cell adhesion, humoral immune response, etc.) (Fig. 6d, e). This is consistent with the coordinated anti-tumor immune response underscored in our results, focusing on the role of plasma cells. And the observation of plasma cells in a highly adhesive state suggests their readiness for further differentiation or homing to specific effector tissues, enabling them to effectively carry out their anti-tumor immune effects. We further conducted an analysis of the correlation between ATM score and the Tertiary Lymphoid Structure (TLS) signature [54,55,56], and in all five datasets, we consistently observed that a low ATM score is associated with a higher TLS signature (Fig. 6f–h).