Adenosine signaling as a prognosis factor of melanoma

To identify effective biomarkers and therapeutic targets that could be used to improve the efficacy of adenosine signaling-targeted treatment, we focused on a set of adenosine signaling targets from ongoing clinical trials to evaluate patient survival (Table S1). In addition to expand the pool of candidate biomarkers, we collected these highly actionable targets and other adenosine signaling-related genes from a literature review (Table S4) [46,47,48,49]. To characterize the roles of adenosine signaling in melanoma, we (1) applied nonnegative matrix factorization (NMF) algorithms to classify each patient into an adenosine signaling subtype group, (2) modified and derived an adenosine signaling score (ADO-score) model, and (3) validated the model in 7 independent public and in-house metastasis melanoma and immunotherapy datasets covering multiple omics datasets. We found that decreased adenosine signaling mediates aberrant interactions between tumor immune cells and potentially reshapes the TME (Fig. 1).

Initially, we analyzed the changes in adenosine signaling at the pan-cancer level (see Supplementary Material). Subsequently, we classified 352 metastatic melanoma samples from the TCGA-SKCM cohort into 3 subtypes based on the expression levels of adenosine signaling-related genes (Fig. 2A, Fig. S4A, with optimal k = 3). Patients in the subtype 1 group (n = 157) exhibited worse survival than did those in the subtype 2 (n = 62) and subtype 3 (n = 133) groups, with no significant differences between Clusters 2 and 3 (Fig. 2B, log-rank test, p = 5.7 × 10–7). Interestingly, the genes significantly overexpressed (labeled diff_UP) in subtype 1 were mainly associated with the metabolic enzymes ATP/ADP (e.g., NME1 and AK1), while those in subtype 2 or subtype 3 were involved in the synthesis, degradation and signal transduction of adenosine (labeled diff_DN; e.g., NT5E, ENTPD2, ADORA1, CD38 and ADA2). Based on these distinct expression patterns, we further developed the adenosine signaling score (ADO-score), calculated as the enrichment score (ES) of diff_UP minus that of diff_DN, to quantify the adenosine signaling status of individual samples (Fig. S4B, see Method). The ADO-score accurately reflected the individual characteristics of the NMF subpopulations and the inherent co-expression patterns of adenosine genes (Fig. S4C and S4D). For prognosis, we found that Cluster 1, which had the worst survival, exhibited a significantly greater ADO-score than did Cluster 2 and Cluster 3 (Fig. S5A, Kruskal‒Wallis test p = 2.9 × 10–36).

Fig. 2

Clinical relevance of adenosine signaling in metastasis melanoma. A The landscape of adenosine signaling genes in metastasis melanoma between clinical factors in three adenosine metabolism subtypes stratified by NMF. B Kaplan–Meier curves of overall survival in TCGA-SKCM metastasis melanoma in three adenosine metabolism subtypes. C–E Kaplan–Meier curves of overall survival in 3 independent metastasis/stage IV melanoma datasets in high and low adenosine signaling group stratified by median ADO-score (D: GSE54467; E: GSE19234; F: GSE22155). F The ADO-score level various in response status in independent melanoma immunotherapy transcriptomic cohort GSE91061. G Overall survival in independent melanoma immunotherapy transcriptomic cohort GSE91061

The Univariate Cox hazard analysis and Kaplan‒Meier curves from the TCGA cohort showed that the ADO-score could be a prognostic factor for metastatic melanoma (Fig. S5B, hazard ratio (HR) = 1.7, 95% CI [1.2–2.2]; log-rank test, p = 6.7 × 10–4). Multivariate Cox regression analysis further revealed the independent prognostic value of the ADO-score in metastatic melanoma (Fig. S5C). The robustness of the ADO-score was validated using 3 independent metastatic/stage IV melanoma microarray datasets (GSE54467 [50], GSE19234 [51], and GSE22155 [52]) generated by microarray analysis, and corresponding clinical prognostic information were available. Consistently, the high ADO-score group presented worse OS that the low ADO-score group in the univariate Cox analysis (Fig. 2C–E; GSE54467: HR = 1.8, 95% CI [1–3.2], p = 4.7 × 10–2; GSE19234: HR = 3, 95% CI [1.2–7], p = 7 × 10–3; GSE22155: HR = 1.7, 95% CI [0.95–3.2], p = 7.3 × 10–2) and the multivariate Cox analysis (Fig. S5C; GSE54467: HR = 1.3, 95% CI [0.7–2.4], p = 0.37; GSE19234: HR = 3.7, 95% CI [1.4–9.6], p = 6 × 10–3; GSE22155: HR = 2.1, 95% CI [0.9–4.8], p = 7 × 10–2).

Drug therapy (both targeted therapy and immunotherapy) is a critical component of melanoma treatment. To investigate the ability of the ADO-score in immune checkpoint blockade (ICB) efficacy, we further evaluated the value of the ADO-score in stratifying immunotherapy efficacy in cohorts with published immunotherapy response and transcriptomic data (Fig. 2G, H) [53, 54]. The responders treated with nivolumab had significantly lower ADO-score values than the non-responders (GSE91061: Wilcoxon test, p = 0.05; Fig. 2F). Patients with high ADO-score value had worse overall survival (OS) (GSE91061: log-rank test, p = 0.014; Fig. 2G). These results indicate that the ADO-score is a prognostic biomarker for melanoma. To further evaluated the effect of ADO-score in melanoma chemotherapeutic response, we selected common targeted therapy for melanoma (including BRAF and MEK inhibitor) for sensitivity analysis. The results show that lower ADO-score patients exhibited significantly better dabrafenib and trametinib sensitivity (Figs. S5D, S5E). Therefore, dabrafenib and trametinib may be more effective in treating patients with low ADO-scores.

The cell-type specificity of adenosine signaling protein expression is associated with the prognosis of melanoma patients receiving PD-1 mAb therapy

To further validate the phenomenon that we observed in transcriptome and to clarify the ability of adenosine signaling components in clinical biopsy sections to act as biomarkers, we performed immunofluorescence staining on the collected pathological sections. We collected clinical biopsies of 19 melanoma patients treated with PD-1 mAb therapy (toripalimab; Fig. 3A, B and Table S4) from Xiangya Hospital and quantified adenosine signaling-related gene expression using a fluorescent multiplex immunofluorescence (multiplex IF) assay (for example, for ADORA1, ENTPD2/CD39, NT5E/CD73, CD38, AK1, AK2 and NME1). We use CD3, CD20, CD68 and SOX10 to identify the T cells, B cells, Macrophages and Mali Analysis of the immunofluorescence results demonstrated that the expression of AK1 and NME1 was significantly positively correlated with the expression of SOX10 (AK1: R = 0.63, p = 0.0046, NME1: R = 0.52, p = 0.025; Fig. 3C), while the expression of CD38 and ENTPD2 was significantly positively correlated with the expression of CD3 (CD38: R = 0.52, p = 0.024, ENTPD2: R = 0.74, p = 0.00043; Fig. 3D), suggesting potential co-expression pattern. Multiplex IF microscopy analysis indicated that AK1, AK2 and NME1 could be detected in SOX10-positive tumor cells, and NT5E, ENTPD2, and ADORA1 were detected in the surrounding CD20+ B-cell compartment. In addition, cell surface-bound CD38, ENTPD2 and ADORA1 were expressed in both tumor-infiltrating CD68+ macrophages and CD3+ T cells (Fig. 3E). We observed an interesting pattern in which adenosine signaling molecules had a spatial distribution bias in terms of their cell type specificity, and this specificity seems to represent the crosstalk between tumor cells and immune cells in the tumor microenvironment mediated by adenosine signaling.

Fig. 3

The immunofluorescence and tissue specificity of adenosine signaling genes is related to the prognosis of melanoma patients with anti-PD1 therapy. A Overview of melanoma patients with anti-PD1 monotherapy IF staining cohorts. All patients were treated with PD1 monoclonal antibody after surgical resection of primary lesions. See also Table S1. B Clinical status statistics of Inhouse melanoma anti-PD1 monotherapy IF staining cohort. C Correlation of adenosine signaling genes AK1 and NME1 with malignant cell marker SOX10 in IF staining slides. D Correlation of adenosine signaling genes CD38 and NNTPD2 with T cell marker CD3 in IF staining slides. E Biopsy specimens obtained from melanoma samples including responders (No.13, No.19) and non-responders (No.4, No.18) were examined and visualized by multiplex IHC. With this staining technique, visible structures include tumor cell panel labeled with DAPI (blue), SOX10 (yellow), AK1 (red), AK2 (green), NME1 (magenta), macrophage panel labeled with DAPI (blue), CD68 (white), CD38 (red), ENTPD2 (magenta), ADORA1 (orange), B cell panel labeled with DAPI (blue), CD20 (cyan), NT5E (red), ENTPD2 (magenta), ADORA1 (orange), and T cell panel labeled with DAPI (blue), CD3 (green), CD38 (red), ENTPD2 (magenta), ADORA1 (orange). Majority of positive cells are located in the yellow dashed box and highlighted with white arrows. DAPI is a nuclear counterstain. Scale bars, 50 μm. F Kaplan–Meier survival curve of melanoma patients’ progression-free survival (PFS) following anti-PD1 therapy grouped by adenosine signaling

After observing the specificity of adenosine-related gene expression in different cells within the cancer immune microenvironment, we further explored the clinical translational value of these genes. Considering the functional relevance of adenosine-related genes in adenosine metabolism, we further assessed whether the ADO-score could represent the celltypes specificity patterns and interaction between tumor and immune cells. (See Materials and Methods). Overall, protein tissue analysis reflected the ADO score derived from gene expression. Meanwhile, we observed that patients receiving PD-1 blockade therapy with low ADO-score values exhibited prolonged progression-free survival (HR = 3.1, 95% CI [1.1–8.8], p = 0.027; Fig. 3F). The above phenomenon suggests that the distinct expression patterns of adenosine-related genes on immune cells and tumor cells can be used to differentiate the response status of patients to immunotherapy. Taken together, our results showed that adenosine signaling might mediate crosstalk between tumor cells and immune cells and could be a prognostic factor in melanoma patients.

Deciphering adenosine signaling at the single-cell and spatial transcriptomics levels

To gain more insight into the utility of the ADO-score model, we investigated in the different adenosine signaling-related genes and ADO-score at the single-cell level. We evaluated the adenosine signaling status of each cell in the melanoma single-cell transcriptomic dataset (GSE115978 and GSE72056) using the same method mentioned above (see Methods). We observed the highest ADO-score value in malignant cells, followed by fibroblasts, and the lowest value was exhibited in immune cells (Fig. 4A–C and D–F). To our surprise, the genes from the diff_UP and diff_DN gene sets were specifically expressed in melanoma and immune cells, respectively (Fig. S6A and S6B). The genes from the diff_UP gene set (AK1, AK2 and NME1) were all generally related to the transformation of adenosine phosphate compounds, such as the ATP-ADP transformation, while the genes from diff_DN mainly participated in adenosine signaling activation and adenosine metabolism through the transformation of adenosine phosphate compounds, including ATP, ADP and AMP, into adenosine in the melanoma microenvironment. In addition, the genes from the diff_DN gene set also exhibited cell type-specific expression (Fig. S6C and S6D). CD38 exhibited the highest expression in macrophages, while NT5E was expressed in B cells and stromal cells, especially in fibroblasts. These results implied that the ADO-score represents the crosstalk between melanoma cells and immune cells centered on energy metabolism.

Fig. 4

Distribution and functional association of the ADO-score in single cell and spatial transcriptomics level. A UMAP embedding of single-cell RNA-seq profiles from GSE115978. B UMAP plot show ADO-score profiles of whole tissue cells from GSE115978. C Boxplot and dotplot show the difference and percentage of ADO-score in different cell types in melanoma form GSE115978. D UMAP embedding of single-cell RNA-seq profiles from GSE72056. E UMAP plot show ADO-score profiles of whole tissue cells from GSE72056. F Boxplot and dotplot show the difference and percentage of ADO-score in different cell types in melanoma form GSE72056. G Spatial enhanced-resolution clustering performed by the BayesSpace algorithm identified four clusters corresponding to the original histopathological annotations. H Spatial Feature plot shows the difference of ADO-score profile among four clusters at the enhanced-resolution condition. I Violin plot show the difference of ADO-score in different cell types among four clusters at the enhanced-resolution condition. J, K The biological functions (J cancer hallmarks and K Reactome biological pathways) divergence for malignant cells with different adenosine signaling states stratified by ADO-score

We further examined these phenomena by integrating spatial information. We first acquired spatial transcriptomic data from melanoma [55] and applied the BayesSpace algorithm to obtain higher resolution images of malignant cells (expressing PMEL), T cells (expressing CD2 and CD3D), B cells (expressing CD19), fibroblasts (expressing COL1A1) and macrophages (expressing C1QB). The annotation clustering results revealed that malignant (red region) and T/B cells (blue region) were separated into stromal (orange region) and macrophage (green region) populations (Fig. 4G). Similarly, the distribution of the ADO-score showed a distinct pattern among the different regions. The ADO-score was greater in the red region, which included a significantly higher proportion of malignant cells than stromal cells, immune cells and macrophages (Fig. 4H, I).

We further characterized the biological functions of malignant cells in different adenosine signaling states stratified by the ADO-score. In the high-ADO-score group, we observed significant enrichment of metabolic and proliferation hallmarks, including oxidative phosphorylation, fatty acid metabolism, the G2M checkpoint, spermatogenesis and MYC targets, while in the low-ADO-score group, we observed enrichment of signatures related to antitumor immunity, including TNF-alpha signaling, promotion of the inflammatory response and interferon alpha/gamma response, activation of the complement cascade and malignant cell apoptosis (Fig. 4J and S6E). According to the biological pathway enrichment results from Reactome, malignant cells with low ADO-score values were enriched in immunoregulatory interactions between a lymphoid cell and a non-lymphoid cell and between the complement cascade and interferon gamma signaling, while tumor cells with high ADO-score values were enriched in several metabolic reprogramming and energy metabolism pathways. The citric acid cycle, respiratory electron transport and ATP synthesis were also evaluated (Fig. 4K and S6F). All of these findings indicate the heterogeneity of adenosine signaling in the tumor microenvironment; in particular, malignant cells with different adenosine signaling states exhibit distinct biological functions.

Biological effects of adenosine signaling in melanoma

To investigate the potential effects of adenosine signaling, we performed a comprehensive analysis of oncogenic and immune features, including cancer hallmarks, immune cell infiltration, immune checkpoints, tumor mutation burden (TMB), the T-cell-inflamed gene expression profile (GEP), cytolytic activity (CYT), PD-L1 protein level, B-cell receptor (BCR) enrichment, T-cell receptor (TCR) enrichment, interferon gamma (IFN-gamma) response and the tumor-infiltrating lymphocyte (TIL) regional fraction, in the TCGA-SKCM and 7 independent melanoma datasets (Figs. 5, S7). Generally, we found that the ADO-score was negatively correlated with immune features and positively correlated with the enrichment of oncogenic pathways. Immune cell infiltration was estimated with multiple algorithms (including ssGSEA estimation for immune cell score, EPIC, ESTIMATE, MCP-counter, quanTIseq and Xcell), and the ADO-score was found to be negatively correlated with immune infiltration cells in multiple independent datasets (Fig. 5A, S7-F). We further collected 40 druggable targets or potential biomarker candidates among immune checkpoint genes (Table S6). The ADO-score was negatively correlated with the levels of nearly all immune checkpoint genes, implying the presence of a potentially inhibitory immune microenvironment (Fig. 5B).

Fig. 5

Biological effects of ADO-score in metastasis melanoma. A Correlation between ADO-score and immune cell infiltration estimated by TIMER in 8 independent melanoma datasets. B Correlation between ADO-score and immune checkpoints in independent melanoma datasets. C–F The correlation between ADO-score with known immune therapy prognosis factors (C Tumor mutation burden, D CYT scores, E GEP score. F PD-L1 expression). G–J Difference between high and low ADO-score in immune repertoire (G BCR richness, H TCR richness, I IFN-γ response, J TIL region fraction)

Immunotherapy is the latest treatment for melanoma patients, while the efficacy exhibited much heterogeneity. Patient stratify using promising biomarkers is critical for therapeutic approaches selection. Several predictive biomarkers for immunotherapy efficacy, including CYT andGEP, also exhibited rather strong significant negative correlations with the ADO-score but not with PD-L1 and TMB (Fig. 5C–F). In addition, the diversity of the immune repertoire was also significantly greater in the low-ADO-score group than in the high-ADO-score group; the IFN-gamma response was significantly enriched in the low-ADO-score group, but the difference in the proportion of TILs was not significant (Fig. 5G–J). For 50 pathways related to cancer hallmarks, we also revealed that the ADO-score was negatively correlated with the enrichment of immune-related pathways (e.g., interferon-alpha and interferon-gamma response) but positively correlated with the enrichment of oncogenic pathways (e.g., MYC target; Fig. S7A). Together, our results demonstrated the relationship between adenosine signaling and the immune microenvironment and implied that patients with lower ADO-score values are more likely to have a “hot” tumors microenvironment, while patients with higher ADO scores may have a “cold” tumor microenvironment. The ADO-score may be a promising indicator of “cold” tumor immune micro-environment (TIME) because of the strong negative association between ADO-score and the levels of both activating and inhibitory immune checkpoints.

Malignant cells with different adenosine signaling pathways exhibit distinct cell communications in the tumor microenvironment

We further explored the differences in cellular crosstalk mediated by malignant cells with high and low adenosine signaling in GSE115978 with CellChat [43]. Overall, disparate communication signals were observed between malignant cells in the high ADO-score group and those in the low ADO-score group (Fig. 6A: upregulated signaling in the high ADO-score group; Fig. 6B: upregulated signaling in the low ADO-score group). In the high ADO-score group, we observed significant upregulation of FN1, which is a ligand that activates integrin (ITG) superfamily NOTCH signaling among nonmalignant cells, especially fibroblasts (Figs. 6A; S8A and B). Tumor cancer cell migration and invasion are potentially involved, especially in generating lymph angiogenesis and tumor cell colonization and triggering invasive protrusions and pro-invasive EMT signaling [56,57,58,59]. The other major high-ADO-score malignant signaling pathways shaping the immunosuppressive microenvironment were MDK and MPZL1 (Fig. 6A). MDK is involved in promoting immunosuppressive macrophage differentiation and endothelial tube formation [60, 61], and MPZL1 has been shown to have an oncogenic function by promoting cell proliferation, migration and invasion but inhibiting cell apoptosis [62, 63]. The relative information flow also revealed several predominant cancer signaling pathways, including the CADM, VEGF and MPZ signaling pathways, in the high ADO-score group (Fig. S8C). For specific ligand–receptor pairs, malignant cells in the high-ADO score group could interact with endothelial cells and fibroblasts via VEGFB–VEGFR1, supporting the use of antiangiogenic drugs (Fig. 6C). The activation of cell‒cell signaling mediated by the COL1A2–ITG family implies a greater risk of tumor invasion [64].

Fig. 6

Distinct cellular communications associated with adenosine signaling in malignant cells. A, B Circos plot of the cellular communications associated with high and low ADO-score. Showing only significantly differentiated signaling pathway using Cellchat. C Volcano plot and lollipop plot show activated and top 10 activated LR-pairs associated with high and low ADO-score using NATMI. D The functional enrichment of top targets inferred by Nichenet for CD8T cells affected by different adenosine signaling

Conversely, the upregulated signals in malignant cells with low ADO-score values were associated mainly with type I HLA signaling between malignant cells and Tex and Tprolif cells (Fig. 6B). For malignant cells in the low ADO-score group, we also found stronger interactions with CD8 Tex cells than with those in the high ADO-score group via HLA-A/B/C-CD8A/B pairs, mediating the CD8+ T-cell-dependent killing of cancer cells by efficient presentation of tumor antigens [65] (Figs. 6B and S8D). Moreover, HLA molecules also communicate with KIRK1, which encodes NKG2D and is constitutively expressed on CD8+ T cells, to authenticate the recognition of a stress-induced target and enhance TCR signaling [66]. Conversely, decreased expression of HLA class I molecules on tumors and impaired signal transduction may facilitate tumor immune escape [67], which to some extent explains why patients with high ADO-score values exhibit poorer survival and immunotherapy resistance.

To further confirm the differences in cellular communication between malignant cells with high and low adenosine signaling and to comprehensively characterize their tumor–immune interactions, we further applied NATMI [44] to identify the unique activated ligand‒receptor interactions between malignant cells with high adenosine signaling and those with low adenosine signaling. Similar to the CellChat results, in the malignant cells with high adenosine signaling, the primary ligand‒receptor interactions predominantly occurred between the tumor and stromal cells (fibroblasts and endothelial cells; Fig. 6C, left upper panel), whereas in the low-adenosine signaling group, more interactions between tumor cells and immune cells (including macrophages and T cells; Fig. 6C, left lower panel) occurred. In the high ADO-score group, the activated pathways, which involved pathways such as ANGPT1, EGFR, and integrins, primarily contributed to angiogenesis, cell adhesion, and tumor invasion and metastasis [68, 69]. In contrast, in the low ADO-score group, significant complement activation-mediated interactions were observed, which may exert antitumor effects [70].

To further investigate the potential functional effects of divergent adenosine signaling states on T cells, we identified the target genes of the 20 most differentially expressed receptors on CD8Tex cells and performed functional analysis of using NanoNet [45] (Figs. 6D, S8E). Multiple pathways involved in T-cell activation, chemotaxis, antigen recognition and cytotoxic functions, such as the TNF signaling pathway, T-cell receptor signaling pathway, leukocyte transendothelial migration pathway and chemokine signaling pathway, were enriched (Fig. 6D). Additionally, recent studies have suggested that the IL-17 signaling pathway is associated with CD8+ cytotoxic T-cell infiltration and cytotoxic function via the ICOS-ICOSL interaction [71].

In summary, malignant cells with high ADO-score values may contribute to cancer progression via communication with endothelial cells and fibroblasts, while malignant cells with low ADO-scores prohibit cancer progression by communicating with different adenosine signaling pathways, which results in more frequent communication with immune cells.