Ethics statement

This study was performed with the approval of the University of Bath ethics committee and in full accordance with the Animals (Scientific Procedures) Act 1986, under Home Office Project Licenses 30/2937 and P87C67227.

Fish husbandry

Embryos were obtained from natural crosses. Staging was performed according to Kimmel et al. [34]. Unless stated otherwise, we used the WIK stock for experiments in WTs, and the following mutant lines: sox10^t3 [26], mitfa^w2[19], ltk^ty82 [9] and tfec^ba6 (Petratou et al., in prep.). The tfec^ba6 allele shows recessive loss of function, and was generated via CRISPR/Cas9 directed mutagenesis. It corresponds to the deletion of 6 nucleotides from the 7^th exon of the gene. This deletion of two amino acids is predicted to interrupt critical spacing in the second alpha-helix of the dimerization domain [35]. Phenotypically, homozygotes show a phenotype (nearly complete loss of iridophores and failure to inflate the swim-bladder) identical to those resulting from frameshift mutations in the DNA binding domain (Petratou et al., in prep.). Embryos were obtained by incrossing heterozygous carriers for each mutant allele, with WT siblings were used as controls.

Transcript detection in whole mount embryos

Detailed information on the preparation of materials and the protocols for performing chromogenic whole mount in situ hybridisation (WISH) as well as multiplex fluorescent RNAscope can be found in Petratou et al. [36]. Probes used for chromogenic WISH were sox10 [26], foxd3 [37], ltk [9], pnp4a [15], mitfa [19]) and tfec (NM_001030105.2). To generate the tfec probe, cDNA prepared from total RNA extracted from 72 hpf zebrafish embryos was amplified with the following primers: forward 5’-AGCCAACAATCACGACAGTG-3’ and reverse 5’-CCAATAGAAACGGGAGGTCA-3’. The product was cloned into pCR II-BluntTOPO vector (Invitrogen) and the orientation assessed by sequencing. The plasmid was linearised with PstI restriction enzyme (NEB) and in vitro transcription was with the SP6 polymerase of the DIG labelling kit (Roche; Cat# 11175025910). For multiplex RNAscope, the following probes were used: ltk (ACD; Cat No. 444641), tfec (ACD; Cat No. 444701), mitfa (ACD; Cat No. 444651), sox10 (ACD; Cat No. 444691) and foxd3 (ACD; Cat No. 444681).

Embryos were imaged using an upright compound Imager 2 microscope (Zeiss). WISH samples were imaged under transmitted light, with an Axiocam 506 colour camera (Zeiss). RNAscope samples were imaged with dsRed, YFP and DAPI filters (supplied by Zeiss), using the Orca Flash 4.0 V2 camera (Zeiss) and Apotome.2 (Zeiss). Images were processed using the ZEN software (Zeiss), the FIJI package and Adobe Photoshop CS6.

We note that mutant and WT embryos subjected to WISH were morphologically indistinguishable after fixation and were usually processed together. The Pearson’s chi-squared test [38–40] was used to test the null hypothesis that in a sample of mixed WT, heterozygous and homozygous mutant embryos, observed alternative gene expression patterns correspond to the expected Mendelian ratios (75% of embryos are expected to show WT phenotype and 25% to potentially show altered gene expression). For this test, degrees of freedom = 1. The chi-squared table [41] was used to calculate the probability that the number of observed embryos with an alternative expression phenotype was consistent with the expected number of homozygous mutants. For p-value > 0.1 the null hypothesis was accepted. For p-value < 0.1 it was assumed that alternative phenotypes in our samples were due to effects independent of the mutant genotype.

RNAscope results were derived from two independent experimental repeats. From each repeat, between 3 and 5 embryos were examined. For each experiment, the number of embryos (2 or more representative individuals for each stage) used to score cells is indicated in S3 Table.

Overexpression by microinjection

For overexpression assays, 50–70 pg of purified mRNA diluted in sterile water were injected in each WT (WIK) one-cell stage embryo using standard methods [19]. Capped mRNA was prepared from plasmid templates using the SP6 mMessage mMachine kit (Ambion) for the overexpression of GFP, Sox10^WT/Sox10^m618 [26] and Mitfa^WT/Mitfa^b692 [19], [42]. For Tfec^WT/Tfec^ba6 overexpression, in vitro capped and polyadenylated mRNA was prepared using the mMessage mMachine T7 Ultra transcription kit (Ambion). Total RNA was isolated using TRI reagent (Sigma) from dissected trunks of 10–15 72 hpf WT (WIK) or homozygous tfec^ba6/ba6 embryos and WT or mutant cDNA was generated using the SuperScript III First Strand Synthesis Supermix kit (Invitrogen). The tfec coding sequence (ENSDART00000164766.1) was amplified from the cDNA templates using the following primers: forward 5’-AGCGAGATCCTCCTGCTTCG-3’, reverse 5’-ATTCTGAGAGTGCGGTCCAG-3’. The T7 promoter was fused to the 5’ end of the amplicons through additional PCR amplification using the same reverse primer and the following forward: 5’-TAATACGACTCACTATAGGGAGAAGCGAGATCCTCCTGCTTCG-3’. The resulting amplicons were used as templates for in vitro transcription.

Quantitative real-time PCR

Total RNA was extracted using TRI reagent (Sigma) from 8 embryos per sample at 6 hours post-injection. cDNA was synthesised from 1 μg of total RNA using the iScript cDNA synthesis kit (Biorad; Cat# 1708890). qRT-PCR was performed in duplicate using Fast SYBR Green Master Mix (Applied Biosystems; Cat# 4385617) and the StepOnePlus Real-Time PCR System (Thermo Scientific; Cat# 4376600). For normalisation, we used expression of the housekeeping gene rlp13 (primers ready-made from Primerdesign Ltd.). Primers for tfec transcript detection: forward 5’-GGAGCTTGGATTGCATGGAG-3’, reverse 5’-TTGATCAGCACCGTACACCT-3’. Primers for pnp4a transcript detection: forward 5’-TGGATGCAGTTGGAATGAGT-3’, reverse 5’-TTGACAGTCTCGTTGTCCTCA-3’. The ΔΔCt method [43] was used to evaluate relative changes in pnp4a expression, whereas absolute levels of tfec transcript were assessed using a standard curve. The null hypothesis that there was no change in the level of gene expression between control samples (injected with GFP or with null transcripts) and overexpression samples was rejected if p-value < 0.05 using a two-sample t-test without assuming equal variances. Unpaired, two-tailed t-tests were performed using Microsoft Excel.

Mathematical modelling

Modelling of gene interactions was done using the approach presented in Greenhill et al. (2011). For each model, gene expression dynamics were described using a system of ordinary differential equations (ODEs; see S1 Text). The equations were solved numerically using the ode45 solver in MatLab software. Solving these equations returns the average concentration of gene output, measured in nM, across a homogeneous cell population. The Monte Carlo sampling algorithm used for randomising the constant parameters and subsequently scoring model outputs was run on MatLab. By random uniform logarithmic draw, we let all parameters vary in the range between a multiple 1/3.5 and 3.5. For each random draw, the system of ODEs of interest is solved and the resulting gene output dynamics are scored for biological relevance. For the scoring criteria refer to the results section, and for the mathematical functions used to calculate the scoring measure (S) for each system of equations, with randomly assigned parameter sets, see S2 Text.

Identification of the iridophore lineage throughout zebrafish embryogenesis

To interpret the expression dynamics of genes of interest during iridophore development in wild-type (WT) embryos, as well as changes of expression patterns in different loss of function contexts, it was crucial to distinguish cell populations comprising the different stages of iridophore development. Gene expression in the zebrafish NC is highly dynamic, reflecting both the rapid fate specification and differentiation of NC derived lineages in zebrafish and the multiphasic expression patterns of many key genes. For example, in previous studies of the ltk marker, we have proposed at least three phases of expression, one in multipotent premigratory progenitors, and two representing iridoblasts and differentiated iridophores respectively [9]. Tfec was first identified as a Mitf-related bHLH-ZIP transcription factor expressed in premigratory NC and later in a pattern reminiscent of iridophores [29]. We have recently shown that Tfec is crucial for fate specification of the iridophore lineage, and that tfec expression labels both early iridoblasts and differentiated iridophores (Petratou et al, in prep.). Building on these previous studies, we assessed the spatio-temporal locations of presumed iridophore progenitors at key stages of embryogenesis by examining tfec expression in whole-mount embryos. WISH on single embryos at 72 hpf confirmed that tfec is a definitive marker of differentiated iridophores (Fig 1A), similar to the established iridophore lineage marker, ltk [9]. Moreover, double labelling of ltk and tfec expression using multiplexed fluorescent RNAscope revealed that tfec is expressed throughout iridophore development (Fig 1C and 1C’; S3 Table). Furthermore, tfec transcripts were first seen in premigratory NC at very early stages, considerably before ltk (Fig 1B). Consequently, we interpret tfec expression as an excellent marker of ‘iridophore potential’ during zebrafish embryogenesis and use it here to produce a working classification of the stages of iridophore development. We note explicitly that this classification is intended to provide a framework for interpretation of mutant phenotypes, and that assessment of tfec alone cannot provide insight into the multipotency of cells at any specific stage.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Detection of developing iridophores using expression of tfec.

(A) Chromogenic WISH identifies tfec transcription in positions matching iridophore (iph) positions at 72 hpf. Images show the posterior trunk of the same individual, live and post-WISH processing. (B) tfec is expressed in sox10-positive NCCs (region in brackets) at 18 hpf, which have downregulated foxd3 but have not yet detectably activated early pigment lineage markers such as ltk, mitfa and pnp4a. (C) RNAscope reveals co-expression of tfec and ltk (arrows) in the anterior region of the posterior trunk (ARPT) during the timecourse of iridoblast specification, both in medially migrating cells at 24 hpf, and in dorsally located specified iridoblasts (ib(sp)) and iph at 30 hpf and 48 hpf, respectively. (C’) Nearly 100% of tfec+ cells co-express ltk during iridophore development (see S3 Table). Error bars indicate standard deviation between embryos. p-values were derived using a two-tailed, unpaired t-test. Using WISH at 24 hpf (D,E) tfec transcription is detectable in most or all cells of the multipotent premigratory NC domain (D, red line) and in a subset of cells of the posterior trunk located dorso-laterally to the spinal cord (D, vertical arrow) and in the medial pathway, between the somites and the notochord (D, horizontal arrow; Ei, arrows). Expression is only more weakly detectable in cells on the lateral pathway along the ARPT, between epidermal keratinocytes and somites (Eii, arrow). At 30 hpf (F,G) tfec is expressed in premigratory NCCs of the tail (red line) and in expected iridoblast positions, specifically dorsally located and medially migrating cells of the posterior trunk, the developing lateral patches and the eye (arrows). At 48 hpf (H) tfec is expressed along the dorsal, ventral and yolk sac stripes (vertical arrows), as well as in the lateral patches (arrowhead) and overlying the eye (horizontal arrow), in a pattern distinctive of differentiated iridophores. (I) RNAscope indicates presence of mitfa transcript in tfec-positive cells (arrows) migrating along the medial pathway in the ARPT at 24 hpf, but not in those located dorsally at 30 hpf, or at 48 hpf. In all three stages, mitfa+;tfec- cells are detectable (asterisks). (I’) Mean percentages of tfec+ cells co-expressing mitfa at each stage of development (refer to S3 Table). Error bars indicate the standard deviation between embryos. A two-tailed, unpaired t-test was used to derive indicated p-values. All panels show lateral views, except dorsal views in (A; D,F insets). Head towards the left. RNAscope panels: single focal planes shown. e, epidermis; ICM, intermediate cell mass; LPs, lateral patches; no, notochord; RL, reflected light. Scale bars: (A,B,D,E,F,G,H) 50 μm; (C,I) 20 μm.

https://doi.org/10.1371/journal.pgen.1007402.g001

At 18 hours post fertilisation (hpf), premigratory NCCs reside along the dorsal trunk and tail (Fig 1B). Towards the posterior tail, these precursors are characterised by WISH as expressing markers such as sox9b, sox10, snai1b and foxd3 [9], [26], [44] and likely correspond to fully multipotent early NCCs (eNCCs). At the same stage, more anteriorly (i.e. posterior trunk and rest of tail), sox9b, snai1b and foxd3 are downregulated, while sox10 is retained. Interestingly, tfec is expressed in premigratory NC cells throughout the trunk and anterior tail (Fig 1B), in a manner similar to sox10, even though fate-mapping of premigratory NCCs shows that only a relatively small subset of NCCs will generate iridoblasts [26]. At this stage neither the melanoblast marker mitfa nor two other early iridoblast markers ltk and pnp4a were detectable by WISH in tfec-positive NCCs of the trunk (Fig 1B) [9], [44]. However, mitfa and ltk are activated widely by 22 hpf [9], [19], with pnp4a following soon after [15]. We consider that these premigratory cells of the Anterior Region of the Posterior Trunk (ARPT; Fig 1B) which express tfec, but no longer foxd3, and which only later detectably upregulate other pigment markers, are multipotent iridophore progenitors corresponding to the proposed partially restricted pigment cell progenitor [9]. Given the spatiotemporal gradient of development that is so pronounced during the stages of NC development, throughout this paper we will largely focus on a readily-defined anatomical zone, the ARPT, when considering expression patterns at different stages, thus minimising the heterogeneity of the examined population of cells. The ARPT lies above the anterior yolk sac extension (YSE) (approximately the region of somites 9–11; bracketed in Fig 1B).

At 24 hpf, NCCs of the trunk have entered the medial and lateral migratory pathways, whereas less developed NCCs of the tail remain in premigratory positions. Although chromogenic WISH revealed maintenance of strong tfec expression in the premigratory NC domain, the majority of trunk NCCs located dorsal to the neural tube in the ARPT downregulate tfec, presumably due to cells becoming specified towards alternative lineages (Fig 1D). In this region, prominent tfec expression was retained in a small subset of precursors scattered over the spinal cord (Fig 1D), as well as in presumed iridoblasts migrating through the medial pathway (Fig 1Ei). tfec was only weakly detectable in cells entering the lateral migratory pathway, consistent with the medial migration pathway bias for iridoblasts noted previously [45] (Fig 1Eii). To further characterise the tfec-positive medially migrating cells, we used multiplexed fluorescent RNAscope to determine co-expression of tfec with the melanocyte lineage marker, mitfa (Fig 1I and 1I’; S3 Table). In all stages we scored dorsally located and medially migrating tfec+ cells of the trunk for mitfa expression. Importantly, the majority of tfec-expressing cells on the medial migration pathway were found to co-express mitfa, leading us to distinguish these as iridophore progenitors which retain at least bipotency. In the context of this paper, we will designate such tfec-expressing cells on the medial pathway as specified iridoblasts (but equally they could also be considered specified melanoblasts, using the definition of [46]). Numerous cells on the medial pathway were positive for mitfa, but displayed very weak or completely lacked expression of tfec (Fig 1I). We interpret these as having lost, or being on the way towards losing, iridophore potential, and are likely definitive melanoblasts. At 24 hpf, for quantitation cells were only scored if the tfec expression signal was elevated (>3 spots of fluorescence surrounding the nuclei) compared to the widespread low level (1–2 spots of fluorescence) expression displayed by numerous cells of the premigratory domain. This threshold was set in order to minimise inclusion of NC-derived cells which were in the process of downregulating tfec expression, thus focusing our analysis on cells most likely to be adopting the iridophore fate. At this stage, 10.3% of the tfec+ cells scored lacked mitfa. Although conceivable that this subpopulation might reflect an alternative pathway to generate iridophores, one that does not require an mitfa+ progenitor, we suspect this small number of cells simply reflects technical limitations of the RNAscope technique.

By 30 hpf, tfec transcript is present in scattered non-melanised cells along the dorsal posterior trunk and in medially migrating cells along the posterior trunk and anterior tail. Bilaterally patterned tfec+ premigratory NCCs were only detectable in the posterior tail (Fig 1F). In the APRT, co-expression analyses via RNAscope revealed that expression of tfec and mitfa has now resolved to be non-overlapping in the majority of cells (76.9%) (Fig 1I and 1I’; S3 Table), which corresponds to a statistically significant drop in the proportion of cells co-expressing the two genes. We distinguish these cells as definitive iridoblasts, whether in the dorsal position characteristic of differentiated iridophores at later stages, or migrating on the medial migration pathway. Finally, we assessed tfec expression at 48 hpf, a stage when in live embryos differentiating, light-reflecting iridophores are distinguishable, interspersed along the dorsal, ventral and yolk sac stripes, occupying the lateral patches and overlying the eye. tfec expression was detectable by chromogenic WISH in all of these positions (Fig 1H). Co-expression analyses using RNAscope showed discrete expression of tfec and mitfa (Fig 1I and 1I’; S3 Table), confirming that tfec expression at this stage definitively marked differentiating iridophores. The small percentage (16.6%) of tfec+;mitfa+ cells is most likely an artefact, due to the very close proximity and overlap between iridophores and melanophores along the dorsal stripe. This can result in expression from one cell being falsely assigned to an adjacent one.

In summary, in addition to the fully multipotent early NCC state (eNCC), these marker studies clearly distinguished four sequential phases of iridophore development in WT zebrafish embryos: 1) premigratory NCCs presenting with widespread expression of tfec, but which have downregulated eNCC markers (for example, foxd3), interpreted as broadly multipotent pigment cell progenitors (chromatoblasts, Cbl); 2) scattered cells strongly maintaining tfec expression dorsal to spinal cord and on migration pathways, but also expressing mitfa and so interpreted as at least bipotent iridoblast progenitors (specified iridoblasts, ib(sp)); 3) scattered undifferentiated cells in iridophore positions, showing rounded morphology and prominent expression of iridophore markers, but not mitfa (definitive iridoblasts, ib(df)); and 4) discrete iridophore marker-expressing (and reflective in live fish) cells in characteristic definitive iridophore pattern (mature iridophores, iph). We note that ib(sp) and ib(df) can only be distinguished in double WISH, and so will not be strictly distinguishable in many experiments, although from the above discussion it can be seen that at 30 hpf most such cells in the ARPT will be ib(df). This characterisation provides a framework for assessment of expression patterns of other genes and for the interpretation of expression patterns seen in mutant embryos.

sox10 expression is maintained throughout development of the iridophore lineage

Although well-known as a key factor in iridophore specification and a key marker of multipotent NCCs [9], [26], [47], sox10 expression has yet to be characterised in the iridophore lineage. We used both WISH and RNAscope to investigate the transcriptional dynamics of sox10 during iridophore development. We imaged iridophores of live embryos at 72 hpf using reflected light, and subsequently detected sox10 transcript in individual fish using chromogenic WISH. sox10 expression was readily detected in all iridophores (e.g. in each of the dorsal, ventral and yolk sac stripes; Fig 2A–2D). We then employed RNAscope to assess whether sox10 expression was maintained throughout all stages of iridophore specification, as opposed to its becoming re-activated in differentiated cells. We found that cells expressing the iridophore lineage marker, ltk, consistently co-expressed sox10 at each of 24 hpf, 30 hpf and 48 hpf (Fig 2E–2T). Therefore, just like tfec, sox10 expression in premigratory multipotent NCCs (Fig 1B) is maintained throughout fate restriction to ib(df) and their subsequent differentiation as iridophores.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. sox10 expression is maintained throughout iridophore development.

(A, B) lateral views of the anterior tail of a single embryo at 72 hpf, imaged live under reflected light (RL) (A) and then post-WISH to detect sox10 transcript (B). (C, D) dorsal views of the trunk and anterior tail of a second individual pre- (C) and post- (D) WISH processing for sox10 expression at 72 hpf. Differently coloured arrowheads point to individual iridophores expressing sox10. sox10 is also detected in developing oligodendrocytes (B, arrow), Schwann cells (B) and in iridophores along the yolk sac stripe (A, B asterisk). ltk-positive cells detected via RNAscope (E, J, O arrowheads) all show sox10 transcript (F, K, P; H, M, R arrowheads), at each of 24 (E-I), 30 (J-N) and 48 hpf (O-S). At 24 hpf, cells on the medial migration pathway are shown (I, boxed region). At 30 hpf and at 48 hpf, cells along the developing dorsal stripe are presented (N, S, boxed regions). (T) Quantification of the proportion of ltk+ cells co-expressing (green), or not co-expressing (red) sox10 by RNAscope, at 24, 30 and 48 hpf. Error bars indicate the corresponding standard deviations. (E-S): lateral views of single focal planes. (A-S): heads positioned towards the left. Sc, Schwann cells; pLLn, posterior lateral line nerve; no, notochord; RL, reflected light. Scale bars: (A-D) 100 μm; (E-H, J-M, O-R) 20 μm; (I, N, S) 50 μm.

https://doi.org/10.1371/journal.pgen.1007402.g002

sox10 maintains tfec expression in iridoblasts undergoing specification

Previous studies with ltk have concluded that iridophore specification fails and that pigment cell progenitors are trapped in a multipotent progenitor state (Cbl) in sox10 mutants [9], [26]. We re-assessed this proposed role of sox10 using both loss and gain of function assays. In sox10 mutants, both at 24 hpf and at 30 hpf (Petratou et al., in prep; Fig 3A and 3B), tfec expression is prominently retained in the premigratory NCC domain, which unlike in WT siblings extends anteriorly throughout the embryo. These observations indicate that Sox10 function is not required for establishment of the tfec-positive multipotent progenitor. Importantly, tfec transcripts were undetectable in ventrally migrating iridophore progenitors, indicating a requirement for sox10 to maintain tfec expression in a subset of cells (ib(sp) and ib(df)). By 48 hpf, we could not detect tfec expression in sox10 mutant embryos (Fig 3E and 3F), consistent with apoptotic elimination of NC derivatives which fail to become specified, including progenitors for all chromatophore lineages [26]. To test the sufficiency of Sox10 for expression of tfec, we overexpressed WT Sox10 or a null mutant version [22], [26] by injection of mRNA into single cell stage WT embryos and assayed absolute tfec transcript levels using qRT-PCR at 6 hours post-injection. Our data clearly showed that functional Sox10, but not the null version of Sox10, ectopically activated expression of endogenous tfec (Fig 3O). Together these data show that Sox10 is not essential for initial activation of tfec in early NCCs, but it is required for maintenance of tfec expression as multipotent progenitors become specified towards an iridophore fate, i.e. for iridoblast fate specification.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. A sox10-dependent tfec/ltk positive feedback loop is required for iridophore specification.

WISH to detect tfec (A-H) and ltk (I-L) expression at 30 hpf (A-D, I, J) and 48 hpf (E-H, K, L) in WT and mutant embryos, with quantitation (M,N). At 30 hpf, sox10 mutants (B) lack tfec expression in medially migrating iridoblasts (A, arrows) and instead show a striking anteriorly expanded multipotent progenitor domain (B, arrowheads). Both ltk (C) and tfec mutants (D) display a reduced number of medially migrating ib(sp) (horizontal arrows) at 30 hpf, while dorsally located ib(sp) of the posterior trunk and tail (vertical arrows) are only significantly reduced in ltk mutants. At 48 hpf, sox10 (F) and ltk (G) mutants lack tfec-labelling in the position of WT iph (E, asterisks), while tfec mutants (H) display a reduced number of cells in these positions expressing tfec (H, asterisks; N). At both 30 hpf (I, J) and 48 hpf (K, L), tfec mutants lack ltk expression in ib(df) (I, J, arrows) and in iph (K, L, asterisks) locations, with the exception of rare escaper cells. (M) Quantitation of tfec-expressing iridophore lineage cells at 30 hpf. Counts for tfec-positive ib(df) along the dorsal posterior trunk, the migration pathways, and the lateral patches are shown from left to right. (N) Quantitation of tfec-expressing cells along the posterior trunk and tail at 48 hpf. At this stage, scored cells are in iph positions. tfec mutants display almost complete lack of ltk-positive cells both along the dorsal and the ventral stripes compared to their siblings, while tfec positive cells in both regions are partially reduced. p-values given in (M,N) were produced using a two-tailed unpaired t-test to determine the statistical significance of differences between means. ns: p > 0.05. Box-and-whisker plots in M and N show mean, 1^st and 3^rd quartiles and range excluding outliers, which are indicated by dots. (O) Overexpression of WT sox10 mRNA results in an increased number of tfec transcripts, compared to overexpression of null sox10 mRNA. Graph represents 2 biological replicates, each with 3 technical replicates. Error bars indicate respective standard deviations. A-L) Lateral views, head positioned towards the left. Scale bars: 100 μm.

https://doi.org/10.1371/journal.pgen.1007402.g003

Tfec and Ltk generate a positive feedback loop

We next asked what roles Tfec and Ltk played in the iridophore GRN. Importantly, at both 18 and 24 hpf, we were unable to distinguish differences in tfec expression between ltk mutants and their WT siblings (S1 Table). Specifically, the premigratory NCC domain as well as specified iridoblasts in the posterior dorsal trunk and on the medial migration pathway of the trunk were unaffected in all examined embryos. Thus, tfec is activated in NCCs and is maintained at early stages of iridoblast specification, independently of Ltk activity. Nevertheless, from 30 hpf we observed a statistically significant decrease in the number of ib(df) located in the dorsal posterior trunk of ltk mutants (Fig 3A, 3C and 3M), and by 48 hpf no cells expressing tfec were identifiable in the ib(df) positions of the embryonic trunk in these mutants (Fig 3E and 3G).

Study of ltk expression in tfec mutants by chromogenic WISH suggests that Tfec function is crucial for ltk expression from the earliest stages onwards. At 24 hpf, approximately 25% of assessed embryos completely lacked ltk expression, with the exception of very rare escaper cells (Petratou et al., in prep.). This phenotype remained clearly identifiable at 30 hpf and persisted until at least 48 hpf (Fig 3I–3L). Moreover, examination of tfec expression in tfec mutant embryos, readily distinguishable owing to lack of melanin pigment in the RPE (Petratou et al., in prep; Fig 3A and 3D insets), revealed a subtle but consistent reduction in the numbers of tfec-positive ib(df) in mutants from 30 hpf (Fig 3A, 3D and 3M). Specifically, tfec mutants displayed a 35% and a 45% decrease in the number of tfec-expressing cells along the migratory pathways and in the developing lateral patches respectively, compared to WT siblings. Similarly, at 48 hpf, tfec mutant embryos (identified by the clear eye phenotype; Fig 3E and 3H insets) showed reductions in tfec expressing cells in the dorsal stripe and the ventral stripe to 58% and 45% of those in WT siblings (Fig 3E, 3H and 3N). Although the remaining tfec-expressing cells show a distribution consistent with their being iph, we note that they lack both ltk expression and visible reflective platelets and hence cannot correspond to ib(df), which are positive for ltk expression, nor to mature iridophores. We speculate that these cells represent an interesting state in which iridoblasts are trapped in a very early stage of their development, where tfec expression, but not other markers, continue to be maintained.

Taken together, our data strongly support the model that Tfec and Ltk function in a positive-feedback loop to maintain each other in specified iridoblasts, although tfec can be activated in this cell type independently of Ltk function.

pnp4a is temporally regulated by variable activators

The gene pnp4a has been defined as an iridophore lineage marker, although it is expressed rather widely in NCCs and long before iridophores differentiate [15]. RNA-seq analysis of gene expression in purified iridophores and melanocytes has shown that it is expressed at high levels not only in differentiated iridophores, but also, albeit at lower levels, in melanocytes [48]. We used chromogenic WISH studies to assess pnp4a expression in iridophore development and in various key mutants.

We first examined the WT expression pattern of pnp4a, and compared it to that of the iridophore and melanocyte lineage markers tfec and mitfa, respectively, at 24 hpf, 30 hpf and 48 hpf (Fig 4A–4I). At 24 hpf, it is notable that the expression pattern of pnp4a strikingly resembled that of mitfa, rather than that of tfec (Fig 4A, 4D and 4G). Specifically, we see clusters of cells just posterior to the otic vesicle and numerous cells on the medial migration pathway in the expression patterns of both mitfa and pnp4a, although both these regions are only sparsely positive for tfec (Fig 4A, 4D and 4G and insets). By 30 hpf, this pattern is still detectable, and indeed now pnp4a transcripts are detectable in differentiating melanocytes clustered behind the otic vesicle, as well as in melanised cells of the head, in a pattern similar to mitfa, but not tfec, expression (Fig 4B, 4E and 4H and insets). In addition, at this stage, the pattern of tfec and pnp4a expression in both the dorsal posterior trunk as well as overlying the RPE showed strong similarities; mitfa transcript is absent from the latter region (Fig 4B, 4E and 4H). By 48 hpf, the pattern of pnp4a strikingly resembled that of tfec, with both transcripts detected in iph positions, consistent with previously reported data [15], whereas mitfa was expressed in melanised cells of the head and of the dorsal, lateral and ventral stripes (Fig 4C, 4F and 4I). Considered together, our data suggest that while at later stages pnp4a is a definitive marker of differentiated iridophores, initially it is expressed widely in specified and differentiating melanoblasts.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. pnp4a is activated by mitfa in premigratory and migrating NC, and by cooperative action of sox10 and tfec in iridoblasts.

At 24 hpf and at 30 hpf, chromogenic WISH reveals strong similarities between the pattern of pnp4a and mitfa in migrating NC in the trunk and in cells posterior to the otic vesicle (A, B, D, E, black arrows), but which is distinct from tfec, which is present in more restricted groups of cells (G, H, blue arrows). At this stage, melanising cells in the head and anterior trunk show distinct expression of both mitfa and pnp4a (B,E, insets, black arrows). At 48 hpf, mitfa is expressed in melanised cells (F, black arrows), but pnp4a is not detectable in these melanocytes (C). From 30 hpf, some aspects of pnp4a expression are similar to those of tfec (B, H, blue arrows) and at 48 hpf both genes are expressed in ib(df) locations (C, I, blue arrows). (J-Q) Mutant analysis. At 30 hpf, mitfa (M) and tfec (L) mutants retain only a subset of the WT pnp4a expression; remaining cells in the former display a ib(df) pattern and in the latter a melanoblast pattern (M, L, arrows). tfec mutants lack pnp4a expression at 48 hpf, with the exception of rare escaper cells in iridophore positions (P, asterisks). Note that embryo in N, but not in O or P, was treated with PTU to inhibit melanisation; dark pigment in P is melanin. In sox10 mutants, pnp4a is largely absent, although weak expression persists in a few premigratory NCCs (K, arrowheads) at 30 hpf, and in rare escaper cells in iridophore positions at 48 hpf (O, asterisks). (Q) qRT-PCR measurement of pnp4a expression after expression of Mitfa or Tfec in early zebrafish embryos. Overexpression of WT Mitfa results in ectopic activation of pnp4a in injected embryos at 6 hours post-injection, whereas mutant Mitfa (null) does not. Interestingly, neither WT nor mutant Tfec is sufficient to drive pnp4a expression at this stage. Fold activation is calculated following normalisation to pnp4a levels upon overexpression of GFP. Graph is representative of 2 biological replicates, each with 3 technical replicates. Error bars indicate respective standard error of the mean. p-values indicate the significance of mean fold change for each sample when compared to the mean of GFP, using a two-tailed, unpaired t-test. Lateral views, head towards the left. Scale bars: 100 μm. Inset scale bars: 50 μm.

https://doi.org/10.1371/journal.pgen.1007402.g004

This suggested that pnp4a expression might be regulated by both Tfec and Mitfa. We began by investigating pnp4a expression in tfec mutants and WT siblings. At 30 hpf, WTs showed prominent pnp4a expression along the dorsal and ventral posterior trunk and the migratory pathways across the trunk and tail (Fig 4J). In contrast, tfec mutants displayed partial loss of pnp4a-positive cells (Fig 4L). Specifically, compared to WT siblings, tfec mutants showed decreased numbers of cells (expressing relatively low levels of pnp4a) along the dorsal trunk, and have comparatively few cells both on the migration pathway and in the ventral trunk, mostly more anterior. This partial reduction was also observed at 24 hpf (Petratou et al., in prep.), and principally affected cells in the ventral trunk and premigratory NC. By 48 hpf, pnp4a expression in iph locations was eliminated in tfec mutants (Fig 4N and 4P; Fig 3N). Thus, pnp4a expression in iridophores and in ib(sp) is dependent upon Tfec, whereas the persistence of pnp4a expression in a subset of developing NC derivatives until 30 hpf suggested that pnp4a expression also depends on additional inputs.

Due to the striking similarity of their expression patterns, we investigated a possible interaction between mitfa and pnp4a. Interestingly, pnp4a-expression was nearly eliminated in homozygous mitfa mutants, compared to their WT siblings at both 24 hpf (S1A–S1B’ Fig) and at 30 hpf (Fig 4J and 4M). Whereas in the former stage very few pnp4a expressing cells persisted along the posterior trunk of homozygous mutants (S1A–S1B’ Fig), in the latter a distinct group of dorsally located cells patterned in a ib(df)-like manner along the posterior trunk and anterior tail region were retained. In contrast, medially migrating cells were almost absent (Fig 4M). These results indicated that mitfa is an important regulator of pnp4a in premigratory and migrating NC cells, but that pnp4a activation in ib(df) is not affected. Thus, pnp4a expression appears to switch from Mitfa to Tfec-dependency during the transition from Cbl to ib(df), and to be detectable transiently in all melanoblasts and early differentiating melanocytes.

We then asked whether Mitfa or Tfec were alone sufficient to drive pnp4a expression. We overexpressed each transcription factor, or a null mutant variant, in 1-cell stage WT embryos and measured pnp4a expression at 6 hours post-injection by qRT-PCR (Fig 4Q). As a negative control, we injected GFP RNA, allowing us to measure the relative (fold) change of expression between samples injected with GFP mRNA, compared to those injected with RNAs encoding WT or mutant. As expected, neither mutant Mitfa nor mutant Tfec altered pnp4a transcript levels, compared to overexpression of GFP. Interestingly, introducing WT Mitfa led to a statistically significant 4-fold increase (Fig 4Q). Surprisingly, however, overexpression of WT Tfec did not result in a statistically significant ectopic activation of pnp4a (Fig 4Q), suggesting that Mitfa, but not Tfec, is sufficient in this ectopic context to upregulate pnp4a.

Loss of function studies using sox10 mutant embryos revealed that pnp4a was completely absent (Fig 4K), and that the gene remained inactive throughout the investigated developmental time-course (Fig 4O). We concluded that sox10 function was directly or indirectly required for all aspects of pnp4a expression, including the tfec-dependent pnp4a upregulation to occur. More broadly, we conclude that pnp4a regulation is more complex than previously assumed and that it should not be considered a definitive marker of the iridophore lineage at early stages.

A preliminary GRN governing iridophore development

Bringing the above interactions together, we propose a preliminary iridophore GRN, comprising model A (Fig 5A–5C). We use solid lines to describe known direct interactions and dashed lines to indicate interactions where their nature is unknown. Sox10 has been previously shown to bind directly to the promoter of mitfa in zebrafish, and to activate its expression [21], we thus include that interaction. Furthermore, Ltk relies on intracellular cascades and effector transcription factors to activate gene expression, thus its input is always indirect. However, for the remainder of the interactions it remains unclear whether Sox10, Tfec or Mitfa bind directly to the promoters of downstream genes.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Mathematical modelling and refinement of the preliminary iridophore GRN.

Graphical representation, simulation outputs and PCA of Monte Carlo results for models A1 (A,D,G), A2 (B,E,H) and A3 (C,F,I)). In the schematics (A-C), green rectangular nodes are used for transcription factor-coding genes, blue rectangles for transmembrane receptors and red circles for enzymes. Transcriptional activation is indicated by red arrows. Solid edges represent known direct interactions, whereas dashed lines indicate known indirect or interactions where their nature is unknown. These diagrams are mathematically described using ODEs, the numerical solutions of which indicate how the gene expression dynamics progress in WT and mutant scenarios (D-F). In the simulation outputs (D-F), lines and arrowheads (where lines are overlapping), coloured according to the key in (D), indicate the change of concentration (nM) for each gene output during the course of iridophore development (hpf). The simulation considers the developing iridophore lineage starting from tfec+;sox10+;ltk-;mitfa-;pnp4a- NCCs of the ARPT at 18 hpf. Qualitative assessment of simulation outputs is indicated as follows for each gene. Green tick: modelling predictions match experimental observations. Crosses: prediction for a specific gene (indicated by colour of the cross, matching gene reference colour as indicated in the legend) is refuted by experimental evidence. (G-I) PCA plots for models A1-A3. Axes plot principal components (PC) 1, 2 and 3. All trials using model A1 fail the set criteria, thus scoring between 0–0.2 (G, turquoise spots). Under model A2, a small proportion of parameter sets scores within the range of 25–55 (H, dark blue and purple spots) and no trials score above S = 55 (magenta spots). Under model A3, performance improves notably, with more trials scoring between 25 and 55 (I, dark blue and purple spots), but still no trials score above S = 55 (magenta spots).

https://doi.org/10.1371/journal.pgen.1007402.g005

We observe that in tfec and ltk mutants, sox10, tfec and ltk are all absent from the iridophore lineage from the ib(sp) stage onwards. Due to the nature of the sox10-dependent Tfec/Ltk positive feedback loop, and its central role in gene regulation within the GRN, loss of function experiments similar to those used to derive the aforementioned interactions cannot clearly indicate which of the candidate genes (Tfec, Ltk or Sox10 itself) is responsible for sox10 maintenance in iridophores. We propose three distinct variants of Model A, distinguished by the mechanism of sox10 maintenance in iridoblasts and iridophores (Models A1, A2 and A3 (Fig 5A–5C)). In A1, sox10 maintenance occurs through an autoregulatory positive feedback loop (Fig 5A), which maintains sox10 expression from premigratory progenitors (Fig 1B). In A2, the presumed iridophore master regulator, tfec, is directly or indirectly responsible for sox10 activation in the context of iridoblasts (Fig 5B). Finally, model A3 proposes that sox10 maintenance is dependent upon Ltk signalling, independent of Ltk’s action in maintenance of Tfec expression (Fig 5C).

Mathematical exploration of the preliminary iridophore GRN

These three Model A variants all share positive feedback loops between Sox10, Tfec and Ltk, and whilst biologically distinct, it is not obvious intuitively how they could be distinguished without detailed investigation of transcriptional regulatory mechanisms. However, like others, we have previously demonstrated the value of simple predictive mathematical modelling of GRNs in developing a robust understanding of their biological implications [22]. Thus, we utilised mathematical modelling of each model A variant, to assess more rigorously whether they could be distinguished. If so, we wished to identify the model offering highest predictive power, i.e. the network that was best able to recapitulate the experimentally observed gene expression dynamics in the iridophore lineage (S3 Fig).

We generated systems of ordinary differential equations (ODEs) describing the interactions in each of the proposed networks (see S1 Text). The changes in the expression of each gene over time were determined using an ODE which incorporated all activatory and repressive influences from other members of the network, as well as a term for degradation of the gene’s protein product. The model aimed to capture the average output (nominally as protein product, assuming direct correlation with transcript production) of each gene in a homogeneous group of cells at a given time. The necessary parameter values characterizing the regulatory dynamics (mRNA maximum production rates (g), protein degradation rates (d), dissociation constants for transcription factors binding (K)) were chosen following exploration of existing literature to identify physiologically relevant values (see S1 Text; S2 Table).

We solved the systems of ODEs numerically in MatLab. To define initial conditions (here at t = 18 hours), we chose the population of premigratory NCCs occupying the dorsal ARPT at 18 hpf. Using chromogenic WISH, sox10 and tfec, but not ltk, mitfa or pnp4a transcripts were detectable in this population of cells (Fig 1B), allowing us to approximate the initial conditions for our simulations. As time proceeds, the simulations were tested for their ability to broadly replicate the changes in gene expression of tfec-expressing cells as they transition through the stages of ib(sp), ib(df) and then iph. In vivo, differentiating iridophores are observable from 42 hpf and prominent by 48 hpf. To account for inaccuracies in our default parameter sets (see S1 Text; S2 Table), we allowed computational simulations to progress until 60 hpf, thus helping to ensure that any biologically meaningful steady state could be successfully reached. Specifically, in the WT context, it was crucial that Sox10 and Tfec stay upregulated in mature pigment cells, i.e. reach a positive steady-state. Similarly, Ltk and Pnp4a concentrations should increase and reach a positive steady-state. Mitfa levels should initially rise rapidly, reflecting the widespread expression of mitfa in all pigment cell progenitors [32], but should then drop to a distinctly lower level at differentiation stages (this work). We note that the Mitfa concentration is not required to attain zero, but simply to drop to a lower steady-state value; given the expectation that in situ hybridisation techniques have a ‘detection threshold’, we consider that this final lower value would reflect expression levels undetectable by our detection methods in differentiated iridophores, although they would still be measurable by microarray in pooled isolated iridophores [48]. As a further test of each of Models A1, A2 and A3, we used MatLab to predict gene expression changes in the context of different mutant scenarios, when function of Sox10, Ltk or Tfec were individually ablated in silico (Fig 5). In the sox10 mutant context we asked that Sox10, Tfec and Ltk acquire positive values, as expression has been identified in trapped chromatoblasts ([9], [26] and this work), however at no point do either Mitfa or Pnp4a become upregulated. Loss of Ltk function was required to predict initial rise of Ltk, Tfec and Pnp4a concentrations (at approximately 24–30 hpf), followed by gradual downregulation to undetectable levels. Similarly, loss of Tfec function should be accompanied by a peak and subsequent decline of Tfec, Sox10 and Pnp4a concentrations within 30–50 hpf. Ltk should never become detectable in the developing iridophore population in this context.

These simulations showed that regardless of the interaction underlying sox10 maintenance, in the WT context all iridophore markers were appropriately upregulated in the course of iridophore development, consistent with biological observations. In all three models, expression of the melanocyte marker, mitfa, a direct target of Sox10 [21], was predicted to be upregulated and then maintained in the iridophore lineage (Fig 5D–5F). These predictions are in contrast to previously published experimental data, showing that maintenance of mitfa is restricted to melanocytes [21], [22], although lower level mitfa expression has been detected by RNA-seq in differentiated iridophores [48]. We used RNAscope to assess directly the relative changes in mitfa expression in the iridophore lineage. Even with this technique, notably more sensitive than conventional chromogenic WISH, we confirmed that co-expression of mitfa with the iridoblast marker, tfec, occurred at 24 hpf, but such overlap was not detectable at 30 hpf and at 48 hpf (Fig 1I). Thus, all three models failed to correctly predict the expected initial peak, followed by downregulation, of mitfa expression in the iridophore lineage. These observations are readily explained by the absence of a mechanism for repression of melanocyte fate in our model; we explore this later.

For all three versions of model A, simulation of loss of Sox10 function appropriately predicted maintenance of tfec (Fig 3B) and of ltk [9], consistent with observations that tfec+;ltk+ progenitors remain trapped in the dorsal trunk and tail. Likewise, they appropriately predict the failure to upregulate both mitfa [21] and pnp4a (Fig 4K). However, it has been previously shown that dorsally trapped progenitors continue to express sox10 upon loss of Sox10 function [26], a feature predicted successfully by models A2 and A3, but not by A1. Similarly, computational implementation of tfec loss of function revealed that model A1 did not generate biologically accurate predictions, whereas models A2 and A3 performed better. Specifically, Model A1 with simulated loss of Tfec function did not result in the experimentally observed lack of ltk and gradual downregulation of both tfec and pnp4a expression (Fig 3D; Fig 4L). In models A2 and A3, ltk expression was correctly predicted to remain undetectable throughout iridophore specification and differentiation, while tfec and pnp4a were gradually diminished. Finally, in silico inhibition of Ltk signalling in models A1 and A2 failed to predict the experimentally observed initial activation, followed by downregulation, of ltk [9], tfec (Fig 3C) and pnp4a (S1C–S1F Fig) expression as iridoblasts differentiate into iridophores. Model A3, however, successfully predicted gradual elimination of iridophore marker gene expression in the lineage, more accurately reflecting the current experimental observations.

Based on the above observations, we conclude that Model A3 has the highest degree of predictive power using the default parameter set. These parameters were chosen based on ranges indicated from the literature as physiologically relevant, nevertheless the exact values were assigned somewhat arbitrarily. We, therefore, conducted an unbiased assessment of whether the experimentally set output requirements, as outlined above, could be achieved using alternative sets of parameter values in any of the models. To that effect, we designed a Monte Carlo algorithm able to randomly assign parameters drawn from a pre-assigned range, spanning two orders of magnitude from 5x lower than the physiological mean value to 5x higher than that value. For each model, the outputs for each of 20,000 combinations of parameters were scored computationally by a suitably designed scoring function, according to our set of qualitative criteria (see S2 Text), which took into account that only qualitative expectations of gene regulatory dynamics could be tested.

The scoring function for a given model output is the multiplication product of the individual scoring measures for each of the gene expression curves. These individual functions could be binary, adopting either 0 or 1 values, if the assessed feature is absolutely required for an output to be considered as biologically relevant (for example successful upregulation of iridogenic genes). Alternatively, individual scores may be quantitative (for example the score of the efficiency of the Mitfa ‘rise and drop’ behaviour), meaning different curve behaviours would result in relatively higher or lower score values (S). This feature implies that the highest score (S_max) achieved by each model may be used an indicator of its ability to produce outputs that closely match experimental observations. Importantly, S values act as relative ranking tools, comparing the capabilities of our models, but note that the exact values bear no biological significance. Furthermore, the frequency by which acceptable and high scores are achieved is crucial in identifying models robustly predicting the experimentally observed gene expression dynamics. In the principal component analysis (PCA) plots (Fig 5G–5I), the frequency of high scores in each model is visualised by the density of dark blue, purple and magenta spots. The three principal components depicted represent linear combinations, each pointing to the direction of maximal variance, with respect to their score-weighted position vectors.

PCA was used to visualise the frequency of scores within parameter space and to compare the three models’ respective capacities to reproduce experimental observations (Fig 5G–5I). Interestingly, all outputs derived from randomly assigning parameter values in the set of equations representing model A1 failed to predict crucial aspects of the biology, thus consistently achieving zero scores (Fig 5G). Model A2 was found able to predict those features correctly, although only limited subsets of parameters achieved admissible outputs (Fig 5H). Model A3 performed similarly to Model A2, except that it generated predictions broadly consistent with the known biology for a wider range of parameter combinations (Fig 5I).

Repression of mitfa in the iridophore lineage requires an unknown Factor R

Although the analysis of the model A variants identified a more favourable model (A3) for most aspects, none of the model A alternatives were able to reproduce the expected Mitfa dynamics (i.e. sufficient downregulation of Mitfa in differentiating iridophores), while simultaneously maintaining relatively high outputs of iridogenic gene products (Fig 5; S6A–S6C Fig). We used model A3 as a starting point to improve this aspect of the iridophore GRN. Since repression of mitfa in the iridophore lineage has not thus far been investigated, we asked which interactions would be able to produce appropriate outputs when mathematically implemented. After testing predictions of alternative models with our default parameter set (S4 Fig), we concluded that upregulation in the iridophore lineage of an unknown mitfa repressor, which we termed factor R, was crucial. The resulting Model B (Fig 6A) incorporated Tfec-dependent activation of factor R, which our implementations suggested should be absent in the multipotent progenitors of the ARPT at t = 18 hours. Manually adjusting the parameters in the system of ODEs describing Model B revealed that the experimentally determined rise and drop of mitfa expression in our group of cells could be achieved using the default parameter set (Fig 6B), and even enriched with alternative parameter values, within the determined physiologically relevant range (S6E and S6G Fig). Random assignment of parameters and algorithmic scoring of respective outputs (see S1 Text) indicated that, of all the tested models, model B best reflected experimental observations regarding gene expression dynamics. Specifically, compared to models A1-A3, a broader range of model B trials achieved high scores, with absolute values higher than those attainable through models A1-A3 (Fig 5; Fig 6). Notably, model B (derived from model A3) consistently scored higher than a designated alternative model B(2), which was derived by introducing factor R into model A2 (S6F and S6G Fig).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. Model B accurately reflects observed gene expression dynamics.

(A) Graphical representation of model B, where blue (blunt-ended) edges indicate transcriptional repression. (B) Simulation outputs. (C) PCA plot for model B. A large proportion of random parameter trials yields scores between 50–100 (dark, blue and purple spots), while several trials achieve scores from 100 to as high as 150 (purple/magenta spots).

https://doi.org/10.1371/journal.pgen.1007402.g006

We considered Foxd3 as a candidate for factor R, as the transcriptional regulator has been previously implicated in mitfa repression [15]. Our modelling predicted that, in the WT context, factor R should be expressed in undetectable levels in the ARPT at 18 hpf, and should then be robustly upregulated in all developing iridophores, before reaching a stable plateau at differentiation. We tested these predictions for foxd3 using RNAscope. Surprisingly, we detected only low levels of foxd3 transcript, and these in only half of ltk-positive iridophore lineage cells at each of 24 hpf, 30 hpf, 36 hpf and 48 hpf (S5 Fig; S3 Table), making it unlikely to be the key factor repressing mitfa expression in the iridophore lineage. As a further test, we examined foxd3 mutants and WT siblings by WISH at 24 hpf, to determine whether absence of Foxd3 activity resulted in derepression of mitfa in cells on the medial pathway (ib(sp) and later stages of iridophore differentiation). Contrary to this prediction, numbers of mitfa positive cells were found to be somewhat reduced in this region of homozygous mutants compared to WT siblings (S5 Fig). These observations were inconsistent with the hypothesis that Foxd3 mediated mitfa repression during iridophore lineage differentiation. Hence, we conclude that whilst previous data indicates a role for foxd3 in pigment cell development, it is unlikely to perform the role of factor R in iridophore GRN.