Melanoma cell states shape spatial tumor-immune ecosystems to dictate the efficacy of anti-PD1 immunotherapy
Authors: Félix K Pham, Marion Dufeu, Valentin Benboubker, Maxime Grimont, Amélie Lhorisson, Justine Berthet, Marie Donzel, Raphaël Schneider, Laurie Tonon, Anne-Claire Doffin, Félix Boivin, Simon Durand, Bertrand Dubois, Jonathan Lopez, Christophe Caux, Jenny Valladeau-Guilemond, Anaïs Eberhardt, Stéphane Dalle, Julie Caramel
Published: 2026

ABSTRACT
AbstractBackgroundMelanoma shows one of the highest response rates to immune checkpoint inhibitors (ICIs), yet nearly half of patients experience primary or acquired resistance. While immune contexture strongly influences therapeutic efficacy, tumor cell-intrinsic features are increasingly recognized as key regulators of antitumor immunity. In particular, intratumoral heterogeneity driven by melanoma cell plasticity underlies diverse immune escape mechanisms. How this plasticity shapes ICI outcomes in patients remains poorly defined.MethodsTumor cell states and immune contexture were assessed in 57 primary cutaneous melanomas from stage III patients, collected prior to adjuvant anti-programmed cell death protein-1 (anti-PD-1) therapy and stratified according to 2-year relapse status (33 relapse-free, 24 relapsed). Whole slide multiplex immunofluorescence was combined with spatial transcriptomics (Visium, n=4) to investigate the spatial architecture of melanoma cell states, T cells, tumor-associated macrophages (TAMs), dendritic cell subsets, and tertiary lymphoid structures.ResultsUnsupervised clustering of melanoma cells identified distinct phenotypic states that formed spatially restricted homotypic patches. From these data, we defined a melanoma plasticity ratio (undifferentiated/differentiated tumor patches), which was significantly associated with reduced relapse-free survival. Integrated immune analyses recapitulated prognostically distinct immunotypes, with macrophage subsets displaying striking spatial compartmentalization. Antitumoral macrophages preferentially infiltrated differentiated melanoma regions, while protumoral macrophages localized to undifferentiated patches. Spatial transcriptomics confirmed that melanoma cell states tightly shape the neighboring immune microenvironment, with macrophages emerging as pivotal players. Their polarization was further influenced by tumor-derived signals in addition to microenvironmental cues (interferon-gamma, hypoxia). Entropy-based integration of melanoma cell states, TAMs, and T cell subsets uncovered two dominant spatial ecosystems with opposing associations to ICI efficacy. Ecosystems enriched in differentiated melanoma cells, programmed death-ligand 1 (PD-L1)+ TAMs, and PD-1+CD8+, and CD4+ T cells correlated with favorable outcomes, whereas ecosystems composed of undifferentiated melanoma cells with PD-L1− protumoral TAMs and PD-1−CD8+ T cells were associated with relapse.ConclusionsOur study uncovers how cancer cell plasticity shapes spatially organized tumor-immune ecosystems that critically modulate adjuvant ICI efficacy in melanoma. These findings highlight melanoma cell plasticity as a key driver of immune evasion via macrophages reprogramming, through targetable interactions that may represent novel therapeutic avenues to enhance ICI efficacy.

WHAT IS ALREADY KNOWN ON THIS TOPIC

Intratumor heterogeneity driven by melanoma cell plasticity contributes to immune evasion but its precise contribution to immune checkpoint inhibitor (ICI) resistance remains unclear. Although an enrichment of undifferentiated melanoma cells has been observed in metastatic samples from non-responding patients after anti-programmed cell death protein-1 (anti-PD1) immunotherapy, the spatial organization of melanoma cells, their interactions with the tumor immune microenvironment, and their clinical relevance remain to be elucidated in early-stage primary cutaneous tumors.

Intratumor heterogeneity driven by melanoma cell plasticity contributes to immune evasion but its precise contribution to immune checkpoint inhibitor (ICI) resistance remains unclear. Although an enrichment of undifferentiated melanoma cells has been observed in metastatic samples from non-responding patients after anti-programmed cell death protein-1 (anti-PD1) immunotherapy, the spatial organization of melanoma cells, their interactions with the tumor immune microenvironment, and their clinical relevance remain to be elucidated in early-stage primary cutaneous tumors.

WHAT THIS STUDY ADDS

This study demonstrates that melanoma cell states form spatially restricted patches, with the balance between undifferentiated and differentiated patches associated with the relapse risk. It further shows that this intratumor heterogeneity dictates the polarization and spatial organization of tumor-associated macrophages and T cells through tumor-derived signaling pathways, thereby defining distinct tumor-immune ecosystems linked to ICI outcomes.

This study demonstrates that melanoma cell states form spatially restricted patches, with the balance between undifferentiated and differentiated patches associated with the relapse risk. It further shows that this intratumor heterogeneity dictates the polarization and spatial organization of tumor-associated macrophages and T cells through tumor-derived signaling pathways, thereby defining distinct tumor-immune ecosystems linked to ICI outcomes.

HOW THIS STUDY MIGHT AFFECT RESEARCH, PRACTICE OR POLICY

These findings establish melanoma cell plasticity as a driver of anti-PD1 immunotherapy efficacy through its crosstalk with immune cells in spatially restricted areas of the tumor. They further underscore the clinical relevance of considering how tumoral and immune compartments crosstalk towards promoting immune evasion, highlighting promising targetable interactions to interfere with resistance.

These findings establish melanoma cell plasticity as a driver of anti-PD1 immunotherapy efficacy through its crosstalk with immune cells in spatially restricted areas of the tumor. They further underscore the clinical relevance of considering how tumoral and immune compartments crosstalk towards promoting immune evasion, highlighting promising targetable interactions to interfere with resistance.

Introduction

Immune checkpoint inhibitors (ICIs) targeting programmed cell death protein-1 (PD-1) have significantly improved the outcome of patients with advanced melanoma; however, primary or acquired resistance remains a clinical challenge for over half of treated patients.1 2 T cell subsets are key mediators of the response to ICI,36 and their activation through cross-presentation by dendritic cells (DCs) is critical.79 The presence of tertiary lymphoid structures (TLSs) in tumors has also been linked to improved clinical outcomes, possibly because of their role in facilitating local immune priming and cellular interactions.10 11 Moreover, recent research has highlighted the central role of tumor-associated macrophages (TAMs) in shaping responses to ICI. Depending on their phenotype and spatial context, TAMs may suppress or enhance anti-tumor immune responses, ultimately tipping the balance towards resistance or responsiveness.1214 Importantly, the spatial organization of the tumor immune microenvironment (TIME) is increasingly being recognized as pivotal for anti-PD1 (aPD1) ICI efficacy. Specific spatial niches composed of CD8+ T cells,15 conventional type 1 DC (cDC1)9 and/or TAM16 17 have recently been described in several tumor contexts, emphasizing the functional relevance of the tissue architecture in the regulation of the immune response. However, most studies have focused on immune-restricted parameters.

Although the genetic determinants of ICI resistance have been described,1822 the majority of melanoma tumors lack resistance-conferring genetic variants.23 Instead, melanoma exploits cancer cell plasticity as a key mechanism of therapeutic resistance,24 25 especially in targeted therapies.26 27 Similar to the epithelial–mesenchymal transition (EMT) that occurs in epithelial cancers, melanoma cell plasticity is related to reversible phenotypic switches between differentiated/melanocytic and undifferentiated/invasive cell states. These transitions, regulated by transcription factors, enable melanoma cell adaptation to the tumor microenvironment (TME)28 and drive intratumoral heterogeneity (ITH). Several immune escape mechanisms linked to melanoma cell plasticity have been described, suggesting a key contribution in shaping antitumor immunity.2934 Nevertheless, the extent to which these mechanisms contribute to ICI resistance remains unclear. While a single-cell RNA sequencing (scRNA-seq) study recently deciphered the dynamic changes in melanoma cell states during aPD1 ICI in the metastatic setting,35 the spatial organization of melanoma cells, their interactions with the TIME, and their clinical impact remain to be precisely addressed in primary cutaneous melanoma.

In this study, we investigated the tumor-immune ecosystems in treatment-naïve primary cutaneous melanoma lesions from 57 stage III patients with lymph node metastasis who underwent surgery and adjuvant aPD1 ICI by integrating RNA sequencing (RNA-seq), spatial multiplex immunofluorescence (mIF), and spatially resolved transcriptomics (SRT). Through spatial analysis of melanoma ITH, we demonstrated that specific tumor-immune ecosystems emerge and are shaped by the crosstalk between melanoma cells and the immune microenvironment. Importantly, these spatially defined ecosystems affect ICI efficacy and metastatic relapse rates. Our findings revealed a spatially integrated tumor-immune architecture tuned by targetable interactions that contribute to therapeutic outcomes in advanced melanoma.

Results

Spatial organization of differentiated and undifferentiated tumor cells within human primary cutaneous melanoma is associated with ICI efficacy

To investigate the heterogeneity and spatial organization of melanoma cells according to their differentiation status in primary cutaneous melanoma lesions, we performed whole-slide mIF analysis, RNA-seq, and SRT on tumor lesions from resectable stage III patients with melanoma (n=57) treated with wide local excision and complete lymph node excision when required, followed by adjuvant aPD1 ICI for 12 months with at least 48 months of follow-up from the time of diagnosis. The main clinical endpoints were relapse status after 2-year follow-up (relapse-free, RF or relapse, REL) and recurrence-free survival (RFS) (figure 1A). The demographic, clinical, histopathological, and genetic information for each patient is presented in table 1. The formalin-fixed, paraffin-embedded (FFPE) primary melanoma lesions analyzed in this study were treatment-naïve before ICI infusion. Importantly, the stained sections were analyzed as whole slides to comprehensively capture the heterogeneity of the lesions with single-cell resolution.

Spatial organization of differentiated and undifferentiated melanoma cells within primary tumors is associated with ICI efficacy. (A) Schematic diagram of the analyses of pre-ICI stage III primary melanoma samples (n=57 patients), showing the clinical follow-up above, and mIF, bulk-RNAseq and spatial transcriptomic analyses performed on FFPE tissues below. Created with BioRender. (B) Representative tissue (of n=45) section showing the mIF markers used to detect differentiated and undifferentiated melanoma cells. DAPI counterstaining appears as blue within each panel. Scale bar=100 µM. (C) Unsupervised k-means clustering of each tumor cell (n=7,644,292) from the sections stained in (B) to define six specific phenotypic clusters based on the MFI of each marker. (D) Boxplots showing the association between the proportion of each individual phenotypic cluster and ICI outcome at 2 years: REL (in orange) and RF (in purple). Mann-Whitney U tests, p values are shown on the plot. (E) Spatial informatic reconstruction of the different phenotypic clusters in one tumor representative of n=45. Scale bar=1.5 mm. (F) mIF image of the C6 cells characterized by the co-expression of SOX10 and ZEB1 (left). These C6 melanoma cells (red) form homotypic patches (border in blue, middle) across the whole tumor (right). Left and middle scale bars=200 µm. Right scale bar=1.5 mm. (G) Coefficients of the features included in a LASSO regularized Cox proportional hazard model. Features with a positive coefficient are associated with worsened RFS, whereas those with a negative one are associated with improved RFS. (H) Definition of the MPR which is the ratio between undifferentiated C4+C6 cells and differentiated C2+C3 cells involved in a homotypic patch (left). Association between the MPR and the 2-year status (middle, boxplots, Mann-Whitney U tests) and RFS (right, Kaplan-Meier curves with MPRhigh patients in red and MPRlow patients in blue, Log-rank test) of patients. (I) mIF image (left) and spatial reconstruction of the indicated clusters (right). First row: a mIF image of an RF-associated MPRlow tumor displaying few C4-NGFRhigh patches (yellow) and numerous C2-MITFhigh patches (green). Second row: a mIF image of a REL-associated MPRhigh tumor displaying a high number of C6-ZEB1high patches (orange) and only few C2-MITFhigh patches (green). CLND, complete lymph node dissection; DAPI, 4′,6-diamidino-2-phenylindole, DC, dendritic cell; FFPE, formalin-fixed, paraffin-embedded; ICI, immune checkpoint inhibitor; LASSO, Least Absolute Shrinkage and Selection Operator; MFI, mean fluorescent intensity; mIF, multiplex immunofluorescence; PD1, programmed cell death protein-1; REL, relapse; RF, relapse-free; RFS, relapse-free survival; RNAseq, RNA sequencing; TAM, tumor-associated macrophage; TLS, tertiary lymphoid structure; TME, tumor microenvironment.

Table 1Clinical characteristics of the cohort2-year statusRFN=33 (57.9%)RELN=24 (42.1%)P valueStatistical testAge at diagnosis0.81Wilcoxon Years, mean±SD63.3±14.463.3±15.8Follow-up duration0.22Wilcoxon Months, median (range)33 (27–42)32.5 (19–42)Sex0.78χ² Female, n (%)1110 Male, n (%)2316BRAF/NRAS mutational status0.64χ² BRAF, n (%)1413 NRAS, n (%)96 Non-mutated,* n(%)105Breslow index0.02Wilcoxon mm, median (range)3.4 (1.0–18.0)5.9 (1.2–20.0)Ulceration0.03χ² Yes, n (%)16 (49)20 (83) No, n (%)18 (51)5 (17)AJCC-eight stages0.0047χ²† IIIA, n(%)5 (15)1 (4) IIIB, n(%)9 (27)1 (4) IIIC, n(%)17 (52)19 (79) IIID, n(%)2 (6)3 (13)Location of metastases at relapseN/AN/A Skin06 Lymph node012 Lung011 Liver07 Surrenal glands04 Spleen01 Bones01 Retroperitoneal01 Brain03Number of anti-PD1 infusions‡<0.0001Wilcoxon Median (range)17 (4–19)7 (2–18)*Or rare mutations: ARID2 (n=1), KIT (n=1), MAP3K1 (n=1), PTEN (n=1), RB1 (n=1), TERT (n=3).†Stages IIIA and IIIB were combined, and stages IIIC and IIID were combined for statistical purposes.‡Until first relapse.AJCC, American Joint Commission on Cancer; N/A, not applicable; PD1, programmed cell death protein-1; REL, relapse; RF, relapse-free.

Clinical characteristics of the cohort

Or rare mutations: ARID2 (n=1), KIT (n=1), MAP3K1 (n=1), PTEN (n=1), RB1 (n=1), TERT (n=3).

Stages IIIA and IIIB were combined, and stages IIIC and IIID were combined for statistical purposes.

Until first relapse.

AJCC, American Joint Commission on Cancer; N/A, not applicable; PD1, programmed cell death protein-1; REL, relapse; RF, relapse-free.

To explore the phenotypic heterogeneity of melanoma cells, a seven-color mIF panel composed of melanoma cell state markers was designed (figure 1B and online supplemental figure S1A). In addition to SOX10,36 37 MITF28 38 39 and ZEB240 41 were selected as markers of a melanocytic/differentiated cell state. Conversely, ZEB1,40 42 NGFR30 43 and SOX943 44 have been included as markers of invasive/neural crest stem-like (NCSC)/undifferentiated cell states.

Overall, 7,644,292 SOX10+ melanoma cells were detected in all the tumor samples. As these markers are not exclusive, we performed k-means clustering of melanoma cells based on the mean fluorescence intensity (MFI) of each marker. This analysis revealed six phenotypic clusters that displayed a gradient in SOX10 expression, with C3 exhibiting the highest SOX10 expression levels (figure 1C and online supplemental figure S1B). ZEB2 was expressed in the C1 (MITF-int) and C2 (MITF-high) clusters, whereas NGFR expression was mostly restricted to C4 (figure 1C and online supplemental figure S1B). A proportion of C6 (ZEB1high) and C5 (SOX9high) clusters displayed very low SOX10 levels (figure 1C and online supplemental figure S1B), yet we were still able to detect these cells, thus identifying at least a subset of undifferentiated tumor cells that may have a downregulated SOX10 expression. Coexpression of SOX9 with MITF in C5 cells (figure 1C, online supplemental figure S1B) may reflect a previously described intermediate cell state.44 45 A cluster presenting very low expression of all six markers was considered undefined and was removed from further analyses. For clarity, we refer to these key phenotypic clusters as C1-ZEB2high, C2-MITFhigh, C3-SOX10high, C4-NGFRhigh, C5-SOX9high and C6-ZEB1high.

The proportion of each cluster was quantified in all patients (online supplemental figure S1C). Despite trends towards enrichment in the baseline frequencies of the C2-MITFhigh cluster in RF tumors and of the C6-ZEB1high cluster in REL tumors, none of these phenotypic clusters was significantly associated with the 2-year relapse status (figure 1D). This finding is similar to that of a recent study, in which no baseline difference in melanoma cell state frequency was found between aPD1 ICI-responding and non-responding stage IV patients with melanoma using scRNA-seq.35 Remarkably, when considering the spatial organization of these tumors, melanoma cells belonging to the same phenotypic cluster tended to be spatially aggregated (figure 1E and online supplemental figure S2A–I). To quantify this spatial organization, we defined homotypic tumor patches46 as regions containing at least five melanoma cells from the same phenotypic cluster within a 25 µm radius (figure 1F). Interestingly, melanoma cells of the C2-MITFhigh cluster predominantly formed homotypic patches, whereas cells of the C6-ZEB1high phenotype were mostly isolated, which may reflect their invasive capacity (online supplemental figure S3A).

Although the overall density of melanoma patches was not correlated with patient outcome after aPD1 treatment, a Cox regression model fitted with Least Absolute Shrinkage and Selection Operator (LASSO) regularization identified the most relevant cluster patches associated with RFS (figure 1G and online supplemental figure S3B,C). Notably, melanoma cells involved in the C2-MITFhigh and C3-SOX10high patches were associated with better outcomes, while cells involved in the C4-NGFRhigh and C6-ZEB1high patches had worse outcomes.

To further investigate the clinical relevance of this melanoma cell state organization, we calculated the melanoma plasticity ratio (MPR), defined as the ratio of C4-NGFRhigh and C6-ZEB1high melanoma cells organized as patches to C2-MITFhigh and C3-SOX10high cells (figure 1H). A higher MPR was found in REL patients and was consistently associated with a worse RFS (figure 1H–K). In contrast, a non-spatial MPR ratio that did not take into account tumor cell spatial organization was not significantly associated with ICI outcomes (online supplemental figure S3D), highlighting the role of pre-treatment undifferentiated melanoma cells’ spatial distribution in reducing the efficacy of adjuvant aPD1 ICI. Interestingly, whereas the Breslow index (the most relevant histopathological prognostic metric), was significantly increased in REL patients (table 1), the prognostic value of the MPR ratio was independent of it (online supplemental figure SE–G).

To further quantify the diversity in melanoma cell clusters, we calculated the Shannon diversity index across every tumor, acknowledging that a higher index reflects a more diverse representation of every melanoma cell cluster.47 Highly undifferentiated tumors exhibited a higher Shannon diversity index, suggesting that acquisition of a more diverse melanoma cell state profile endows tumors with a survival advantage (online supplemental figure S3H).

These findings revealed that melanoma cells are spatially organized in homotypic patches and that the ratio between undifferentiated and differentiated melanoma patches in the pretreatment primary tumor influences the efficacy of aPD1 ICI.

Integrative analysis of immune contexture highlights spatial organization of TAM subsets as a crucial determinant of ICI efficacy

We then performed a comprehensive characterization of the immune context using mIF across multiple tumor sections. This included assessing the densities and spatial organization of T cells (online supplemental figure S4A, analyzed with the CD3, CD8, PD1 and KI67 markers), DCs subpopulations (online supplemental figure S4B)—(cDC1 (CLEC9A+), conventional type 2 DCs (cDC2, CLEC10A+), plasmacytoid DCs (pDCs, BDCA2+), and migratory DCs (LAMP3+))—TLS (online supplemental figure S4C, assessed by the presence of at least one segregated T cell and B cell zones across the whole slide using the CD3, CD20 and DC-LAMP markers), and TAM (figure 2A and online supplemental figure S4D, analyzed with the pan-TAM marker CD68, pSTAT1 as a readout of proinflammatory TAM, CD163 as a marker of protumoral TAM and programmed death-ligand 1 (PD-L1) expression).

TAMs density and spatial location dictate response to immune checkpoint inhibitors. (A) mIF antibodies used to study TAM subsets. Scale bar=40 µm. (B) Density and percentage of TAMs in all tissue sections (n=48). Graphs represent all TAMs, PDL1+, pSTAT1+CD163−, and CD163+pSTAT1− TAMs, as well as the TAM ratio according to the 2-year status. Mann-Whitney U tests, p values are indicated on the boxplots. (C) Proximity score illustrating the intratumoral infiltration by TAM subsets in a perimeter of 25 µm around a tumor cell. A higher score indicates more colocalization between both cell types independently of their frequencies. Friedman test and paired Wilcoxon tests, with p values adjusted by the Benjamini-Hochberg method. (D) Reconstructed tissue section displaying the specific spatial pattern of PDL1−pSTAT1 TAM and PDL1−CD163 TAM subsets within the same tumor. Scale bar=4 mm. (E) Corresponding mIF images from the tumor reconstructed in (D). Right panels show enlarged views of inserts (scale bar=50 µm) in the images on the left (scale bar=4 mm). (F) Proximity scores between PDL1+pSTAT1 TAMs and target cells in a radius of 25 µm. Mann-Whitney U tests between RF (purple) and REL (orange) patients. P values are indicated in the boxplot. (G) Heatmap of immune cell frequencies from the mIF analyses, in addition to the IGR hallmark established by ssGSEA on RNA-seq data. Each column represents a patient (n=48). Patients are clustered with partition around medoids clustering into three immunotypes. (H) Kaplan-Meier curves of RFS according to the immunotypes obtained in G. Log-rank test, the p value is indicated on the graph. cDC1, conventional type 1 dendritic cell; cDC2, conventional type 2 dendritic cell; DAPI, 4′,6-diamidino-2-phenylindole, DC, dendritic cell; IGR, interferon-gamma response; IT, intratumoral; mIF, multiplex immunofluorescence; NA, not available; pDC, plasmacytoid dendritic cell; PD-L1, programmed death-ligand 1; REL, relapse; RF, relapse-free; RFS, relapse-free survival; RNA-seq, RNA sequencing; ssGSEA, gene-centric single sample gene set enrichment analysis; TAM, tumor-associated macrophage; TLS, tertiary lymphoid structure.

Within the T cell compartment, we observed a trend towards increased CD8+ T cell density in RF patients, whereas REL patients exhibited a higher proportion of PD1+CD8+ T cells (online supplemental figure S5A). We found that higher densities of cDC1, pDCs, and LAMP3+ DCs were associated with improved aPD1 efficacy (online supplemental figure S5B), as previously reported in stage IV melanoma.9 TLS was detected in only four of the tumors in our primary melanoma cohort, precluding any conclusion regarding their prognostic role. While total TAM (CD68+ cell) densities were comparable between RF and REL patients, we observed an increased proportion of pSTAT1+CD163− TAM (pSTAT1 TAM) in RF patients (figure 2B and online supplemental figure S5C,D). Although no significant difference was detected in the relative proportion of CD163+pSTAT1− TAM (CD163 TAM), the pSTAT1/CD163 TAM ratio was increased in RF patients, suggesting that proinflammatory TAM are associated with improved aPD1 efficacy (figure 2B).

By integrating the immune cell densities, we identified three distinct immunotypes at the patient level (figure 2G,H). Immunotypes 1 and 3 correspond to previously described inflamed (“hot”) and non-inflamed (“cold”) TIMEs, respectively.48 More precisely, immunotype 1 was characterized by high intratumoral infiltration of CD8+ and CD4+ T cells with DCs, along with numerous pSTAT1 TAM and PDL1+ TAM. As expected, this density pattern was fueled by high interferon-gamma (IFN-γ) signaling as evidenced by gene set enrichment analysis (GSEA) of RNA-seq data (figure 2G),49 with the upregulated genes from immunotype 1 (vs immunotype 3) being associated with immune cell recruitment and activation pathways (online supplemental figure S5E,F). Immunotype 2 may represent an intermediate immune profile, as proposed in the “intermediate immunoscore” model.48 Quantitative analyses of immune cell densities alone did not allow for accurate stratification of patients according to their ICI outcome status (online supplemental figure S5G), prompting us to further characterize the spatial organization of immune cells.

Accordingly, we performed spatial analyses using a normalized mixing score that quantifies the proximity between two given cell types, normalized by their respective proportions.50 These analyses revealed that PD1+CD8+ T cells were predominantly localized within the tumor core compared with other T cell subsets (online supplemental figure S6A), and that their increased infiltration into the tumor correlated with increased aPD1 ICI efficacy (online supplemental figure S6B). DCs tended to spatially cluster into homotypic and heterotypic aggregates (online supplemental figure S6C–E), with cDC1 cells being closer to the tumor core than other DC subsets (online supplemental figure S6F). Comparisons between RF and REL patients demonstrated that the proximity of LAMP3+ DCs to CD8+ T cells or cDC2 cells was significantly correlated with favorable prognosis (online supplemental figure S6G).

The description of TAM heterogeneity and its association with ICI efficacy has emerged as a major focus of recent research; however, their spatial distribution still requires further analysis. Paired analyses within each tumor revealed that CD163 TAM and pSTAT1 TAM were spatially aggregated in different locations, with PDL1−CD163 TAM being primarily located within the tumor core, whereas pSTAT1 TAM and PDL1+ TAM were more frequently found at the tumor margin (figure 2C–E). Interestingly, increased tumor infiltration by PDL1+pSTAT1 TAM was associated with favorable ICI outcomes, and their colocalization with CD163 TAM was similarly associated with enhanced aPD1 efficacy (figure 2F). No other spatial interaction was significantly different between the RF and REL patients.

Altogether, these integrated analyses of the immune landscape highlight novel patterns of immune cell spatial organization and raise important questions about how the spatial interactions between these immune cells and melanoma cell subsets affect the efficacy of aPD1 ICI.

Differentiated and undifferentiated melanoma cells colocalize with functionally distinct macrophage populations

Given the highly organized spatial distribution of melanoma cells and the immune contexture within the tumor, we hypothesized that the composition of the TIME may differ within these melanoma cell-state-defined areas. To test this hypothesis, we integrated data from the three mIF panels—melanoma cell states, T cells, and TAM—across eight tumors containing at least one large patch of differentiated (C2-MITFhigh and/or C3-SOX10high) and undifferentiated (C4-NGFRhigh and/or C6-ZEB1high) areas (figure 3A,B). Differentiated and undifferentiated tumor cell-related areas were manually contoured on the fused images, and the differential MPR ratio was confirmed before quantifying TAM and T cell densities within each area (figure 3B,C). T cell densities did not significantly differ between the differentiated and undifferentiated areas (online supplemental figure S7A). However, TAM polarization exhibited a consistent shift, with a reduction in CD163 TAM density and a concomitant increase in pSTAT1 TAM density in differentiated areas as compared with undifferentiated areas within the same tumor, leading to an elevated pSTAT1/CD163 TAM ratio in differentiated areas (figure 3C and online supplemental figure S7B).

Melanoma cell states colocalize with specific TAM subsets. (A) Fused mIF images of n=8 tumors showing TAM, T cells and melanoma cells, fused on a single image using HALO AI software. (B) Differentiated and undifferentiated tumor areas within the same tumor were manually contoured and the densities of T cell and TAM subsets within each area were calculated. (C) Paired comparisons of TAM subsets between differentiated and undifferentiated areas. Paired Wilcoxon test. P values are indicated in the graphs. (D) Identification of tumor-enriched spots. First row: tumor-enriched spots labeled using the MLG signature from Pozniak et al, 2024 (left, binarized, with tumor-enriched spots in green), as compared with the pathological annotation of the tumor area (right, in dashed green). Second row: identification of MITFhigh (left, in blue) and MITFlow (right, in red) spots based on the expression of MITF. (E) Radar plots characterizing tumor-enriched MITFhigh and MITFlow spots compared with previously published gene signatures. (F) Spatial plot of GOBP pigmentation (blue) and peripheral nervous system development (red) in tumor-enriched spots. (G) Schematic diagram of the selection process of the tumor-enriched differentiated and undifferentiated spots and their neighboring spots (two adjacent rows corresponding to an approximate perimeter of 200 µm) in tumor samples (n=4). Deconvolution of TAM subsets within spots using the myeloid annotations from an in-house snRNAseq dataset from two treatment-naïve cutaneous metastases. UMAP showing the annotations of the myeloid compartment in this dataset. The selected spots were then clustered to test for a differential abundance between differentiated and undifferentiated tumor-enriched spots and their neighboring spots. Created with BioRender. (H) Milo differential abundance of TAM clusters between differentiated and undifferentiated tumor-enriched spots and their neighboring spots. FDR (false discovery rate) threshold set at 0.1, with correction for multiple comparisons. Gray dots represent non-significant neighborhoods according to the FDR threshold. GOBP, Gene Ontology Biological Process; IFN, interferon; mIF, multiplex immunofluorescence; MLG, melanoma lineage gene; MPR, melanoma plasticity ratio; snRNAseq, single-nuclei RNA sequencing; SRT, spatially resolved transcriptomics; TAM, tumor-associated macrophage; UMAP, uniform manifold approximation and projection.

To explore the TIME in differentiated versus undifferentiated areas in more detail, we performed SRT analyses with the Visium platform on four out of the eight previously described samples. Tumor-enriched spots were identified based on the expression of a previously established minimal lineage gene (MLG) signature,35 compared with manual pathological annotation, and subsequently stratified into MITFhigh and MITFlow spots according to normalized MITF expression (figure 3D and online supplemental figure S8A). Differential pathway analyses confirmed that MITFhigh spots largely overlapped with known differentiated or melanocytic signatures,51 thus defining differentiated tumor spots (figure 3E). In contrast, MITFlow spots exhibited gene expression profiles consistent with undifferentiated states, including the NCSC, invasive, and EMT signatures, thereby delineating undifferentiated tumor spots (figure 3E and online supplemental figure S8B). Gene Ontology Biological Process (GOBP) enrichment analyses further supported these distinctions, with pigmentation-related pathways being enriched in differentiated tumor spots, peripheral nervous system development, and neuron projection guidance pathways characterizing undifferentiated tumor spots (figure 3F and online supplemental figure S8C).

To more precisely delineate TAM transcriptomic profiles within the spatial context of differentiated versus undifferentiated tumor spots, we performed cellular deconvolution of the spots using the CARD algorithm52 (figure 3G and online supplemental figure S8D,E), leveraging an in-house single-nucleus RNAseq (snRNAseq) dataset from two cutaneous melanoma metastases, which is finely annotated for the myeloid compartment.53 Consistent with the increase in pSTAT1 TAM in differentiated areas quantified by mIF, analyses of TAM transcriptomic profile within tumor spots and their neighboring spots showed that differentiated tumor spots were enriched in C1QH TAM and IFN-primed TAM (IFN-TAM), a proinflammatory TAM subset overlapping with pSTAT1 TAM (figure 3H and online supplemental figure S8F–H). Conversely, undifferentiated tumor spots were primarily adjacent to proangiogenic TAM (angio-TAMs),53 a subtype of protumoral TAM that expresses hypoxia-related genes and promotes metastatic dissemination and EMT (figure 3H and online supplemental figure S8F–H).5355

Thus, these findings suggest that TAM polarization differs locally according to melanoma cell state, suggesting a preferential and spatially regulated crosstalk between TAM and melanoma cells.

Tumoral and microenvironmental cues govern TAM polarization

To investigate the molecular cues originating specifically from the tumor cells, we performed ligand-receptor interaction analyses in the SRT data with CellChat V.256 focusing on cell-cell communications inferred between tumor-enriched and adjacent macrophage-enriched spots (figure 4A,B and online supplemental figure S9A,B). Deconvolution analyses confirmed the enrichment of TAM in these spots compared with other stroma-enriched spots (online supplemental figure S9B). Differential analyses led to a list of ligand–receptor interactions, specifically originating from differentiated or undifferentiated tumor spots (figure 4B). We next analyzed in the scRNA-seq dataset from Pozniak et al35 the expression pattern of the ligands in the distinct melanoma cell clusters, as well as the level and specificity of expression of the receptors in TAM versus other cells from the TME (online supplemental figure S10). Among the top interactions enriched in the differentiated spots was interleukin-16 (IL16)-CD4 (figure 4B). While IL-16 was described to stimulate the production of proinflammatory cytokines,57 58 its function in the antitumor immune response remains poorly known.59 IL-16 was highly expressed in differentiated spots and predominantly expressed by the “interferon αß response” melanoma cell subset (figure 4C,D). IL-16 receptor, CD4, is not only expressed by T cells but also by a proportion of TAM (figure 4E and online supplemental figure S10A). IFNγ-induced T-cell attracting chemokines CXCL9, 10 and 11, which signal through the CXCR3 receptor, were also enriched in differentiated spots, and specifically expressed by the “Antigen presentation” melanoma cell subset (figure 4C,D). This is reminiscent of our previous results showing that a ZEB1-dependant decrease in CXCL10 melanoma cell-specific expression was associated with immune escape.32 Similarly, CCL5-CCR5 T cell-attracting signaling was enriched in differentiated spots, with CCL5 being specifically enriched in the “Antigen presentation” melanoma cell subset (figure 4C–E). Importantly, tumor cell production of CCL5 and CXCL9 coexpression was previously associated with prolonged survival and better response to ICI.60 Since CXCL9, 10 and 11 are also strongly produced by TAM on exposure to T cell-derived IFNγ, an amplification loop may further sustain this antitumoral immune ecosystem. Consistent with this, differentiated tumor spots exhibited a high expression of pSTAT1 signature and IFN-γ response pathways (online supplemental figure S11A), reminiscent of the strong IFN-γ signaling observed in Immunotype 1 tumors enriched in pSTAT1 TAM (figure 2A).

Differential cell–cell communication between melanoma cells and TAM according to the melanoma cell state. (A) Labeling of cell spots of patient #1, including adjacent TAM-enriched spots in yellow (TAM-enriched spots located in a radius of ≤250 µm to the nearest tumor-enriched spot) and distant TAM-enriched spots in light yellow (TAM-enriched spots located in a radius of >250 µm to the nearest tumor-enriched spot). (B) Bubble plot diagrams illustrating cell–cell communications, inferred by CellChat, between ligands from differentiated/undifferentiated tumor-enriched spots and receptors from adjacent TAM-enriched spots. Results from the four tumors. (C) Heatmap showing the gene expression as a z-score of the ligand-coding genes among malignant clusters in the Pozniak et al (2024) dataset. (D) Spatial plot of ligand-coding genes IL16 and CCL5 (differentiated spots) and SEMA3B and MIF (undifferentiated spots) on patient #1 sample. (E) Spatial plot of receptor-coding genes CD4 (receptor of IL16), CCR5 (CCL5), CD74 (MIF), and NRP1 (SEMA3B) on patient #1 sample. IL16, interleukin-16; MIF, macrophage migration inhibitory factor; TAM, tumor-associated macrophage.

On the other hand, in undifferentiated spots, various ligands already known to promote both melanoma cell invasive state and TAM polarization towards protumoral states, as well as tumor angiogenesis and resistance to ICI,61 were upregulated, with on top of the list, MIF (macrophage migration inhibitory factor) and its receptor CD74 (figure 4B–E). MIF is strongly expressed by “stress hypoxia” melanoma cells, and the targeting of MIF-CD74 signaling on TAM was shown to restore the antitumor immune response in melanoma.62 63 Semaphorins and their receptors, neuropilins and plexins may also play a central role at the crossroads between TAM attraction and angiogenesis. The SEMA3B isoform was indeed specifically increased in the “NCSC” and “mesenchymal-like” melanoma cell subsets (figure 4C,D), while expression of its receptor NRP1 by TAM (figure 4E and online supplemental figure S10A) was previously shown to regulate their entry into hypoxic niches.64 Consistently, undifferentiated tumor spots displayed a hypoxic environment65 (online supplemental figure S11A), along with enrichment of the angiogenic VEGFA-VEGFR1R2 pathway66 67 (online supplemental figure S11B,C). Notably, VEGFA is specifically expressed by “stress hypoxia” melanoma cells (figure 4C). Other ligands enriched in “NCSC” and “mesenchymal-like” cells, the signaling of which was known to directly regulate TAM polarization, notably included Midkine (MDK) and its receptor LRP168 69 (online supplemental figure S11B,C).

Overall, integrating spatially resolved protein-based, transcriptomic, and cell–cell communication analyses revealed that distinct melanoma cell states differentially engage with TAM within spatially conserved areas to modulate their polarization. Beyond TAM, these melanoma cell subsets also likely interact with other immune cells, particularly T cells, thereby potentially shaping functionally divergent tumor-immune ecosystems.

Distinct immune ecosystems emerge around melanoma cell state patches

To validate that distinct tumor-immune ecosystems coexist within primary melanoma tumors, we integrated data from our three mIF panels, namely melanoma cell states, T cells, and TAM, as described above (figure 3A), but across 45 tumors. To investigate the immediate surroundings of melanoma subsets, the borders of the C2-MITFhigh, C3-SOX10high, C4-NGFRhigh and C6-ZEB1high patches were computationally extended by 50 µm, defining an immune cellular neighborhood (CN) (figure 5A). Louvain clustering of 24,082 CNs based on their immune cell composition identified seven specific and recurrent CNs (figure 5B–D and online supplemental figure S12A,B). Interestingly, differentiated (C2-MITFhigh, C3-SOX10high) and undifferentiated (C4-NGFRhigh and C6-ZEB1high) patches displayed preferential immune CNs; differentiated patches were enriched in CN3 and CN7 (figure 5E) containing PDL1+ TAM co-occurring with CD8+PD1+ T cells, undifferentiated patches were enriched in CN5 characterized by the colocalization of PDL1−CD163+ TAM with CD8+PD1− T cells, and CN6 was enriched in PDL1− TAM (figure 5E and online supplemental figure S12C).

Crosstalk between melanoma cell subsets and immune cell neighborhoods defines different ecosystems. (A) Schematic diagram of the analysis pipeline for the identification of immune CNs around each C2/C3/C4/C6 tumor patch in fused mIF melanoma cells, T cells, and TAM images from 45 tumor samples. Created with BioRender. (B) UMAP of the Louvain clustering of the n=24,082 CNs across the 45 samples. Each dot represents one CN. (C) Heatmap of the most variable cell types across the seven identified CNs. Colors depend on the z-score of the most variable cell types between CNs. (D) Violin boxplots of the normalized counts of T cell and TAM subsets in the different CNs. The counts were normalized against the area of each CN and against the proportion of the given cell type across the whole slide. (E) UMAP (left) and Milo beeswarm plot (right) of the differential abundance of CNs between differentiated (C2/C3, in blue) and undifferentiated (C4/C6, in red) patches. FDR (false discovery rate) threshold set at 0.1, with correction for multiple comparisons. PD1, programmed cell death protein-1; PDL1, programmed death-ligand 1; TAM, tumor-associated macrophage; UMAP, uniform manifold approximation and projection.

Overall, our results suggest that patches of distinct melanoma cell states are associated with specific immune neighborhoods, thus defining preferred ecosystems.

Spatial tumor-immune ecosystems dictate the efficacy of ICI in melanoma

Next, we investigated whether these ecosystems influenced the efficacy of aPD1 ICI. Supporting this hypothesis, the seven CNs displayed differential enrichment between REL and RF patients, with CN3 and CN7 being associated with favorable outcomes, whereas CN5 and CN6 were enriched in REL patients (figure 6A). While CN analyses highlighted distinct microenvironmental ecosystems associated with patient outcomes, this approach was guided by predefined hypotheses centered on melanoma cell-centric neighborhoods. To gain a broader and more comprehensive understanding of the spatial determinants of prognosis, we complemented this analysis with an entropy-based framework to systematically address the prognostic impact of every possible pairwise colocalization between melanoma cell states, T cells, and TAM subsets (figure 6B). For this, we used a previously published analytical pipeline50 that quantified the colocalization between these cell types, weighted for their respective densities, to generate a unique colocalization score that accounts for inter-sample heterogeneity (See Methods and figure 6B). Using this approach, we identified specific pairwise colocalizations associated with either improved or worsened RFS following adjuvant aPD1 ICI (online supplemental figure S13A). Strikingly, PDL1+pSTAT1 TAM was the cell type for which the colocalization score was the most frequently linked to better RFS. Interestingly, despite similar overall densities of CD163 TAM in RF and REL patients (figure 2B), high colocalization between PDL1−CD163 TAM and either CD8+PD1− T cells or C3-SOX10high melanoma cells was significantly associated with poorer survival, suggesting an immunosuppressive role played by this TAM subset.

Tumor-immune spatial co-organization defines functional ecosystems driving ICI efficacy. (A) UMAP (left) and Milo beeswarm plot (right) of the differential abundance of CNs between RF-associated (in purple) and REL-associated (in orange) patches. FDR (false discovery rate) threshold set at 0.1, with correction for multiple comparisons. (B) Schematic diagram of the pipeline used to calculate the colocalization score of given cell types. Fishnet approach to divide each tumor into a grid containing 20×20 tiles, excluding empty squares to calculate the localized entropies between TAM, T cells and melanoma subsets. As an example, localized entropies of a pattern of two cell types (TAM and tumor cells) are displayed and compared with the global entropy calculated on the whole tumor. The colocalization score is the percentage of tiles with a localized entropy above the global entropy for a given pattern. (C) Ecosystem U1 composition and prognostic impact (n=45). In the upper row, the left panel represents one point per cell, the middle panel is the binned representation of the localized entropy of U1 cell types, and the right panel represents in red the tiles with a localized entropy above that of the global entropy of the given tumor (hotspots). In the lower row, an enlarged view of one area of high colocalization and boxplots of the colocalization score according to the 2-year status REL (relapse) or RF (relapse-free) of patients. Mann-Whitney U test: the p value is indicated on the graph. (D) Ecosystem D1 composition and its prognostic impact, as described in (B) for the U1 ecosystem. Mann-Whitney U test. (E) Description of the four main scenarios according to the proportion of U1 and D1 ecosystems across the 45 tumors (left) and their Kaplan-Meier curves (right). Log-rank test, the p value is indicated in the graph. CNs, cellular neighborhoods; D1, differentiated-1; ICI, immune checkpoint inhibitor; PD1, programmed cell death protein-1; PDL1, programmed death-ligand 1; REL, relapse; RF, relapse-free; RFS, relapse-free survival; TAM, tumor-associated macrophage; U1, undifferentiated-1; UMAP, uniform manifold approximation and projection.

Given the frequent recurrence of specific pairwise colocalization, we hypothesized that these interactions likely represented larger spatial ecosystems involving multiple heterotypic cell types. To explore this, we focused on cell types involved in detrimental colocalization, that is, C6-ZEB1high, C4-NGFRhigh, and C3-SOX10high melanoma cells, mostly corresponding to undifferentiated tumor areas, CD8+PD1− T cells, and PDL1−CD163 TAM. Entropy-based colocalization analyses revealed the presence of this undifferentiated-1 (U1) ecosystem across 43 (out of 45) tumors, which were enriched in REL-associated tumors (figure 6C). Conversely, by identifying recurrent pairwise interactions associated with favorable outcomes across the tumors, we identified the differentiated-1 (D1) ecosystem that was associated with prolonged RFS and was composed of C2-MITFhigh, C3-SOX10high, PDL1+pSTAT1, PDL1+CD163 TAM, CD8+PD1+ and CD4+PD1− T cells (figure 6D and online supplemental figure S13B,C).

Most patients displayed both U1 and D1 ecosystems; thus, we hypothesized that the relative balance between these ecosystems could influence aPD1 efficacy. We identified four distinct prognostic scenarios based on their frequency (figure 6E). The best-case scenario implied only a few hotspots of the U1 ecosystem associated with a high prevalence of the D1 ecosystem (n=13, n=1 REL), whereas the worst-case scenario consisted of a high number of hotspots in the U1 ecosystem associated with only a few to no hotspots of the D1 ecosystem (n=12, n=9 REL). In contrast, two opposite scenarios converged to the same intermediate prognosis: either the tumor was roughly equally enriched in the U1 and D1 ecosystems (n=11, n=3 REL) or both ecosystems were scarcely present across the tumor (n=7, n=3 REL). The clinical relevance of these ecosystems was further highlighted by their ability to predict the efficacy of ICI with higher accuracy than the Breslow index, American Joint Commission on Cancer, eighth edition (AJCC-8) staging, MPR alone, or immunotypes (online supplemental figure S13C). Interestingly, through regularized Cox regression on the MPR, immunotypes, ecosystems, and their interaction terms, we computed a composite score that yielded a high predictive value in discerning REL and RF patients, underscoring the clinical relevance of combining tumoral, immune, and spatial parameters (online supplemental figure S13D,E).

Overall, these findings suggest that distinct tumor-immune ecosystems coexist in spatially specific areas of melanoma tumors, and that their balance is a key determinant of aPD1 ICI efficacy.

Discussion

In this study, we leveraged integrative RNA-seq, mIF and SRT analyses from treatment-naïve stage III patients with melanoma treated with ICI to explore the spatial organization of ITH and immune cell context and its impact on patient outcome. Our findings revealed non-random spatially arranged tumor-immune ecosystems driven by heterotypic tumor-immune crosstalk, which correlated with aPD1 ICI efficacy.

Although prior studies have demonstrated the implication of melanoma cell plasticity in immune escape,293234 undifferentiated melanoma cells found in metastases have only recently been reported to predict resistance to aPD1 ICI following a single treatment cycle.35 Whether such undifferentiated melanoma populations contribute to ICI efficacy in the pretreatment setting and in primary cutaneous melanoma remains unclear. Whole-slide spatial profiling of treatment-naïve primary cutaneous melanomas revealed that melanoma cells sharing phenotypic identity colocalize together into homotypic spatial patches, mirroring the architecture observed in lymph node metastases.35 Notably, an enrichment in undifferentiated NGFRhigh and ZEB1high patches, relative to differentiated MITFhigh and SOX10high patches prior to treatment, was associated with decreased ICI efficacy. Interestingly, the same ratio of melanoma cells without accounting for these homotypic patches was not predictive of treatment outcomes, highlighting the importance of the spatial context. These findings suggest that the formation of homotypic tumor cell patches, rather than cellular composition alone, modulates the therapeutic response, with the balance between distinct melanoma cell states serving as a potential determinant of ICI efficacy. The non-random arrangement of melanoma cell states implies that signals originating from the TME, such as chemokines or cytokines, namely IFN signaling, orchestrate the emergence of these tumoral homotypic patches.70 Conversely, the organization of these patches may, in turn, influence the phenotypic states of the surrounding immune and stromal populations, suggesting crosstalk within the TME that may critically shape therapeutic outcomes. This is illustrated through unsupervised clustering of the immune cell composition surrounding these homotypic melanoma patches that unveiled recurrent immune CNs. Remarkably, differentiated melanoma patches were more frequently surrounded by CNs composed of PDL1+ TAM and PD1+CD8+ T cells than undifferentiated melanoma patches. This finding is reminiscent of the colocalization of IFN-responsive tumor cell states with macrophages and T cells observed in a pan-cancer analysis.71 Additionally, intratumoral analyses of TAM phenotypes identified the colocalization between CD163 TAM and undifferentiated melanoma niches that mirrored observations in primary acral lentiginous melanoma where M2-like APOE+CD163+ TAM promotes an undifferentiated melanoma cell state through the IGF1-IGF1R signaling pathway in spatially restricted regions.33

Here, our newly generated spatial transcriptomic data highlighted how spatially constrained microenvironmental cues, such as IFNγ and hypoxia, might engage both tumoral and non-tumoral compartments. Beyond these signals, cell–cell communication analyses performed on these spatial transcriptomic data uncovered several pathways through which differentiated and undifferentiated melanoma cell subsets might shape the phenotype of neighboring TAM. In hypoxic regions enriched in undifferentiated melanoma cells, increased production of a range of well-known proinvasive and proangiogenic molecules such as MIF, SEMA3-B, MDK, or VEGFA, may favor TAM polarization towards protumor subsets. In contrast, in regions enriched in differentiated melanoma cells, associated with a strong IFN-γ signaling, our results highlighted, besides IFN-dependent chemokines, increased production of poorly studied cytokines such as IL-16. Interestingly, decreased IL-16 production was recently associated with resistance to ICI in patients with breast and lung cancer, by altering CD4+ T cell function and TAM phenotypes in the TME.59 Further investigations will be needed to determine whether and how it may play a role in the antitumor immune response in melanoma. The precise network of cytokines and chemokines that may favor proinflammatory TAM and their dialogue with T cells will deserve further characterization.

Collectively, these findings support the growing paradigm of tumor ecosystems,72 whereby heterotypic interactions between melanoma cell states and immune populations modulate the efficacy of ICI. First, entropy-based analyses identified PDL1+ pSTAT1 TAM as the cell population most consistently associated with favorable ICI outcome. This result goes beyond prior reports implicating the proinflammatory, “M1-like” TAM in enhancing ICI response within spatially restricted tumor regions.1517 Indeed, the D1 ecosystem, enriched in RF patients, and marked by the previously reported colocalization of CD8+PD1+ T cells with PDL1+ pSTAT1 TAM,16 was found in preferential interactions with differentiated MITFhigh and SOX10high melanoma cells. Conversely, the U1 ecosystem, associated with metastatic relapse, was characterized by undifferentiated NGFRhigh, ZEB1high, and SOX10high—melanoma cells, colocalized with PDL1− CD163 “M2-like” TAM and CD8+PD1− T cells. These findings add to the growing body of evidence supporting the coevolution of ITH and immune cell polarization.17 32 33 71

Importantly, the association of the U1 ecosystem in the primary tumor with macroscopic metastatic relapse suggests a role for this ecosystem in fostering metastatic founder cells (MFCs). This hypothesis is supported by recent data in mice, which demonstrated that in primary tumors, melanoma cells with an undifferentiated/mesenchymal-like cell state can act as a reservoir of MFCs.73 Interestingly, recent data on human sentinel lymph node biopsies revealed that early disseminated cancer cells (DCC) entering the lymph node undergo a switch from a transitory phenotype to a more NCSC state on exposure to IFN.34 Thus, the heterotypic tumor-immune ecosystems identified here might also foster the outgrowth of DCC towards overt metastases.

Overall, while our study identifies novel spatial ecosystems strongly associated with clinical outcome, evaluation of their translational potential will require further validation in larger cohorts to determine whether these spatial signatures could serve as clinical biomarkers for stratifying patients for aPD1 therapy. Furthermore, longitudinal post-treatment analyses would provide additional insights into the predictive value of such ecosystems on aPD1 treatment. Moreover, although our spatial transcriptomic analyses offer valuable insights into the molecular pathways underlying these ecosystems and highlight promising targetable interactions to interfere with resistance, we acknowledge that the small number of samples used for these analyses is a putative limitation of this study. Mechanistic studies will be essential to elucidate their developmental origins and functional interdependencies. Future analyses will also enable us to assess the interest of these potential targets as novel therapeutic approaches targeting tumor-immune communications, which could be combined with aPD1 immunotherapy to improve efficacy.

In summary, our study suggests that melanoma and immune cells co-evolve together in spatially restricted areas through cell–cell communications and integration of microenvironmental signals, forming diverse heterotypic ecosystems that can either enhance or limit the efficacy of adjuvant aPD1 ICI, ultimately influencing macroscopic relapse.

Materials and methods

Clinical data of the patient cohorts

The cohort consisted of 57 pretreatment (FFPE) sections of the entire primary melanoma from 57 patients with resected stage III (AJCC-8) melanoma from Lyon Sud Hospital (Hospices Civils de Lyon, Lyon, France). Availability of a FFPE primary melanoma lesion was mandatory for including the patients in the cohort. Pathological review by a pathologist of the H&E slide allowed us to validate the presence of tumor tissue, but no additional criteria were applied. All patients were treated with a wide local excision of the tumor bed and complete lymph node dissection (if needed), followed by adjuvant aPD1 ICI with either nivolumab (480 mg intravenously every 4 weeks) or pembrolizumab (400 mg intravenously every 6 weeks) for 12 months (between 2019 and 2020). Patients should not receive any ICI or treatment with MAPK inhibitors prior to surgery. Patients were followed-up with cerebral, chest, and abdominal CT scans every 3 months for 3 years, every 6 months for 2 years, or until relapse after the first aPD1 infusion. All patients were followed for at least 2 years after the first immunotherapy infusion. Patients were classified according to their 24-month relapse status. Relapse was defined as either radiological recurrence, clinical progression according to the investigator, or death, whichever came first. A summary of the patients’ clinical and demographic (including age and sex) data is presented in table 1.

Multiplexed immunofluorescence staining

3 µm tissue sections were cut from FFPE human melanoma specimens. Sections were stained with five 7-color panels using OPAL technology (Akoya Biosciences, Marlborough, Massachusetts, USA) on a Leica Bond RX (Leica Biosystems, Wetzlar, Germany) on the Research pathology platform East (ANAPATH RECHERCHE). The panels included a melanoma cell plasticity panel targeting ZEB1 (RRID:AB_1844977), ZEB2 (RRID:AB_10603840), MITF (RRID:AB_10603840), SOX10 (RRID:AB_10844002), SOX9 (RRID:AB_1080067), NGFR (RRID:AB_10839265), a T cells panel targeting CD3 (RRID:AB_2631163), CD8 (RRID:AB_2075537), SOX10, ZEB1, PD1 (RRID:AB_2894867), KI67 (RRID:AB_2142367), a TAM panel targeting CD68 (RRID:AB_2074844), pSTAT1 (RRID:AB_10863745), CD163 (RRID:AB_2756375), PD-L1 (RRID:AB_2887952), SOX10, a DC panel targeting CLEC9A (RRID:AB_2884022), CLEC10A (RRID:AB_3714545), LAMP3 (RRID:AB_1148761), BDCA2 (RRID:AB_354762), CD8, SOX10, and a TLS panel targeting CD3, CD20 (RRID:AB_2282030), IgA (RRID:AB_3094763), IgG (RRID:AB_2335700), LAMP3, SOX10. After staining, the slides were digitized using Vectra Polaris (Perkin Elmer, Waltham, Massachusetts, USA).

Multiplex immunofluorescence image processing and cell segmentation

Visualization and annotation of the regions of interest of each scan were performed with Phenochart V.1.0.9 (RRID:SCR_019156, Akoya Biosciences). Quality control of the images was performed using slides with poor quality (n=4), and heavily pigmented lesions precluding immunostaining (n=2) were excluded. Segmentation of the scans was performed with inForm version 2.4.8 (Akoya Biosciences). The scans were then analyzed using the HALO-AI image analysis software V.3.5 (Indica Labs, London, UK).

Under HALO-AI, tissue and nuclear segmentation were performed using the MiniNet module and the NucSeg FL V.1.0.0 module, respectively. The slides were manually annotated: the first class was the “tumor”, and the second class the “stroma”. The “Analysis” module “FISH-IF 1.0.0” was used with the positivity threshold for the MFI of each OPAL channel set, depending on its nuclear or cytoplasmic/membrane affinity and visual control. The resulting output was a .csv file per slide containing all recognized cells.

Fusion and informatic reconstruction of the three panels under HALO-AI

The three immunofluorescence images were reconstructed using HALO-AI. The images were then registered and fused using the “serial qptiff fusion” module. Next, fusion was performed using HALO-AI, which applied deformation and translation to match the three images. The .csv files were extracted as the outputs.

RNA sequencing analyses

RNA extraction was performed using the CRCL Cancer Genomics platform on an Illumina NovaSeq machine with a paired-end protocol (64M reads). Libraries were prepared using an RNAexome FFPE kit (Illumina). Raw sequencing reads were aligned to the human genome (GRCh38) using STAR (V.2.7.8a, RRID:SCR_005622), with the annotation of known genes from gencode V.37. Gene expression was quantified using Salmon (V.1.4.0) and the annotation of protein-coding genes from gencode V.37. Differential gene expression between immunotypes 1 and 3 was analyzed using DESeq2 (V.1.42.1, RRID:SCR_000154). Genes with a log2FoldChange >1 and an adjusted p value (padj)<0.05 were kept for further analysis. A volcano plot was constructed using ggplot2 (SCR_014601). GSEA was performed using the enrichGO function of the clusterProfiler package (V.4.10.1, RRID:SCR_016884) to identify enriched GO terms based on the identified differentially expressed genes (DEGs). Finally, pathway enrichment results were sorted based on gene ratio and visualized with bar plots using ggplot2.

Visium data generation and preprocessing

Selection of cases relied on the eight tumors in which inter-area analyses were performed (figure 3A–C). Out of these eight tumors, only four could be further used for SRT following review of the H&E slide by the pathologist. Spatial transcriptomic data were generated by the CRCL (Cancer Research Center Of Lyon) Cancer Genomics Platform using Visium CytAssist Spatial gene expression for FFPE, Human Transcriptome, kit (10X Genomics), according to the manufacturer’s guidelines. Briefly, five micron FFPE sections were placed on histological slides (TOMO), deparaffinized, and stained with H&E by the Research Pathology Platform at CRCL. After imaging, 6.5×6.5 mm areas were selected in the tissue section for each patient. After unstaining and decrosslinking, human Visium probes were hybridized on the samples, ligated, and transferred onto the Visium capture slides using the CytAssist instrument (10X Genomics). Libraries were then prepared and sequenced on a NovaSeq 6000 sequencer (Illumina), targeting 50,000 reads/spot. Visium data were then processed using SpaceRanger from 10X Genomics (V.3.0.0), with the “count” algorithm and the human genome reference refdata-gex-GRCh38-2020-A for quantification. Transcript capture was performed using Visium_Human_Transcriptome_Probe_Set_v2.0_GRCh38-2020-A.csv probeset from 10X Genomics.

Identification of tumor-enriched spots

Spot labeling was performed by (1) histopathological reviews by our pathologist specializing in dermatopathology (MDo) and (2) using the minimal lineage gene (MLG) signature described by Pozniak et al,35 which is designed to discriminate melanoma cells from cancer-associated fibroblasts. The AddModuleScore (Seurat V.5.1.0) function was applied with a threshold of 4 to classify spots; an MLG score >4 indicated tumor enrichment of the spot, whereas an MLG score <4 indicated stromal enrichment of the spot. To further characterize the tumoral spots, the normalized expression of MITF was used to label them as MITFhigh (MITF-normalized expression >1.8) and MITFlow (MITF-normalized expression <1.5). To access the DEGs across spots, the FindAllMarkers function was applied with the following thresholds: avg_log2FC >1 and padj <0.05. The tool “What is my melanoma signature” (https://wimms.tanlab.org/)51 was used to generate the radar plots, based on the identified DEGs, allowing to characterize the spots based on public gene signatures. GSEA was also performed on the same DEGs using the fgsea function from the fgsea package (V.1.28.0) and the hallmark and GOBP database from the package msigdbr (V.7.5.1). GOBP pathways were visualized on Visium slides using AddModuleScore, the SpatialFeaturePlot function from the Seurat package, and ggplot2 (V.3.5.1).

snRNAseq data generation and preprocessing

The generation of the snRNA-seq dataset was performed following the snPATHOseq protocol from,74 74 using FFPE tissue samples from two treatment-naïve cutaneous metastases. The libraries were prepared using the GEM-Flex protocol. SnRNA-seq data were processed using CellRanger from 10X Genomics (V.7.1.0), using the “multi” algorithm and the human reference genome refdata-gex-GRCh38-2020-A for quantification. Transcript capture was performed using ChromiumHumanTranscriptomeProbeSetv1.0.1_GRCh38-2020-A.csv probeset from 10X Genomics.

Cellular deconvolution of the spots

To constitute an in-house reference dataset, snRNA-seq data from two treatment-naïve cutaneous metastases were merged. Quality control filtering was performed using the following thresholds: nFeature_RNA>150, nCount_RNA<1 50 000. Data were normalized using SCTransform (V.0.4.1) and regressed on the cell cycle and nFeature_RNA. Principal component analysis (PCA) was performed using runPCA. The batch effect was corrected using harmony (V.1.2.1). Clusters were defined using FindNeighbors and FindClusters. Finally, a uniform manifold approximation and projection (UMAP) is computed using RunUMAP. The doublets were removed using Scrublet (V.0.2.3) in Python (SCR_008394). The entire TME was largely annotated using signatures and gene markers from the literature. The myeloid compartment was then isolated and reclustered using the same pipeline and TAM signatures53 were used to annotate this compartment. The CARD (V.1.1) package52 was used to perform cellular deconvolution of the spots using the snRNAseq dataset composed of two cutaneous metastases as a reference.

Spots neighborhood analyses

Spot neighborhoods were defined based on the specific hexagonal geometry of the visium spots. A spot was considered to be within the neighborhood of a tumor spot if it was located within the first two rows (corresponding to an approximate radius of 200 µm) around the spots of interest and labeled as stroma-enriched. Spots that were simultaneously present in the list of differentiated and undifferentiated spot neighborhoods were labeled as mixed-neighboring non-tumoral spots (figure 3G). The spots were CLR-normalized based on the deconvoluted cellular fraction of TAM clusters, clustered using Seurat, and annotated according to the predominant TAM subtypes. The “IFN-TAM & C1Qh-TAM”, “Angio-TAM” and “MoMacs” spots were further kept for abundance analyses between differentiated versus undifferentiated spots using the miloR package (V.2.0.0), using a false discovery rate (FDR) threshold of 0.1 to account for multiple comparisons.75

CellChat analysis between tumor and TAM-enriched spots

TAM-enriched spots were identified using the AddmoduleScore function and the macrophage gene signature derived from MCPcounter signatures. Stroma-enriched spots with a macrophage score>0.5 were classified as TAM-enriched. Adjacent TAM-enriched spots were defined as spots located within the neighborhood of differentiated or undifferentiated tumor spots (as described above), in contrast to distant TAM-enriched spots. Cell-cell interactions were inferred using the CellChat package (V.2.1.2) with Contact=FALSE and a fixed distance of 200 µm between differentiated/undifferentiated spots and TAM-enriched spots within this spatial range. Conserved interactions were those that were either unique to one condition (differentiated or undifferentiated), or with a significantly higher probability in one condition compared with the other. Chord plots representing these interactions were assessed using CellChat and circlize (V.0.4.16) packages.

pSTAT1, interferon gamma response and hypoxia scoring

The pSTAT1 gene signature was scored among macrophage subsets using AddModuleScore. ViolinPlot was visualized using VlnPlot2 from the SeuratExtend package. The pSTAT1 signature was defined as a combination of three signatures7678 (online supplemental table S1). The IFN-γ response signature is extracted from the hallmark database. The hypoxic signature was described by Di Giovannantonio et al.65

The slides of all the patients were annotated to identify differentiated and undifferentiated spots, as described above. Then, all patient slides were merged, and the SCTransform was applied to normalize the data. The AddModuleScore was used to compute the score for each signature or pathway. Then, a Mann-Whitney U test was performed on each score, between differentiated and undifferentiated spots, and box-violin plots were used for visualization purposes, using the ggplot2 package.

Cut-off’s settlement

Every quantitative variable analyzed with the log-rank test was converted into a binary variable (high/low) according to a cut-off value automatically defined using Cutoff Finder V.1.79 Euclidean distance was used for cut-off determination.

Multi-IF spatial analyses

K-means clustering of tumor cells

All SOX10+ tumor cells across all samples were retrieved as a single .csv file containing the MFI of each phenotypic marker. Subsequently, k-means clustering was performed. The optimal number of clusters is determined using the silhouette score.

Melanoma cell states patches detection

For each tumor, the .csv file containing the spatial coordinates and phenotypes of every tumor cell was converted into a SpatialExperiment object. Next, a spatial interaction graph was built (“buildSpatialGraph” function) with a radius of 25 µm using the imcRtools R package (V.1.10.0).46 The “patchDetection” function was then run with a minimal number of cells of 5. This operation returns the SpatialExperiment object with a new column containing the patch identifier for each cell for a given melanoma cell state.

Spatial analyses of immune cells

For the clustering of DCs, the “FindNeighbor” function from the SPIAT (Spatal Image Analysis of Tissues) package (V.1.6.2) was used setting the minimal number of cells to form a cluster to 2 and the radius to 50 µm. Aggregates were defined as at least 10 DCs within a radius of 50 µm. The proximity score corresponded to the normalized mixing score (NMS) from the SPIAT package and was the main spatial metric used to account for the baseline density of both the tested immune cells.50 A radius of 25 µm was chosen for the NMS analyses.

Immunotypes identification

First, the patients for whom TAM and T cell parameters were available were retained (48 of 57 patients). Next, every TAM and T cell density significantly associated with RFS was included in the partitioning around medoid clustering of patients (non-parametric k-means clustering). All features were transformed using the following operation: x=>log10(1e-3+x) and z-normalized across the dataset prior to inclusion in the model. A silhouette score is computed to select the optimal number of clusters.

Comparisons of the immune cell densities between melanoma cell states areas

From the fused images of the eight selected tumors, homotypic differentiated (C2-MITFhigh, C3-SOX10high) and undifferentiated (C6-ZEB1high, C4-NGFRhigh) patches were manually contoured from the same tumors. Then, the surface area was calculated, the counts of TAM and T cells were extracted, and the density of each immune subset was computed within the different patches.

Immune cellular neighborhoods

From the fused images of the 45 tumors, melanoma patches of clusters C2-MITFhigh, C3-SOX10high, C4-NGFRhigh and C6-ZEB1high were retrieved, and their borders were extended by 50 µm with the function “findMilieu” from the imcRtools R package46 to constitute immune CN. Briefly, the number of TAM and T cell subsets was retrieved for each CN. Only CNs with at least 10 TAM and T cells were retained. The area of each CN was calculated by considering them as concave polygons with the functions “concaveman” from the concaveman package (V.1.1.0) and “st_polygon” from the sf package (V.1.0.19). Subsequently, the count of each immune cell was divided by the area of the given CN and then by its corresponding frequency across the entire section of the given tumor to account for interpatient heterogeneity. Next, a Seurat object was built with every CNs from the 45 tumors, and these CNs were clustered using the Louvain function (resolution=0.2) in the Seurat package (V.5.2.0, RRID:SCR_007322). To study the enrichment of the generated CN clusters within specific melanoma cell states and to link it with ICI efficacy, the miloR package (V.2.0.0)75 was used (k=30, d=9) with a subsampling rate of 10% and a FDR threshold set at 0.1.

Entropy-based spatial analyses

The entropy-based analysis, as proposed by Feng et al50 with the SPIAT, automatically applied a fishnet approach to divide the tumor into n*n tiles (“grid_metrics” function, n=20). Each individual tile returned a localized entropy that could then be compared with a threshold value. Next, a colocalization score was computed, consisting of the percentage of tiles with an entropy value above the threshold across the entire tumor. To dynamically determine this threshold for each cell types combination in order to account for the probability of these cell types to be colocalized by chance alone, and the intratumoral and intertumoral variation of the cell types’ density, we calculated the global entropy of a given tumor using the “calculate_entropy” function and the returned value would then become the threshold value. Empty tiles corresponding to the voids around the tumors were excluded.

Identification of the D1 ecosystem

The prevalence score of pairwise entropy-based interactions associated with an improved RFS (HR <1) with a p value <0.10 were included in a co-occurrence matrix across each individual tumor to determine recurrent patterns of cell type interactions using Jaccard’s similarity index from the “CooccurrenceAffinity” package (V.1.0.0).80 These interactions were then clustered based on this index using Leiden clustering to identify highly recurrent pairwise colocalization patterns that might represent broader and more complex ecosystems. This led to the identification of Ecosystem D1.

Regularized cox regression analyses

The Cox proportional hazard models were L1 (LASSO) regularized logistic regression models fitted using the glmnet R package (V.4.1.8, SCR_015505) (Friedman et al, 2010). Briefly, this approach enables the selection of variables with the strongest relationship to the outcome and decreases the risk of overfitting. For MPR, features were computed under the transformation x -=>log(1e-3+ x) and z-normalized across the dataset prior to inclusion in any model.

To calculate the composite score, MPR, immunotypes, and ecosystems were included in the Cox regression analyses. Regarding variable preprocessing, the MPR was binarized according to its cut-off value (“low” vs “high”). The immunotypes corresponded to the cluster number as described before (“1,” “2” or “3,” figure 2B). The ecosystem variable was the scenario number (figure 6E) depending on the U1 or D1 ecosystems’ prevalence (“1”, “2”, “3”, “4”).

Next, multicollinearity was measured with a variance inflation factor (VIF) on a Cox proportional hazard model including MPR, immunotypes, ecosystems, and their interaction terms. Every VIF value was equal to one, indicating that no collinearity was detected. Subsequently, the Cox model was fitted using the L1 penalty. A cross-validated error plot based on Harrell’s C index can then help in selecting the optimal λ value where the curve hits its minimum. This λ value was used in the finalized fitted Cox model. From this model, we extracted the coefficients (w) associated with each of the n included variables (v), and calculated the composite score for each patient as follows:

Composite score = w1*v1 + w2*v2 + … + wn*vn

Notably, in such models, positive coefficients were associated with worse RFS, whereas negative coefficients were associated with better RFS.

Comparisons between Cox proportional hazard models

For each clinical, histological, or biological prognostic variable, a univariate Cox proportional hazard model was computed, and Harrell’s C index was compared with that of a univariate model based on the Breslow index, which is the best prognostic variable in clinical practice. The “compareC” function from the compareC R package (V.1.3.2) was used.

ROC curve analyses

When performed, receiver operating characteristic (ROC) curves were generated with the pROC R package (V.1.18.5, SCR_024286) using the function “roc” with respect to predict the 2-year recurrence status. The area under the curve (AUC) was calculated with the function “auc”.

Survival analyses

A univariate Cox proportional hazards model was fitted for the tested parameters, and a log-rank test was performed to compare the different groups with respect to RFS. The time to event was computed as the difference in months between the date of diagnosis of the primary melanoma and the date of the first recurrence or death, whichever came first. For patients with RF, this date corresponded to the date of the last survival update.

Statistical analysis and plotting

An α risk of 5% was used for all the statistical analyses. Statistical analyses were performed using RStudio (SCR_000432) with R. Statistical tests used included non-parametric Mann-Whitney U tests for unpaired two-group comparisons, paired Wilcoxon tests for paired two-group comparisons, log-rank tests for survival analyses, non-parametric Kruskal-Wallis tests for three or more unpaired groups, and Friedman test for three or more paired groups with post-test between-group comparisons performed with paired Wilcoxon tests, in which p values were adjusted for multiple comparisons with the “Benjamini-Hojberg” method. For correlation tests, a non-parametric Spearman test was used. Regarding boxplots, the whiskers extend to the minimum and maximum values within 1.5 times the IQR from the 25th and 75th percentiles.

Supplementary material