Both male and female H. latifasciatus in adult stage are characterized by four (in some strains five) vertical melanic bars. Individuals of our breeding stock had consistently four bars (n > 30; Figure 1E): one anterior bar above the operculum (vertical bar 1; vb1), two vertical bars in the trunk area that cover the whole dorso-ventral axis (vb2 and vb3) and a more posterior vertical bar (vb4) at the anterior caudal peduncle that only covers the dorsal part up to the horizontal myoseptum. The regions between the bars (referred to as interbars, ib) are yellow to beige with dominant males often having nuptial colors with stronger yellow but anteriorly also red to orange hues (Figures 1Q,R). In contrast to the closely related Lake Victoria cichlids bars are thicker and more pronounced in H. latifasciatus with a clearer demarcation, lower number and without as pronounced physiological color change (Greenwood, 1974; Seehausen et al., 1999).

To investigate the formation of vertical bars (Figure 1E), we described the development of H. latifasciatus larvae between 7 and 21 days post-fertilization (dpf; Figures 1A–D,F–P and Supplementary Figure S1). The developmental progression of vertical bar formation is fully consistent in all three individuals examined (Supplementary Figures S1, S2). At 7 dpf, melanophores are present in the dorsal head region as well as on the dorsal part of the yolk sac. In the trunk area only a few melanophores have formed without any obvious indication of a bar-like pattern (Figure 1A). Starting at 8 dpf vertical bars form in an anterio-posterior sequence (Figures 1B–D,F–P and Supplementary Figure S1). Trunk melanophore number has increased considerably within an anterior dorsal patch (dc1) forming at 8 dpf (data not shown), followed by a more posterior one at 9 dpf (dc2; Figure 1B). The melanophore patches anticipate the position of vb1 and vb2, respectively (Figures 1B,C). At 10 dpf a third and fourth dorsal patch (dc3 and dc4) are appearing at the positions where vb3 and vb4 will later form, respectively (Figures 1C,G and Supplementary Figures S1, S2). After the appearance of the dorsal patches they expand dorsally into the dorsal fin and ventrally forming the four bars. One exception is the posterior bar vb4, where the bar only extends to the horizontal myoseptum. The formation of vb3 is furthermore contributed by a second melanophore cluster (lateral cluster; lc) that forms in a more posterior-ventral position and merges with the developing bar between 12 and 13 dpf (Figures 1C,I,J and Supplementary Figures S1, S2). After this time, at around 2 weeks post-fertilization, the complete adult vertical pattern is already fully formed, but the contrast of the bars further increases until 21 dpf (Figures 1D,P and Supplementary Figure S1) and beyond.

In order to understand the formation of the characteristic bar pattern of H. latifasciatus, we analyzed the progression of the bar pattern formation over time. Specifically, we analyzed how the darkening of the bar regions is generated at the cellular level. We hypothesized that three processes might contribute to the darker appearance of the bar regions: (a) melanosome dispersal (the aggregation and dispersion of melanosomes, the melanin-containing organelles, within melanophores) as it increases the fraction of the tissue covered by melanin, (b) the density of melanophores, and (c) the darkness of the melanophores, i.e., the concentration of melanin.

To investigate the development and the importance of these factors we followed the development of three individuals over twelve days of development between 9 and 21 dpf focusing on the dorsal portion of the two bars in the trunk region, vb2 and vb3 and the yellow interbar region in between (ib2). To assess melanosome dispersal, we calculated the diameter of the melanin-covered part of all melanophores (see section Materials and Methods). We found no strong spatial difference in melanosome aggregation (Figure 2A and Supplementary Figure S1), yet dispersal diameters increased with age, likely because cells are still growing during these early developmental stages (Supplementary Figure S3). Cell density was evenly distributed at 9 dpf, but during the formation of the bars, cells became more densely packed in the bar regions, while cell density decreased in the interbar regions (Figure 2B and Supplementary Figure S1). To measure the darkness of individual melanophores we measured relative gray values (see section Materials and Methods). Here the difference between bar and interbar regions continuously increased suggesting stronger accumulation of melanin in bar melanophores (Figure 2C and Supplementary Figure S1). Indeed, closer observation of melanophores shows that interbar melanophores are often not fully filled with melanin, with parts of the chromatophore, often the center of the cell being poorly pigmented (Figures 1N–P and Supplementary Figure S1).

Next, we investigated the cellular behavior of single melanophores between 9 and 13 dpf. To do so, we photographed the vb2 region of the same individual on five consecutive days and tracked cellular migration and the formation of new melanophores (Figures 2D–M). The data suggests that melanophores mainly move indirectly through the general expansion of the skin. Newly differentiating melanophores could be found across the whole examined area, but at an increased rate in the forming bars. They grow to the full size (diameter: ∼0.06 mm) within 1–2 days. We found no evidence of proliferating melanin-containing cells.

In summary, these results suggest that the bars of H. latifasciatus form through spatial variation in melanophore properties (i.e., melanin content) and melanophore cell density that mainly arise in the second week after hatching (standard length 6–8 mm). The increase in melanophore density is most likely caused by an increased differentiation of melanoblasts within the bar regions.

To investigate the distribution of chromatophores in different integument regions, we estimated cell density and size of both, melanophores and xanthophores in the two interbars (ib1, ib2) and the two bars (vb2, vb3) of the trunk region (Figure 1E). Here, we only quantified pigment cells in female individuals as the vertical bar pattern does not differ between sexes, but substantial interindividual variation in red coloration of males would have complicated quantifications (Figures 1Q,R). Scales and dermis without scales (from hereon called “scales” and “skin”, respectively) were analyzed separately. Three types of chromatophores could be found in both the dark bars and the light interbars: melanophores with black/dark brown pigments, xanthophores with yellow to orange and reddish pigments (we did not differentiate between xanthophores and erythrophores; see discussion in section Materials and Methods), and iridophores with iridescent/reflective properties (Figures 3I,L). To quantify pigment cell quantity and properties we used three measurements: (a) pigment cell coverage, (b) pigment cell dispersion, and (c) pigment cell density.

Pigment cell coverage is influenced by both cell number and size (or intracellular dispersal of pigments in the cell) and was measured by estimating the percentual coverage of the tissue with pigment using light microscopy for melanophores and fluorescence microscopy to detect the autofluorescence of xanthophores. Consistent with the visual impression, melanophore coverage in bars was significantly higher than in the yellowish interbars in both scales (26.4% in bars, 2.0% in interbars) and skin (75.9% in bars, 13.7% in interbars) (Figures 3, 4A,B, Supplementary Figure S4, and Supplementary Table S1). Both melanophore density and melanosome dispersal contributed to the difference in melanophore coverage. The average melanosome dispersal diameter is larger in dark bars (Ø 0.058 mm in scale, Ø 0.082 mm in skin) than light interbars (Ø 0.022 mm in scale, Ø 0.033 mm in skin) (Figures 4C,D, Supplementary Figure S5, and Supplementary Table S1). However, variation within the same skin region was quite high suggesting that both dispersed and aggregated melanophores are widely distributed (Supplementary Figure S5 and Supplementary Table S1). The density of melanophores in both scale and skin were significantly higher in the bars (158 cells/mm² in scale, 333 cells/mm² in skin) than in interbars (52 cells/mm² in scale, 115 cells/mm² in skin) (Figures 4E,F, Supplementary Figure S6, and Supplementary Table S1).

Xanthophore coverage in bars and interbars is similar on scales (15.7% in bars, 18.4% in interbars), but significantly different in the skin (10.9% in bars, 50.4% in interbars) (Figures 3D,H,J,M, 4G,H, Supplementary Figure S4, and Supplementary Table S1). However, xanthophore coverage might be underestimated in the skin of the bar regions due to the high density of melanophores that could aggravate detection of the xanthophore autofluorescence. The dispersal diameter of xanthophores did not differ between interbars (Ø 0.048 mm in scale, Ø 0.025 mm in skin) and bars (Ø 0.044 mm in scale, Ø 0.021 mm in skin), neither in scales nor in skin preparations (Figures 4I,J, Supplementary Figure S5, and Supplementary Table S1). Cell densities were comparable in scales of interbars (157 cell/mm²) and bars (155 cell/mm²) (Figure 4K, Supplementary Figure S6, and Supplementary Table S1).

Next, we analyzed the molecular correlates of the observed differences in pigment cell density and pigment synthesis. Molecular markers for iridophore and xanthophore were also used as we could not analyze differences in these cells’ number due to the high melanophore density in the bar regions. Therefore, we compared expression levels of twelve candidate genes across the same two bar and interbar regions for the chromatophore measurement (ib1, vb2, ib2, vb3) using quantitative real-time PCR (qPCR; Figure 5 and Table 1). The selected genes are including marker genes for chromatophores, melanin synthesis genes and genes involved in the melanocortin signaling pathway. The latter was a particular focus as they have been previously implicated in color pattern formation of cichlids, teleosts and vertebrates more generally (Table 1).

Sox10 is a marker gene for chromatophore progenitors (Dutton et al., 2001). Expression was slightly higher in bar than in interbar regions (Figure 5A) and significantly differed between regions (ANOVA: P < 0.05). However, we found no significant sox10 expression variation between bars and interbars (Tukey HSD: P = 0.067–0.37; Supplementary Table S2). Expression of the melanophores marker, mitfa (Lister et al., 1999; Béjar et al., 2003) was higher within the dark vertical bars than in the adjacent light interbars (Figure 5B and Supplementary Table S2), yet differences were not significant (ANOVA: P = 0.19). Similar expression profiles can be also found in pmel, a melanophore specific gene important for melanin deposition in melanosomes (Schonthaler et al., 2005). Also here, we find that pmel is expressed at a higher yet not significantly higher (ANOVA: P = 0.056) level in dark bars compared to interbars (Figure 5C and Supplementary Table S2). The gene slc24a5, a melanosome-specific cation exchanger (Lamason et al., 2005), showed differential expression in some pair-wise comparisons of bar and interbar regions (ANOVA: P < 0.01; Tukey HSD: P = 0.002–0.18; Figure 5D and Supplementary Table S2). The two melanophore-specific genes that express melanogenic enzymes and are essential for the production of melanin, tyr (Hidehito et al., 1994; Camp and Lardelli, 2001) and tyrp1a (Braasch et al., 2009; Krauss et al., 2014) showed significantly higher expression levels in bars (ANOVA: both P < 0.001; Tukey HSD: all P < 0.01; Figures 5E,F and Supplementary Table S2).

In order to compare the distribution of iridophores and xanthophores, we used the iridophore lineage-specific marker ltk (Lopes et al., 2008) and the xanthophore marker csf1ra (Parichy and Turner, 2003). Notably, both ltk and csf1ra were expressed at similar levels across the differently pigmented bar and interbar regions (Figures 5G,H and Supplementary Table S2). This is in support of a rather homogenous distribution of iridophores and xanthophores across the trunk.

We also examined the gene expression of Agouti family genes (namely Asip/asip1 across all vertebrates and agrp2 in cichlid fishes) and melanocortin receptors, as they have been previously implicated in pigmentation. Two proteins of the Agouti family, Asip1 and Agrp2, likely act as the antagonists for the melanocortin receptors Mc1r and/or Mc5r in teleost skin (Cal et al., 2017). As previous studies suggest, asip1 is a key regulator of dorso-ventral countershading, presumably by regulating melanophore number (Ceinos et al., 2015; Cal et al., 2017, 2019). Interestingly we also found significant variation along the anterior-posterior axis (ANOVA: P < 0.01) with asip1 being expressed significantly higher in some pair-wise post hoc comparisons between bars and interbars (Tukey HSD: P = 0.0016–0.8222; Figure 5I and Supplementary Table S2).

Previous work in cichlids showed that agrp2 regulates presence/absence of stripe patterns while not contributing to shaping the pigment pattern itself through expression variation across the skin (Kratochwil et al., 2018). In contrast to results from the Lake Victoria species Pundamilia nyererei, we found significant differences between skin regions (ANOVA: P = 0.041; Figure 5J and Supplementary Table S2). However, expression mainly differed between anterior and posterior regions (Tukey HSD between ib1 and vb2: P < 0.05) and not consistently between bars and interbars (Tukey HSD: P = 0.046–0.895) (Supplementary Table S2).

Mc1r and Mc5r are two receptors antagonized by Agrp2 and/or Asip1 signaling and mc1r and mc5r have been shown to be expressed in chromatophores. The expression of mc1r has been reported in skin melanophores of zebrafish, barfin, and Japanese flounder (Kobayashi et al., 2010, 2012b; Higdon et al., 2013), xanthophores of goldfish (Kobayashi et al., 2011) and iridophores of zebrafish (Higdon et al., 2013). The expression of mc5r was detected in melanophores and xanthophores of flatfish (Kobayashi et al., 2010, 2012b). Here we found no significant expression differences for mc1r (ANOVA: P = 0.296) and mc5r (ANOVA: P = 0.764) which suggest that the melanocortin receptors may not contribute to shaping the vertical bar pattern (Figures 5K,L and Supplementary Table S2).

Most previous investigations of teleost pigment pattern formation focused on the horizontal stripe patterns of the model teleost Danio rerio, the zebrafish. In zebrafish, melanoblasts, the progenitors of melanophores, migrate from the dorsal neuroectodermal margin along nerve fibers between the myotomes (Dooley et al., 2013). After they settled at their final position, melanoblasts differentiate into melanophores and accumulate melanin to form the dark horizontal stripes that gave the zebrafish its name.

In the haplochromine cichlid H. latifasciatus [divergence time with zebrafish approximately 220 million years (Hughes et al., 2018)] color pattern formation starts with the development of four melanophore clusters that arise in the dorsal rim of the trunk (Figures 1B–D). The dorsal clusters initiate the formation of bar patterns later in development, by spreading dorsally into the dorsal fin and ventrally forming bars on the trunk. The third vertical bar (vb3) fuses with the lateral melanophore cluster (Figures 1C,D,I–P and Supplementary Figure S1). In contrast to the zebrafish where adult color pattern develops indirectly during a post-embryonic metamorphic phase (Parichy et al., 2009), the vertical bars in H. latifasciatus develop directly and are already visible at a time when the larvae are still feeding from their yolk. This is in line with previous reports on the direct development of other morphological traits of African cichlid fishes (Woltering et al., 2018).

By tracking the behavior of individual melanophores we could observe how new melanophores form over a period of 1–2 days (Figure 2 and Supplementary Figure S1). It is unclear how the melanoblasts reach this position. However, it is likely that they migrate dorsally from the neuroectodermal margin and then ventrally within the skin, possibly explaining the gradual expansion of bars from dorsal to ventral. This is interesting, because in zebrafish clones of single pigment cell progenitor cells have been shown to mainly spread dorso-ventrally (Singh and Nusslein-Volhard, 2015; Nüsslein-Volhard and Singh, 2017). It is possible that vertical bars are prepatterned at the progenitor level, potentially already at the level of the dorsal neuroectodermal ridge. Position and size of bars could be driven by anterior-posterior differences in progenitor number or transcriptional identity that later influence proliferation and differentiation within the dermis. It is also possible that migratory routes, environmental factors along the migratory pathways and in the skin as well as cell-cell interactions including reaction-diffusion systems contribute to or constitute the basis of this process. The anterior-posterior sequence of bar formation could be a consequence of the anterior-posterior sequence of somitogenesis. Anterior melanoblast would hereby first receive guidance cues from dermamyotome, sclerotome, and ermerging dermis. Anterior neural-crest cells including pigment cell precursors would therefore migrate first as previously described (Bronner and LeDouarin, 2012).

As soon as melanophores are visible in the skin they do not migrate anymore and grow to the full size within 1–2 days. Hereby the cell-size increases, new dendritic arbors form and melanin levels visibly increase (Figure 2 and Supplementary Figure S1). Therefore, the appearance of the bar pattern is mainly driven by an increase in cell density in the developing bar region as well as changing cellular characteristics including increased melanosome dispersal and melanin production (Figures 2, 6 and Supplementary Figure S1).

As one of the most common pigment patterns in haplochromine cichlids, vertical melanic bars can vary in contrast and number, both between as well as within species (Witte et al., 1976; Kocher et al., 1993; Seehausen et al., 1999). However, compared to most other bared haplochromine cichlids, the vertical bars in H. latifasciatus are fixed in number (mostly four and in some strains five) and are relatively thick. Our study shows that vertical bar number remains constant over the course of development in H. latifasciatus (Figures 1F–R, 2D–M and Supplementary Figure S1), which contrasts with a recent study on the Lake Malawi cichlid Copadichromis azureus, a species with unfixed bar number and thinner bars (Hendrick et al., 2019). In C. azureus, vertical bars split during development resulting in an increase of bar number (Hendrick et al., 2019). However, both processes are not driven by rearrangement of excisting pigment cells, but both broadening of bars in H. latifasciatus and increase in bar number in C. azureus is driven by the formation of new melanophores.

Comparisons of juvenile and adult patterns largely demonstrate that the juvenile patterns bear already most cellular and morphological characteristics of adult patterns. Vertical bars have a higher density of melanophores and melanophores are also evidently darker in the melanic regions (Figures 4E,F, 6). In contrast to juveniles, pigments are more dispersed in adult bar melanophores. Xanthophores show no clear differences, yet reliable quantification was only partially possible due to the high density of melanophores.

Gene expression analyses of known pigmentation genes gave complementary information as well as mechanistic insights on how the spatial variation of chromatophore properties and densities are archieved (Figures 5, 6). Clearly, melanophore-specific genes such as tyr, tyrp1a, but also slc24a5 are differentially expressed between bars and interbars (Figures 5D–F and Supplementary Table S2). All three genes are essential for melanin production. The two tyrosinases are directly involved in the synthesis of melanin, slc24a5, a potassium-dependent sodium/calcium exchanger is thought to modulate melanosomal pH, which is a crucial parameter for melanin synthesis (Ginger et al., 2008). Also pmel, the premelanosome protein, another melanophore-specific gene is slightly, yet not significantly, upregulated in bar regions compared to interbar regions (Figure 5C and Supplementary Table S2). Overall, expression data show a sixfold higher expression of tyrosinases in bars providing an explanation for the darker appearance of bar melanophores (Figure 6 and Supplementary Table S2). Even, if we standardize the differential expression by the increased cell-number or the melanophore-lineage marker mitfa we still find a 2–4.5-fold (Supplementary Figure S7 and Supplementary Tables S3, S4) increase of tyrosinase expression and therefore likely melanin synthesis.

Although we did not find significant differential expression of the transcription factors sox10 and mitfa (Figures 5A,B), both crucial factors for melanophore differentiation, expression is slightly elevated in bar regions. As mitfa expression has been reported to be only weak in differentiated melanophores (Johnson et al., 2011), this result might indicate an elevated number of melanoblasts. Although this would have to be analyzed in detail, it would explain the increased number of newly forming melanophores in the bars of juveniles (Figure 2 and Supplementary Figure S1) and consequently an increase in melanophore density.

Additionally, we used a marker gene for the xanthophore lineage (csf1ra) and the iridophore lineage (ltk). None of them showed any expression differences (Figures 5G,H), suggesting that there is no obvious difference in xanthophore and iridophore number. Therefore, melanophores seem to be the main cell type that clearly differs both in number and expression of marker genes between bar and interbar regions. Although there is some indication of a decreased number and coverage of xanthophores in the bar regions, this analysis was greatly hindered by the dense melanophore coverage. Use of mutant lines as for example of the oca2 gene that will likely not affect the pattering of xanthophores and iridophores (Kratochwil et al., 2019b) would be a possibility to better assess this question. Still, this result is somewhat surprising as cell-cell interactions between xanthophores, melanophores and iridophores seem to have a completely different dynamic in cichlids compared to zebrafish, where these cell types are spatially segregated (Mahalwar et al., 2014; Patterson et al., 2014; Singh et al., 2014; Eom and Parichy, 2017). In zebrafish, especially xanthophores and melanophores are almost mutually exclusive, with a decreased cell number of xanthophores in the melanic stripes. No evidence of such a strong antagonistic interaction could be found in H. latifasciatus. In xanthophore morphology though, zebrafish and H. latifasciatus have some shared featured: xanthophores in melanic regions are more irregular and seem smaller and have fewer and not as evenly distributed dendritic arbors (Figures 3, 4, Supplementary Figure S5, and Supplementary Table S1), although we did not find significant differences in xanthophore diameter (Figures 4I,J and Supplementary Table S1).

Our result on the expression of mc1r and mc5r shows that the melanocortin receptors likely play no role in shaping the vertical bar patterns in H. latifasciatus (Figures 5K,L). The antagonists asip1 and agrp2 show significant differences in expression (Figures 5I,J and Supplementary Table S2). The gene agrp2 has been shown to inhibit the formation of horizontal stripe patterns: species with high agrp2 expression lack any horizontal stripe patterns, while species with horizontal stripes have generally low expression of agrp2 in skin. Yet, spatial variation in agrp2 expression does not seem to be necessary for the formation of stripes (Kratochwil et al., 2018). On the other hand, a gradient in agrp2 expression along anterio-posterior axis could be found in both H. latifasciatus with significant difference from this study (Figure 5J and Supplementary Table S2) and P. nyererei without significant difference from our previous study (Kratochwil et al., 2018). However, based on our previous work that showed no changes in bar patterns in a knockout of agrp2 in P. nyererei (Kratochwil et al., 2018) or association with bar presence and absence (Kratochwil et al., 2019a), a role in vertical bar formation seems rather unlikely.

The variation in asip1 expression is intriguing. The melanocortin signaling antagonist asip1 has higher expression in interbars (twofold difference; Supplementary Table S2) with a particularly high expression in the first interbar. This is surprising as asip1 (and also the tetrapod homolog Asip) was previously known to vary along the dorso-ventral axis where it is involved in generating the dorso-ventral countershading that can be seen in many vertebrates (Manceau et al., 2011; Cal et al., 2017, 2019; Haupaix et al., 2018; Kratochwil, 2019). Here we observed variation along the anterio-posterior axis (Figure 5I and Supplementary Table S2), suggesting that gene expression differences of asip1 might contribute to variation in pigmentation on both axes. As Asip1 acts as Mc1r/Mc5r antagonist, one of the responses of decreased melanocortin signaling would be an increased aggregation of melanosomes. This would be in line with the significant difference in dispersion diameter we observe between bars and interbars. In interbars the diameter is on average 2.8 times smaller (Figures 4C,D and Supplementary Table S1), contributing to the lighter appearance of the interbar region (Figure 6). Yet the variation of asip1 expression serves only as a partial explanation as (a) the expression level in int2 is similar to the adjacent bars and (b) the dispersion/aggregation states of melanophores vary greatly within bars and interbars as well as between individuals. Still, based on these findings, it might be an important take-home-message to consider functions of asip1 beyond dorso-ventral patterning (i.e., countershading).

H. latifasciatus were kept in groups of 10–25 individuals. Fertilized eggs were removed from the gravid females as early as possible. The larvae and juvenile fish were then raised in egg incubators (ZET-E55, Ziss aqua) at room temperature. At 5 days post-fertilization (dpf) three larvae were separated from their siblings and each of them was raised in a single incubator to be able to keep track of a single individual in the following weeks. The hatching fish larvae were provided with nutrients by their yolk sac which lasts for ∼14 days. After 14 dpf the larvae were fed on Artemia nauplii twice a day. Experiments were performed in accordance with the rules of the animal research facility of the University of Konstanz, Germany and with the permission of the animal care committee (Regierungspräsidium) Freiburg, Germany (G18/60 and T16/13).

Photographs of H. latifasciatus embryos and larvae were captured with a stereomicroscope (Leica MZ10F) with a Leica DMC2900 color camera. Fish were first anesthetized with 0.04% tricaine (MS-222, Sigma-Aldrich). Images were taken as previously described (Kratochwil et al., 2015; Kratochwil et al., 2017). In order to capture the color pattern development of H. latifasciatus, photographs were taken from 7 to 21 dpf of three individuals.

For quantitative analyses, a comparable area including vb2 and vb3 of photographs of three individuals from stages between 9 and 21 dpf were put together, aligned (using the anterior dorsal fin as landmark), transformed into a black-and-white image and quantified in three different ways using Fiji (Schindelin et al., 2012). Melanosome dispersal was estimated by measuring the diameter of the minimally sized circle that encloses all melanic parts of a melanophore (minimal enclosing circle). It is therefore affected both by the melanophore size and the dispersal/aggregation state of melanosomes within the melanophores. The relative gray value was measured by taking the mean gray value of the same circle and dividing it by the mean of all melanophore gray value measurements of the same individual and stage. This relative value was used to account for differences between images and stages. Melanophore density was calculated from manual melanophore counts. As the position of the bars would shift as the fish is growing, we corrected the position values by the growth of the individual (size of individual at X dph/size of individual at 9 dph). The values were plotted using non-parametric regression (locally weighted scatterplot smoothing; LOWESS) in R (f parameter: 0.1). For the density calculation the area was split into ten equally-sized zones. To track individual melanophores we took images of the same individual on five consecutive days (from 9 to 13 dpf; Figures 2D–H). The identity of single melanophores across stages was estimated using overlays of the images from multiple days.

Three female individuals (standard length ∼12 cm) were examined for pigment quantification. We only quantified the chromatophores in females as the vertical bar pattern does not differ between sexes, we found more variation in red coloration in males (Figures 1Q,R). For each fish, in total 156 photos were obtained from scales and skin with scales removed (referred to as skin). Firstly, we separated the two flanks with scales still attached to skin. Scales from study regions were then carefully removed from the left flank and kept separately in Ringer’s solution (6.5 g/L NaCl, 0.25 g/L KCl, 0.3 g/L CaCl₂ and 0.2 g/L NaHCO₃) at 4°C before imaging. Five light microscope photographs were taken on each melanic and non-melanic region (the two melanic bar regions vb2 and vb3, and the interbar ib1 which anterior of vb2 and ib3 which between vb2 and vb3) from dorsal to ventral with a Leica MZ10F stereomicroscope equipped with a Leica DMC2900 color camera. Four photos were taken also along the dorsal-ventral axis on both region using Lecia MZ10F equipped with a Leica DFC 3000G black-and-white camera with GFP filter to image the auto-fluoresce of adult xanthophores (Kelsh et al., 1996; Guyader and Jesuthasan, 2002). For each pigmented region, ten scales were imaged using a Leica DMC2900 camera on Leica DM6B upright microscope. Leica Application Suite X software was used to capture the photos using the same setting. To count the number of melanophores, the right flank including both scales and skin tissue was treated with 10 mg/ml epinephrine (SIGMA-ALDRICH) for 20 min at room temperature to aggregate the melanosomes and thereby permit robust cell number quantifications (Figures 3B,F,K,N). After epinephrine treatment, tissue was washed by Phosphate Buffered Saline (PBS, pH7.4) three times and then kept in 4% Paraformaldehyde (PFA) in PBS. Epinephrine-treated scales were removed from skin and photos were taken for both skin and scales as described above. Fluorescence images were not taken as this treatment caused a high autofluorescence background.

All photos were analyzed with the image analysis software Fiji (Schindelin et al., 2012). At the beginning of the analysis we manually adjusted the color threshold to obtain reliable quantification of the bright field photos from non-treated skin and scales. Using this setting we executed the “Analyze Particles” function to obtain the melanophore coverage. For 20 randomly selected melanophores from each skin image and 10 melanophores from each scale image we measured the dispersal diameter of melanin covered parts of the melanophores (as described above). For xanthophore coverage we used the same approach. The dispersal diameter of the pigment filled part of 20 xanthophores of each skin image and 10 xanthophores of each scale image was measured from each photograph. The number of melanophore and xanthophore was counted from epinephrine-treated skin and scale specimens. Single xanthophores on scales could be observed easily by fluorescence microscopy, while not all boundaries of xanthophores in the skin could be easily identified as the cells often overlapped. Therefore, we were able to count the xanthophore number in scales but not in the skin dissections, while dispersal diameter measurements for xanthophores was possible on both scales and skin. Xanthophores and (what likely are erythrophores) was treated as the same cell type. Several studies show that erythrophores also exist in cichlids (Chen et al., 2014, 2015). However, vesicles containing pteridine and carotenoids could be found in the same cells, in which case the overall color depends on the ratio of red and yellow pigments (Matsumoto, 1965; Bagnara, 1966). Hence, the distinction between xanthophores and erythrophores is not always clear. Therefore, we classified yellow/orange/red colored cells all as xanthophores. Although we could identify iridophores in both dark bars and light interbars before and after epinephrine treatment (Figures 3I,K,L,N), iridophores were mostly below or above melanophores hindering reliable measurements. We, therefore, used gene expression of the iridophore lineage marker gene ltk (Lopes et al., 2008).

To compare the expression of coloration and pigment genes between melanic bars and interbars regions, we sampled the whole melanic and non-melanic skin region. Skin tissue was dissected and kept in RNAlater (Invitrogen) at 4°C overnight and then transferred to −20°C for long-term storage. RNAlater was removed prior to homogenization. Skin samples and the appropriate amount of TRIzol (Invitrogen) (1 ml TRIzol per 100 mg sample) were homogenized in 2 ml Lysing Matrix A tube (MP Biomedicals) using FastPrep-24 Classic Instrument (MP Biomedicals). RNA was extracted according to the manufacturer’s recommendations (Invitrogen) with an additional wash step by 75% Ethanol. Subsequent purification and on-column DNase treatment were performed with the RNeasy Mini Kit (Qiagen) and RNase-Free DNase Set (Qiagen). Following extraction and purification, RNA was quantified using Qubit RNA BR Assay Kit (Invitrogen) with Qubit Fluorometer (Life Technologies).

Gene expression analyses were performed on two melanic (vb2 and vb3) and two non-melanic (ib1 and ib2) skin regions. First strand cDNA was synthesized by using 1μg of total RNA with a GoScript Reverse Transcription System (Promega). qPCRs were performed with 2 μl of 100 μl synthesized first strand cDNA that was diluted ten times from 20 μl of initial reaction volume as a template, 10pmol of each forward primer and reverse primer, and GoTaq qPCR Master Mix (Promega) with nuclease-free water to make the final volume of 20 μl in a 96-well plate. Twelve genes were processed to examine the expression level including sox10, mitfa, csf1ra, ltk, pmel, slc24a5, tyr, tyrp1a, asip1, agrp2, mc1r, mc5r (Table 1). Primers are listed in Supplementary Table S5. We used 40 cycles of amplification on a CFX96 Real-Time PCR Detection System (Bio-Rad). The amplification program was: initial denaturation at 95°C for 10 min, 40 cycles of 95°C for 20 s, 60°C for 60 s. At the end of the cycles, melting curve of the products was verified for the specificity of PCR products. Only samples with one peak in the melting curves were processed to analyses. We assayed gene expression in triplicate for each sample and normalized the data using the reference genes β-actin and gapdh. Ct values were defined as the point at which fluorescence crossed a threshold (R_Ct) adjusted manually to be the point at which fluorescence rose above the background level. Next, we compared the relative expression between samples using the 2^–ΔΔCT method (Nolan et al., 2006). For group comparisons, we used ANOVA followed by Tukey’s HSD. All statistical tests were performed in R (R Development Core Team, 2019).