1. Cancer Biology

Restraint of melanoma progression by cells in the local skin environment

  1. Yilun Ma
  2. Mohita Tagore
  3. Miranda V Hunter
  4. Ting-Hsiang Huang
  5. Emily Montal
  6. Joshua M Weiss
  7. Yingxiao Shi
  8. Tuulia Vallius
  9. Peter K Sorger
  10. Richard M White  Is a corresponding author
  1. Weill Cornell/Rockefeller/Sloan Kettering Tri-Institutional MD-PhD Program, United States
  2. Department of Cancer Biology and Genetics, Memorial Sloan Kettering Cancer Center, United States
  3. Cell and Developmental Biology Program, Weill Cornell Graduate School of Medical Sciences, United States
  4. Laboratory of Systems Pharmacology, Harvard Program in Therapeutic Science, Harvard Medical School, United States
  5. Department of Clinical Computational Oncology, Dana-Farber Cancer Institute, United States
  6. Ludwig Center, Harvard Medical School, United States
  7. Department of Systems Biology, Harvard Medical School, United States
  8. Nuffield Department of Medicine, Ludwig Institute for Cancer Research, University of Oxford, United Kingdom
https://doi.org/10.7554/eLife.101974.3
Share this article
Cite this article
  1. Yilun Ma
  2. Mohita Tagore
  3. Miranda V Hunter
  4. Ting-Hsiang Huang
  5. Emily Montal
  6. Joshua M Weiss
  7. Yingxiao Shi
  8. Tuulia Vallius
  9. Peter K Sorger
  10. Richard M White
(2026)
Restraint of melanoma progression by cells in the local skin environment
eLife 13:RP101974.
https://doi.org/10.7554/eLife.101974.3
  2. Download BibTeX
  3. Download .RIS

eLife Assessment

In this important study, the authors used a zebrafish model and scRNAseq analysis to show that a subset of keratinocytes within melanoma microenvironment highly up-regulate Twist and undergo Epithelial-Mesenchymal Transition (EMT). Surprisingly, when overexpressing Twist in keratinocytes, the resulting alteration in keratinocytes is inhibitory for melanoma invasion in both zebrafish and human cell culture models. The results are supported by convincing experimental data that provide new insights into the interactions between melanoma cells and their environment.

https://doi.org/10.7554/eLife.101974.3.sa0

Abstract

Keratinocytes, the dominant cell type in the melanoma microenvironment during tumor initiation, exhibit diverse effects on melanoma progression. Using a zebrafish model of melanoma and human cell co-cultures, we observed that keratinocytes undergo an epithelial-mesenchymal transition (EMT)-like transformation in the presence of melanoma, reminiscent of their behavior during wound healing. Surprisingly, overexpression of the EMT-transcription factor Twist in keratinocytes led to improved overall survival in zebrafish melanoma models, despite no change in tumor initiation rates. This survival benefit was attributed to reduced melanoma invasion, as confirmed by human cell co-culture assays. Single-cell RNA-sequencing revealed a unique melanoma cell cluster in the Twist-overexpressing condition, exhibiting a more differentiated, less invasive phenotype. Further analysis nominated homotypic jam3b–jam3b and pgrn–sort1a interactions between Twist-overexpressing keratinocytes and melanoma cells as potential mediators of the invasive restraint. Our findings suggest that EMT in the tumor microenvironment may paradoxically limit melanoma invasion through altered cell–cell interactions.

Introduction

The complex interplay between cancer cells and their microenvironment has emerged as a critical determinant of tumor progression and therapeutic response. In melanoma, the tumor microenvironment (TME) encompasses diverse cell types, including immune cells, fibroblasts, and endothelial cells (Lee and Herlyn, 2007). However, during melanoma initiation, the dominant cell type in the TME is the keratinocyte, an epithelial cell which makes up the majority of our skin surface. In normal homeostasis, each melanocyte reciprocally interacts with 30–40 keratinocytes (Fitzpatrick and Breathnach, 1963), and this interaction is essential for skin and hair color (Kunisada et al., 1998; Tanimura et al., 2011). Despite decades of research, our understanding of keratinocytes in the context of melanoma remains incomplete.

Keratinocytes have been shown to both promote and inhibit tumor initiation. They are tightly adherent to melanocytes, and this can inhibit tumor development because nascent melanoma cells cannot escape the epidermis (Fukunaga-Kalabis et al., 2008). Previous work has demonstrated that newly initiated melanoma cells must escape from this adhesive network by loss of proteins such as PAR3 (Mescher et al., 2017). In contrast, keratinocytes can also promote tumor development through secretion of growth factors such as endothelins or via GABAergic crosstalk between the two cell types (Jamal and Schneider, 2002; Tagore et al., 2023).

These conflicting data highlight that interactions between keratinocytes and nascent melanoma cells are likely dynamic and change rapidly during tumor initiation. Studying the nature of these interactions in human samples is challenging because biopsies are taken after the patient has come to the clinic, meaning that the earliest interactions in tumor initiation will be missed. This necessitates models which faithfully recapitulate the earliest stages of tumor initiation, yet have the cellular resolution to measure interactions between melanoma cells and keratinocytes.

In this study, we utilized a zebrafish model of melanoma to investigate the earliest interactions between melanoma cells and their neighboring keratinocytes (Callahan et al., 2018). Zebrafish have emerged as a powerful tool for cancer research due to their genetic tractability, conserved biology, and the ability to visualize tumor development and progression in real-time within the context of an intact organism (Kaufman et al., 2016; Li and Uitto, 2014). Using a combination of cell-type-specific genetic manipulations, in vivo imaging, and single-cell transcriptomics, we found that tumor-associated keratinocytes (TAKs) undergo changes associated with EMT, similar to what is found in wounded skin. Unexpectedly, we found that this keratinocyte EMT suppresses melanoma progression. This change in the keratinocytes occurs shortly after melanoma initiation and results in keratinocytes which are more adhesive to these nascent tumor cells and prevents their movement out of the epidermis. Our data suggests that melanoma initiation revises an evolutionarily conserved wounding response in the nearby skin environment, which acts as a cell-extrinsic tumor suppressor to prevent newly transformed cells from becoming fully tumorigenic.

Results

Melanoma initiation is associated with EMT in keratinocytes

To investigate the relationship between keratinocytes and melanoma cells in vivo, we created a transgenic zebrafish line in which GFP is expressed under the krt4 promoter (Gong et al., 2002). This line faithfully marks all adult keratinocytes present throughout the fish epidermis and scales, similar to previous lines using this promoter (Figure 1A). We then initiated melanomas in this background using the TEAZ method (Transgene Electroporation of Adult Zebrafish) (Callahan et al., 2018; Montal et al., 2024), in which plasmids containing oncogenes or sgRNAs against tumor suppressors can be introduced directly into the skin (Figure 1B). The major advantage of this method is that we can visualize melanoma initiation when the tumor is in its early stages, consisting of a small number of cells. We initiated tumors with a combination of BRAFV600E, sgRNAs against PTEN (ptena/b), and germline loss of p53 (Figure 1C). To account for skin wounding from electroporation in assessing changes in TAKs, we also performed TEAZ using a control vector that labels melanocyte-precursors but does not induce melanoma formation. Fluorescent imaging 8 weeks post-electroporation with the control vector demonstrated an injury-free epidermis, in contrast to the pronounced melanoma development in zebrafish administered with the oncogenic vectors (Figure 1D).

Figure 1
Download asset Open asset
Keratinocytes in the melanoma microenvironment undergo EMT-like changes.

(A) Generation of transparent zebrafish with GFP-labeling of keratinocytes. Casper Triples (mitfa−/−;mitfa:BRAFV600E;tp53−/−;mpv17−/−) were injected with the Tol2Kit 394 vector containing a krt4:eGFP cassette and tol2 mRNA. Brightfield shows a transparent zebrafish while fluorescence imaging shows eGFP-labeling of keratinocytes (scale bar = 5 mm). Panel A was created with BioRender.com. (B) Schematic of TEAZ (Transgene Electroporation of Adult Zebrafish). Plasmid mix containing miniCoopR:tdTomato, mitfa:Cas9, zU6:sgptena, zU6:sgptenb, and the tol2 plasmid was injected superficially in the flank of the zebrafish. Electroporation of the injection site results in rescue of melanocyte precursors and the generation of a localized melanoma that could be analyzed by microscopy and FACS. Panel B was created with BioRender.com. (C) Brightfield and immunofluorescence of zebrafish 8 weeks post-TEAZ with localized and fluorescently labeled melanoma (scale bar = 5 mm). (D) Immunofluorescence imaging of TEAZ region after 8 weeks, comparing empty vector control vs. miniCoopR:tdT conditions, with yellow dotted circles indicating general area of dissection for FACS (scale bar = 1 mm). (E) Confocal imaging of zebrafish epidermis. Normal epidermis of Tg(krt4:eGFP) Casper Triple post-TEAZ with empty vector control shows eGFP-labeled, polygonal shaped keratinocytes regularly connected while epidermis with melanoma generated with miniCoopR-driven melanocyte rescue shows disrupted epidermis and irregularly shaped keratinocytes. Inset highlights melanoma cells with adjacent keratinocytes (scale bar = 50 µm). (F) qPCR of FACS sorted zebrafish epidermis with or without melanoma (n = 3 per condition). tdTomato-labeled melanoma cells, eGFP-labeled keratinocytes, and non-fluorescently labeled tumor microenvironment (TME) cells were isolated by dissection (as indicated in G) and FACS. Comparison of mitfa and krt4 expression of samples normalized to non-fluorescent cells, either TME-Other in tumor samples or Other in non-tumor samples, shows enrichment of mitfa in melanoma samples and krt4 in keratinocyte samples. Mean ± SE. ns = no significance, *p < 0.05, ***p < 0.001, ****p < 0.0001 by Tukey’s multiple comparisons test. (G) Comparison of the EMT markers vim (vimentin) and cdh2 (N-cadherin) shows enrichment in TME keratinocytes vs. keratinocytes from epidermis without melanoma. Mean ± SE. *p < 0.05 by Welch’s t-test. (H) Schematic of keratinocyte–melanoma co-culture experiment. HaCaTs were cultured in monoculture or co-culture with A375 melanoma cells in triplicates for 21 days, followed by FACS isolation of keratinocytes for RNA-sequencing comparing co-culture vs. monoculture keratinocytes (n = 3 per condition). Panel H was created with BioRender.com. (I) Top 5 enriched Hallmark pathways in HaCaTs co-cultured with A375 melanoma cells compared with HaCaTs in monoculture. (J) Normalized counts of EMT biomarkers vimentin (VIM) and N-cadherin (CDH2). Mean ± SE.

To better investigate the changes in the keratinocytes, we used confocal microscopy on fish 8 weeks after TEAZ. This revealed a marked disruption of keratinocyte morphology in the tumor-bearing fish, which was not seen in control fish. We specifically noted disrupted cell–cell junctions, a disorganized pattern of keratinocytes, and loss of the normal hexagonal cell layer (Figure 1E). These changes were reminiscent of keratinocyte EMT, which has been previously noted to occur in wounded epidermis (Haensel and Dai, 2018; Leopold et al., 2012; Moreno-Bueno et al., 2009). To further assess this possibility, we excised tissues from both tumor and control skin and used FACS to isolate keratinocytes (GFP+) and melanoma cells (tdTomato+) and performed qPCR. As expected, we found enrichment of mitfa in melanoma cells and krt4 in keratinocytes, validating the successful cell-type isolation (Figure 1F). Comparative analyses between TAKs and normal keratinocytes (NKCs) from tissue without melanoma revealed upregulation of EMT markers vimentin and N-cadherin, consistent with our imaging results (Figure 1G).

We next asked whether these changes were also seen in human samples. To address this, we performed co-culture experiments between keratinocytes and melanoma cells. We grew GFP-labeled HaCaT keratinocytes either alone or with A375 melanoma cells for 21 days, followed by isolation by FACS for bulk RNA-sequencing (Figure 1H; Tagore et al., 2023). Consistent with our in vivo results in the fish, the top pathway altered in the co-cultured HaCaT cells was enrichment of EMT (Figure 1I). Differential gene expression analysis showed notable upregulation of the mesenchymal markers vimentin and N-cadherin in co-cultured keratinocytes compared to monocultured control (Figure 1J), similar to what was found in vivo. Collectively, our data indicate that melanoma cells induce morphological and molecular markers of EMT in nearby keratinocytes.

EMT-transcription factors are upregulated in TAKs

EMT is usually driven by upstream transcription factors, which then act on downstream targets to repress adhesion molecules such as E-cadherin or activate other adhesion molecules such as N-cadherin (Moreno-Bueno et al., 2009). We next wanted to understand which of these transcription factors was responsible for the EMT-like behavior in our keratinocytes. To address this, we utilized an existing scRNA-sequencing dataset of a BRAFV600E-driven zebrafish melanoma (Figure 2A). Dimensionality reduction with UMAP and subsequent clustering revealed two keratinocyte populations as indicated by module scoring for genes enriched in keratinocyte populations (Figure 2B, C). Subsequent differential gene expression and Gene Set Enrichment Analysis (GSEA) of the two keratinocyte clusters revealed one cluster with enrichment for EMT, similar to what we observed above (Figure 2D, E). We refer to this EMT cluster as TAKs, and the other cluster as NKCs. We focused on three EMT-transcription factors expressed in zebrafish keratinocytes, snai1a, snai2, twist1a, zebrafish homologs of human SNAIL, SLUG, and TWIST. Differential expression showed significant enrichment in snai1a and twist1a in TAK vs. NKC clusters (Figure 2F). The significant enrichment of twist1a in the TAK cluster, coupled with its rare expression in NKCs, positioned twist1a as a promising candidate for further investigation into its potential role in driving EMT-like changes in keratinocytes and, consequently, its impact on melanoma progression.

Figure 2
Download asset Open asset
Zebrafish scRNA-sequencing shows upregulation of EMT-TFs in tumor-associated keratinocytes.

(A) Schematic of scRNA-sequencing experiment. Embryo injection with miniCoopR:eGFP and tol2 mRNA in Casper Triple results in melanocyte rescue and subsequent melanoma formation. Melanoma was dissected and dissociated to single-cell suspension for FACS isolation of eGFP+ melanoma cells and non-fluorescent tumor microenvironment (TME) cells for single-cell RNA-sequencing. Panel A was created with BioRender.com. (B) UMAP highlighting two keratinocyte clusters, melanoma, and other TME cells. (C) Violin plot of keratinocyte module scores comparing between keratinocyte clusters, TAK and NKC, vs. melanoma and other TME cells. (D) Top 6 GSEA Hallmark analysis comparing TAK vs. NKC (NES = Normalized Enrichment Score). (E) Hallmark EMT pathway enrichment in TAK vs NKC. (F) Enrichment of EMT-transcription factors in TAK vs. NKC, melanoma, and other TME cells.

Keratinocyte TWIST restrains melanoma invasion in zebrafish

Having identified Twist as a potential driver of EMT-like changes in TAKs, we next asked how it affected melanoma phenotypes. We created new transgenic zebrafish in which the krt4 promoter was used to drive the two zebrafish TWIST paralogs (twist1a and twist1b) in the context of BRAF-driven melanomas. In this experiment, we used standard 1-cell injection of the plasmids rather than TEAZ to ensure that the results were more generalizable to different melanoma initiation assays. Injected embryos were sorted at 5 days post-fertilization for tdTomato (melanocyte) and GFP (keratinocyte) expression, indicating successful melanocyte and keratinocyte transformation (Figure 3A). These were mosaic animals, rather than stable lines, so we could directly assess the effect of TWIST on melanoma growth (although this precluded us from examining the effect of TWIST in keratinocytes without a nearby melanoma). We monitored the fish for tumor-free survival as well as overall survival over the ensuing 26 weeks. Interestingly, we found no difference in melanoma initiation rate in the TWIST overexpression condition compared to empty vector (Figure 3B). Unexpectedly, however, we noted that overall survival was improved in the transgenic animals expressing TWIST in the keratinocytes (Figure 3B).

Figure 3 with 1 supplement see all
Download asset Open asset
Overexpression of twist1a/b results in improved survival of fish with melanoma.

(A) Schematic of zebrafish melanoma model with labeling and perturbation of keratinocytes. Twist1a and twist1b are overexpressed under the keratinocyte-specific, krt4, promoter in the TWIST condition, and an empty vector control was used in the CTRL condition. Fish were sorted at 5 days for eGFP and tdTomato positivity as markers of successful keratinocyte labeling and melanocyte rescue (n = 135 CTRL, n = 118 TWIST). Panel A was created with BioRender.com. (B) Tumor-free survival and overall survival of A. ns = no significance, **p < 0.005 by log-rank (Mantel–Cox) test. (C) Sample images of zebrafish with melanoma at 26 weeks post-injection. Keratinocytes are labeled by GFP, melanomas are labeled by tdTomato. Melanomas are pigmented. Scale bar = 5 mm. (D) H&E and IHC of cross-sections through zebrafish body and melanoma. Dotted yellow line demarcates the border of the body. % of tumor in body is calculated as tumor area within body border divided by total tumor area. Scale bar = 500 µm.

This discrepancy between melanoma-free and overall survival suggested that the tumors in the TWIST condition should be phenotypically distinct from the control tumors. To assess this, we performed immunohistochemistry on the tumors (n = 3 in each condition) and surrounding tissues from the control and TWIST conditions (Figure 3C). The oncogenic driver in this melanoma model is hBRAFV600E and serves as an IHC marker for the tumor cells. Comparison of hBRAFV600E staining revealed significant melanoma infiltration into the zebrafish body in the CTRL condition, as opposed to a nearly non-invasive tumor in the TWIST condition (Figure 3D). The lack of melanoma invasion was also observed when the melanoma developed in other anatomical locations as well (Figure 3—figure supplement 1), although because survival started to decrease around 26 weeks, we could not let these fish grow longer to see if they eventually became invasive at very late time points. These findings suggest that twist1a/twist1b overexpression in keratinocytes does not impair tumor initiation, but instead impairs melanoma invasion and improves survival when expressed in the microenvironment.

Because this was a somewhat unexpected finding, we wanted to ensure that this was not an artifact of the zebrafish system. We therefore developed an assay to measure interactions between human keratinocytes and melanoma cells. To simulate melanoma invasion into neighboring cells, we developed a novel cell culture-based invasion assay using HaCaT lines and multiple melanoma cell lines. Melanoma cells cultured on poly-l-lysine coated coverslips were placed atop keratinocytes, allowing assessment of melanoma cell infiltration after 20 hr (Figure 4A). Due to the migratory nature of melanoma cells, they will migrate off the coverslip and infiltrate the layer of keratinocytes, allowing us to assess relative differences between culture conditions.

Figure 4 with 1 supplement see all
Download asset Open asset
Zebrafish findings are recapitulated in human cell lines.

(A) Schematic of coverslip cell infiltration assay. Melanoma cells are plated on a coverslip and allowed to attach overnight. The coverslip is then transferred into a well of keratinocytes to assess melanoma infiltration into keratinocytes. Panel A was created with BioRender.com. (B) Generation of a HaCaT cell line overexpressing TWIST1. HaCaT-BFP was infected with lentivirus containing a cassette with nuclear localized GFP and CMV-driven TWIST1 or no ORF. Infected cell lines were allowed to grow for a week before sorting for nuclear GFP as a marker of successful integration. Panel B was created with BioRender.com. (C) Western blot for Twist expression in HaCaT-CTRL and HaCaT-TWIST. (D) Immunofluorescence imaging for Twist localization in HaCaT-CTRL and HaCaT-TWIST. TWIST staining is pseudo-colored in magenta, DAPI in blue, phalloidin in white. Scale bar = 50 µm. (E) Immunofluorescence imaging of coverslip cell infiltration assay with HS294T-tdT (orange) melanoma cells in co-culture with either HaCaT-CTRL or HaCaT-TWIST (cyan). Insets highlight areas of melanoma cell infiltration into keratinocyte monolayers. Scale bar = 500 μm. (F) Quantification of E. Infiltrating HS294T melanoma cells from each image were counted and averaged across four images per well. Resulting cell counts were normalized to average cell counts of HaCaT-CTRL from each set. N = 9, 3 sets, 3 replicates/wells per set. Mean ± SE. *p < 0.05 by t-test. (G) Immunofluorescence imaging of coverslip cell infiltration assay with SKMEL2-tdT (orange) melanoma cells in co-culture with either HaCaT-CTRL or HaCaT-TWIST (cyan). Insets highlight areas of melanoma cell infiltration into keratinocyte monolayers. Scale bar = 500 μm. (H) Quantification of G. Infiltrating SKMEL2 melanoma cells from each image were counted and averaged across four images per well. Resulting cell counts were normalized to average cell counts of HaCaT-CTRL from each set. N = 9, 3 sets, 3 replicates/wells per set. Mean ± SE. *p < 0.05 by t-test.

Figure 4—source data 1

Labeled scans of western blot films for data presented in Figure 4C.

https://cdn.elifesciences.org/articles/101974/elife-101974-fig4-data1-v1.zip
Figure 4—source data 2

Original scans of western blot films for data presented in Figure 4C.

https://cdn.elifesciences.org/articles/101974/elife-101974-fig4-data2-v1.zip

HaCaT keratinocytes were transformed via lentiviral infection to overexpress TWIST1 (HaCaT-TWIST) or an empty vector control (HaCaT-CTRL) (Figure 4B). Western blot analysis confirmed robust Twist protein overexpression in HaCaT-TWIST compared to HaCaT-CTRL (Figure 4C), and immunofluorescence imaging revealed nuclear localization of Twist (Figure 4D). Co-culture of HS294T melanoma cells with HaCaT-TWIST resulted in significantly reduced melanoma cell invasion into keratinocytes compared to HaCaT-CTRL (Figure 4E, F), similar to what we had observed above in our zebrafish system. This finding was recapitulated using SKMEL2, demonstrating the inhibitory effect of TWIST1 overexpression in keratinocytes on melanoma invasion across different cell lines (Figure 4G, H). However, it was also possible that it was something lacking in the TWIST keratinocytes that led to less melanoma invasion (rather than actively suppressing invasion). To test this, we also performed a control assay in which we allowed tdTomato melanoma cells to invade into unlabeled melanomas of the same line. Using both HS294T and SKMEL2 cells, these cells do have some migratory capacity in this setting (Figure 4—figure supplement 1). This assay cannot discern whether the TWIST keratinocytes actively suppress invasion vs. lack of something that promotes invasion, which will be an area for future studies. Collectively, our results demonstrate that induction of EMT in keratinocytes is associated with reduced melanoma invasion and improvement in animal survival.

Keratinocyte EMT promotes aberrant adhesion to nascent melanoma cells

While EMT within the tumor cell is well recognized to promote invasion, our data suggest that EMT in the microenvironment paradoxically restrains tumor invasion. We wanted to better understand the downstream mechanisms accounting for this result. We therefore analyzed our control vs. TWIST tumors using single-cell RNA-sequencing, which would allow us to understand potential mechanisms by which melanoma cells were interacting with these keratinocytes. Melanoma and adjacent skin from three fish per condition (CTRL or TWIST) were dissociated into single-cell suspensions and FACS-sorted for eGFP+ keratinocytes and tdTomato+ melanoma cells. To enrich the keratinocyte population, keratinocyte and melanoma cell suspensions were combined at a 7:3 ratio for each condition and prepared for scRNA-sequencing (Figure 5A). The resulting CTRL and TWIST datasets were filtered for quality control as described (Figure 5—figure supplement 1A). UMAP dimensional reduction of the integrated dataset shows distinct clusters of melanoma and keratinocytes, with high eGFP expression in keratinocyte clusters and high tdTomato expression in melanoma clusters (Figure 5B). Furthermore, we identified two KC clusters in the present in both the CTRL and TWIST conditions and a third cluster of KCs in the TWIST condition. We performed differential gene expression analysis between clusters and GSEA on differentially expressed genes between the two clusters present in both conditions showed enrichment in EMT as seen in Figure 2, identifying EMT-enriched cluster as TAK and other as NKC (Figure 5—figure supplement 1B). As expected, we also identified a unique cluster of KCs only present in the TWIST dataset that highly expresses twist1a/b (Figure 5—figure supplement 1C), the genes that we overexpressed in the TWIST condition, and labeled this cluster Twist-High (Figure 5C).

Figure 5 with 1 supplement see all
Download asset Open asset
scRNA-sequencing shows unique keratinocyte–melanoma communication with twist1a/b overexpression in keratinocytes.

(A) Schematic of scRNA-sequencing protocol. Melanoma and surrounding tissue were dissected from 26 weeks old zebrafish from either CTRL or TWIST conditions as shown in Figure 3. Samples were dissociated to single-cell suspensions for FACS isolation of keratinocytes (GFP) and melanoma (tdTomato). Keratinocytes and melanoma were recombined per condition at a ratio of 7:3 for enrichment of keratinocytes for scRNA-sequencing. Panel A was created with BioRender.com. (B) UMAP dimensional reduction and feature plots of scRNA-sequencing dataset. CTRL and TWIST samples were sequenced, and cell types were identified using eGFP+ for keratinocytes and tdTomato+ for melanoma, after which the datasets were integrated. (C) UMAP of keratinocyte clusters from CTRL and TWIST conditions combined, showing normal keratinocyte cluster (NKC), tumor-associated keratinocyte (TAK), and Twist-High cluster (TWIST-specific). CTRL: 32.9% TAK. TWIST: 14.3% TAK, 56.3% Twist-High. (D) UMAP highlighting keratinocyte clusters, with a TAK cluster, an NKC, and a Twist-High cluster unique to the TWIST condition. (E) Melanoma cell state analysis of the melanoma cluster unique to the Twist condition indicated by arrow in C. (F) Schematic overview of CellChat analysis. In CTRL condition, we analyzed Ligand-Receptor pairs with NKC as sender and CTRL-Melanoma as receiver. In TWIST condition, we analyzed L–R pairs with both NKC and Twist-High as sender and TWIST-Melanoma as receiver. Panel F was created with BioRender.com. (G) CellChat analysis results. L–R pairs shown at p < 0.01, with color scale indicating communication probability of L–R pair.

We first analyzed the melanomas that arose in the control vs. TWIST animals by comparing their gene signature to well-annotated signatures of human melanoma cell states. It is now widely recognized that melanomas exist in at least four different transcriptional states, ranging from undifferentiated/invasive to neural crest, to intermediate, to melanocytic/proliferative (Tsoi et al., 2018). These distinct states are mutually exclusive, yet exhibit a high degree of plasticity, with interconversion between the states (Tsoi et al., 2018). Interestingly, we found a cluster of melanoma cells that developed only in the TWIST condition was enriched for the melanocytic/proliferative state but not undifferentiated/invasive gene markers (Figure 5E). Consistent with this, GSEA of the TWIST condition also showed significant enrichment for oxidative phosphorylation (Figure 5—figure supplement 1D–F), which we previously showed is representative of the melanocytic state (Lumaquin-Yin et al., 2023). Because this state is among the least invasive ones, this is consistent with our in vivo observations that these melanomas are phenotypically less invasive.

We hypothesized that this change in cell state might be induced by physical interactions between the TWIST keratinocytes and the melanoma cells. To address this, we analyzed potential cell–cell interactions using CellChat, a software tool that allows us to quantitatively characterize and visualize cell–cell communications using a curated zebrafish ligand–receptor interaction database (Jin et al., 2021; Figure 5F). Two unique ligand–receptor pairings were identified that only occur between TWIST keratinocytes and the melanomas that arose in these animals: a homotypic jam3b–jam3b interaction and a pgrn–sort1a (progranulin–sortilin) interaction (Figure 5G).

The jam3b interaction was of particular interest to us, as this protein has been recently identified as one required for melanophore survival in zebrafish (Eom et al., 2021) and for human melanoma metastasis (Arcangeli et al., 2012; Langer et al., 2011). In normal human skin, JAM3 is expressed in keratinocytes of the superficial epidermis, whereas its heterotypic partner JAM1 is only expressed in basal keratinocytes (Uhlén et al., 2015). Since melanomas largely arise in this basal area, this suggests that TWIST expression in the keratinocytes was leading to aberrant expression of jam3b, forming stronger homotypic jam3b–jam3b attachment between these keratinocytes and the melanoma cells and preventing their invasion.

Conservation in human studies

Finally, we wished to determine if the TWISThi state in our zebrafish models was also present in human melanomas. We use the GeoMX spatial transcriptomics platform to interrogate a series of early melanoma precursor lesions, a subset of a larger study using this method to look at spatial organization of melanoma stages (Vallius et al., 2026). We first stained sections with antibodies to cytokeratin (to mark keratinocytes) and SOX10/MART1 (to mark melanocytes and melanoma cells). This allows for selection of microregions of interest (MR) which are then subject to RNA-seq (Figure 6A–D). From this, we then queried for the TAK/TWIST gene sets found in the zebrafish. Whereas normal skin melanocytes showed no such enrichment, there was an enrichment of these signatures as the lesions progressed from early to late precursor (Figure 6E). This suggests that EMT-like alterations in keratinocytes are an early event in both zebrafish and human melanoma.

Figure 6
Download asset Open asset
Enrichment of keratinocyte signatures in microregions with melanocytic atypia.

(A–D) Visualization of a subset of reanalyzed GeoMx microregions (MRs); GeoMx data (n = 8 MRs) is reanalyzed from Vallius et al., 2026. Immunofluorescence staining is shown for morphology markers DNA (gray), panCK (green), and SOX10/MART1 (red; middle and bottom rows). The selected MRs are illustrated with white polygons (middle), and the shaded area within the MRs represent the area for oligonucleotide barcode collection (referred to as a ‘segmentation mask’ per NanoString; bottom). Serial sections were H&E-stained and the regions mapping to the sequenced MRs are shown (top row). The MRs were annotated as epidermis containing morphologically normal melanocytes (panel A), and regions of melanocytic atypia (precursor; panels B–D). Scale bars = 50 μm. (E–G) Melanocyte-adjusted TAK and TWIST gene-signature activity across reanalyzed MRs (n = 8 GeoMx MRs). (E) Heatmap showing gene-signature scores derived from the top 100 upregulated genes defining the TAK and TWIST programs across MRs. Each column represents an MR, and the MRs are annotated for sample ID and histopathological stage, with melanocyte fraction (melanocyte count divided by total nuclei count within the MR) per GeoMx MR displayed as a top annotation.

Discussion

In this study (summarized in Figure 7), we observed that keratinocytes in a zebrafish of melanoma undergo an EMT-like transformation in the presence of melanoma. This alteration is reminiscent of keratinocyte behavior during wound healing, in which keratinocytes exhibit markers and morphological changes associated with EMT in development (Leopold et al., 2012; Koike et al., 2020). Interestingly, we also observed an increase in N-cadherin expression in KC, which is usually attributed to melanoma as it becomes more aggressive and invades into the dermis to associate with fibroblasts. Our findings would suggest that if keratinocytes are also upregulating N-cadherin expression, their normal contact-based regulatory controls on melanoma may still be relevant, though not acting through E-cadherin. Additionally, reanalysis of a published zebrafish melanoma scRNA-sequencing dataset showed distinct populations of keratinocytes that expressed markers of EMT, demonstrating the feasibility of studying this KC population using our zebrafish model and nominating Twist1 as a potent EMT-TF in this cell (Hunter et al., 2021). Twist expression has been found to be upregulated at the edge of wounded skin upon treatment with bFGF, a well-characterized growth factor in melanoma (Koike et al., 2020; Halaban et al., 1988). If melanoma acts as an open-wound in the skin, then it is possible that Twist1 might be upregulated in TAKs to close this wound.

Figure 7
Download asset Open asset
Keratinocyte EMT and TWIST upregulation restrains melanoma invasion.

Schematic model of melanoma–keratinocyte interactions in the skin. During early melanoma development, melanocytes reside within an epidermis composed of normal keratinocytes. In the context of low TWIST expression (TWISTlo keratinocytes), melanoma cells adopt an invasive, undifferentiated phenotype and are able to breach the epidermal compartment. In the presence of EMT/TWISThi keratinocytes, melanoma cells are restrained within the epidermis and adopt a more differentiated, less invasive cell state. These findings suggest that EMT-like changes in tumor-associated keratinocytes, rather than promoting tumor spread, act as a cell-extrinsic barrier to melanoma invasion through altered cell–cell interactions in the local skin microenvironment. Created with BioRender.com.

Our zebrafish melanoma survival experiment showed improvements in overall survival in zebrafish with KCs overexpressing Twist compared to those that received an empty vector, despite the fish forming tumors at the same rate. This survival improvement was shown to be caused by a decrease in melanoma invasion, raising the possibility that Twist-overexpressing KCs could restrain melanoma invasion. Human cell co-cultures with a HaCaT cell line overexpressing Twist showed a similar finding to our in vivo zebrafish model, with reduced melanoma cell infiltration into keratinocytes. However, an important caveat to these cell culture experiments is that the keratinocyte densities vary when co-cultured with different melanoma lines (i.e. HS294T vs. SKMEL2), the reasons for which are not known. However, because we saw a consistent pattern across these two lines, it is less likely that our results are due to these differences in keratinocyte density alone.

To learn more about the dynamics of melanoma and TME KCs, we performed scRNA-sequencing on our zebrafish melanoma model to account for both KCs in contact with melanoma and those in the periphery. This yielded three KC clusters in fish overexpressing Twist, while two clusters in fish that received the empty vector; accounting for KCs overexpressing Twist and the KC populations we had observed previously in our reanalysis of a zebrafish melanoma scRNA-sequencing dataset (Hunter et al., 2021).

Interestingly, we found a cluster of melanoma cells unique to the TWIST condition, which shared gene signatures similar to that of the genes observed in both transitory and melanocytic states as published by Tsoi et al., 2018. These cell states were defined to be more differentiated with higher MITF expression and correlated to a more proliferative but less invasive cohort from Hoek et al., 2006. Thus, it is likely that this melanoma cluster unique to the TWIST condition could be responsible for the reduced overall melanoma invasion.

Further analysis using CellChat nominated pgrn–sort1a and jam3b–jam3b as unique interactions between Twist-High KCs and TWIST-Melanoma cells. As previously described, jam3b has been identified as a critical protein in zebrafish melanophore survival with known involvement in melanoma metastasis. The aberrant expression of jam3b on Twist-overexpressing KCs could indicate strong homotypic interaction with melanoma jam3b that retains the melanoma in the epidermis. Although the progranulin–sortilin interaction has not been characterized in melanoma, sortilin has been identified as a key regulator of progranulin levels (Hu et al., 2010). Progranulin is known to be constitutively expressed by KCs, which could be cleaved to epithelins that promote KC proliferation (Daniel et al., 2000; Shoyab et al., 1990). Progranulin is also a potent mediator of the wound response produced by dermal fibroblasts in addition to epidermal keratinocytes (He et al., 2003). Perhaps the increase in progranulin is cleared by melanoma cells through sortilin, resulting in the endocytosis and lysosomal transport of sortilin (Tanimoto et al., 2017). The degradation of sortilin could be responsible for decreased cell migration and invasion, as sortilin is required for the interaction of proNGF, a neurotrophin produced by melanoma, with NGFR in promoting melanoma migration (Truzzi et al., 2008). Further studies are needed to elucidate the precise mechanisms underlying the nominated interactions, jam3b–jam3b and pgrn–sort1a, and to explore their potential as therapeutic targets in melanoma.

Methods

Zebrafish husbandry

Zebrafish were maintained in a dedicated facility with controlled temperature (28.5°C) and salinity. The fish were kept on a 14-hr light/10-hr dark cycle and fed a standard zebrafish diet consisting of brine shrimp followed by Zeigler pellets. Embryos were obtained through natural mating and incubated in E3 buffer (5 mM NaCl, 0.17 mM KCl, 0.33 mM CaCl2, 0.33 mM MgSO4) at 28.5°C. For procedures requiring immobilization, zebrafish were anesthetized using Tricaine-S (MS-222, Syndel) prepared as a 4 g/l stock solution with a pH of 7.0. The stock solution was protected from light exposure and diluted to the appropriate concentration to achieve fish immobilization. All experimental procedures and animal protocols described in this manuscript were conducted in compliance with the Institutional Animal Care and Use Committee (IACUC) protocol #12-05-008, approved by the Memorial Sloan Kettering Cancer Center (MSKCC).

Generation of zebrafish line with fluorophore labeled keratinocytes

Embryos at the one-cell stage from the Casper Triple zebrafish line (mitfa:BRAFV600E;p53−/−;mitfa−/−;mpv17−/−) (Callahan et al., 2018; White et al., 2008) were injected with a krt4:eGFP expression cassette in the 394 vector of the Tol2Kit (Kwan et al., 2007) with tol2 mRNA. Larvae were sorted for positive GFP fluorescence at day 3 and raised to adult for breeding. F0 fish were in-crossed and resulting F1 were outcrossed with Casper Triple zebrafish for consistent GFP expression. Starting from F2, the krt4:eGFP zebrafish line was maintained by out-crossing with Casper Triple zebrafish and sorting for GFP expression.

Transgene Electroporation in Adult Zebrafish

TEAZ was utilized to generate melanoma as previously described (Callahan et al., 2018; Montal et al., 2023). Krt4:eGFP zebrafish (krt4:eGFP Casper Triple) were anesthetized with tricaine and injected with a plasmid solution containing miniCoopR-tdT (250 ng/μl), mitfa:Cas9 (250 ng/μl); zU6:sgptena (23 ng/μl) with target sequence GAATAAGCGGAGGTACCAGG, zU6:sgptenb (23 ng/μl) with target sequence GAGACAGTGCCTATGTTCAG, and the tol2 plasmid (55 ng/μl). For control electroporation without generating melanoma, zebrafish were injected with mitfa:tdTomato (250 ng/μl), mitfa:Cas9 (250 ng/μl), zU6:non-targeting (46 ng/μl), and tol2 plasmid (55 ng/μl). All fish were injected on the left flank below the dorsal fin and electroporated with the BTX ECM 830 electroporator using 3 mm platinum Tweezertrodes (BTX Harvard Apparatus; #45-0487). Electroporator settings used: LV Mode, 40 V, 5 pulses, 60-ms pulse length, and 1-s pulse interval. Electroporated zebrafish were screened for successful electroporation 7 days post-electroporation by tdTomato expression using fluorescence microscopy and melanoma tracked by imaging once per week. All live zebrafish imaging were performed with the Zeiss AxioZoom V16 fluorescence microscope.

sgRNA sequences for TEAZ listed below:

  • Nontargeting: AACCTACGGGCTACGATACG

  • ptena: GAATAAGCGGAGGTACCAGG

  • ptenb: GAGACAGTGCCTATGTTCAG

Confocal imaging of zebrafish epidermis

Zebrafish with or without melanoma were anesthetized in tricaine (MS-222, Syndel) as described above. Site of injection is visually identified by the presence of melanoma or the area below the dorsal fin. Scales were removed with tweezers and fixed in 4% PFA in PBS (Santa Cruz 281692) in a 96-well plate for 15 min. Fixed scales were washed three times with PBS and permeabilized with 0.1% Triton X-100 (Thermo Scientific 85111) in PBS, then blocked with 10% goat serum (Thermo Fisher 50062Z). Scales were incubated with 1:250 GFP polyclonal antibody, Alexa Fluor 488 (Thermo Fisher A21311) overnight at 4°C. The next day, scales were washed three times with PBS, incubated with 1:1000 Hoechst 33342 (Thermo Fisher H3570) for 1 hr and mounted onto slides with VECTASHIELD Vibrance Antifade Mounting Media (Vector Laboratories H-1700). Samples were imaged on the Zeiss LSM880 inverted confocal microscope, and images were processed using FIJI v1.53.

Flow cytometry of adult zebrafish cells

Zebrafish were euthanized using ice-cold water. Melanoma and adjacent skin were dissected from fish with melanoma, and skin alone was dissected from below the dorsal fin. Subsequently, samples were cut into 1-mm strips using a clean scalpel and placed into 15-ml conical tubes (Falcon 352099) with 3 ml of DPBS (Gibco 14190250) and 187.5 μl of 2.5 mg/ml Liberase TL (Roche 5401020001). Samples were incubated in dissociation solution at room temperature for 30 min on a shaker with gentle movement to prevent tissue from settling at the bottom of the tube. At 15 min of incubation, a wide bore p1000 pipette tip (Thermo Scientific 2079G) was used to gently pipette the sample up and down for 90 s. After 30 min, 250 μl of FBS (Gemini Bio) was added to stop the enzymatic activity of Liberase TL, and samples were pipetted up and down using a wide bore p1000 pipette tip for 90 s. Dissociated cells were then filtered through a 70-μm cell strainer (Falcon 352350) into a 50-ml conical tube (Falcon 352098) placed on ice. Samples were centrifuged at 500 × g at 4°C for 5 min and supernatant was removed by pipetting. The cell pellet was resuspended in 500 μl of PBS with 5% FBS and filtered again through 40 μm tip filters (Bel-Art H136800040) into 5 ml polypropylene tubes (Falcon 352063). For subsequent FACS analysis, 0.5 μl of 1000x DAPI (Sigma-Aldrich D9542) was added to each sample. Samples were FACS sorted (BD FACSAria) at 4°C for GFP-positive keratinocytes and tdTomato-positive melanoma gated using fluorophore-negative zebrafish controls.

Zebrafish tissue RNA extraction and real-time quantitative PCR

FACS sorted zebrafish cells were deposited directly into 750 μl TRIzol LS Reagent (Invitrogen 10296010) in Eppendorf DNA LoBind Tubes (Eppendorf 022431021). After collection, samples were snap-frozen using dry ice and stored at −80°C. RNA extraction was performed per TRIzol LS manufacturer protocols. For precipitation of RNA, 10 µg supplemental glycogen (Roche 10901393001) was used per sample to account for low cell numbers. Resulting RNA was resuspended in Nuclease-free water (Fisher Scientific AM9937). 25 ng RNA per sample was transcribed to cDNA using Superscript III First-Strand Synthesis System (Invitrogen 18080051). cDNA mix was diluted 1:10 with Nuclease-free water for real-time quantitative PCR using Power SYBR Green PCR Master Mix (Applied Biosystems 4368708) and the Bio-Rad CFX384 Touch Real-Time PCR System (Bio-Rad 1855484). Resulting Cq values were normalized to hatn10 as previously described.

qPCR primer sequences:

  • hatn10 fwd: 5′-TGAAGACAGCAGAAGTCAATG-3′

  • hatn10 rev: 5′-CAGTAAACATGTCAGGCTAAATAA-3′

  • mitfa fwd: 5′-GGCACCATCAGCTACAATGA-3′

  • mitfa rev: 5′-GAGACAGGGTGTTGTCCATAAG-3′

  • krt4 fwd: 5′-GGAGGTGTTTCCTCTGGTTATG-3′

  • krt4 rev: 5′-GAACCGAATCCTGATCCACTAC-3′

  • vim fwd: 5′-GGATATTGAGATCGCCACCTAC-3′

  • vim rev: 5′-GACTCTCGCAGGCTTAATGAT-3′

  • cdh2 fwd: 5′-GAGCCATCATCGCCATACTT-3′

  • cdh2 rev: 5′-CTTGGCCTGTCTCTCTTTATCC-3′

scRNA-sequencing analysis of Miranda et al.

Zebrafish scRNA-seq data from Hunter et al., 2021 was reanalyzed using R 4.2.0 and Seurat 4.3.0 (Butler et al., 2018; Hao et al., 2021). Cluster identities were maintained as published. Keratinocyte Module Scores were calculated using the AddModuleScore function with default parameters using published gene lists. Differential gene expression analyses between clusters were performed using FindMarkers. Differentially expressed gene lists were converted from zebrafish genes to human orthologs using DIOPT as previously described (Hu et al., 2011; Campbell et al., 2021). GSEA analysis on differentially expressed genes between keratinocyte clusters was performed using fgsea 1.22.0 and the Hallmark pathways set from MSigDB (Subramanian et al., 2005; Liberzon et al., 2015).

Twist overexpression in zebrafish keratinocytes

Twist1a (ENSDART00000043595.5) and twist1b (ENSDART00000052927.7) were TOPO cloned into the attL1-L2 Gateway pME vector and LR cloned into the pDestTol2pA2 vector (Tol2Kit 394) with p5E-krt4 promoter and p3E-polyA (Tol2Kit 302). To generate the zebrafish melanoma model as previously described, one-cell stage Casper Triple zebrafish embryos (mitfa:BRAFV600E;p53−/−;mitfa−/−;mpv17−/−) were injected with miniCoopR-tdTomato, krt4-eGFP, either krt4-twist1a and krt4-twist1b for Twist overexpression condition or empty vector for control condition, tol2 mRNA, and phenyl red. Injections were performed three times on different days with parents from the same clutch. Embryos were grown at standard conditions and sorted at 5 days post-injection for eGFP and tdTomato expression using the Zeiss AxioZoom V16 fluorescence microscope. eGFP+/tdTomato+ fish in the CTRL (n = 135) and TWIST (n = 118) conditions were maintained to adulthood.

Zebrafish imaging and tumor-free survival tracking

Zebrafish were regularly monitored for melanoma formation and survival every 4 weeks, beginning at 10 weeks post-fertilization. Melanoma formation was screened visually using the Zeiss AxioZoom V16 fluorescence microscope under ×20 magnification. Kaplan–Meier curves and corresponding statistics were generated using GraphPad Prism 9. Statistical differences in survival between conditions were determined by the Mantel–Cox log-rank test.

Histology of zebrafish samples

Zebrafish were euthanized in tricaine (MS222, Syndel). Each fish was dissected in three sections consisting of head, body, and tail. Samples were placed in 4% PFA in PBS (Santa Cruz 281692) for 72 hr on a shaker at 4°C, then paraffin embedded. Histology was performed by HistoWiz Inc (histowiz.com) using a Standard Operating Procedure and fully automated workflow. Samples were processed, embedded in paraffin, and sectioned at 5 μm. Immunohistochemistry was performed on a Bond Rx autostainer (Leica Biosystems) with enzyme treatment (1:1000) using standard protocols. Sections were stained with H&E or IHC with antibodies including BRAFV600E (ab228461) and GFP (ab183734). Bond Polymer Refine Detection (Leica Biosystems) was used according to the manufacturer’s protocol. After staining, sections were dehydrated and film coverslipped using a TissueTek-Prisma and Coverslipper (Sakura). Whole slide scanning (40x) was performed on an Aperio AT2 (Leica Biosystems).

Cell culture

Human melanoma lines A375, HS294T, and SKMEL2 were obtained from ATCC (CRL-1619, HTB-140, HTB-68). Human keratinocyte line HaCaT was obtained from AddexBio. All cells were routinely tested and confirmed to be free from mycoplasma. Cells were maintained in a humidified incubator at 37°C and 5% CO2. Cells were maintained in DMEM (Gibco 11965) supplemented with 10% FBS (Gemini Bio) and split when confluent, approximately 2–3 times per week.

Twist overexpression in HaCaT

The HaCaT cell line was labeled with eBFP to allow for identification during co-culture. 293T (ATCC CRL-11268) was transfected with the pLV-Azurite plasmid (Addgene 36086) with pMD2.5 (Addgene 12259) and psPAX2 (Addgene 12260) using Invitrogen Lipofectamine 3000 Transfection Reagent (Invitrogen L3000015) according to the manufacturer’s protocol. HaCaTs were infected with lentivirus containing CMV:eBFP and selected for eBFP positivity using ampicillin and FACS. Subsequently, HaCaT-eBFP was infected with lentivirus containing CMV:TWIST1 (Horizon Precision LentiORF Human TWIST1 OHS5898-202622685) or CMV:empty control created by removing ORF of CMV:TWIST1 plasmid. HaCaT line overexpressing TWIST1 was labeled as HaCaT-TWIST and HaCaT line with empty vector was labeled as HaCaT-CTRL. HaCaT lines were subsequently sorted for nuclear eGFP expression present in the plasmid as part of the Precision LentiORF system. HaCaT-CTRL and HaCaT-TWIST were cultured as previously described.

Western blot

Cells were washed with DPBS (Gibco 14190250) and lysed in RIPA buffer (Thermo Scientific 89901) with the addition of protease and phosphatase inhibitors (Thermo Scientific 78440). Lysates were centrifuged at 13,000 × g at 4°C and quantified using the Pierce BCA Protein Assay Kit (Pierce 23227). Samples were reduced with the Laemmli SDS-sample buffer (Boston BioProducts BP111R) and boiled for 10 min. For gel electrophoresis, samples were loaded into 4–15% precast protein gels (Bio-Rad 4561084), then transferred to 0.2-µm nitrocellulose membranes (Bio-Rad 1704158). Membranes were washed in TBST and blocked with 5% milk in TBST (Boston BioProducts P1400) for 1 hr at RT. Membranes were washed and incubated with primary antibodies overnight. Antibodies used include Twist (Abcam ab50887) and beta-actin (CST 3700S). On the next day, membranes were washed with TBST and incubated with appropriate secondary antibodies for 1 hr. Blots were incubated with the Immobilon Western Chemiluminescent HRP Substrate (Millipore WBKLS0500) and imaged with the Amersham ImageQuant 800.

Immunofluorescence

Cells were cultured on chamber slides (Thermo Scientific 154739) overnight at standard cell culture conditions. Culture media was washed with DPBS (Gibco 14190250) and fixed with 2% PFA in PBS (Santa Cruz 281692) for 15 min at RT. Cells were washed with DPBS and permeabilized with 0.1% Triton X-100 (Thermo Scientific 85111) in PBS, then blocked with 10% goat serum (Thermo Fisher 50062Z) for 1 hr at RT. Primary antibodies used include Twist (Abcam ab50887). Cells were incubated with primary antibody in 10% goat serum overnight at 4°C and washed with DPBS the next day, before incubation with 1:1000 Hoechst 33342 (Thermo Fisher H3570) for 1 hr and mounted with VECTASHIELD PLUS Antifade Mounting Medium (Vector Laboratories H1900). Slides were imaged on the Zeiss LSM880 inverted confocal microscope and images were processed using FIJI v1.53.

Melanoma infiltration assay

RFP-labeled melanoma cell lines, including A375, SKMEL2, HS294T, were plated on poly-L-lysine coated round glass coverslips (Corning 354085) placed in 24-well plates at 150–200k cells per coverslip. HaCaT cell lines were plated in 6-well plates at 250–300k cells per well. Cells were allowed to attach overnight and the coverslip containing melanoma cells was transferred to 6-wells containing HaCaT cell lines using tweezers. All coverslips were placed in the center of the well. KC–melanoma co-cultures were incubated in standard cell culture conditions for 24 hr. Co-cultures were imaged by fluorescence microscopy at four locations of each coverslip: top, right, bottom, and left, to capture variations in melanoma cell infiltration into KC lines. FIJI v1.53 was used to count the number of infiltrating melanoma cells per image and average infiltrating melanoma cells were calculated per well. All experiments were performed in 3 sets, with 3 replicates per set per condition. Average infiltrating cell numbers per well were normalized to the average infiltrating cell number per well in the HaCaT-CTRL condition. scRNA-sequencing analysis of zebrafish melanoma.

Six zebrafish, three each from CTRL and TWIST conditions at 26 weeks post-injection were selected for scRNA-sequencing of melanoma tumors. To account for set differences, one fish from each of three injection sets was chosen in each condition. Melanoma and adjacent skin were dissected from the fish, then enzymatically and mechanically dissociated into single-cell solutions as described above. The samples were FACS sorted for GFP and RFP positivity, corresponding to eGFP expressed by keratinocytes and tdTomato expressed by melanoma. The sorted cells were placed in (Gibco 11965) supplemented with 10% FBS (Gemini Bio) and 1% penicillin–streptomycin–glutamine. To enrich for KCs, sorted KCs and melanoma cells from each fish were recombined at a 7:3 KC:melanoma ratio. Sorted cells were pelleted and resuspended in DPBS + 0.1% BSA. Samples were also combined based on their genetic perturbation condition. Droplet-based scRNA-seq was performed using the Chromium Single Cell 3′ Library and Gel Bead Kit v3 (10X Genomics) and Chromium Single Cell 3′ Chip G (10X Genomics). 10,000 cells were targeted for encapsulation. GEM generation and library preparation were performed according to kit instructions. Libraries were sequenced on a NovaSeq S4 flow cell. Resulting reads were aligned to the GRCz11 reference genome with the addition of eGFP and tdTomato sequences using CellRanger v5.0.1 (10x Genomics). scRNA-sequencing analysis was performed as detailed above. In addition, melanoma cells were scored using AddModuleScore to assess their enrichment of genes associated with the four main melanoma cell states and intermediate states (Tsoi et al., 2018). The highest scoring gene module for each cell was annotated as its cell state. CellChat (Jin et al., 2021) was used to analyze cell–cell communication between KC and melanoma clusters using its zebrafish L–R database.

Spatial transcriptomics (GeoMx) data reanalysis

The GeoMx DSP spatial transcriptomics data was reanalyzed from Vallius et al., 2026. First, the reanalyzed MRs were selected based on the availability of the keratinocyte compartment within the sequenced MRs. Only the MRs with keratinocytes present within the profiled region (referred to as the ‘segmented area’ per Bruker Spatial Biology (formerly NanoString); examples shown in Figure 6A–D) were reanalyzed and included in this manuscript. The H&E-stained serial sections were used for histopathological annotation of the samples, as described in Vallius et al. The biospecimen metadata and GeoMx MR annotations are shown in Supplementary file 1.

To account for the variable cell-type composition across MRs, we first computed the melanocyte fraction for each MR as the ratio of melanocyte count to total nuclei count. The melanocyte count was determined by the SOX10/MART1 immunofluorescence staining, by visually counting the nuclei staining positively for SOX10 and/or MART1. A non-melanocyte scaling factor was then defined as 1 − melanocyte fraction. With the top 100 upregulated genes in the TAK and TWIST programs, the raw gene-signature scores were calculated using Gene Set Variation Analysis and subsequently normalized by dividing each score by the corresponding non-melanocyte scaling factor, yielding keratinocyte-adjusted signature values. For visualization, keratinocyte-adjusted signature scores were assembled into a matrix (signatures × MRs) and standardized by row (z-score scaling across MRs) to emphasize the relative differences between MRs. Heatmaps were generated using the ComplexHeatmap package in R, with columns representing samples and rows representing gene signatures. Rows were hierarchically clustered, while columns were ordered using optimal leaf ordering based on pairwise sample distances.

Statistics and reproducibility

Statistical analysis and figures were generated by GraphPad Prism 9, R Studio 4.2.0, and Biorender.com. Image processing was performed in FIJI v1.53. Statistical tests are described in figure legends and methods. Experiments were repeated at least three times unless otherwise noted. All animal and cell experiments were performed with a reasonable number of replicates by power calculations or feasibility of the experimental method.

Data availability

The sequencing data generated in this study are available via the Gene Expression Omnibus (GEO) database under accession code GSE333353. GeoMx gene expression data (Vallius et al., 2026) will be released via GEO coincident with the final publication of that work.

The following data sets were generated
    1. Ma Y
    (2026) NCBI Gene Expression Omnibus
    ID GSE333353. Restraint of melanoma progression by cells in the local skin environment.
    https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE333353

References

    1. Fitzpatrick TB
    2. Breathnach AS
    (1963)
    The epidermal melanin unit system
    Dermatologische Wochenschrift 147:481–489.
    1. Halaban R
    2. Kwon BS
    3. Ghosh S
    4. Delli Bovi P
    5. Baird A
    (1988)
    bFGF as an autocrine growth factor for human melanomas
    Oncogene Research 3:177–186.
  32. Book
    1. Montal E
    2. Suresh S
    3. Ma Y
    4. Tagore MM
    5. White RM
    (2024) Cancer modeling by transgene electroporation in adult zebrafish (TEAZ)
    In: Amatruda JF, Houart C, Kawakami K, Poss KD, editors. Zebrafish: Methods and Protocols. Springer US. pp. 83–97.
    https://doi.org/10.1007/978-1-0716-3401-1_5

Article and author information

Author details

  1. Yilun Ma

    1. Weill Cornell/Rockefeller/Sloan Kettering Tri-Institutional MD-PhD Program, New York, United States
    2. Department of Cancer Biology and Genetics, Memorial Sloan Kettering Cancer Center, New York, United States
    3. Cell and Developmental Biology Program, Weill Cornell Graduate School of Medical Sciences, New York, United States
    Contribution
    Conceptualization, Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Writing – original draft, Writing – review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-2297-9980

  2. Mohita Tagore

    Department of Cancer Biology and Genetics, Memorial Sloan Kettering Cancer Center, New York, United States
    Contribution
    Conceptualization, Formal analysis, Investigation, Methodology
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-2586-1427

  3. Miranda V Hunter

    Department of Cancer Biology and Genetics, Memorial Sloan Kettering Cancer Center, New York, United States
    Contribution
    Conceptualization, Data curation, Formal analysis, Visualization
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-6971-1738

  4. Ting-Hsiang Huang

    Department of Cancer Biology and Genetics, Memorial Sloan Kettering Cancer Center, New York, United States
    Contribution
    Conceptualization, Data curation, Visualization, Methodology
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-8739-0571

  5. Emily Montal

    Department of Cancer Biology and Genetics, Memorial Sloan Kettering Cancer Center, New York, United States
    Contribution
    Data curation, Formal analysis, Investigation
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-7855-3631

  6. Joshua M Weiss

    1. Weill Cornell/Rockefeller/Sloan Kettering Tri-Institutional MD-PhD Program, New York, United States
    2. Department of Cancer Biology and Genetics, Memorial Sloan Kettering Cancer Center, New York, United States
    3. Cell and Developmental Biology Program, Weill Cornell Graduate School of Medical Sciences, New York, United States
    Contribution
    Data curation, Formal analysis, Investigation, Methodology
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-7340-4004

  7. Yingxiao Shi

    1. Laboratory of Systems Pharmacology, Harvard Program in Therapeutic Science, Harvard Medical School, Boston, United States
    2. Department of Clinical Computational Oncology, Dana-Farber Cancer Institute, Boston, United States
    Contribution
    Data curation, Formal analysis, Investigation, Visualization, Methodology
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-2783-9642

  8. Tuulia Vallius

    1. Laboratory of Systems Pharmacology, Harvard Program in Therapeutic Science, Harvard Medical School, Boston, United States
    2. Ludwig Center, Harvard Medical School, Boston, United States
    3. Department of Systems Biology, Harvard Medical School, Boston, United States
    Contribution
    Conceptualization, Data curation, Formal analysis, Validation, Investigation, Visualization, Methodology
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-3006-4887

  9. Peter K Sorger

    1. Laboratory of Systems Pharmacology, Harvard Program in Therapeutic Science, Harvard Medical School, Boston, United States
    2. Ludwig Center, Harvard Medical School, Boston, United States
    3. Department of Systems Biology, Harvard Medical School, Boston, United States
    Contribution
    Resources, Supervision, Project administration
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-3364-1838

  10. Richard M White

    1. Department of Cancer Biology and Genetics, Memorial Sloan Kettering Cancer Center, New York, United States
    2. Nuffield Department of Medicine, Ludwig Institute for Cancer Research, University of Oxford, Oxford, United Kingdom
    Contribution
    Conceptualization, Formal analysis, Supervision, Funding acquisition, Writing – original draft, Project administration, Writing – review and editing
    For correspondence
    richard.white@ludwig.ox.ac.uk
    Competing interests
    Senior editor, eLife
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-9099-9169

Funding

National Institutes of Health (F30CA265124)

  • Yilun Ma

National Institutes of Health (R01CA238317)

  • Richard M White

National Institutes of Health (T32GM007739)

  • Yilun Ma

National Cancer Institute (1K99CA266931)

  • Miranda V Hunter

National Institute of Health (K00CA223016)

  • Emily Montal

National Institute of Health (F30CA236442)

  • Joshua M Weiss

National Institute of Health (T32GM008539)

  • Joshua M Weiss

American Cancer Society (PF-24-1316850-01-CD)

  • Tuulia Vallius

NIH/NCI Cancer Center (P30 CA008748)

  • Richard M White

National Institute of Health (R01CA229215)

  • Richard M White

Memorial Sloan Kettering Cancer Center (Marie-Josée Kravis Women in Science Endeavor Postdoctoral Fellowship)

  • Mohita Tagore

Memorial Sloan Kettering Cancer Center (Experimental Immuno-Oncology Scholars Fellowship)

  • Mohita Tagore

Melanoma Research Alliance

  • Richard M White

Ludwig Institute for Cancer Research

  • Richard M White

Ludwig Center at Harvard

  • Tuulia Vallius
  • Peter K Sorger

The Mark Foundation for Cancer Research (ASPIRE Award)

  • Tuulia Vallius
  • Peter K Sorger
  • Yingxiao Shi

The funders had no role in study design, data collection, and interpretation, or the decision to submit the work for publication.

Acknowledgements

We thank members of the White lab for valuable discussions and feedback on this project. We thank the MSKCC aquatics core facility for their support of zebrafish maintenance, the MSKCC flow cytometry core for assistance with cell sorting, the Single-Cell Analytics Innovation Lab (SAIL) for their help with scRNA-sequencing sample preparation and data processing, and the Molecular Cytology Core for their help with imaging. YM was supported by a Medical Scientist Training Program grant from the NIH under award number T32GM007739 to the Weill Cornell/Rockefeller/Sloan Kettering Tri-Institutional MD-PhD Program and a Kirschstein-NRSA predoctoral fellowship under award number F30CA265124. MT was supported by the Marie-Josée Kravis Women in Science Endeavor Postdoctoral Fellowship and the Experimental Immuno-Oncology Scholars Fellowship from MSKCC. MVH was supported by a K99/R00 Pathway to Independence Award from the National Cancer Institute under award number 1K99CA266931. EM was supported by the Predoctoral to Postdoctoral Fellow Transition Award (F99/K00) from the NIH under award number K00CA223016. JMW was supported by a Kirschstein-NRSA predoctoral fellowship from the NIH under award number F30CA236442 and a predoctoral fellowship from the Cell and Developmental Biology Program at Weill Cornell Graduate School under NIH award number T32GM008539. TV was supported in part by a Research Scholar Grant (PF-24-1316850-01-CD) from the American Cancer Society. PKS, TV, and YS are supported by the Ludwig Center at Harvard, and an ASPIRE Award from The Mark Foundation for Cancer Research. R.M.W. was funded through the NIH/NCI Cancer Center Support Grant P30 CA008748, the Melanoma Research Alliance, NIH Research Program Grants R01CA229215 and R01CA238317, and the Ludwig Institute for Cancer Research.

Ethics

The GeoMX data were performed on de-identified biopsy specimens.

This study was performed in strict accordance with the recommendations in the Guide for Care and Use of Laboratory Animals of the NIH. All animals were handled according to approved MSKCC IACUC protocol number 12-05-008.

Version history

  1. Sent for peer review:
  2. Preprint posted:
  3. Reviewed Preprint version 1:
  4. Reviewed Preprint version 2:
  5. Version of Record published:

Cite all versions

You can cite all versions using the DOI https://doi.org/10.7554/eLife.101974. This DOI represents all versions, and will always resolve to the latest one.

Copyright

© 2024, Ma et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 781
    views
  • 42
    downloads
  • 0
    citations

Views, downloads and citations are aggregated across all versions of this paper published by eLife.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

  1. Yilun Ma
  2. Mohita Tagore
  3. Miranda V Hunter
  4. Ting-Hsiang Huang
  5. Emily Montal
  6. Joshua M Weiss
  7. Yingxiao Shi
  8. Tuulia Vallius
  9. Peter K Sorger
  10. Richard M White
(2026)
Restraint of melanoma progression by cells in the local skin environment
eLife 13:RP101974.
https://doi.org/10.7554/eLife.101974.3

Share this article

https://doi.org/10.7554/eLife.101974

Categories and tags

Research organism