ZEB1 controls a lineage-specific transcriptional program essential for melanoma cell state transitions

Cell plasticity sustains intra-tumor heterogeneity and treatment resistance in melanoma. Deciphering the transcriptional mechanisms governing reversible phenotypic transitions between proliferative/differentiated and invasive/stem-like states is required. Expression of the ZEB1 transcription factor is frequently activated in melanoma, where it fosters adaptive resistance to targeted therapies. Here, we performed a genome-wide characterization of ZEB1 transcriptional targets, by combining ChIP-sequencing and RNA-sequencing, upon phenotype switching in melanoma models. We identified and validated ZEB1 binding peaks in the promoter of key lineage-specific genes crucial for melanoma cell identity. Mechanistically, ZEB1 negatively regulates SOX10-MITF dependent proliferative/melanocytic programs and positively regulates AP-1 driven invasive and stem-like programs. Comparative analyses with breast carcinoma cells revealed lineage-specific ZEB1 binding, leading to the design of a more reliable melanoma-specific ZEB1 regulon. We then developed single-cell spatial multiplexed analyses to characterize melanoma cell states intra-tumoral heterogeneity in human melanoma samples. Combined with scRNA-Seq analyses, our findings confirmed increased ZEB1 expression in Neural-Crest-like cells and mesenchymal cells, underscoring its significance in vivo in both populations. Overall, our results define ZEB1 as a major transcriptional regulator of cell states transitions and provide a better understanding of lineage-specific transcriptional programs sustaining intra-tumor heterogeneity in melanoma.


INTRODUCTION
Cell plasticity contributes to intra-tumor heterogeneity and sustains tumor adaptation and treatment resistance [1].Cutaneous malignant melanoma is an aggressive form of skin cancer arising from melanocytes.Despite recent advances in targeted therapies and immunotherapies for the treatment of metastatic melanoma, nearly 60% of patients still develop resistance.A major mechanism of resistance to treatments relies on the ability of cancer cells to adapt to their environment and change phenotypes, i.e. to display cellular plasticity [2].One of these non-genetic adaptive mechanisms, the epithelial-mesenchymal transition (EMT), is a reversible trans-differentiation process finely regulated by a network of transcription factors (EMT-TFs) belonging to the SNAIL, TWIST and ZEB families [3].EMT-TFs are aberrantly reactivated in many cancers [4][5][6], particularly in carcinomas, where they play an oncogenic role by fostering metastasis and endowing cells with stem-like features [7,8].Although EMT cannot be formally defined in non-epithelial cancers, a related process of cellular plasticity contributes to intratumor heterogeneity (ITH) in melanoma and relies on reversible phenotypic transitions between proliferative/differentiated and invasive/stem-like states [9].Loss of MIcrophthalmia-associated Transcription Factor (MITF), the master regulator of melanocyte differentiation, induces a reprogramming towards an invasive and stem-like phenotype in melanoma cells [10][11][12].Gene expression analyses of tumors at the single-cell level refined this phenotype-switching model, by including the description of intermediate states and major molecular regulators [13][14][15][16].Reprogramming towards a Neural Crest Stem Cell-like (NCSC) phenotype was proposed as an adaptive response to targeted therapy, accounting for therapy resilience [17].However, a deeper understanding of cellular and molecular mechanisms underlying phenotypic adaptations and thus, the exceptional capacity of melanoma cells to develop resistance to current therapeutic strategies, is still needed.
We previously showed that a switch from ZEB2 to ZEB1 expression is a poor prognostic factor in melanoma [18,19].ZEB2 is expressed in normal melanocytes and its expression progressively decreases during transformation to melanoma, while ZEB1 expression increases.ZEB2 supports melanoma cell proliferation and differentiation by activating MITF expression [20].ZEB1, on the contrary, inhibits MITF expression, promotes transition to an invasive phenotype and resistance to targeted therapies [21].Though direct target genes of ZEB1 have been characterized in carcinoma models [22,23], they remain unknown in melanoma.Nonetheless, cell type-specific targets are particularly expected, given the antagonistic functions of ZEB1/ZEB2 in melanoma.
Herein, in order to characterize ZEB1 function and provide a comprehensive view of its transcriptional target genes in a genome-wide manner, we performed ChIP-sequencing combined with RNA-sequencing upon phenotype switching in melanoma cells.We define ZEB1 as a major transcriptional regulator of genes associated with phenotypic transitions in melanoma.Specific markers were validated as ZEB1 direct target genes, upon ZEB1 gain-or loss-of-function.Their relevance in human samples was further addressed through single cell multiplexed spatial analyses and analyses of public single-cell RNA-Seq datasets.Intra-tumor heterogeneity of markers of melanoma cell states according to ZEB1 expression in human samples demonstrated co-expression of ZEB1 with both stem-like and invasive markers, highlighting the relevance of ZEB1 in these two sub-populations.

RESULTS
Modeling phenotype switching towards ZEB1 high /MITF low NCSC/invasive state in vitro In order to study the role of endogenous ZEB1 during phenotypic transitions in melanoma cells, we used two BRAF V600 patientderived short-term cultures, established with a low number of passages after culture (GLO and C-09.10).These two short-term cultures display a ZEB1 low /MITF high 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 ZEB1 high /MITF low 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 ZEB1 high /MITF low state in GLO cells (Fig. 1A, B).ZEB1 protein expression increased, while ZEB2 expression decreased upon treatment (Fig. 1A).Drastic downregulation 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).
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 MITF high Fig. 1 Modeling phenotype switching towards ZEB1 high /MITF low 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 RNAseq 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 (nonsignificant) means P > 0.05.or a MITF low state.TNFα treatment decreased the proportion of the MITF high population and led to the emergence of a MITF low / NGFR high phenotype (representing about 70% of the population at day 7), before a transition towards a MITF low /NGFR low 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 MITF high population and a transient increase in MITF low /NGFR int/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 ZEB1 high 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 (ZEB1 low ) and TNFα-treated (ZEB1 high ) GLO cells at day 14, but also in the A375 melanoma cell line, which displays a ZEB1 high NCSClike expression pattern (MITF low , NGFR high , 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 ZEB1 low (untreated) and ZEB1 high (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).
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 ZEB1 low 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 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.

Hoek Proliferative Hoek Invasive
Tsoi Melanocytic Tsoi NCL 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 ZEB1 high melanoma models.

Tsoi Transitory Tsoi Undifferentiated
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 downregulation 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 ZEB1 high 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 ZEB1 high 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).
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 ZEB1 high 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 overexpression 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 overexpression 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.
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.melregulon 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.melregulon in single-cell RNAseq data.ZEB1.mel regulon displayed an increased activity in  ns

Relative mRNA expression
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.
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].
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 ZEB1 high tumors and evidenced cell populations with decreased SOX10 intensity.As illustrated in the MM10 tumor, which bore the presence of a well-defined ZEB1 high clone, increased ZEB1 expression was not only associated with low ZEB2 and MITF expression, but also with decreased SOX10 levels (Figs.7C and 8A).Quantitative analyses of SOX10 intensity according to ZEB1 expression level (high, intermediate, low and negative) confirmed a significant decrease in SOX10 levels when ZEB1 expression increases in melanoma cells (Fig. 8B).
We then investigated a putative gain in the expression of NCSC and mesenchymal markers, NGFR and SOX9, in tumors presenting a high ZEB1 expression.Although high SOX9 or NGFR expression was only detected in a limited number of tumors (4 and 5 out of 30; respectively), these tumors were enriched in ZEB1 expression (Fig. 8C).In most cases, NGFR expression was observed in scarce sub-populations of cells.Interestingly, the ZEB1 high clone from MM10 displayed NGFR positivity (Figs.7C and 8D), and NGFR levels were significantly higher in ZEB1 high melanoma cells (Fig. 8E).As regards to SOX9 intra-tumoral expression, antagonistic expression with SOX10 could clearly be evidenced at the singlecell level (Fig. 7A).Moreover, SOX9 levels were correlated with ZEB1 expression as exemplified in the ZEB1 high MM10 tumor (Figs.7D and 8F-G).Interestingly, only a few cells displayed triple positivity for ZEB1, NGFR and SOX9 (Fig. 8D).
Overall, single-cell and spatial characterization of melanoma intra-tumor heterogeneity demonstrate that increased ZEB1 expression is correlated with decreased MITF and SOX10 levels and increased NGFR or SOX9 expression, thus highlighting the importance of ZEB1 in vivo in both the NCSC and mesenchymal cell populations.

DISCUSSION
Our study reports a genome-wide characterization of the transcriptional functions of ZEB1 in melanoma, providing a better understanding of the molecular mechanisms underlying phenotype plasticity and intra-tumor heterogeneity in melanoma.We identified and validated the direct binding of the ZEB1 transcription factor to the promoter of genes specific to the melanocytic lineage or driving melanoma cell identity.Gainor loss-of-function of ZEB1, combined with function analyses, further demonstrates that ZEB1 negatively regulates proliferative-melanocytic programs and up-regulates invasive/ stem-like programs.
Overall, this study defines ZEB1 as a major regulator of melanoma cell identity and phenotype switching.Although initially described as a transcriptional repressor, our data confirm previous ChIP-seq analyses in carcinoma models [23] showing the capacity of ZEB1 to mediate both transcriptional activation and repression in similar proportions.Interestingly, at the basal level, ZEB1 mostly represses the expression of melanocytic genes, while increased ZEB1 expression upon phenotype switching, is associated with de novo binding, driving the up-regulation of invasive and NCSC genes.Although ZEB1 binding peaks are observed, ZEB1 may not increase NCSC and mesenchymal markers in the same cells, nor at the same time.This is consistent with models proposed in carcinoma, where ZEB1 may promote stemness features (partial EMT state associated with tumor initiation) but not necessarily invasive/EMT features, these two features being uncoupled.Overall, ZEB1 may not be associated with the 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.melmelanoma 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.acquisition of a given cell state but may regulate reversible cell state transitions in a dynamic manner.
ZEB1 not only represses MITF expression and subsequently the MITF-transcriptional program, but also directly regulates known MITF targets.Previous MITF ChIP-seq analyses demonstrated that MITF directly and positively regulates genes involved in DNA replication, repair and mitosis, while repressing genes involved in melanoma invasion [40].ZEB1 and MITF may thus bind to the same genes but with different consequences.AP-1 motif was found enriched at ZEB1 binding sites and the AP-1 subunit FOSL2 was shown to co-occupy similar loci with ZEB1, suggesting a cooperation of ZEB1 with AP-1 in melanoma, in line with data in breast cancer cells [23].Furthermore, the TF activity of both JUN and FOS was enriched upon ZEB1 activation, suggesting a positive feedback loop on AP-1 activity.
Importantly, even if some ZEB1 target genes are shared between melanoma and breast cancer cell lines [23], such as CDH1 and other EMT genes, our study reveals cell type-specific effects of ZEB1, through the regulation of melanocytic lineagespecific genes.ZEB1 notably binds to the promoter of ZEB2 and represses its expression in melanoma cells, while these two factors are co-expressed in mesenchymal cells.We further designed a melanoma-specific ZEB1 regulon which can be used by the scientific community to accurately study ZEB1 transcription factor activity in a melanoma context.We validated the increased activity of the ZEB1-melanoma-specific regulon in single-cell-RNA-seq data in both Neural-Crest-like and Mesenchymal populations.The precise characterization of ZEB1 cofactors, as well as its relationships with other key TFs regulating cell states, such as the recently described PRRX1 [41], would be required in order to better comprehend their interplay and hierarchy during gene regulation in a context and cell type-specific manner.
Importantly, our quantitative spatial analyses at the single-cell resolution of markers of melanoma cell states in human samples provided further validation of the co-expression of ZEB1 with NCSC (NGFR) or undifferentiated (SOX9) markers and its inverse correlation with the melanocytic markers MITF and SOX10.Although NGFR and SOX9 are not sufficient to define the NCSC and undifferentiated states respectively, ZEB1 expression may be found in these two subpopulations of cells.Importantly, ZEB1 is not only expressed in the invasive front but also in the bulk, where it may sustain stem-like features and tumor-initiating properties.SOX10 intra-tumor heterogeneity in melanoma samples was consistent with recent reports [42].
TNFα and TGFβ mimic only part of the signals emanating from the tumor microenvironment.Indeed, we recently highlighted the major crosstalk existing between melanoma cells and their immune microenvironment [39,43] and recent work from the Marine lab highlighted the major role of endothelial cells in promoting a mesenchymal state [41].ZEB1 and other melanoma markers of intra-tumor heterogeneity may thus significantly be modified by the tumor microenvironment in specific niches that will deserve further characterization.Moreover, if we focused on markers of melanocyte cell identity, this ChIP-seq approach also revealed additional targets related to inflammatory or interferon responses, consistent with the pleiotropic roles of ZEB1, that extend beyond invasion, including immune escape [39].
Hence, this study provides important insights into the way ZEB1 orchestrates gene expression, with a precise combination of both down-regulation and up-regulation of key genes of melanoma cell state, which in turn mediate reversible phenotypic plasticity, known to foster the acquisition of resistance to treatment in melanoma.Although targeting ZEB1 remains challenging, this work highlights new candidates/pathways that could represent interesting targets to dampen melanoma cell plasticity as a strategy to overcome treatment resistance.

MATERIALS AND METHODS Human tumor samples
Melanoma tumor samples were obtained through the Biological Resource Center of the Lyon Sud Hospital (Hospices Civils de Lyon) and were used with the patient's written informed consent.This study was approved by a regional review board (Comitéde Protection des Personnes Ile de France XI, Saint-Germain-en-Laye, France, number 12027) and is registered in ClinicalTrial.gov (MelBase, NCT02828202).30 cutaneous melanoma patients were used for multi-immunofluorescence analyses.All melanoma biopsies were cutaneous, either primary melanoma or cutaneous metastases.

Cell culture and treatments
The A375 human melanoma cell line was purchased from ATCC and cultured in DMEM complemented with 10% fetal bovine serum (FBS) (Cambrex) and 100 U/ml penicillin-streptomycin (Invitrogen).In order to authenticate the cell lines, the expected major genetic alterations were verified by NGS sequencing.The absence of Mycoplasma contamination was verified every 3 weeks with the MycoAlert detection kit (Lonza).Previously described patient-derived short-term cultures ( < 10), GLO and C-09.10, established from BRAF V600 metastatic melanomas [21], were grown in RPMI complemented with 10% FBS and 100 U/ml penicillinstreptomycin.TNFα (100 ng/mL) and TGFβ (20 ng/mL) (Peprotech) were replaced in the culture medium every 3 days.The BRAF inhibitor PLX4032/ vemurafenib was purchased from Selleck Chemicals (Houston, TX, USA) and reconstituted in DMSO.

Immunoblot analyses
Cells were washed twice with phosphate buffered saline (PBS) containing CaCl 2 and then lysed in a 100 mM NaCl, 1% NP40, 0.1% SDS, 50 mM Tris pH   8.0 RIPA buffer supplemented with a complete protease inhibitor cocktail (Roche, Mannheim, Germany) and phosphatase inhibitors (Sigma-Aldrich).Loading was controlled using anti-GAPDH.Horseradish peroxidaseconjugated goat anti-rabbit polyclonal antibodies (Glostrup) were used as secondary antibodies.Western blot detections were conducted using the Luminol reagent (Santa Cruz).Western Blot Digital Imaging was performed with the ChemiDoc™ MP Imager (Bio-Rad).The list of antibodies used is detailed (Supplementary Table 1).

RT-Q-PCR
Total RNA was isolated using the RNeasy Kit (QIAGEN) and reversetranscribed using a high cDNA capacity reverse transcription kit (Maxima First Strand cDNA synthesis Kit, Thermoscientific) following the manufacturer's instructions using 1000 ng of RNA as a reverse transcription template in a 20 μL final volume.The samples were incubated for 10 min at 25 °C, followed by 15 min at 50 °C and 5 min at 85 °C in T100 Thermal Cycler (1861096, Bio-Rad).Real-time qPCR reactions were performed using SsoAdvanced Universal SYBR® Green Supermix (1725274, Bio-Rad) according to the manufacturers protocol.Reaction were done using 15 ng of cDNA template and 1 µM of each primer.All reactions, including notemplate controls and RT controls were performed in triplicate on a CFX96 (Bio-Rad) with 40 cycles at 95 °C for 5 s followed by 10 s at 60 °C.Results were analyzed with the Bio-Rad CFX manager software.Human GAPDH was used for normalization.Exon-spanning probes were designed using the ProbeFinder software (Roche).The list of primers used is detailed (Supplementary Table 2).

RNA-seq analyses
RNA libraries were prepared with the TrueSeq poly-A+ kit from Illumina and sequenced on the genomic platform of the CRCL, on an Illumina NovaSeq 6000 sequencing machine with a paired-end protocol (2 x 75 bp, 32 Mp reads).Raw sequencing reads were aligned on the human genome (GRCh38) with STAR (v2.7.8a), with the annotation of known genes from gencode v37.Gene expression was quantified using Salmon (1.4.0) and the annotation of protein-coding genes from gencode v37.
Bulk RNA-seq data of melanoma cell lines from Tsoi et al. [29].was retrieved from GEO, with accession number GSE80829.Single-cell RNAseq data of patient derived xenograft from Rambow et al. [15].were obtained from GEO with accession number GSE116237.Single-cell RNAseq data of melanoma cell lines from Wouters et al. [37].were retrieved from GEO (GSE134432) and of Pozniak et al. [38] from the KU Leuven Research Data Repository.Single-cell RNAseq data were analysed and visualized using Seurat (4.3.0) and SCpubr (1.1.2) packages.

Chromatin immunoprecipitation and ChIP-sequencing analyses
The ChIP assay was carried out according to the protocol from the iDeal ChIP-Seq Kit for Transcription Factors (Diagenode).Briefly, cells from one 15-cm dish were cross-linked with 1% formaldehyde at RT for 10 min and quenched in 125 mM glycine for 5 min.Cross-linked chromatin was isolated and sonicated to generate DNA fragments averaging 200-500 bp in length by Bioruptor plus sonication device (Diagenode).Chromatin fragments were immunoprecipitated with antibodies directed against ZEB1 (1 μg, Genetex, GTX105278, RRID:AB_11162905), or IgG (1 µg, Bio-Rad, PRABP01, RRID:AB_321631) as negative control.Immunoprecipitated DNA was purified and dissolved into 25 µl of H 2 O.To build the Illumina library, 5 ng of input and 30 pg of IP were used.Sequencing was performed on the CLB genomic platform, on an Illumina NextSeq machine with a paired-end protocol (2 x 75 bp, 64 Mp reads).To assess the efficacy of the ChIP before sequencing, MITF (positive control) was analyzed by qPCR.Primers were specified to amplify genomic DNA using the sequence from the peak (ChIPSeq data from MDAMB231, A375 or GLO cells) (Supplementary Table 3).Relative promoter enrichment was normalized against chromatin inputs.
DNA libraries were prepared with the Diagenode MicroPlex Library Preparation Kit v2, and sequenced on a NextSeq sequencing machine (paired-end protocol, 75 bp, 80 M reads) on the genomic platform of the Centre Léon Bérard (CLB).

Flow cytometry
To analyze the expression of the NGFR/CD271 cell surface marker, cells 1 × 10 6 cells per condition were incubated with an AlexaFluor647-conjugated anti-CD271 antibody (BD Pharmingen, RRID: AB_1645403) for 1 h in the dark before being counted on a BD LSRFortessa™ Flow Cytometer (BD Biosciences-IN).Data were analyzed using the FlowJo_V10 software.

Transwell migration assay
Falcon® Cell Culture Inserts for 24-well plates were placed into a 24-well plate, and 20% FBS RPMI was added into the well and 0.5% FBS RPMI was added into inserts.The plate was then incubated for 1 h for an initial equilibrium period.200,000 GLO or 250,000 A375 cells were trypsinated and rinsed with 0.5% FBS RPMI medium.100 µl 0.5% FBS RPMI containing 200,000 GLO cells was seeded into the insert, then 600 µl 20% FBS RPMI was added to the well.After 24 h incubation, cells inside the insert were removed carefully and the migrated cells on the membrane were fixed and colored using 4% PFA and Brilliant Blue, then rinsed with PBS.When inserts were completely dry, cells were viewed using phase-contrast microscopy.Photos were taken and analyzes by ZEISSZEN Microscope software.Migrated cells were counted in 3 different fields.
Incucyte® Live Cell mortality measurement 30,000 cells were seeded onto a 24-well plate and treated with TNFα ± TGFβ.After 24 h, cell medium was renewed with indicated treatments, as well as PLX4032 and propidium iodide (Sigma-Aldrich®, 1/3000).The EssenBioScience IncuCyte Zoom Live-Cell Analysis System was used to measure and analyze real-time cell mortality every 2 h.Dead cells were marked with propidium iodide.Data were then converted into Excel in order to draw graphs.

7-color immunofluorescence multiplex analyses
3-µm tissue sections were cut from formalin-fixed paraffin-embedded human melanoma specimens.The sections underwent immunofluorescence staining using the OPAL™ technology (Akoya Biosciences) on a Leica Bond RX.A 7-color panel was designed (Supplementary Table 4).DAPI was used for nuclei detection.Sections were digitized with a Vectra Polaris scanner (Perkin Elmer, USA).Using the Inform software (Perkin Elmer), an autofluorescence treatment of images was carried out and tissue segmentation was performed to identify epidermis, stroma and tumor.Cell segmentation was then applied to analyze the expression of each marker in each cell.The matrix of phenotype containing the X-and Ypositions of each cell as well as the mean nuclear-, cytoplasmic-and membrane-intensities of each fluorescence staining was then further analyzed using the R software.Tumors were spatially reconstructed using the R plot() function.
For quantitative analyses, melanoma cells were classified following the expression level of each marker, different cut-off values for each marker were defined (Supplementary Table 5).

Statistical analyses
To ensure adequate power and decreased estimation error, we performed multiple independent repeats and experiments were conducted at least in triplicate.Data are presented as mean ± s.d. or ± s.e.m as specified in the figure legends.Statistical analyses were performed using GraphPad Prism 8 software (GraphPad Software, Inc., San Diego, USA) or R software (v4.1.0)and plots were generated with ggplot2 (v3.3.5).All statistical tests were two-tailed and p-values were corrected, when indicated, with the Benjamini-Hochberg method.Paired student's t tests were used to compare the means of two groups.To determine significant differences between two groups, student's t tests or Mann Whitney tests were used as indicated in the figure legends.

Fig. 3
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.

Fig. 4
Fig.4Validation 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.

Fig. 6
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 nontargeting control (NTC).F Violin plot of the transcription factor activity of ZEB1.mel regulon in Pozniak et al. single cell RNA-seq data.

Fig. 7
Fig.7Single-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 ZEB1 low (top) and a ZEB1 high (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 ZEB1 low B, and two ZEB1 high 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.

Fig. 8
Fig. 8 ZEB1 antagonistic expression with SOX10 and co-expression with NGFR and SOX9 within melanoma lesions.A Representative pictures of ZEB1 and SOX10 staining showing antagonistic expression of ZEB1 and SOX10 in the ZEB1 high clone from Patient MM10.B Violin plots showing the expression levels of SOX10 for cells grouped with respect to their ZEB1 expression status (high, intermediate, low or not expressed) in two representative ZEB1 high melanoma cases.C ZEB1, NGFR and SOX9 status annotated in high, intermediate and low expression based on the protein expression level in the multi-IF analysis.D Representative pictures of ZEB1, NGFR and/or SOX9 staining showing co-expression of ZEB1 with NGFR and/or SOX9.E Violin plots showing the expression levels of NGFR for cells grouped with respect to their ZEB1 expression status (high, intermediate, low or not expressed) in two representative ZEB1 high melanoma cases.F Representative pictures of ZEB1 and SOX9 staining showing co-expression of ZEB1 with SOX9.G Violin plots showing the expression levels of SOX9 for cells grouped with respect to their ZEB1 expression status (high, intermediate, low or not expressed) in one representative ZEB1 high melanoma.The median is shown with a red dot.P values were determined by Mann-Whitney test with Bonferroni correction B, E, G. Differences were considered statistically significant at *P ≤ 0.05, **P < 0.01 and ***P < 0.001.ns (non-significant) means P > 0.05.