Modeling phenotype switching towards ZEB1high/MITFlow NCSC/invasive state in vitro

In order to study the role of endogenous ZEB1 during phenotypic transitions in melanoma cells, we used two BRAFV600 patient-derived short-term cultures, established with a low number of passages after culture (GLO and C-09.10). These two short-term cultures display a ZEB1low/MITFhigh proliferative phenotype. As previously described, C-09.10 cells are highly melanocytic, while GLO cells tend towards a transitory state, with intermediate MITF expression [21]. To induce phenotype switching towards a ZEB1high/MITFlow state, cells were treated every 3 days, for up to 14 days, with the inflammatory cytokine TNFα, a known inducer of dedifferentiation in melanoma cells [24]. As expected, TNFα treatment decreased proliferation, but no significant cell death was observed (Supplementary Fig. 1A). Treatment with TNFα induced a switch towards a ZEB1high/MITFlow state in GLO cells (Fig. 1A, B). ZEB1 protein expression increased, while ZEB2 expression decreased upon treatment (Fig. 1A). Drastic down-regulation of MITF expression was observed after treatment. Interestingly, the expression of the NCSC marker NGFR (Nerve Growth Factor Receptor) [25, 26] and of the receptor tyrosine kinase AXL [27] was increased. Analysis of SOX10 and SOX9 expression [28] (Fig. 1A) confirmed a progressive switch, towards a putative undifferentiated state, losing SOX10, in favor of SOX9, according to the four melanoma cell states as per the nomenclature proposed by Tsoi et al. [29] (melanocytic, transitory, neural-crest like NCL and undifferentiated). mRNA expression levels of these markers of melanoma cell states were consistently modified with the protein (Fig. 1B), except for ZEB2 mRNA which was not modified, consistent with previous reports, suggesting an additional post-translational regulatory mechanism [20]. We confirmed in the RNA-seq dataset from Tsoi et al. the progressive increase in ZEB1 expression in the NCL and undifferentiated states, while ZEB2 was strongly expressed in the melanocytic state and gradually decreased during dedifferentiation (Supplementary Fig. 1B). This switch appeared to be reversible, since the withdrawal of TNFα promoted the return to baseline expression levels (Supplementary Fig. 1C).

Fig. 1: Modeling phenotype switching towards ZEB1high/MITFlow neural crest stem cell /invasive state in vitro.

Western blot A and RT-qPCR B analyses of ZEB1, ZEB2, MITF, NGFR, AXL, SOX10 and SOX9 expression after 7 (D7) and 14 (D14) days of TNFα (100 ng/mL) treatment in GLO cells. GAPDH was used as loading control. Histograms represent quantitative analyses of relative expression (n = 4 independent experiments). C Longitudinal intra-tumor heterogeneity characterization of MITF and NGFR expression by flow cytometry in GLO pMITF-GFP cells, upon TNFα treatment during 7 (D7) or 14 (D14) days. NGFR was marked by anti-NGFR antibody coupled with APC. The proportion of cells with MITF high, intermediate or low and with NGFR high, intermediate or low statuses is indicated. D Transwell migration assays in GLO cells upon TNFα treatment after 7 (D7) or 14 (D14) days. Cells were fixed after 24 h, the number of migrating cells is plotted (n = 3). E Incucyte assay showing the relative increase in cell death upon PLX4032 (500 nM) treatment over time, in cells previously treated with TNFα for 7 or 14 days. F RNA-seq analyses of GLO cells after 7 (D7) or 14 days (D14) of TNFα treatment. Heatmap of ssGSEA scores of the most relevant hallmarks and of melanoma states signatures from Hoek, Tsoi, Verfaillie and Widmer. Clustering Ward.D2 / distance: Euclidean. G ssGSEA scores of the melanoma signatures (Tsoi and Verfaillie) in cells treated for 7 (D7) or 14 days (D14) with TNFα (GLO) or TNFα + TGFβ (C-09.10). H Inference of transcription factors (TF) activity in gene expression data using VIPER algorithm. Barplot of DoRothEA TF Normalized Enrichment Score (NES) comparing untreated versus TNF treated (D14) GLO cells. Data are shown as the mean ± SEM. P values were determined by a two-tailed paired student t test B, D and ANOVA test E. Differences were considered statistically significant at *P ≤ 0.05, **P < 0.01 and ***P < 0.001. ns (non-significant) means P > 0.05.

In C-09.10 cells, TNFα was combined with TGFβ in order to ensure an efficient phenotype switching, evidenced by decreased MITF expression (Supplementary Fig. 1D). Drastic up-regulation of ZEB1 expression was associated with progressive ZEB2 protein level down-regulation upon TNFα + TGFβ treatment (Supplementary Fig. 2A, B). NGFR and AXL expression were also induced, however, the SOX10/SOX9 switch was not observed in this model.

Precise monitoring of intra-tumor heterogeneity during phenotype transitions over time was achieved by flow cytometry using stable cell lines established with a MITF promoter-GFP reporter construct and combined with NGFR membrane staining. As previously mentioned, GLO cells exhibit a transitory phenotype, with about half of the cell population harboring either a MITFhigh or a MITFlow state. TNFα treatment decreased the proportion of the MITFhigh population and led to the emergence of a MITFlow/NGFRhigh phenotype (representing about 70% of the population at day 7), before a transition towards a MITFlow/NGFRlow population occurred (representing 22% of the population at day 14) (Fig. 1C). Similarly, in C-09.10 cells, FACS analyses of NGFR in a MITF-GFP reporter line, showed a decrease in the MITFhigh population and a transient increase in MITFlow/NGFRint/high cell population upon TNFα + TGFβ treatment (Supplementary Fig. 2C).

Finally, we performed functional assays to validate the transition towards a more invasive state. Transwell migration assays validated that TNFα treatment progressively increased the migratory capacity of GLO cells (Fig. 1D). Consistent with their increased migratory capacity, the sensitivity of TNFα-treated GLO cells to the BRAF inhibitor (BRAFi) PLX4032 was also decreased compared to control cells, as assessed by Incucyte live-cell analysis (Fig. 1E). TNFα + TGFβ treated C-09.10 cells also exhibited increased resistance to BRAFi (Supplementary Fig. 2D). These data confirmed the transition towards a more invasive and targeted-therapy resistant state.

In order to further characterize dysregulated pathways, we performed RNA-seq at day 7 and day 14 after TNFα ± TGFβ treatment in GLO and C-09.10 cells. Pathway analyses of the 4531 differentially expressed (DE) genes at day 14 compared to untreated GLO cells (p < 0.001 & |lFC | > 1) (2490 up and 2041 down), confirmed a decrease in proliferation hallmarks, as well as an enrichment in TNFα signaling, inflammatory response, and EMT hallmarks (Fig. 1F and Supplementary Fig. 3A, C). We next analyzed previously described proliferative and invasive melanoma signatures from Hoek et al. and Verfaillie et al. [30, 31]. Gene Set Enrichment Analysis (GSEA) confirmed a progressive switch from a proliferative/melanocytic (untreated), towards a more invasive state upon TNFα treatment (Fig. 1F–G and Supplementary Fig. 3C). The NCSC (i.e. NCL) and undifferentiated state signatures from Tsoi et al. [29]. were also activated at day 14, while the transitory state signature decreased in this model. RNA-seq analyses confirmed that C-09.10 cells display a more melanocytic phenotype than GLO cells, but similar pathways and signatures were consistently altered in this model (Fig. 1G and Supplementary Fig. 3B, D, E). Computational inference of transcription factor (TF) activity, with the VIPER algorithm, confirmed MITF TF decreased activity (Fig. 1H and Supplementary Fig. 3F). Moreover, the activity of the Activator Protein-1 (AP-1) complex members JUN and FOS, the major regulators of the mesenchymal state [32], was induced upon TNFα treatment in both GLO and C-09.10 cells, as well as that of the NF-κB subunits (RELA and NFKB1) (Fig. 1H and Supplementary Fig. 3F). Of note, ZEB1 and ZEB2 TF activity could not be reliably assessed based on DoRothEA pan-cancer database, because of melanoma cell-type specificities.

Overall, we developed two suitable in vitro models of phenotypic transitions of melanoma cells towards ZEB1high NCSC-like and invasive states.

Determination of ZEB1 direct target genes during phenotype switching

In order to obtain a comprehensive view of endogenous ZEB1 direct target genes in a genome-wide manner, we performed chromatin immunoprecipitation coupled to deep sequencing (ChIP-seq) analyses with an anti-ZEB1 antibody, in untreated (ZEB1low) and TNFα-treated (ZEB1high) GLO cells at day 14, but also in the A375 melanoma cell line, which displays a ZEB1high NCSC-like expression pattern (MITFlow, NGFRhigh, SOX10+, SOX9−) (Fig. 2). Consistent with increased ZEB1 expression, twofold greater ZEB1 peaks were found in TNFα-treated GLO cells compared to untreated cells (Fig. 2A). 33% of ZEB1 peaks were conserved between ZEB1low (untreated) and ZEB1high (TNFα-treated) cells, while 67% were acquired upon TNFα-induced ZEB1 expression (Fig. 2B). Interestingly, 44% of ZEB1 peaks observed in TNFα-treated GLO cells were conserved in A375 cells (Fig. 2A, B). A large proportion of ZEB1 peaks (62% in TNFα-treated GLO cells) were found in promoter regions (−1000, +0 bp), predominantly centered on the Transcription Start Site (Fig. 2C, D). Motif enrichment analyses of the top 1000 peaks in TNFα-treated GLO cells, when considering a 50 bp region centered on the ZEB1 peak, confirmed an enrichment in the ZEB1 binding motif, which ranked first (Fig. 2E), sustaining the notion of a direct binding of ZEB1. Motif enrichment for GFX and AP-1 complex were also evidenced. Consistently, the clustering of the peak density of ZEB1 with publicly available ChIP-Seq data of FOSL2 [33], a member of AP-1 complex, revealed a co-occupancy of FOSL2 at ZEB1 binding loci (Supplementary Fig. 4A).

Fig. 2: ZEB1 ChIP-sequencing analyses in two melanoma cell lines.

ZEB1 ChIP sequencing was performed in GLO cells, untreated or after 14 days (D14) of TNFα treatment and in A375-AS3 (control) cells. A The number and size of peaks are indicated. B Venn diagram showing the overlap between peaks found in GLO cells untreated (green) or after 14 days of TNFα treatment (D14) (red) and in A375-AS3 cells (blue). A hypergeometric test confirmed that overlaps between the three conditions are significant and not by chance alone (P < 0.001). Localization of the peaks C and distance to the transcription start site D. E Top3 HOMER-identified enriched motifs in GLO after 14 days of TNFα treatment. The associated p-values, the percentages of motif representation on target and background are indicated. F Heatmap of all genes presenting a ZEB1 binding peak in TNFα-treated cells at day 14. The most significantly enriched hallmarks and melanoma state signatures are indicated on the right. Clustering Ward.D2 / distance: Euclidean. G Integration of ChIP-Seq with RNA-seq data in GLO cells. Heatmap of DE genes presenting a ZEB1 peak in GLO cells after 14 days of TNFα treatment. Presence of a ZEB1 peak in the gene is indicated by a green line (untreated condition) or a red line (TNFα D14 condition) on the right. The most significantly enriched hallmarks and melanoma state signatures within down- and up-regulated genes in D14 versus untreated cells are indicated on the right.

We initially focused on all genes presenting a ZEB1 binding peak in GLO cells after TNFα treatment (Fig. 2F). Pathway analyses on all these genes showed, aside from proliferation (E2F targets, G2M checkpoint) and inflammation (TNF, interferon signaling) hallmarks, a striking enrichment in melanoma signatures (Fig. 2F). Next, integration with RNA-seq data (TNFα differentially expressed genes, p < 0.001 and |lFC | > 1) allowed us to correlate binding of ZEB1 with transcriptional up- or down-regulation (Supplementary Fig. 4B). Interestingly, a significant enrichment in the proportion of genes bound by ZEB1 was found, with 45% of dysregulated genes during phenotype switching exhibited a ZEB1 peak, 50% among down-regulated genes, and 40% among up-regulated genes (Supplementary Fig. 4C). Pathway analyses on differentially expressed genes displaying a ZEB1 peak demonstrated that ZEB1 directly binds to down-regulated genes involved in proliferation hallmarks and to up-regulated genes involved in TNF signaling and invasion/EMT hallmarks (Fig. 2G). Highly significant enrichment in invasive, undifferentiated and NCSC melanoma signatures was unveiled among up-regulated genes presenting a ZEB1 binding peak. Importantly, ZEB1 was more frequently bound in baseline conditions, in untreated ZEB1low cells, to genes that are down-regulated upon phenotype switching, while it was more frequently recruited de novo to genes that are activated during the transition to the NCSC-like and invasive states (Fig. 2G).

A further analysis of melanoma phenotype signatures from Hoek et al. [31]. demonstrated that ZEB1 binding peaks were present in 38% of genes of the proliferative signature, which are down-regulated upon phenotype switching (among which MITF, PMEL, CDH1, RAB38, or ASAH1) (Fig. 3A); 44% of invasive signature genes (which are activated upon phenotype switching) also displayed a ZEB1 binding peak (including ZEB1 itself, AXL, EGFR, BIRC3, THBS1, ITGA2, and ITGA3) (Fig. 3A). With respect to Tsoi et al. signatures, ZEB1 binding peaks were found in the promoter of 43% of genes of the melanocytic signature (including TSPAN10), 35% of NCL signature genes (among which TGFA and TGFBI) and 33% of undifferentiated signature genes (including EGFR again, CITED2, KRT7, KRT18, KRT80, AJUBA) (Fig. 3B). Only 24% of transitory signature genes displayed ZEB1 binding peaks (Fig. 3B). ZEB1 peaks were also identified in 44% and 34% of genes from the proliferative and invasive signatures from Verfaillie et al., respectively (Supplementary Fig. 4D). Importantly in the A375 cell line, the percentages of genes presenting a ZEB1 binding peak were largely similar to those described in GLO cells, more specifically in the invasive, NCL and undifferentiated signatures (Fig. 3A, B and Supplementary Fig. 4D), highlighting the overall conservation in ZEB1 binding specificity in these two ZEB1high melanoma models.

Fig. 3: Specific analyses of ZEB1 binding on genes from melanoma cell state signatures.

Heatmap of genes from the melanoma signatures from Hoek et al. A, Tsoi et al. B in untreated or TNFα-treated GLO cells at day 14. The presence of a ZEB1 peak is indicated by a green (untreated), red (TNFα D14) or blue (A375-AS3 cells) square. The percentage of genes of the corresponding signature presenting a ZEB1 peak for each condition is indicated. Clustering Ward.D2 / distance: Euclidean.

Overall, combined RNA-seq and ChIP-seq analyses performed in two models, led to the identification of novel ZEB1 direct target genes, specific to the melanocytic lineage, including down-regulation of proliferative/melanocytic genes and up-regulation of NCSC and undifferentiated genes.

ZEB1 directly regulates the expression of lineage-specific major markers of melanoma cell states

We then focused on major markers of melanoma cell states. ZEB1 was already bound, in untreated GLO cells, to the promoters of ZEB2, MITF and SOX10, the expression of which is down-regulated upon phenotype switching towards a ZEB1high state (Fig. 4A). In contrast, a ZEB1 peak was acquired de novo during phenotype switching in genes that are activated, such as ZEB1 itself, NGFR and AXL. Although no statistically significant peak was identified at the SOX9 locus, a ZEB1 binding signal seems to be increased in TNFα-treated cells. ZEB1 binding peaks were also found in other major markers of melanoma cell identity, namely BIRC3, ITGA2 and EGFR, which are up-regulated upon phenotype switching towards the ZEB1high state as confirmed by RT-qPCR (Supplementary Fig. 5A, B). Importantly, most ZEB1 binding peaks were conserved in A375 cells (Fig. 4A and Supplementary Fig. 5A).

Fig. 4: Validation by ChIP-qPCR of ZEB1 binding on the promoters of lineage-specific major markers of melanoma cell states.

A UCSC genome browser captures showing ZEB1 binding peaks in ZEB1, ZEB2, MITF, NGFR, AXL, SOX10 and SOX9 promoters in untreated or TNFα-treated GLO cells at day 14 (TNF D14), A375-AS3 cells and MDA-MB-231 cells. Significantly enriched peaks are marked by a grey square and red line. ZEB1 ChIP-qPCR on ZEB2, MITF, SOX10, NGFR and AXL promoters in GLO cells treated with TNFα for 14 days (D14) B and A375 cells C. Anti-ZEB1 (α ZEB1) or control IgG were used for the IP. Relative promoter enrichment was normalized against chromatin inputs (n = 3). Data are shown as the mean ± SEM. P values were determined by a two-tailed student t test.

To analyze the lineage specificity of ZEB1 binding compared to carcinoma models, we performed a comparative analysis with a previously published ZEB1 ChIP-seq dataset performed in the MDA-MB-231 breast cancer cell line [23]. Interestingly, while some ZEB1 peaks were conserved in MDA-MB-231, peaks in ZEB2, SOX10 and NGFR were only found in melanoma cells (Fig. 4A). Other genes involved in oncogenesis bore melanoma-specific ZEB1 binding sites which were either lost in MDA-MB-231 cells, such as for the WNT regulators TLE4 (Groucho) SFRP1 or FOXC2 [34] (Supplementary Fig. 5C). Furthermore, several melanocyte differentiation-related genes, namely the anti-apoptotic gene BCL2, a known MITF target [35], and BLOC1S5, the mutation of which is associated with defects in pigmentation [36], also displayed melanoma-specific ZEB1 binding peaks, further supporting lineage specificity of ZEB1 binding in melanoma cells compared to carcinoma cells.

Binding of ZEB1 to the sites defined by ChIP-seq was then validated by ChIP-qPCR in both GLO and C-09.10 models, upon phenotype switching towards a ZEB1high state and in A375 cells (Fig. 4B, C and Supplementary Fig. 5D). Consistent with increased ZEB1 expression, an enrichment in the binding of this TF to the promoters of ZEB2, MITF, NGFR, and SOX10 was observed in A375 cells, and in the GLO and C-09.10 models upon phenotype switching.

In order to reinforce the driving role of ZEB1 in the regulation of these genes, their expression was analyzed upon ZEB1 over-expression in C-09.10 and ZEB1 knock-out in A375 cells (Fig. 5A, B). A pair of A375 ZEB1 control (AS3) and knocked-out (AZ1) clones was analyzed, which did not show any defect in proliferation, albeit ZEB1 knock-out in A375 cells was confirmed to decrease cell migration (Fig. 5C), further validating the key role of ZEB1 in this process. MITF and NGFR expression following ZEB1 dysregulation has previously been described [21]. Here, we further witnessed that AXL is upregulated by ZEB1, following the same expression profile as NGFR. AXL expression decreased upon ZEB1 knock-out in A375 cells. We further demonstrated that SOX10 expression decreased upon ZEB1 over-expression in C-09.10 although its expression was not further increased upon ZEB1 knock-out in A375 cells. Inversely, SOX9 expression increased upon ZEB1 over-expression in C-09.10, and remained low upon ZEB1 knock-out in A375 cells. To go beyond the analyses of these markers, RNA-Seq analyses upon ZEB1 over-expression were performed in C-09.10, further demonstrating the increase in EMT/invasion pathways (Fig. 5D). The most significantly enriched melanoma signatures among activated genes were the invasive and the NCSC, while the undifferentiated signature was only partly induced, suggesting that ZEB1 is able to promote the expression of some but not all undifferentiated markers (Fig. 5D, E). We indeed validated that ZEB1 ectopic expression in C-09.10 was associated with an increase in additional NCSC/invasive markers, such as ITGA2, BIRC3 and AQP1, but also to a strong upregulation of TNFAIP2 and SECTM1, which are both markers of the undifferentiated phenotype [29] (Fig. 5F). Consistently, the knock-out of ZEB1 in A375 cells lead to a significant decrease in the expression of some of these undifferentiated markers, notably TNFAIP2, as well as KRT7, which is not detected in C-09.10 but strongly decreased in A375 (Fig. 5G). SECTM1 was not down-regulated in A375 suggesting that other factors may compensate for the loss of ZEB1 in this model.

Fig. 5: ZEB1-dependent regulation of markers of melanoma cell states in gain or loss of function models.

Western blot analyses of phenotype markers (ZEB1, MITF, NGFR, AXL, SOX10, SOX9) in C-09.10 cells with ZEB1 over-expression (ZEB1) A and in A375 control (AS3) or ZEB1 knocked-out (AZ1) clones B. GAPDH was used as loading control. (n = 3). C Transwell migration assays in A375-AS3 and A375-AZ1 ZEB1 knocked-out cells. Cells were fixed after 24 h, the number of migrating cells is plotted (n = 4). RNA-seq analyses of C-09.10 cells upon ZEB1 over-expression (ZEB1 OE). D The most significant hallmarks and melanoma signatures enriched in up-regulated genes in ZEB1 OE are indicated. E ssGSEA scores of invasive, NCL and undifferentiated melanoma signatures are plotted in C-09.10 cells upon ZEB1 over-expression or upon TNFα + TGFβ treatment at day 7 and day 14 for comparison. F RT-qPCR analyses of ITGA2, AQP1, BIRC3, TNFAIP2 and SECTM1 expression in C-09.10 cells with ZEB1 over-expression. Histograms represent quantitative analyses of relative expression (n = 3 independent experiments with two technical replicates for each). G RT-qPCR analyses of TNFAIP2, KRT7 and SECTM1 expression in A375-AS3 and A375-AZ1 ZEB1 knocked-out cells. Histograms represent quantitative analyses of relative expression (n = 3 independent experiments with two technical replicates for each). H Gene filtering strategy used to define the ZEB1.mel melanoma specific regulon. I Inference of transcription factors (TF) activity in gene expression data using VIPER algorithm with ZEB1.mel added to the list of regulons. Barplot of DoRothEA TF Normalized Enrichment Score (NES) comparing control versus ZEB1 OE C-09.10 cells. Data are shown as the mean ± SEM. P values were determined by a two-tailed paired student t test (C), and by a two-tailed unpaired student t test (F) and (G). Differences were considered statistically significant at *P ≤ 0.05, **P < 0.01 and ***P < 0.001. ns (non-significant) means P > 0.05.

Given the lack of relevance of the ZEB1 regulon from the DoRothEA database in the context of melanoma (Fig. 1H and Supplementary Fig. 3F), we processed our data to define a melanoma specific ZEB1 regulon (referred to as ZEB1.mel), that would be a useful tool for the scientific community. To achieve this, we selected the intersection of genes bound by ZEB1 in two cell lines, GLO cells upon TNFα treatment and A375_AS3 cell lines, with the highly differentially expressed genes in GLO upon TNFα treatment (Fig. 5F, Supplementary Table 6). Interestingly, there are no common genes between the ZEB1.mel and the ZEB1 pan cancer regulon. Subsequently, we tested the ZEB1.mel regulon in TNFα-treated and ZEB1 overexpressing C-09.10 cells. While no enrichment in the activity of the conventional ZEB1 regulon was observed, we could show an increase in the melanoma-specific ZEB1.mel regulon in these models (Fig. 5G and Supplementary Figs. 3F and 6A). Regarding other transcription factors, similarly to TNFα + TGFβ treatment, overexpression of ZEB1 was sufficient to induce decreased activity of MITF TF, as well as increased activity of RELA and JUN (Fig. 5G). Altogether, these results further support the conclusion that a major part of TNFα-mediated transition towards a more invasive/NCSC-like state is regulated by ZEB1.

Single-cell and spatial analyses of ZEB1 intra-tumor heterogeneity in melanoma patient samples

To further investigate the correlation between the expression of ZEB1 and markers of melanoma cell states at the single-cell level, we first used publicly available single-cell RNA-seq datasets. In the scRNA-seq dataset of melanoma PDX tumors from Rambow et al. [15], ZEB1 expression was significantly increased in the NCSC and the invasive populations compared to the pigmented state (Fig. 6A). Importantly, ZEB1 is part of the top 200 genes enriched in the NCSC signature. ZEB1 expression was also found in intermediate states, including the SMC (Starved Melanoma Cell) phenotype (Fig. 6A). Additionally, in the panel of 10 melanoma cell lines from Wouters et al. [37], including the A375 cell line, ZEB1 was preferentially expressed in mesenchymal-like cells, while SOX10 and MITF were found in melanocyte-like cells (Fig. 6B). Of note, both ZEB1 and SOX10 were expressed in the A375 cell line, which displays a neural-crest-like phenotype. We then investigated the relevance of the ZEB1.mel regulon in single-cell RNA-seq data. ZEB1.mel regulon displayed an increased activity in mesenchymal-like cell lines and in the neural-crest-like A375 cell line while the currently used Dorothea ZEB1 regulon did not show any significant variation between states (Fig. 6C, D). Consistently, ZEB1.mel regulon displays an antagonistic pattern when compared to MITF regulon. Moreover, ZEB1.mel regulon activity is increased in three short-term cultures of invasive-like switching induced by SOX10 knock-down (Fig. 6E). We further confirmed the specificity of ZEB1.mel regulon in patient single-cell RNAseq data from Pozniak et al. [38], we observed an enhanced activity of ZEB1.mel regulon in mesenchymal cells (Fig. 6E and Supplementary Fig. 6B, C) which was not observed with the Dorothea ZEB1 regulon. Altogether, we were able to create and validate the melanoma-specific ZEB1 regulon and confirm the activity of ZEB1 in mesenchymal and neural-crest-like cells in melanoma at single cell level.

Fig. 6: ZEB1 expression and melanoma specific ZEB1 regulon activity in public single cell RNA-seq dataset of melanoma models.

A ZEB1 expression levels in the different cell states defined by Rambow et al., 2018 in single cell RNA-seq data of melanoma patient-derived xenografts (PDXs). P values were determined by Mann-Whitney test. Differences were considered statistically significant at *P ≤ 0.05, **P < 0.01 and ***P < 0.001. ns (non-significant) means P > 0.05. B UMAP visualizations of single-cell RNA-seq data of 10 melanoma cell lines from Wouters et al. The expression levels of ZEB1, MITF and SOX10 genes are indicated. C UMAPS visualizations with transcription factor activity of ZEB1 regulon given by Dorothea collection (pancancer), the melanoma-specific ZEB1.mel and MITF regulons (Wouters et al.). D Violin plots showing transcription factor activity of ZEB1 and ZEB1.mel regulons (Wouters et al.). E. Violin plot of the transcription factor activity of ZEB1.mel regulon in Wouters et al. single cell RNA-seq data from 3 melanoma short term cultures transfected with SOX10 siRNA or non-targeting control (NTC). F Violin plot of the transcription factor activity of ZEB1.mel regulon in Pozniak et al. single cell RNA-seq data.

In order to further investigate ZEB1 co-expression or antagonistic expression with markers of melanoma cell states in patient samples, we performed spatial multi-immunofluorescence analyses (7 colors, OPAL, Perkin-Elmer) in a cohort of 30 cutaneous melanomas, previously annotated for ZEB1 expression as low, int or high. We analyzed at single-cell resolution the protein expression levels of ZEB1, ZEB2, MITF, NGFR, SOX10 and SOX9 to precisely define the frequency and spatial organization of the different phenotypes (Fig. 7A). This technique enables the specific quantification of the level of expression of markers of melanoma cell state in ZEB1-expressing melanoma cells, and excludes other cells from the microenvironment [39].

Fig. 7: Single-cell spatial analyses of markers of melanoma cell states according to ZEB1 status in melanoma samples.

A 7-color multiplexed immunofluorescence analyses of human melanoma samples with ZEB1 (red), ZEB2 (white), SOX10 (blue), SOX9 (yellow), NGFR (orange), MITF (green) and DAPI. Representative pictures of a ZEB1low (top) and a ZEB1high (bottom) cutaneous melanoma showing antagonistic expression of ZEB1 and ZEB2, MITF and NGFR, and SOX10 and SOX9. B–D. Reconstruction of three representative heterogeneous tumors as whole slides: a ZEB1lowB, and two ZEB1high tumors C, D. Each dot represents a cell. The expression levels of ZEB1, ZEB2, MITF and SOX9 are indicated in red (high), yellow (intermediate), blue (low) and grey (not expressed). SOX10 display 3 levels of expression defined as high, low and not expressed; and NGFR display only 2 levels of expression defined as high and not expressed.

Spatial reconstitution of the intensity of each marker at whole tumor level, revealed differential patterns of expression (Fig. 7B–D). Thin primary melanomas (such as MM28, Breslow = 0.8 mm) (Fig. 7A) or thick primitive melanomas (such as MM25, Breslow = 16 mm) (Fig. 7B) displayed a proliferative/differentiated ZEB2+ MITF+ SOX10+ phenotype with no or low ZEB1, SOX9 and NGFR expression (Fig. 7A, B). A gain in ZEB1 expression could be observed not only at the invasive front of primary melanomas, but also in the bulk, either in specific clones (as in MM10) (Fig. 7C) or in the whole tumor (as in MM14) (Fig. 7D). As shown in previous cohorts by immunohistochemistry [18, 21], antagonistic expression patterns of ZEB1 and ZEB2, as well as ZEB1 and MITF were confirmed at the intra-tumoral level by immunofluorescence (Fig. 7A, C, D).

We next analyzed SOX10 intra-tumor heterogeneity in these ZEB1high tumors and evidenced cell populations with decreased SOX10 intensity. As illustrated in the MM10 tumor, which bore the presence of a well-defined ZEB1high clone, increased ZEB1 expression was not only associated with low ZEB2 and MITF expression, but also with decreased SOX10 lev