Generation of a double binary transgenic zebrafish model to study myeloid gene regulation in response to oncogene activation in melanocytes

Immune cells and their associated molecules are responsible for detecting and eliminating cancer cells and their precursors at early stages of cancer development (Williams et al., 2016; de Visser et al., 2006). A complex network of inflammatory genes is closely linked to somatic cell transformation and malignant disease. Importantly, analysis of the initial response of immune cells to cancer might lead to the discovery of targets for both prevention and treatment. As with pathogens, cancer cells can evade the immune system and as such tumour growth and ultimately the emergence of clinically detectable cancer is largely dependent on the capacity of cancerous cells to evade and manipulate the immune response (Blair and Cook, 2008). Myeloid cells are equipped with sensors for damage and inflammation as well as subsequent effector mechanisms for resolving inflammation (Feng et al., 2010). However, transformed cells are able to shape the nature of the myeloid cells and evoke an immunosuppressive response, enabling progression of the disease (Engblom et al., 2016). Further studies are needed to understand the different inflammatory signalling pathways in tumour initiation and progression, and how they may be modulated for immunotherapy. Because we cannot predict when and where transformed cells may arise in the human host, little is known about the sterile inflammatory signalling pathways that closely follow somatic cell transformation (Feng et al., 2012).

The zebrafish is an accepted model in cancer research, with many aspects of carcinogenesis being conserved between teleosts and humans (White et al., 2013; Stern and Zon, 2003). Unlike mammalian models, the zebrafish is easily amenable to techniques that allow the study of both cancer initiation and progression (Feng et al., 2010). Zebrafish are particularly useful in studying melanoma because their melanocytes are externally visible and large single cells can be directly visualised live (Ceol et al., 2008).

Previous transgenic zebrafish melanoma models used oncogenic BRAF^V600E (Patton et al., 2005), NRas^Q61K (Dovey et al., 2009) and HRas^G12V (Michailidou et al., 2009; Anelli et al., 2009) under the melanocyte-specific microphthalmia-associated transcription factor a (mitfa) promoter, coupled to an additional mutation in the tumour suppressor gene p53. In the mitfa-driven HRas^G12V transgenic line, melanoma development is rare but secondary mutations in genes involved in the PI3K signalling pathway contribute to its occurrence (Patton et al., 2005; Dovey et al., 2009; Anelli et al., 2009). Santoriello et al. subsequently developed a conditional transgenic zebrafish line that uses the Gal4 system and the melanocyte-specific kita promoter to drive HRas^G12V expression (Vallone et al., 2007; Santoriello et al., 2010). kita:HRas^G12V induced early-onset melanoma without additional mutations in tumour suppressor genes. However, some reports suggest that high levels of Gal4 expression can be toxic and may result in developmental defects (Emelyanov and Parinov, 2008). Moreover, upstream activating sequences (UAS) have been shown to be susceptible to DNA methylation leading to transcriptional silencing of the transgene, which is minimal in the first generation but exacerbated upon propagation through later generations (Goll et al., 2009; Akitake et al., 2011).

The study of the myeloid response to the earliest precursors of melanoma would require temporal control of carcinogenesis. The current zebrafish models for benign naevus and cutaneous melanoma rely on the targeted expression of a human oncogene in melanocytes, limiting the ability to control melanoma initiation. Additionally, constitutive expression of the oncogene renders maintenance of stable transgenic lines difficult because these models may develop severe tumours before the fish reach reproductive age (Nguyen et al., 2012).

To overcome these limitations, we set out to develop an inducible system of melanoma initiation. Here, we report the generation of a new inducible zebrafish model, specifically designed to study the onset of melanoma using the mifepristone-inducible LexPR system developed by Emelyanov and Parinov (Emelyanov and Parinov, 2008). Additionally, we have developed an innovative transgenic zebrafish model system for in vivo biotinylation and subsequent isolation of neutrophil (previously described in Kenyon et al., 2017) and macrophage nuclei based on the biotagging method (Trinh et al., 2017). Using this approach, it is possible to perform genome-wide analysis of the active nuclear transcriptomes of neutrophils and macrophages. By combining the biotagging approach that enables isolation and analysis of either neutrophil or macrophage nuclei with the inducible model for melanoma initiation, we have created a powerful double binary system for regulatory profiling of myeloid cells that respond to the precursors of melanoma. In a proof-of-concept analysis, we reveal a number of interesting factors upregulated by neutrophils responding to the oncogene-activated melanocytes that may promote melanoma progression.

To generate an inducible model for the initiation of melanoma tumorigenesis, we adapted a hormone-inducible binary system for targeted gene expression in zebrafish, previously used to drive liver carcinogenesis (Emelyanov and Parinov, 2008; Nguyen et al., 2012). In this system, the transcriptional activator expressed by the transgenic driver is a ligand-dependent chimeric transcription factor, termed the LexPR transactivator. The LexPR transactivator is a fusion protein containing the DNA-binding domain (DBD) of the bacterial LexA repressor, a truncated ligand-binding domain (LBD) of the human progesterone receptor and the activation domain of NF-κB/p65 protein. The effector transgenic fish line harbours the operator-promoter sequence, consisting of the synthetic LexA operator fused to a minimal 35S promoter sequence from the Cauliflower mosaic virus, upstream of the reporter gene. In the presence of the progesterone agonist, mifepristone, which binds within the LexPR transactivator, LexPR activates the transcription of the target reporter gene by binding to the LexA-binding sites within the LexA operon (LexOP) positioned upstream (Fig. 1A).

In our system, we placed the LexPR transactivator, followed by Cerulean reporter, under the control of the melanocyte-specific promoters kita and mitfa using bacterial artificial chromosome (BAC)-mediated transgenesis (Fig. S1). Two melanocyte-specific driver lines [TgBAC(kita:LexPR-Cerulean)^ox123 (hereafter referred to as kita:LexPR-Cerulean) and TgBAC(mitfa:LexPR-Cerulean)^ox124 (hereafter referred to as mitfa:LexPR-Cerulean)] were generated (Fig. 1B,C).

To generate the effector lines, different Ras oncogenes (NRas^Q61K, HRas^G12V and KRas^G12V) fused to an mCherry fluorescent reporter were cloned downstream of the LexOP. The effector cassettes were placed within a non-autonomous maize dissociation (Ds) element, which allows for effective activator (Ac)-mediated genomic integration using the transposase Ac/Ds (AcDs) system (Emelyanov et al., 2006) (Fig. S1). This system is highly efficient and results in multiple independent integrations into the zebrafish genome, thus allowing for high-level expression of the Ras effectors activated by LexPR in the context of oncogenic transformation. Three Ras-effector lines were generated: Tg(Crysβ:ECFP-LexOP:mCherry-NRas^Q61K)^ox125 (hereafter referred to as LexOP:mCherry-NRas^Q61K), Tg(Crysβ:ECFP-LexOP:mCherry-HRas^G12V)^ox126 and Tg(Crysβ:ECFP-LexOP:mCherry-KRas^G12V)^ox127. The Ras-effector lines all contain the selection marker enhanced cyan fluorescent protein (ECFP) under the control of the lens-specific zebrafish crysβ promoter (crystallin betaβ), which results in blue fluorescence in the lens. This can be used to screen for transgene integration, without driving the oncogene (Fig. 1D).

When a driver line is crossed to an effector line, in the presence of the synthetic progesterone ligand, mifepristone, mCherry–Ras-oncogene is transcribed only in melanocytes. Two different melanocyte-specific driver lines (kita and mitfa) and three different mCherry–Ras-oncogene lines were generated in this study, paving the way for multiple combinations of the system. We characterised two of these lines, Tg(kita:LexPR-Cerulean;Crysβ:ECFP-LexOP:mCherry-NRas^Q61K)^ox129 (Fig. 1E) and Tg(mitfa:LexPR-Cerulean; Crysβ:ECFP-LexOP:mCherry-NRas^Q61K)^ox130 (Fig. 1F) (hereafter referred to as kita:LexPR;LexOP:NRas^Q61K and mitfa:LexPR;LexOP:NRas^Q61K, respectively) and showed that, when embryos positive for both alleles were reared in the presence of mifepristone, melanocyte-specific expression of the mCherry-NRas fusion protein was detected. Transformed melanocytes are seen in pairs or clumps of oncogene-activated red cells (as shown by the white arrowheads in Fig. 1G-G‴), indicative of their proliferative potential. In the mitfa:mCherry-NRas^Q61K genotype, low levels of melanin were occasionally observed in melanocytes expressing mCherry (data not shown).

To assess and characterise the effect of mitfa- and kita-targeted oncogene activation, we compared the proliferative potential and morphology of control and activated melanocytes. We performed this analysis at two different time points: at 5 days post-fertilisation (dpf) (Fig. 2A-C), when the melanocytes are likely to be post-mitotic as the embryonic pigmentation pattern is thought to be largely complete at 48 h post-fertilisation (Kelsh et al., 2009, 1996) and at 12 dpf (Fig. 2D-I), just before the adult pigmentation patterning commences at 2 weeks post-fertilisation (Rawls et al., 2001).

To ensure optimal visualisation of the mCherry-NRas^Q61K-expressing melanocytes at 5 dpf, embryos were partially bleached to ensure that excessive pigment in melanoblasts did not mask the mCherry signal and the signal amplified with an anti-Cherry antibody (Fig. 2A-B‴). A total of 5-9 z-stacks were acquired for the dorsal cranial region of each larvae and the melanocytes were counted within each corresponding plane (Fig. 2A-B‴). A nuclear stain (magenta) was used to aid in the visualisation of individual melanocytes (Fig. 2B). Strikingly, we observed that oncogene-activated melanocytes were present in clusters spanning multiple cell layers, in contrast to controls, where melanocytes lay within a single basal skin layer (Fig. 2A′-B″). We sought to quantify the degree of melanocyte proliferation in kita:mCherry-NRas^Q61K and mitfa:mCherry-NRas^Q61K embryos by precisely counting the total number of melanocytes within the outlined area in the presence and absence of mifepristone (+/− mif). Both kita and mitfa transgenic lines showed a statistically significant increase (P<0.01 and P<0.0001, respectively) in the number of melanocytes in the region of interest at 5 dpf (Fig. 2C). To ensure that the increase in melanocytes was a result of oncogene activation rather than an artefact of mifepristone application, we quantified the number of melanocytes in 5-dpf embryos reared as previously stated in the presence or absence of mifepristone at 5 dpf and could detect no difference (Fig. S2). To characterise the progression of the oncogene-activated cells at 12 dpf, we similarly examined both morphology and the proliferative state of the cranial melanocytes. As previously described, kita:mCherry-NRas^Q61K embryos were reared in the presence or absence of mifepristone. Whereas, in control embryos, the melanocytes appear to round up at 12 dpf, mCherry-NRas^Q61K-expressing melanocytes exhibited an elongated fibroblast-like shape (Fig. 2D,D′, Fig. S3). Although larval melanophores are mostly replaced during larval-to-adult metamorphosis, some melanophores leave their initial position to form the stripes of the adult. Interestingly, these cells are characterised as small round black cells and are increasingly melanised (Parichy, 2006). To quantify the melanocytes, nuclear stain (green) was used to aid in the counting and visualisation of individual melanocytes cells (Fig. 2E,E′). Strikingly, cranial melanocytes were clearly present as multinucleated clusters compared to controls (Fig. 2E,E′). Melanocytes of control embryos were present in a single plane, in contrast to the oncogene-activated melanocytes, which were present in several overlapping layers of interlinked cells (Fig. 2F). Quantification of the melanocytes within the outlined area in Fig. 2D in the presence and absence of mifepristone at 12 dpf (+/− mif) (Fig. 2G-H‴) revealed a statistically significant increase (P<0.0001) in the number of melanocytes in the region of interest (Fig. 2I). All experiments are representative of two or three individual experiments.

Genome-wide analysis of the regulatory networks that govern the inflammatory response in vivo is complicated by the fact that relatively small numbers of responding cells are present within complex tissues. The current zebrafish lines for studying macrophage and neutrophil function make use of fluorescent reporters under the control of the mpeg1 and mpx promoters, respectively (Ellett et al., 2011; Renshaw et al., 2006). These lines have largely been used for functional studies at the cellular level, allowing for a detailed examination of host-pathogen interactions, comparative analysis of macrophage and neutrophil behaviours responding to inflammation, and the dynamic interaction between the two cell types (Renshaw et al., 2006; Ellett et al., 2011). Based on the biotagging method (Trinh et al., 2017), we have developed a zebrafish model to study the gene regulatory landscape of macrophages and neutrophils by cell-specific in vivo biotagging of nuclei, thereby enabling us to obtain the active transcriptome of these cell types (Fig. 3A). In our system, the biotagging driver line harbours cell-specifically (macrophage or neutrophil) expressed E. coli biotin ligase BirA. Myeloid BirA driver lines TgBAC(mpeg1:BirA-Citrine)^ox122 [hereafter referred to as mpeg1:BirA-Citrine (Fig. 3B) and TgBAC(mpx:BirA-Citrine)^ox121 (hereafter referred to as mpx:BirA-Citrine (Fig. 3C) (described in Kenyon et al., 2017)] were generated using BAC-mediated transgenesis (Fig. S1). The previously characterised effector line Tg(bactin:Avi-Cerulean-RanGap)^ct700a (hereafter referred to as bactin:Avi-Cerulean-Rangap) (Trinh et al., 2017) features ubiquitous expression of a zebrafish-compatible version of Avi-tagged nuclear-envelope-associated fusion protein. The Avi-tag is a 14-amino-acid minimal biotinylation substrate, specifically biotinylated by biotin holoenzyme synthetase, BirA, which adds a single biotin moiety onto a lysine residue within this peptide (de Boer et al., 2003). In addition to the biotinylatable Avi-tag, this effector protein consists of the carboxyl-terminal domain of Ran GTPase-activating protein1 (RanGap1), which acts as an outer-nuclear-envelope-targeting sequence, and a fluorescent variant Cerulean, allowing visualisation of the fusion protein (Fig. 3D,E). When the driver line is crossed to the effector line, in the specific cells (neutrophils or macrophages) that carry both biotagging alleles (driver and effector), the Avi-tagged fusion protein gets biotinylated and localises to the outer nuclear envelope. This process effectively results in the biotinylation of nuclei in a cell-specific (macrophage or neutrophil) fashion, allowing isolation of nuclei from these cell types using streptavidin-coated magnetic beads. Subsequently, the content of these nuclei can be analysed in a genome-wide fashion by next-generation sequencing (NGS) techniques. The use of mpeg1:BirA-Citrine and mpx:BirA-Citrine drivers combined with the Avi-Cerulean-Rangap effector offers the unique possibility to interrogate initial, often crucial, regulatory responses in specific cell types. Cell isolation protocols that are based on expression of fluorescent reporters always feature a delay between the time of expression of the cell-type-specific marker and the downstream analysis, because the fluorophore requires a certain amount of time to reach steady-state levels for visualisation or isolation by fluorescence-activated cell sorting (FACS). Given that cell-specific biotinylation begins very shortly after the cells start expressing BirA, the use of the biotagging system minimises such delays and enables faithful representation and readout of the relevant cell states. This is particularly crucial for analysis of cell populations in rapidly developing zebrafish embryos.

To show that the mpeg1/mpx driver lines generated in this study faithfully recapitulate larval inflammatory phenotypes, tailfin transections (Renshaw et al., 2006) were carried out at 3 dpf in mpeg1:BirA-Citrine and mpx:BirA-Citrine embryos, respectively. Injured specimens were fixed at 5 h post-injury (hpi) (mpeg1:BirA-Citrine) and 1 hpi (mpx:BirA-Citrine), and stained using anti-GFP antibody to amplify the Citrine signal and thus detect macrophages and neutrophils expressing BirA driver transgenes and their associated fluorescent reporter. Both mpeg1:BirA-Citrine (Fig. 3F) and mpx:BirA-Citrine (Fig. 3G) showed infiltration of macrophages and neutrophils to the site of injury.

To study the initial myeloid response to oncogene-transformed melanocytes, we have generated a double binary system in zebrafish by coupling the inducible LexPR driver/effector to the biotagging system for nuclear transcriptional profiling of macrophages and neutrophils (Fig. 4A, Fig. S4). Since NRas^Q61K is the most relevant oncogene to human melanoma, with 30% of melanomas having contributing mutations in NRas (Ceol et al., 2008), it was prioritised for the purpose of this study. We hypothesised that, if NRas^Q61K-transformed melanocytes generated by the inducible LexPR binary system were immunogenic, an increase in neutrophils in proximity to melanocytes following oncogene activation may be expected. kita:LexPR;LexOP:NRas^Q61K fish were in-crossed and embryos were grown to 24 hpf, at which point mifepristone is added to the embryo rearing medium in experimental embryos or omitted in non-treatment controls. Larvae were reared until 5 dpf, fixed and stained with an anti-mpx antibody to detect neutrophils. We chose to focus on the cranial region as this would be the focus of future experiments and capture only the neutrophils in the skin of the larvae in the vicinity of the melanocytes as compared to the controls. Images were captured as 11 z-stack projections to a depth of 100 μm from the surface of the embryo in kita:LexPR;LexOP:NRas^Q61K (+mif) and in genotype and stage-matched controls (−mif) larvae, and neutrophils were counted within the outlined area (Fig. 4B,B′). We observed a significant increase in the number of neutrophils when compared to stage- and genotype-matched controls with no oncogene activation (P<0.01) (Fig. 4C). Furthermore, high-resolution images show that neutrophils were detected in the vicinity as well as directly interacting with the melanocytes in the tail (Fig. 4D) and cranial region (Fig. 4D′) following oncogene activation.

To test the efficacy of this model, we carried out a proof-of-concept study and performed transcriptional profiling of neutrophils responding to NRas^Q61K-activated melanocytes, with the aim to identify components of the cellular response that might contribute to tumour progression. All subsequent experiments were carried out by crossing the kita:LexPR;LexOP:NRas^Q61K inducible melanoma system to the mpx:BirA;bactin:Avi-Rangap biotagging system, allowing isolation of neutrophil nuclei responding to oncogene activation. NRas^Q61K expression in the melanocytes of resulting embryos from this cross was driven by addition of mifepristone to the water at 24 hpf. Embryos expressing mCherry-NRas^Q61K in the melanocytes were selected for the profiling experiments, allowed to develop until 5 dpf and used for downstream analysis. Approximately 45 larvae from each condition – experimental (+mif) and control (−mif) – were used per experiment. A total of three biological replicates were collected. In each experiment, prior to nuclei isolation and analysis, cranial regions containing major foci of melanocyte oncogenic transformation were dissected to eliminate the anterior yolk sac and intermediate cell mass (ICM), thus enriching for responding neutrophils (Fig. 5A) (Bennett et al., 2001). By focussing our analysis on the isolated oncogene-activated mCherry-expressing regions in the head, we enriched the samples with environing neutrophils responding to transformed melanocytes. Moreover, mCherry-NRas^Q61K-expressing goblet cells found in the trunk region were excluded from the analysis.

‘Biotagged’ nuclei were isolated from dissected tissue by nuclei pulldown. Nuclear RNA pools were extracted and libraries prepared using SMART-seq™ v4 technology for cDNA synthesis and amplification for small cell numbers followed by Illumina Nextera library preparation. Libraries were sequenced on an NGS platform.

Scatterplots of technical replicates show reproducibility between the three independent experiments (Fig. S5). Differential expression analysis comparing active transcription in neutrophil nuclei of embryos with NRas-transformed melanocytes versus control embryos identified 194 upregulated genes and 144 downregulated genes with a statistical significance cut off of P<0.05 (Fig. 5B,C, Fig. S3, Table S1).

The analysis of RNA sequencing (RNA-seq) datasets showed that a number of genes that may contribute to melanoma progression were upregulated in neutrophils upon engaging with transformed melanocytes. Those include mediators of tumour growth, angiogenesis and extracellular matrix (ECM) remodelling (Fig. 5D, Fig. S6), suggesting their potential pro-tumorigenic role in promoting progression to melanoma. In addition, we found significant enrichment for a number of factors involved in neutrophil migration, adhesion and cytoskeletal rearrangement as well as key factors involved in neutrophil biological processes (Fig. S6). To confirm the results of the RNA-seq analysis, we performed quantitative real-time PCR (qPCR) for a number of genes that were upregulated in the RNA-seq analysis. We detected an increase in serpinb1l4 transcripts (Fig. 5D). SerpinB1 is an inhibitor of neutrophil serine proteases, including neutrophil granule proteins, elastase and cathepsin G (Benarafa et al., 2011). One potential mechanism by which neutrophils could eliminate transformed cells is by proteases stored in neutrophil granules. This data suggests that, instead, this mechanism is inhibited and neutrophils exhibit an immunosuppressive phenotype, allowing for melanocyte proliferation. We also detected an increase in members of the fibroblast growth factor (FGF) family, fgf1 and fgf6 (Fig. 5E,F), known to control cellular proliferation, survival, migration and differentiation. Uncontrolled FGF signalling, implicated in diverse tumour types, can drive tumorigenesis (Turner and Grose, 2010). Notably, we detect a statistically significant increase in cathepsin H (Fig. 5G). Cathepsins are a group of lysosomal proteinases or endopeptidases, of which different members can play a role in different tumorigenic processes, including proliferation, angiogenesis and metastasis (Tan et al., 2013). Furthermore, we see a statistically significant upregulation in members of the galectin family – galectin 3 and galectin 1 (Fig. 5H,I) – a group of proteins that bind β-galactoside and whose expression correlates with tumour development and invasiveness (Danguy et al., 2002). Other upregulated transcripts as confirmed by qPCR include tgm, dusp23 and calpain b (Fig. 5J-L).

The prolonged exposure (21 days) to the progesterone agonist mifepristone has biological consequences at the level of gene expression in zebrafish larvae and adults (Bluthgen et al., 2013). To confirm that the transcriptional changes exhibited by the neutrophils were a direct result of oncogene-activated melanocytes rather than the 4-day exposure to mifepristone, nuclei pulldowns were carried out in mpx:BirA;bactin:Avi-Rangap fish. Embryos were grown to 24 hpf and mifepristone added to the E3 medium of the experimental batch, whereas control embryos were reared without mifepristone. Larvae were reared until 5 dpf, nuclear pulldowns performed and corresponding nuclear RNA pools extracted. Reverse transcription (RT)-qPCR assaying a number of candidate genes found to be upregulated in neutrophils upon the response to NRas activation in melanocytes was carried out. We found little or no change in transcript levels of tgm1, fgf6a, fgf1, cathepsin L, gstp2, calpain b, galectin 1, galectin 3, dusp23, cathepsin H and serpinb14 upon the addition of mifepristone (Fig. S7A-K), suggesting that the results obtained after 4 days of exposure to mifepristone were a result of activation of NRas^Q61K rather than mifepristone exposure.

Although progress has been made in understanding melanoma biology, the incidence of melanoma continues to increase and the prognosis for melanoma patients still remains poor (Anelli et al., 2009). Here, we report an inducible model to study the molecular mechanisms following oncogene activation in melanocytes that allows for the spatiotemporal control of oncogene expression in zebrafish melanocytes and thus the potential to study the earliest precursors of melanoma. Additionally, we have developed a binary model designed for regulatory nuclear profiling of myeloid cells. Coupled together, these systems provide the potential to understand the mechanisms by which macrophages and neutrophils may promote melanoma progression or conversely eliminate cancer cells as they establish themselves as the precursors of a tumour. Although not described in this study, the biotagging model can be extended to purification of ribosomes using the translating ribosomal affinity purification method (TRAP) based on the in vivo biotinylation of the Avi-tagged Rpl10 protein in zebrafish (Heiman et al., 2008; Trinh et al., 2017). A full understanding of the RNA landscape of neutrophils at the onset of melanoma would require both ribosomal and nuclear profiling to obtain the translated and active transcriptome, respectively (Trinh et al., 2017).

Importantly, in combination or as individual systems, both binary models presented in this study provide the potential to study different aspects of melanoma biology, including the transformation of melanocytes, myeloid immunity and immunotherapy. For example, the larval phenotype observed in this zebrafish melanoma model can provide a powerful platform for testing drug targets for preclinical testing of therapeutics. Because of the larval phenotype that we see, the mCherry-positive transformed cells can provide a direct biological readout in fast, easy-to-score chemical screens, aimed at finding compounds or drugs that may revert the overproliferation phenotype observed in melanocytes. Previous studies have shown that activated components of the Ras signalling pathway drive a tumour-promoting inflammatory response (Feng et al., 2010). However, these findings do not address the difference in the inflammatory response that different driver mutations may elicit. In our model, the same driver line can be crossed to multiple oncogene effector lines. This provides further opportunities to study mechanistic details of melanoma initiation and progression with different driver mutations.

In this proof-of-concept study, we chose to use this model system to profile the neutrophil response to NRas^Q61K-oncogene-activated melanocytes using the kita driver. Neoplastic cells may only initiate tumour formation in the context of a supportive tumour microenvironment, largely influenced by myeloid cells. Although an immunosuppressive mechanism for neutrophils has been demonstrated in solid tumours, the role for neutrophils in tumour initiation largely remains to be explored. Although it is now accepted that tumour-associated neutrophils develop a pro-tumorigenic phenotype, largely driven by the presence of TGF-β, inhibition of TGF-β can modulate the cells to be tumouridical (Fridlender et al., 2009). In untreated tumours, neutrophils were reported to contribute to tumorigenesis by secreting factors that promote tumour growth, ECM remodelling and suppress the immune system (Galdiero et al., 2013). A study using a Tet-on inducible transgenic zebrafish line expressing oncogenic KRas^G12V in hepatocytes resulted in accelerated hepatocellular carcinoma (HCC) (Chew et al., 2014; Yan et al., 2017). The authors attributed the gender disparity of HCC, which occurs more frequently and aggressively in men than in women, to differences in cortisol levels, and reported a higher number of tumour-infiltrating neutrophils and macrophages in male zebrafish. After 7 days of doxycycline exposure in 3-month-old zebrafish, both infiltrating macrophages and neutrophils had a pro-tumour gene expression profile, although the effect was higher in male zebrafish than in their female counterparts (Yan et al., 2017).

We show here that neutrophils responding to early melanoma onset may provide an early source of Fgf1 and Fgf6. Given the multiple functional roles for FGF signalling in different tumours, the strong enrichment in FGF components in the neutrophils responding to oncogene-transformed melanocytes suggests a pro-tumorigenic role in promoting progression to melanoma. For example, uncontrolled FGF signalling in tumours can lead to both expansive tumour growth and progression to metastasis (Korc and Friesel, 2009). Furthermore, as FGF signalling can also promote angiogenesis, it is possible that the early role of FGFs in melanoma may be linked to angiogenesis, which not only promotes tumour growth but also the advancement from a pre-malignant to a malignant phenotype (Massi et al., 2010; Yang et al., 2015). In line with this interpretation, excessive angiogenesis, a hallmark of melanoma, and the progression from a radial growth phase to a vertical growth phase has been shown to require high angiogenic activity dependent on FGF1, FGF2 and VEGFA (Massi et al., 2010). In particular, FGF1 is a direct activator of phosphatidyl inositol 3-kinase (PI3K)-AKT signalling in endothelial cells known to initiate migration and invasion (Yang et al., 2015). These findings are particularly meaningful given that naturally targeting paracrine signalling in the form of growth and angiogenic factors is poised to provide exciting targets for cancer therapy.

Traditionally, cysteine cathepsin proteases are largely accepted as degradative enzymes of the lysosome. More recently, secreted cathepsins have emerged as potent effectors of multiple processes during tumour development, including the turnover and degradation of the ECM as well as processing or degradation of various growth factors, cytokines and chemokines. Consequently, cathepsins may promote tumour growth, tissue invasion and metastasis (Olson and Joyce, 2015). In a model of pancreatic islet carcinogenesis, deletion of cathepsin H significantly altered angiogenic switching of pre-malignant hyperplastic islets and ultimately a reduction in the number of tumours formed (Gocheva et al., 2010). Cathepsin H activity was largely attributed to macrophages in close proximity to the blood vessels (Gocheva et al., 2010). Here, we provide evidence that neutrophils may provide a source of cathepsin H in melanoma initiation.

Another interesting class of molecules implicated in tumour progression are the galectins, a family of lectin carbohydrate-binding proteins that function both intracellularly and extracellularly (Fortuna-Costa et al., 2014). Multiple examples exist in the literature for the role of galectins in tumour progression, including facilitating neoplastic transformation, tumour cell survival, angiogenesis and tumour metastasis (Liu and Rabinovich, 2005). Galectin 1, for example, increases cellular growth and exerts its effects by binding to the ECM as well as other cells (Liu and Rabinovich, 2005). In melanoma, galectin 1 has been thought to cause resistance of melanoma cells to cytotoxic stimuli and to enhance angiogenesis (Mathieu et al., 2012). Extracellular galectin 3 may influence tumour progression by impeding the endocytosis of key receptors, including TGFβ and EGF, while at the same time inducing endocytosis of β1 integrins (Furtak et al., 2001; Partridge et al., 2004). Extracellular galectin 3 secreted by tumour cells has also been shown to induce angiogenesis (Liu and Rabinovich, 2005). Furthermore, since galectins are expressed by immune cells, both intracellular and extracellular galectins may serve to regulate immune cell function (Liu and Rabinovich, 2005). Galectin 1 inhibits the release of inflammatory mediators by neutrophils and reduces transendothelial migration of the cells in response to inflammatory stimuli (Liu and Rabinovich, 2005; La et al., 2003). Therefore, targeting galectins derived from neutrophils in benign neoplasms may serve to prevent malignant transformation by promoting a pro-inflammatory environment. Galectin 3, on the other hand, can function as a chemokine attracting both monocytes and macrophages (Sano et al., 2000). In this capacity, galectin 3 may facilitate tumour progression driven by macrophages.

The zebrafish embryo has previously been used as a model organism to study the interactions between macrophages and neutrophils and oncogene-transformed melanocytes and goblet cells using a kita:HRas^G12V transgenic line (Feng et al., 2010). This study revealed a recruitment of both neutrophils and macrophages to oncogene-transformed cells. The authors reported frequent cytoplasmic tethers between the immune cells and transformed cells, including neutrophils, as well as phagocytosis of transformed cells by macrophages (Feng et al., 2010). In agreement with this data, we show, at the transcriptional level, an upregulation in signalling molecules involved in cell migration as well as cytoskeletal reorganisation in neutrophils.

In conclusion, we developed a double binary transgenic system for studying the myeloid response to the earliest precursors of melanoma and demonstrated it as a powerful resource, offering exciting potential to reveal inflammatory mediators that may contribute to tumour promotion and progression. This unbiased, genome-wide approach, carried out at the earliest stages of somatic cell transformation in vivo, is currently the only study of its kind. The ensuing gene set provides a list of interesting pathways and specific candidates to explore further in human melanoma biology. The results described here contribute to the growing body of evidence that suggests a tumour-promoting role of neutrophils in cancer biology and provides insight into the mechanisms by which neutrophils can be harnessed for immunotherapy.

This study was carried out in compliance with local ethical approval from the University of Oxford (UK) and using procedures authorised by the UK Home Office in accordance with UK law (Animals Scientific Procedures Act 1986). Zebrafish were maintained as described previously (Westerfield, 2000). Wild-type embryos for transgenesis were obtained from AB or AB/TL mix strains.

The pDs(cry:ECFP-LexOP:Cherry) construct was a kind gift from Dr Alexander Emelyanov and Dr Sergey Parinov (Emelyanov and Parinov, 2008). pCrysβ:ECFP-LexOP:mCherry-NRas^Q61K, pCrysβ:ECFP-LexOP:mCherry-HRas^G12V and Crysβ:ECFP-LexOP:Cherry-KRas^G12V were cloned by In-Fusion cloning (Clontech, cat. no. 638909) using pBabe-NRas 61K (plasmid #12543; Addgene), mEGFP-HRas^G12V (plasmid #18666; Addgene) and pEFm.6 HA-KRas^v12 (a gift from Prof. Xin Lu, Ludwig Institute for Cancer Research Oxford, Oxford, UK) were used as donor vectors. The optimised open reading frame (ORF) of the LexPR transactivator, which eliminates all putative donor and acceptor sites that could interfere with proper production of the LexPR transcript, was synthesised based on the coding sequence obtained from Dr Parinov (Temasek Life Sciences Laboratory, National Institute of Singapore, Singapore) and assembled using G-Blocks (Integrated DNA Technologies, BVBA, Leuven, Belgium). pGEM LexPR-2A-Cerulean-SV40pA-FRT-Kan-FRT donor plasmid for BAC transgenesis was cloned by In-Fusion cloning. pGEM BirA-2A-Citrine-SV40pA-FRT-Kan-FRT donor plasmid for BAC transgenesis was cloned by fusion PCR of HA-BirA-2A amplified from PMT-HA-BirA-2A-mCherryRas (Trinh et al., 2017) and inserted into pGEM-GFP-SV40pA-FRT-Kan-FRT where the GFP ORF was been replaced with Citrine. All plasmids generated for the purposes of this study are listed in Table S2 and available from Addgene (https://www.addgene.org/Tatjana_Sauka-Spengler/).

The TgBAC(mpeg1:BirA-Citrine)^ox122 and TgBAC(kita/mitfa:LexPR-Cerulean)^ox123/ox124 transgenic lines were generated by BAC-mediated transgenesis as previously described (Trinh et al., 2017). Briefly, PCR-amplified BirA-2A-Citrine-SV40pA-FRT-Kan-FRT and LexPR-2A-Cerulean-SV40pA-FRT-Kan-FRT cassettes were recombined into the first coding exon of a driver gene within the corresponding BAC clone (macrophage-specific mpeg1 gene, and melanocyte-specific mitfa or kita, respectively). In a second recombination step, an iTol2-Ampicillin cassette (provided by Prof. Kawakami, National Institute of Genetics, Mishima, Japan) was introduced into the BAC backbone as previously published (Suster et al., 2011; Bussmann and Schulte-Merker, 2011). Wild-type embryos were injected at the one-cell stage with 200 ng/μl of purified BAC DNA and 100 ng/μl tol2 transposase mRNA. Putative founders were outcrossed to wild-type fish and offspring screened for Citrine or Cerulean expression, in combination with PCR amplification of the transgene when expression levels were low. The full list of zebrafish transgenic lines generated and used in this study can be found in Table S3.

In transgenic embryos harbouring both a transactivation driver and oncogene effector allele, melanocyte-specific expression of mCherry–Ras-oncogene was activated by addition of 1 μM mifepristone (Sigma, cat. no. M8046) to the E3 medium and changed daily.

Embryos/larvae were fixed in 4% PFA in phosphate buffered saline (PBS) (0.1 M, pH 7.4) for 1 h at room temperature (RT). Embryos were bleached (3% H₂O₂, 0.8% KOH and 0.5% Tween20) and post-fixed in 4% PFA in PBS for 1 h at RT. Embryos were then rinsed 3×15 min in PBT (2% DMSO, 0.5% Triton in PBS) and incubated in block solution (10% donkey serum in PBT) for 2 h at RT, followed by incubation with primary antibody (1:300 in block solution) overnight (O/N) at 4°C. Embryos/larvae were then washed 3-5 times for 1 h at RT, followed by an O/N wash at 4°C in PBT and then incubated with secondary antibody (1:1000 in PBT). Embryos were then washed 6-8 times for 1 h at RT+O/N at 4°C in PBT. The primary antibody used to detect neutrophils was rabbit anti-mpx (GeneTex, cat. no. GTX128379) and the secondary antibody used was Alexa-Fluor-488-conjugated anti-rabbit IgG (Life Technologies, cat. no. R37116). For amplification of mCherry signal we used Living Colors^® mCherry Monoclonal Antibody (1:200) (Clontech, cat. no. 632543) as primary and Alexa-Fluor-568-conjugated goat anti-mouse IgG (H+L) (Life Technologies, cat. no. A-11004) as secondary antibody. For quantification of melanocytes, the far-red nuclear counterstain RedDot1 (Cambridge Biosciences, BT40060) was used. Images of stained embryos/larvae were taken on an inverted Zeiss 780 confocal microscope or captured using the Zeiss Z1 Light Sheet microscope.

At 3 dpf, larvae were anaesthetised with 0.02% 3-amino benzoic acidethylester (Tricaine) in E3 medium and tail transection was performed with a sterile scalpel.

For nuclei isolation, zebrafish embryos were anaesthetised with 0.01% Tricaine and cranial regions dissected away from anterior yolk sac and ICM. Embryos were washed in hypotonic buffer H [20 mM HEPES (pH 7.9), 15 mM MgCl₂, 10 mM KCl, 1 mM DTT and 1× complete protease inhibitor (Roche)] and subsequently re-suspended in 500 μl of buffer H. From here, nuclei isolation was carried out as previously described (Deal and Henikoff, 2010; Trinh et al., 2017). Following isolation, nuclei-beads were re-suspended in RNA lysis buffer from RNAqueous^®-Micro Kit (Life Technologies, cat. no. AM1931).

Total nuclear RNA was extracted using RNAqueous^®-Micro Kit according to the manufacturer's instructions. RNA integrity was checked with an RNA Pico chip (Agilent Technologies, cat. no. 5067-1513) using the Agilent 2100 Bioanalyzer.

cDNA was synthesised and amplified from 100 pg to 1 ng of input RNA using SMART-seq™ v4 ultra-low input kit for RNA (Clontech Laboratories, cat. nos 634888, 634889, 634890, 634891, 634892, 634893 and 634894). Sequencing libraries were prepared using the Nextera XT DNA library preparation kit and NGS was performed on a NextSeq500 platform using a NextSeq™500 150 cycle High Output Kit (Illumina, cat. no. FC-404-1002) to generate 80-basepair paired end reads.

Read quality was evaluated using FastQC (Andrews, 2011). Reads were mapped to the ‘Jul. 2014 Zv10/danRer10 assembly’ version of the zebrafish genome using STAR (v.2.4.2a) splice-aware aligner (Dobin et al., 2013). Count tables were generated using subread FeatureCounts (v1.4.5-p1q), with standard parameters (Liao et al., 2014). Differential expression was carried out using DESeq2 R package (Love et al., 2014). Data generated in this study were submitted to GEO (GSE96534).

RNA extraction and cDNA synthesis were carried out using the RNAqueous^®-Micro Kit (Life Technologies, cat. no. AM1931) and Superscript III Reverse Transcriptase (Invitrogen, cat. no. 18080093), respectively. Quantitative PCR was performed using Fast SYBR Green Master Mix (Applied Biosystems, cat. no. 4385612) on a 7500 Fast Real-Time PCR system (Applied Biosystems). Gene-specific primers were used for tgm (5′-TGCGATTATGACTGCAAGCA-3′ and 3′-TCCAGCACTTCATCAGGTCA), fgf6a (5′-GCCAGACGGAAGGATAAACG-3′ and 3′-GACAAACAGCCCGCTTTTCA-5′), fgf1 (5′-AAAGCAGTGAAAGCAGGAGT-3′ and 3′-AGAAACACTCTTCATTCAGCACA-5′), cathepsin l (5′-CTGCTGCTTTTCTGGCTGTT-3′ and 3′-ACGGTGTGACTCTTCCTCTG-5′), gstp2 (5′-TCTCTTTGGACAGCTGCCTA-3′ and 3′-CGTTCTTTCCATACGCACCA-5′), calpain b (5′-AGCTCAGGACTCCTCTCTCT-3′ and 3′-CCTTGGCAGATGTCAGTCCT-5′), galectin 3 (5′-GACAGCCTGGACAACTGACT-3′ and 3′-AGAGGGAAATCAAGTGGCAC-5′), galectin 1 (5′-TGTGCTTATACAGAATATGGCCT-3′ and 3′-GCGATGTCTTCAGAGCTGTT-5′), serpinb1l4 (5′-CTCGGTGCCAAAGGAAACAC-3′ and 3′-CCATTTGGGCTTCTGGGTTG-5′), cathepsin H (5′-CGTTGTACACCGAGGAAGATG-3′ and 3′-GGTCAATCCTCTTCTTGTTCTCC-5′) and dusp23 (5′-CGAGTCGCAGTCAAATCCTG-3′ and 3′-CGTCCTGCCGTGACCATG-5′). Expression levels were normalised to β-actin.

Graphs and statistical analyses were generated by GraphPad prism, version 6. For quantification of the melanocytes and neutrophils, statistical significance was determined by two-tailed unpaired Student's t-test with Welch's correction. The following P-values are indicated in the figures: not significant, P>0.05; **P<0.01; ****P<0.0001.

We thank Isabella Watts for her assistance on the project and Profs Tudor Fulga and Ahmed Ahmed for their pertinent comments on the manuscript.