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¹. 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², and this interaction is essential for skin and hair color^3,4. Despite decades of research, our understanding of keratinocytes in the context of melanoma remains incomplete.

Keratinocytes have been shown to both inhibit and promote melanoma initiation. They are tightly adherent to melanocytes through homophilic interactions of the cell adhesion molecule E-cadherin⁵. Through E-cadherin, keratinocytes can control melanocyte growth and behavior^6–9. As melanoma progresses, they undergo a cadherin subtype switch in which they downregulate E-cadherin and upregulate N-cadherin, thereby escaping keratinocyte-mediated controls while promoting migration and survival^10–12. Beyond cadherins, invasive melanoma can downregulate matricellular proteins such as CCN3, which usually facilitates melanoma attachment to the basement membrane^13,14. In addition to melanoma intrinsic changes, alterations in epidermal keratinocytes such as the loss of PAR3 expression results in a local environment that facilitates melanoma invasion and metastasis¹⁵. 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^16–18.

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¹⁹. 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^20,21. Using a combination of cell-type specific genetic manipulations, in vivo imaging, and single-cell transcriptomics, we found that tumor-associated keratinocytes 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 clinically meaningful.

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²². 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 in Adult Zebrafish)^19,23, 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 BRAF^V600E, sgRNAs against PTEN (zebrafish ptena/b) and germline loss of p53, using mitfa-Cas9 to ensure that PTEN were inactivated only in melanocytes and not in surrounding skin cells (Figure 1C). To account for skin wounding from electroporation in assessing changes in tumor-associated keratinocytes, 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.

After initiating tumorigenesis, we examined the morphology of tumor-adjacent keratinocytes to investigate whether the growing tumor could be influencing keratinocyte behavior. We performed confocal microscopy on the scales of fish 8 weeks post-electroporation. 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^24–26. 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 (Figure 1B).

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 tumor-associated keratinocytes (TAKs) and normal keratinocytes (NKCs) from tissue without melanoma revealed upregulation of EMT markers vimentin and N-cadherin in TAKs, confirming the EMT-like morphological changes of keratinocytes in 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 as published in Tagore et al.¹⁸ (Figure 1H). 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 keratinocytes (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 tumor-associated keratinocytes

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²⁶. We next wanted to understand which of these transcription factors was responsible for the EMT-like behavior in tumor-associated keratinocytes. To address this, we utilized an existing scRNA-sequencing dataset of a BRAF^V600E-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 zebrafish keratinocyte populations (Figure 2B). Subsequent differential gene expression and GSEA analysis of the two keratinocyte clusters revealed one cluster with enrichment for EMT, similar to what we observed in Figure 1 (Figure 2C-D). We refer to this EMT cluster as Tumor Associated Keratinocytes (TAKs), and the other cluster as a Normal Keratinocyte Cluster (NKCs). We focused on three EMT-transcription factors expressed by zebrafish keratinocytes in this dataset, 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.

Keratinocyte TWIST restrains melanoma invasion in zebrafish

Having identified Twist as a potential driver of EMT-like changes in tumor-associated keratinocytes, we next asked how it affected melanoma phenotypes. To investigate the role of Twist in keratinocytes, we created new transgenic zebrafish in which the keratinocyte-specific krt4 promoter was used to drive the two zebrafish TWIST paralogs (twist1a and twist1b) in the presence of BRAF^V600E-driven melanomas. We injected the plasmids for melanoma initiation to zebrafish embryos and then monitored the fish for tumor-free survival as well as overall survival over the ensuing 26 weeks (Figure 3A). 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.

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 and surrounding tissues from the control and TWIST conditions (Figure 3C). The oncogenic driver in this melanoma model is hBRAF^V600E and serves as an IHC marker for the tumor cells. Comparison of hBRAF^V600E 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 (Supp. Figure 1). 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.

To further validate this finding, we developed a cell culture-based invasion assay by pairing the HaCaT keratinocyte cell line with multiple melanoma cell lines that were pre-cultured on coverslips (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. HaCaT keratinocytes were transformed via lentiviral transduction to overexpress human 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 in vivo. This finding was recapitulated using the SKMEL2 human melanoma cell line, demonstrating the inhibitory effect of TWIST1 overexpression in keratinocytes on melanoma invasion across different cell lines (Figure 4G-H). Collectively, our results demonstrate that induction of EMT in keratinocytes is associated with reduced melanoma invasion and improvement in animal survival.

Figure 4.

Keratinocyte EMT promotes aberrant adhesion to nascent melanoma cells

While EMT of tumor cells 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 CTRL vs. TWIST overexpression tumors and adjacent microenvironment 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 (Figure 5A and Supp Figure 2A). In the resultant scRNA-seq dataset, we identified distinct clusters of eGFP+ keratinocytes and tdTomato+ melanoma cells (Figure 5B). Of the keratinocyte clusters identified, two were present in both the CTRL and TWIST conditions, while a third was only present in samples overexpressing TWIST in keratinocytes (Figure 5C). To identify biological processes that may be active within the keratinocyte, we performed GSEA, comparing differentially expressed genes between the two clusters present in both conditions. GSEA identified an enrichment in the EMT pathway as seen in Figure 2, allowing us to label the EMT-enriched cluster as TAK and other as NKC (Figure 5D and Supp Figure 2B). As expected, we also identified a unique cluster of keratinocytes only present in the TWIST dataset that highly expresses twist1a/b, the genes that we overexpressed in the TWIST condition, and labeled this cluster Twist-High (Figure 5D and Supp Figure 2C).

Figure 5.

To understand how the Twist-high keratinocytes may be restraining melanoma progression, we compared gene expression in melanomas that arose in the presence of control vs. Twist overexpression keratinocytes by comparing their gene signatures to published melanoma gene signatures representing a range of differentiation states²⁸. It is now widely recognized that melanoma cells exist along a trajectory of differentiation, ranging from undifferentiated/invasive, to neural crest, to intermediate, to melanocytic/proliferative²⁸. Interestingly, we found a cluster of melanoma cells that developed in the TWIST condition was enriched for the melanocytic/proliferative state but not undifferentiated/invasive gene markers (Figure 5E). 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-high 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²⁹ (Figure 5F). Two unique ligand-receptor pairings were identified that only occur between Twist-high keratinocytes and the melanomas that arose in these animals: a homophilic 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³⁰ and for human melanoma metastasis^31,32. In normal human skin, JAM1 (or F11R) is expressed in keratinocytes of the superficial epidermis, whereas its heterophilic partner JAM3 is exclusively found in basal keratinocytes^33,34. This distribution is significant because melanomas predominantly originate in the basal area of the skin. Interestingly, basal keratinocytes, but not superficial keratinocytes, have been shown to inhibit melanocyte growth⁷. This observation, combined with our findings, suggests that Twist expression in the keratinocytes was leading to aberrant expression of jam3b, resulting in stronger homophilic jam3b-jam3b attachments between these keratinocytes and melanoma cells, potentially inhibiting melanoma invasion.

Discussion

In this study, we observed that keratinocytes in both zebrafish and human models 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^25,35. Interestingly, we 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 a subpopulation of keratinocytes could maintain contact with melanoma through upregulation of N-cadherin. Additionally, re-analysis 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 transcription factor in this cell type²⁷. Twist expression has been found to be upregulated at the edge of wounded skin upon treatment with bFGF, a well-characterized growth factor produced by melanoma^35,36. If melanoma acts as an open-wound in the skin, then Twist1 might be upregulated in tumor-associated keratinocytes as an adaptive response to close this wound.

Our zebrafish melanoma survival experiment showed improvements in overall survival in zebrafish with keratinocytes 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 keratinocytes 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. To learn more about the dynamics of melanoma and TME keratinocytes, we performed scRNA-sequencing on our zebrafish melanoma model to account for both keratinocytes in contact with melanoma and those in the periphery. Compared to the fish that received an empty vector control that had two subpopulations of keratinocytes, fish with Twist overexpression contained an additional novel keratinocyte subpopulation with overexpression of twist1a/b²⁷.

Interestingly, we also 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.²⁸. 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.³⁷. The enrichment for specific melanoma cell states when surrounded by a Twist-overexpressing keratinocyte TME could be responsible for the reduced overall melanoma invasion.

Further analysis using CellChat nominated jam3b-jam3b and pgrn-sort1a as unique interactions between Twist-High keratinocytes 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 keratinocytes could indicate strong homophilic interactions with melanoma jam3b that retains the melanoma in the epidermis. Whether jam3b acts in concert or independent of cadherin-based adhesion remains to be determined in future studies. Although the progranulin-sortilin interaction has not been characterized in melanoma, sortilin has been identified as a key regulator of progranulin levels³⁸. Progranulin is known to be constitutively expressed by keratinocytes, which could be cleaved to epithelins that promote KC proliferation^39,40. Progranulin is also a potent mediator of the wound response produced by dermal fibroblasts in addition to epidermal keratinocytes⁴¹. Perhaps the increase in progranulin is cleared by melanoma cells through sortilin, resulting in the endocytosis and lysosomal transport of sortilin⁴². 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⁴³. 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-hour light/10-hour 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 CaCl₂, 0.33 mM MgSO₄) 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:BRAF^V600E;p53-/-;mitfa-/-;mpv17-/-)^19,44 were injected with a krt4:eGFP expression cassette in the 394 vector of the Tol2Kit⁴⁵ 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)

TEAZ was utilized to generate melanoma as previously described^19,46. Krt4:eGFP zebrafish (krt4:eGFP Casper Triple) were anesthetized with tricaine and injected with a plasmid solution containing miniCoopR-tdT (250ng/μl), mitfa:Cas9 (250 ng/μl); zU6:sgptena (23 ng/μl), zU6:sgptenb (23ng/μl), and the tol2 plasmid (55ng/μl). For control electroporation without generating melanoma, zebrafish were injected with mitfa:tdTomato (250ng/μl), mitfa:Cas9 (250 ng/μl), zU6:non-targeting (46ng/μl) and tol2 plasmid (55ng/μl). All fish were injected on the left flank below the dorsal fin and electroporated with the BTX ECM 830 electroporator using 3mm platinum Tweezertrodes (BTX Harvard Apparatus; #45-0487). Electroporator settings used: LV Mode, 40V, 5 pulses, 60ms pulse length, and 1s 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: 5’-AACCTACGGGCTACGATACG-3’

• ptena: 5’-GAATAAGCGGAGGTACCAGG-3’

• ptenb: 5’-GAGACAGTGCCTATGTTCAG-3’

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 minutes. 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. Next day, scales were washed three times with PBS, incubated with 1:1000 Hoechst 33342 (Thermo Fisher H3570) for 1 hour 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 1mm 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 minutes on a shaker with gentle movement to prevent tissue from settling at the bottom of the tube. At 15 minutes of incubation, a wide bore p1000 pipette tip (Thermo Scientific 2079G) was used to gently pipette the sample up and down for 90 seconds. After 30 minutes, 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 seconds. 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 500g at 4°C for 5 minutes 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 (RT-qPCR)

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, 10ug 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). 25ng 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 RT-qPCR 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’

Re-analysis of zebrafish scRNA-sequencing data from Hunter et al

Zebrafish scRNA-seq data from ref²⁷ was re-analyzed using R 4.2.0. and Seurat 4.3.0^47,48. 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 (DGE) analyses between clusters were performed using FindMarkers. Differentially expressed gene lists were converted from zebrafish genes to human orthologs using DIOPT as previously described^49,50. GSEA analysis on differentially expressed genes between keratinocyte clusters was performed using fgsea 1.22.0 and the Hallmark pathways set from MSigDB^51,52.

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:BRAF^V600E;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 20X 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 72h on a shaker at 4°C, then paraffin embedded. Histology was performed by HistoWiz Inc. (histowiz.com). 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. 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) 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 manufacturer 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,000g 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 minutes. 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 hour at RT. Membranes were washed and incubated with primary antibodies overnight. Antibodies used includes 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 hour. 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 minutes 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 hour 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 hour 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 is 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 hours. Co-cultures were imaged by fluorescence microscopy at 4 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 were 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 keratinocytes, sorted keratinocytes and melanoma cells from each fish was 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 was 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²⁸. The highest scoring gene module for each cell was annotated as its cell state. CellChat²⁹ was used to analyze cell-cell communication between KC and melanoma clusters using its zebrafish L-R database.

Statistics and reproducibility

Statistical analysis and figures were generated by GraphPad Prism 9, R Studio 4.2.0 and Biorender.com. Image processing were 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.

Acknowledgements

We thank the MSKCC Single Cell Analytics Innovation Lab for their consult with single cell RNA-sequencing, the MSKCC Molecular Cytology Core Facility for assistance with imaging, the MSKCC Aquatics Core Facility for zebrafish care and maintenance, and the MSKCC Flow Cytometry Core Facility for assistance with FACS. We also thank the members of the White Lab for helpful discussions. YM was supported by a Medical Scientist Training Program (MSTP) grant from the National Institute of General Medical Sciences of the National Institutes of Health (NIH) under award number T32GM152349 to the Weill Cornell/Rockefeller/Sloan Kettering Tri-Institutional MD-PhD Program, the Kirschstein-NRSA predoctoral fellowship (F30) from the National Cancer Institute (NCI) of the NIH under award number F30CA265124, and the Barbara and Stephen Friedman Pre-doctoral Fellowship from MSKCC. MT was support by the Marie-Josée Kravis Women in Science Endeavor Postdoctoral Fellowship and the Experimental Immuno-Oncology Scholars Fellowship. MVH was funded by a K99/R00 Pathway to Independence Award from the NCI(1K99CA266931), the Scholarship for the Next Generation of Scientists from the Cancer Research Society, and a postdoctoral fellowship from the Canadian Institutes of Health Research. EM was supported by the NIH Individual Predoctoral to Postdoctoral Fellow Transition Award under award number 5K00CA223016. JMW was supported by the F30 fellowship from the NCI of the NIH under award number F30CA236442, a predoctoral fellowship (T32) by the Cell and Developmental Biology Program at Weill Cornell Graduate School supported by the NIH under award number T32GM008539, and a MSTP grant from the NIH under award number T32GM007739. RMW was funded through the NIH/NCI Cancer Center Support Grant P30 CA008748, the Melanoma Research Alliance, NIH Research Program Grants R01CA229215 and R01CA238317, NIH Director’s New Innovator Award DP2CA186572, The American Cancer Society, The Alan and Sandra Gerry Metastasis Research Initiative at the Memorial Sloan Kettering Cancer Center, The Harry J. Lloyd Foundation, and Oxford Ludwig Cancer Research.