kcnj13 regulates pigment cell shapes in zebrafish and has diverged by cis-regulatory evolution between Danio species

Teleost fish produce some of the most diverse pigment patterns in nature, which are of great evolutionary importance as direct targets of natural and sexual selection. Closely related species of the genus Danio, including the widely used model organism zebrafish, Danio rerio, develop amazingly different patterns and are therefore excellent models to investigate the evolution of pigment pattern diversification in vertebrates (Irion and Nusslein-Volhard, 2019; Parichy, 2015, 2021; Patterson and Parichy, 2019; Singh and Nusslein-Volhard, 2015). Recently, the phylogenetic relationships in the Danio genus have been resolved, which led to the insight that a complex evolutionary history underlies their speciation and morphological diversification (McCluskey and Postlethwait, 2015).

The horizontally striped pattern in D. rerio emerges during metamorphosis when multipotent pigment cell progenitors derived from stem cells located at the dorsal root ganglia (DRGs) migrate into the skin (Dooley et al., 2013a; Singh et al., 2016, 2014). Here, they differentiate and form the pattern, presumably by a self-organizing process dependent on multiple cell-cell interactions (Frohnhofer et al., 2013; Patterson and Parichy, 2013). These interactions lead to the acquisition of location-dependent cell shapes, compact/yellow xanthophores and dense/reflective iridophores in the light stripes, and stellate xanthophores and loose/blue iridophores in the dark stripes (Gur et al., 2020; Mahalwar et al., 2014; Singh et al., 2014). Melanophores are restricted to the dark stripes. Precise superimposition of the differentially shaped pigment cells is required for colour and contrast of the pattern. The cellular interactions are, at least partially, mediated by direct cell-cell contacts through gap junctions, adhesion molecules and ion channels. Gap junctions are formed by two connexins (Gja4 and Gja5b) (Irion et al., 2014a; Watanabe et al., 2006, 2016), Igsf11 and Jam3b regulate adhesion (Eom et al., 2012, 2021) and Kcnj13 is an inwardly rectifying potassium channel (Iwashita et al., 2006). The diverse patterns in other Danio fish are produced by the same three types of pigment cells; however, the genetic and cell biological basis of the pattern variation is still largely unexplored. So far, the evolution in two separate cell differentiation pathways, xanthophore-specific Csf1 signalling in D. albolineatus and iridophore-specific Edn signalling in D. nigrofasciatus, has been linked to patterning differences (Parichy and Johnson, 2001; Patterson et al., 2014; Spiewak et al., 2018). This mode of evolution might partly cause changes in the timing and strength of the interactions between pigment cells, with cascading effects on their final distribution within the skin.

In this study, we focus on the diversification of pigment patterns between the two sister species D. rerio and D. aesculapii. Whereas D. rerio develop a very stereotypical pattern of sharp horizontal dark and light stripes on the flanks and in the anal and tail fins (Fig. 1A), in D. aesculapii a more variable pattern of vertical bars with lower contrast is formed anteriorly on the flank that dissolves into irregular spots posteriorly; the fins are not patterned, except for one dark stripe in the anal fin (Fig. 1B). We have shown that the potassium channel gene kcnj13 evolved to contribute to these patterning differences between the two species (Podobnik et al., 2020).

In D. rerio kcnj13 mutants fewer, wider and interrupted stripes develop, and melanophores and compact xanthophores fail to separate completely (Fig. 1C,E-G) (Haffter et al., 1996; Henke et al., 2017; Inaba et al., 2012; Irion et al., 2014a; Iwashita et al., 2006; Maderspacher and Nusslein-Volhard, 2003; Podobnik et al., 2020; Silic et al., 2020). A CRISPR/Cas9-mediated loss-of-function allele of kcnj13 in D. aesculapii showed that the gene is also required for the formation of vertical bars in this species. This null allele leads to a complete loss of any pattern with uniform distribution of mixed pigment cells in the skin (Fig. 1D) (Podobnik et al., 2020). Hybrids between the two species display stripes similar to the pattern in D. rerio. The evolutionary divergence of kcnj13 between D. rerio and D. aesculapii was demonstrated by reciprocal hybrids between wild-type and mutant fish (Podobnik et al., 2020). This genetic test is used to identify evolved genes by comparing the phenotypes of reciprocal hemizygotes; that is hybrids that carry a null allele from either one of the parental species in an otherwise identical genetic background (Stern, 2014). It depends on the ability to generate null alleles in a given species pair, which is possible in several Danio species since the introduction of the CRISPR/Cas9 system. Hemizygous hybrids between D. rerio kcnj13 mutant and D. aesculapii wild type display a spotted phenotype indicating that the D. aesculapii allele fails to complement the D. rerio null-allele, whereas the reciprocal hybrid in which the D. aesculapii allele was mutant displayed the striped phenotype of hybrids between the wild-type species. The different phenotypes demonstrated that the wild-type alleles from the two species are functionally no longer equivalent. Mutations in gja4, gja5b and igsf11 in D. aesculapii revealed functions for all these genes in the formation of the bar pattern. However, all hemizygous hybrids showed patterns indistinguishable from patterns of wild-type hybrids, ruling out functional evolution of these loci. Hybrids between D. rerio kcnj13 mutants and seven additional Danio species suggest that kcnj13 evolved independently several times in the genus, as the wild-type alleles from three different species, Danio aesculapii, Danio tinwini and Danio choprae, do not complement a D. rerio kcnj13 loss-of-function allele in hemizygous fish (Podobnik et al., 2020).

In chimeras produced by blastula transplantations, we corroborate previous studies (Iwashita et al., 2006; Maderspacher and Nusslein-Volhard, 2003) showing that kcnj13 function is cell-autonomously required in melanophores but not in xanthophores for normal stripe formation. In addition, we show that the gene function is also not required in iridophores, the third pigment cell type. In vitro experiments have shown that the function of kcnj13 is required for the depolarization of melanophore membranes upon contact with xanthophores (Inaba et al., 2012). This form of contact-dependent depolarization might underlie the repulsive interactions between melanophores and xanthophores during the establishment of the striped pattern. To test the effects of kcnj13 loss-of-function on the shapes of pigment cells in vivo we performed further blastula transplantations, fluorescence imaging of labelled pigment cells and cell-lineage tracing of marked clones. We find that the shapes of all three types of pigment cells are altered in the mutants, suggesting that cell-cell interactions responsible for the location-dependent acquisition of cell shapes are dependent on kcnj13 function and defective in the mutants. Using a newly generated CRISPR/Cas9-mediated knock-in reporter line, we detect kcnj13 expression in only very few differentiated melanophores in the skin, suggesting that kcnj13 function might be required only during a short period or in a subset of cells for a longer time during pattern formation.

The coding sequence for kcnj13 is highly conserved within the Danio genus, with very few non-synonymous changes between the species. However, it was not clear whether these changes between D. rerio and D. aesculapii are functionally relevant, or whether cis-regulatory evolution underlies kcnj13 divergence (Podobnik et al., 2020). We show that transgenic rescue of the kcnj13 mutant phenotype is possible with the wild-type coding sequences of both D. rerio and D. aesculapii, suggesting that both proteins are functionally equivalent. Strikingly, we observe a much higher expression of the D. rerio allele compared with the D. aesculapii allele in the skin of wild-type hybrids. We conclude that regulatory rather than protein changes underlie the evolution of the gene between D. rerio and D. aesculapii. The differences in the two patterns might result in part from the lower expression of kcnj13 in D. aesculapii leading to variation in pigment cell distribution and shapes reminiscent of those in D. rerio mutants deprived of kcnj13 activity.

To understand the function of kcnj13 during pattern formation, we focused on its role during stripe formation in D. rerio. Multiple dominant alleles of kcnj13 have been found in several independent genetic screens (Haffter et al., 1996; Henke et al., 2017; Inaba et al., 2012; Irion et al., 2014a; Iwashita et al., 2006; Maderspacher and Nusslein-Volhard, 2003; Podobnik et al., 2020; Silic et al., 2020). Fish homozygous for two dominant alleles (Fig. 1E,F) and homozygotes for a recessive loss-of-function allele (Fig. 1C) develop similar but variable phenotypes with fewer, wider and interrupted stripes. To test whether this variability in our stocks is attributable to the nature of the allele (dominant or recessive) or the genetic background, we compared different allelic combinations in F2 fish with the same genetic background and found that all of them led to indistinguishable phenotypes. This indicates that dominant and recessive alleles cause the same developmental effects in homozygous mutants (Fig. 1G,H), showing that the dominant alleles are dominant-negative in heterozygotes.

We followed the development of the mutant pattern during metamorphosis. As previously described (Maderspacher and Nusslein-Volhard, 2003), and comparable with wild type, melanophores in the mutants are cleared from the region of the first light stripe, where compact iridophores and xanthophores develop (Fig. 1I-M). However, unlike in wild-type fish, iridophores later fail to initiate the consecutive light stripes, which leads to a phenotype of fewer and broader stripes in the mutants with occasional interruptions (Fig. 1N-O′).

Melanophores, but not xanthophores, require kcnj13 function for stripe formation, as shown in chimeras created by blastula transplantations (Maderspacher and Nusslein-Volhard, 2003). We confirmed these findings and also tested the requirement of kcnj13 in iridophores. In these experiments, the donor embryos were mutant for kcnj13 and genetically able to provide only one of the three pigment cell types. Hosts were wild type for kcnj13 but lacking this pigment cell type. Thus, in three sets of transplantations, the resulting chimeras had one mutant pigment cell type placed adjacent to the other two wild-type cell types. In contrast to mutant xanthophores and iridophores, only mutant melanophores could not contribute to wild-type patterns in chimeras (Fig. 2A-C), leading to the conclusion that kcnj13 is cell-autonomously required in melanophores but not in xanthophores or iridophores. By transplanting kcnj13 mutant cells into slc45a2 mutant hosts (also known as albino) we further tested whether mutant melanophores can integrate into a normal pattern with wild-type melanophores in the chimeric animals. We observed disruptions in the striped pattern wherever mutant (pigmented) melanophores were present (Fig. 2D-D″). Similar severe pattern defects are never observed in chimeras that have received either wild-type melanophores (Dooley et al., 2013a) or no melanophores at all, suggesting the absence of any functional requirement in non-pigment cells. These results indicate that stripe formation requires kcnj13 function autonomously only in melanophores or their progenitors.

To investigate when kcnj13 functions in the melanophore lineage, we used CRISPR/Cas9-mediated homology-directed repair to produce a KalTA4::Venus knock-in line (for details see Materials and Methods) as a reporter for endogenous kcnj13 expression in D. rerio. In early larvae we observed expression in the pronephros, hindbrain and melanophores, a pattern very similar to previously published results obtained by in situ hybridization (Silic et al., 2020), suggesting that our reporter line faithfully recapitulates endogenous kcnj13 expression (Fig. 3A). During later stages, at the onset of metamorphosis, expression is detected in patches of cells in the spinal cord along the entire anterior-posterior (A-P) axis of the fish (Fig. 3B). These positions do not overlap with the DRGs, where the neural crest-derived stem cells for the pigment cells are located (Dooley et al., 2013a; Singh et al., 2016, 2014) (Fig. 3C). We conclude that kcnj13 does not provide a function for stripe formation in these cells, as our transplantation experiments indicate no functional requirement in non-pigment cells (Fig. 2D). Although the signals in the kidney and spinal cord persist throughout metamorphosis, we do not find expression of the reporter in pigment cell progenitors, but instead in a few xanthophores and melanized melanophores in the skin during the time of pattern formation (Fig. 3D-D″). These results show that kcnj13 is expressed at detectable levels only in a small subset of melanophores at any given time during pattern formation.

A key aspect of pigment pattern formation in D. rerio is the location-specific acquisition of different pigment cell shapes. In the dark stripes of wild-type D. rerio, melanophores are densely packed and compact, only cells located at the boundaries to the light stripes form long protrusions, possibly interacting directly with xanthophores and iridophores (Frohnhofer et al., 2013; Hamada et al., 2014). To investigate cell shapes in kcnj13 mutants we observed fish carrying Tg(kita::mCherry), which labels both xanthophores and melanophores. Some cells were unlabelled due to the variegation of the transgene, which allows us to visualize the shapes of the tightly packed melanophores. Similar to previous findings (Iwashita et al., 2006) we observed that in the dark regions in the mutants melanophores were less compact and less tightly packed compared with wild-type cells. We also found that the melanophores bordering the light stripes were less polarized, generally lacking the very long protrusions present in wild-type cells (Fig. 4A,B). We confirmed this observation by measuring the polarity and length of the fluorescent signal for each melanophore (Fig. S1). We then applied clustering to distinguish polarized from non-polarized cells and found a significantly higher proportion of polarized melanophores in the stripes of wild-type compared with mutant fish (24.5% and 5.3%, respectively; P<0.001, Chi-square test). This suggests that kcnj13 mutant melanophores do not interact with one another and with xanthophores and iridophores in the same way that wild-type melanophores do.

Next, we investigated the effect of kcnj13 mutations on xanthophore behaviour during stripe formation. Upon transplanting wild-type xanthophores, labelled with Tg(sox10:mRFP), into kcnj13 mutants, these cells acquire compact shapes in the dark stripe regions, where they normally appear stellate (Fig. 4C,D). Similar to findings from in vitro studies (Inaba et al., 2012), these results suggest that wild-type xanthophores are not always able to interact with mutant melanophores, which causes patterning defects in vivo.

To assess the effects of mutations in kcnj13 on iridophores, we induced fluorescently labelled clones in the mutants using a Tg(sox10:cre-ERt2) line (Singh et al., 2014) and followed labelled iridophores during metamorphosis. We found clones of dense iridophores, which are characteristic for light stripes, in the dark stripe area (Fig. 4E,F). This result suggests that iridophores require the presence of and interaction with melanophores to acquire the loose form; and that this interaction depends on kcnj13 function. Thus, iridophores might not be able to recognise mutant melanophores and therefore develop ectopically in the dense form in the dark stripe regions. We conclude that kcnj13 function, required in melanophores, is important for homotypic and heterotypic pigment cell interactions, which control the location-dependent cell shape acquisition of all three pigment cell types during pattern formation. These cumulative effects might inhibit the reiteration of dark and light stripes in the mutant fish.

Melanophores in D. rerio produce pronounced polarized protrusions towards compact xanthophores and both cell types are strictly separated between the light and dark stripes. The polarity of the protrusions is lost in kcnj13 mutants, where both cell types also mix occasionally (Fig. 4A,B,G,H). This mutant phenotype is similar to the situation in wild-type D. aesculapii, where we found a mixing of cells and no pronounced polarity of melanophores towards xanthophores (Fig. 4I). The contrast of the bar pattern is therefore reduced; there is no contrast in D. aesculapii kcnj13 mutants, where all pigment cells mix and no bars are formed (Podobnik et al., 2020). Our observations suggest that the divergence of the pigment patterns between D. rerio and D. aesculapii could partially be due to evolutionary changes in the interactions between all three pigment cell types, which influence the cell shapes.

To investigate the channel structure of Kcnj13 (Kir7.1), we expressed the D. rerio protein fused to mCherry using a Multibac-derived baculovirus/insect cell expression system (Altmannova et al., 2021; Bieniossek et al., 2012), purified the recombinant protein by affinity and size-exclusion chromatography, and measured the molecular mass with mass photometry (Young et al., 2018) (Fig. S2). The results suggest that Kcnj13 exists as a homo-tetramer, which can explain the dominant-negative effects observed in alleles carrying point mutations affecting the selectivity filter or the second transmembrane helix (Fig. 1A,B) as caused by mutant proteins negatively interfering with wild-type copies in the complex in heterozygous fish (Haffter et al., 1996; Henke et al., 2017; Irion et al., 2014a; Iwashita et al., 2006; Podobnik et al., 2020; Silic et al., 2020). We constructed homology-based and AlphaFold-multimer models of the homo-tetrameric Kcnj13 channel (the models are available in ModelArchive at https://modelarchive.org/doi/10.5452/ma-xpcgr). These models agree with published structures of similar potassium channels. The protein sequences of D. rerio and D. aesculapii differ only by two amino acid residues (Q23L and D180G; Fig. 5E, magenta) in the cytoplasmic domain (Podobnik et al., 2020); structure modelling of the two alleles is insensitive to these differences.

Reciprocal hemizygosity tests showed that the divergence of kcnj13 must reside within the locus, either in the protein-coding region or in cis-regulatory elements, but cannot be due to trans-acting factors (Podobnik et al., 2020). To test whether the amino acid changes identified between the two species contribute to the evolution of kcnj13, we used Tol2 transgenesis to express the coding regions from D. rerio or D. aesculapii under the control of the melanophore-specific mitfa promoter in kcnj13 null-mutant D. rerio (Fig. 5A-D). In both cases the transgenes were able to restore the striped pattern in the trunk of the fish, indicating that the protein from D. aesculapii can function in a similar manner to the D. rerio protein (Fig. 5E). We observed some differences in the rescue capabilities of the transgenes among the lines we established, possibly due to copy number variations and expression differences of the randomly inserted transgenes. The striped pattern of the caudal fin was never restored in the transgenic lines, most likely due to the inactivity of the promoter at the appropriate time points in this tissue, corroborating the finding of fundamental mechanistic differences in pigment pattern formation between the trunk and fin (Frohnhöfer et al., 2013). Our results suggest that the coding regions from both species function similarly and that the protein-coding changes are irrelevant for kcnj13 divergence.

Therefore cis-regulatory changes likely underlie kcnj13 evolution and patterning differences between the two species. To test this prediction, we produced hybrids between the two species and performed allele-specific expression analysis in the skin and posterior trunk, which includes the skin, of adult fish. We found significantly higher levels of the D. rerio allele compared with the D. aesculapii allele in skin and trunk (Fig. 5F; Fig. S3), indicating species-specific regulation of the locus and thereby confirming cis-regulatory evolution. Quantitative differences in expression levels might cause differences in pigment cell interactions and shapes observed between D. rerio and D. aesculapii. Based on the repeated and independent evolution of the ancestral kcnj13 function in the Danio genus (Podobnik et al., 2020) we speculate that similar cis-regulatory changes might also have occurred in D. tinwini and D. choprae. Our results highlight the Danio genus as an excellent model system to study the molecular, genetic and cellular basis of pigment pattern diversification in vertebrates.

Teleost fish produce some of the most intricate pigmentation patterns in nature. However, in only a few species have the pattern forming mechanisms been studied in detail. D. rerio, an excellent vertebrate model organism widely used in research, shows a conspicuous pattern of horizontal stripes on the flank and on the anal and tail fins. This pattern is produced by three types of pigment cells interacting in complex ways to self-organize into dark and light stripes. During pattern formation the horizontal myoseptum serves as an anatomical pre-pattern for the orientation of the stripes. The stripes in the anal and tail fins are contiguous with the stripes in the body. However, the fin pattern is formed by a different, possibly somewhat simpler, mechanism that involves only two cell types, melanophores and xanthophores. Cellular interactions mediated by direct cell-cell contacts depending on gap junctions and adhesion molecules are essential for stripe formation as demonstrated by the spotted phenotypes of gja4, gja5b, igsf11 and jam3b mutants (Eom et al., 2012, 2021; Irion et al., 2014a; Watanabe et al., 2006). In addition, mutations in kcnj13 lead to defects in the pattern with fewer, wider and interrupted stripes and occasional mixing of compact xanthophores with melanophores (Haffter et al., 1996; Henke et al., 2017; Inaba et al., 2012; Irion et al., 2014a; Iwashita et al., 2006; Podobnik et al., 2020; Silic et al., 2020). Kcnj13 regulates the membrane potential of melanophores (Inaba et al., 2012), which might be important for the repulsion between xanthophores and melanophores. By interspecies complementation tests in Danio hybrids it was previously shown that of these four genes, only the function of kcnj13 diverged within the Danio genus, probably several times independently (Podobnik et al., 2020).

To better understand the role of kcnj13 in pattern formation and diversification, we examined its function in D. rerio in more detail. All kcnj13 alleles isolated in genetic screens are dominant with a relatively weak heterozygous and considerably stronger homozygous phenotype. We previously produced a loss-of-function allele, which is completely recessive (Podobnik et al., 2020). The phenotypes of homozygous fish for a dominant or the recessive allele in the same genetic background are indistinguishable (Fig. 1G,H). This demonstrates that the dominant alleles are in fact dominant-negatives and not neomorphs. The variability we observe in our mutant strains is dependent on the genetic background.

Phenotypic analysis of chimeras obtained by blastula transplantations had already demonstrated the autonomous requirement of kcnj13 function in melanophores but not in xanthophores (Maderspacher and Nusslein-Volhard, 2003). We repeated these transplantation experiments including the third pigment cell type, iridophores. Our results show that kcnj13 function is required only in melanophores for stripe formation in D. rerio, but not in any other cell type (Fig. 2A-C). In addition, we find that mutant melanophores lead to strong patterning defects when transplanted into wild-type fish (Fig. 2D). This shows that the mutant cells are not guided by their wild-type neighbours but influence the patterning process cell-autonomously, possibly failing to instruct neighbouring xanthophores and iridophores.

Our results support the previous observation that a kcnj13 transgene expressed under the control of the mitfa promoter, which is known to be active in melanophores and their stem cells (Dooley et al., 2013a), can rescue the mutant phenotype in the trunk (Inaba et al., 2012). As these experiments were conducted in the presence of a dominant-negative kcnj13 allele, which impedes the wild-type channel function, a complete rescue could not be expected. In our transgenic rescue experiments using the recessive mutant, expression of kcnj13 using the mitfa promoter restores the stripes on the flank of the fish to a pattern very similar to the one observed in wild types (Fig. 5B), which further supports the notion that kcnj13 is required in melanophores. The striped pattern in the anal and tail fins is not restored by the transgenes, suggesting that expression under the melanophore-specific mitfa promoter does not recapitulate all aspects of the endogenous expression pattern of kcnj13, and mechanisms that form stripes in the fins are fundamentally different from those that form stripes in the trunk (Frohnhofer et al., 2013).

To visualize the expression pattern of kcnj13 in D. rerio, we made a reporter line by homology directed knock-in of an optimized GAL4 coding sequence (KalTA4) into the endogenous locus. In combination with a UAS:Venus transgene, this reporter line shows expression in early larvae in the pronephros and melanophores (Fig. 3A,B), very similar to published data from in situ hybridizations (Silic et al., 2020), indicating that our line faithfully recapitulates kcnj13 expression. Later, during metamorphosis when the pigment pattern is formed and also in adult fish, we detected expression in neurons of the spinal cord (Fig. 4C). During these stages in situ hybridizations are difficult in D. rerio and we rely on the reporter to indicate expression of the gene. As our transplantation experiments clearly show a cell-autonomous requirement of kcnj13 in melanophores or their precursors (Fig. 2D) we can rule out a function of the gene for pattern formation in these neuronal cells. We also found expression of the reporter line during later stages in few xanthophores and, unexpectedly, only in a small subset of melanophores (Fig. 4D). Expression of the reporter in xanthophores might reflect earlier activation in a common precursor for melanophores and xanthophores and the long persistence of the proteins (KalTA4 and Venus). Alternatively, kcnj13 could genuinely be expressed in xanthophores but without any obvious function in stripe formation. Our observation that we cannot detect kcnj13 expression in all melanophores at any given time point suggests that it is either required only very transiently or that only a few cells depend on kcnj13 function and then influence the behaviours of all the pigment cells. Alternatively, our reporter might not be sensitive enough to allow the detection of very low expression levels, which could nevertheless be relevant for pattern formation. A different possibility is that the channel protein might be very stable and present in the cell membrane for prolonged times even after transcription has ceased and also the reporter is no longer detectable. In any case, our data is consistent with published data from single-cell RNA-sequencing (Saunders et al., 2019), which also show expression of kcnj13 to be low and limited to a very minor fraction of pigment cell progenitors as well as differentiated melanophores and xanthophores.

We conclude that kcnj13 is only required in melanophores during pattern development. Mutant melanophores are less compact and less tightly packed, affecting the tiling within the dark stripe. Mutant melanophores at the stripe boundaries also do not form polarized protrusions towards the light stripes (Fig. 4A,B). The significance of these protrusions is unclear, they could be used for direct repulsive interactions with xanthophores or iridophores to delineate the boundary between light and dark stripe (Frohnhofer et al., 2013; Hamada et al., 2014). In kcnj13 mutants homotypic and heterotypic interactions, among melanophores and between melanophores and the other two pigment cell types, are affected, as seen, for example, by the mixing of the cells. We find that the shapes of both cell types are affected in kcnj13 mutants, with dense iridophores and compact xanthophores, which are limited to the light stripes in wild type, also appearing in dark stripe regions. Therefore, we conclude that melanophores play a crucial kcnj13-dependent role in directing dark stripe-specific cell shape transitions in both iridophores and xanthophores. In the absence of Kcnj13, all three types of pigment cells may lose their dark stripe-specific shapes, which might indicate that the default shapes for xanthophores and iridophores are the ones these cells acquire in the light stripe region.

The same types of pigment cells that are found in D. rerio form a range of very different patterns in closely related Danio species. The specification and differentiation of pigment cells are similar in D. rerio and D. aesculapii. They both require Mitfa and Kit signalling in melanophores and Csf1 and Ltk signalling in xanthophores and iridophores, respectively (McCluskey et al., 2021a; Podobnik et al., 2020). Mutants indicate that iridophores do not emerge along the horizontal myoseptum, are lower in number and dispensable for bar formation in D. aesculapii, whereas they guide stripe formation in D. rerio (Podobnik et al., 2020). Whether genes required for iridophore development have evolved between these two species is not known. However, for another species, D. nigrofasciatus, it was shown that reduced iridophore proliferation contributes to a reduction in stripe number and integrity (Spiewak et al., 2018). In addition, species-specific differences in the developmental timing of pigment cell proliferation and differentiation can lead to patterning differences as observed for xanthophores, which differentiate precociously in D. albolineatus resulting in a loss of the striped pattern (Patterson et al., 2014). We find that melanophores in D. aesculapii do not form long protrusions towards the light regions (Fig. 4G-I), which is similar to kcnj13 mutants in D. rerio (Fig. 4A,B). In D. rerio these protrusions might partly regulate melanophore survival (Hamada et al., 2014) and the overall stability of the boundary between dark and light stripes. Similar to the D. rerio mutant, the lack of such protrusions in D. aesculapii might indicate a less robust mechanism for the consolidation of the boundary between dark bars and light regions (Fig. 4I), where melanophores and xanthophores frequently mix. The fact that D. rerio mutants and D. aesculapii wild types still develop melanophore protrusions suggest that other factors than kcnj13 are also important for their establishment.

When tested in D. aesculapii the four genes (kcnj13, gja4, gja5b and igsf11), known to function in cell-cell interactions during stripe formation in D. rerio, were found to be also required to form the bar pattern (Podobnik et al., 2020). Although residual patterns of spots or wider and interrupted stripes still form in D. rerio mutants, the bar pattern is completely lost in D. aesculapii mutants and all pigment cells intermingle and distribute evenly in the skin, a phenotype only seen in double mutants in D. rerio. This indicates that cellular interactions in both species occur but are more complex in D. rerio, which could lead to a higher robustness of the patterning mechanism in this species. Reciprocal hemizygosity tests for all four genes lead to the conclusion that there is functional conservation in three cases, gja4, gja5b and igsf11, whereas only kcnj13 diverged between the two species (Podobnik et al., 2020). Thus, the formation of the very different patterns of horizontal stripes and vertical bars involves the same players. Three of these, Kcnj13 and the two gap junction proteins, might be involved in an electric coupling of pigment cells, which could allow coordinated tissue-scale patterning (Harris, 2021). Evolution in kcnj13 between the two species might influence the conditions for these interactions, with the consequence of evolutionary change in patterning.

In our rescue experiments the coding sequences from both species, D. rerio and D. aesculapii, were equally able to restore stripe formation in D. rerio kcnj13 mutants, indicating functional equivalency. However, the use of a non-native promoter and possible position effects due to random integration of the transgenes might obscure subtle functional differences between the two proteins. This question could be addressed in the future by precise exchanges in the coding sequence of the endogenous locus in D. rerio. However, we found allele-specific differences of kcnj13 expression in hybrids with much higher levels of expression from the D. rerio allele (Fig. 5F) clearly indicating regulatory differences between the loci from the two species. Therefore, the functional divergence of kcnj13 between D. rerio and D. aesculapii is most likely caused by evolution of cis-regulatory elements affecting the levels of expression of the gene. Cis-regulatory evolution has been implicated in other cases of pattern diversification of Danio fish. In D. albolineatus, the increased expression of Csf1 causes early differentiation of xanthophores leading to a loss of the striped pattern and the mixing of pigment cells (Patterson et al., 2014). In D. nigrofasciatus iridophore development is reduced due to cis-regulatory changes in the edn3b gene leading to an attenuated pattern with fewer melanophores and stripes, similar to hypomorphic D. rerio mutants (Spiewak et al., 2018). In the rare case of D. kyathit and D. quagga, hybrids between the two species are fertile, which allows for quantitative trait locus (QTL) mapping. QTL analysis for differences between the spotted D. kyathit and the striped D. quagga led to the identification of a complex genetic basis for the pattern differences with multiple candidate loci, probably involving changes in a number of regulatory regions (McCluskey et al., 2021b). In the more distantly related cichlids, bars and stripes evolved repeatedly in species endemic to the Great African Lakes. Here, QTL mapping identified regulatory changes in the gene agouti-related peptide 2 (agrp2) that underly these patterning differences (Kratochwil et al., 2018).

In three-spine sticklebacks, genome-wide association studies identified loci underlying repeated ecological adaptations in independent pairs of fresh- and saltwater populations (Jones et al., 2012). These adaptive loci are predominantly affected by cis-regulatory changes leading to differences in gene expression in the gills (Verta and Jones, 2019). In contrast, divergent development of the pharyngeal tooth plate in sticklebacks is shaped primarily by evolution of trans-acting factors (Hart et al., 2018). It was speculated that the genetic architecture of teeth formation is less complex than the adaptations to salt handling; evolution of trans-acting factors might therefore be less pleiotropic in dental tissue compared with multifunctional gills.

Dominant mutations in kcnj13 in D. rerio cause pigment pattern defects but also late-onset retinal degeneration (Toms et al., 2019a, b), similar to mutations in the human orthologue that are known to cause two rare retinal diseases (Hejtmancik et al., 2008; Sergouniotis et al., 2011). Mutations in mice lead to lethal defects in tracheal development (Yin et al., 2018). Owing to this observed pleiotropy, protein evolution might be highly constrained, favouring regulatory evolution. In general pigment patterns appear to evolve often by regulatory mutations, whereas pigmentation frequently diverges by protein changes (Orteu and Jiggins, 2020). However, constraints on regulatory evolution also exist; ectopic expression of kcnj13 in the dermomyotome leads to a long-finned phenotype (Silic et al., 2020). Our results suggest differential regulation of kcnj13 in the skin, where its expression appears to be restricted to pigment cells. Cis-regulatory evolution in kcnj13 specifically affecting expression in the skin is presumably non-pleiotropic and might therefore be more permissive for evolutionary change influencing pigment cell behaviour.

A basic colour-forming unit in cold-blooded vertebrates, fish, amphibians and reptiles consists of xanthophores in the top layer, iridophores in the middle layer and melanophores in the bottom layer. Melanophores appear black in the absence of shiny iridophores and yellow-orange xanthophores on top, as in D. rerio ltk (shady) or csf1ra (pfeffer) mutants. Modifications of this basic arrangement of pigment cells can yield diverse colourations. By varying the mechanisms that regulate pigment cell shape and layering, differences in colour, brightness and contrast can be achieved. In this regard our study points towards kcnj13 as a key node for evolutionary tinkering that underlies colour pattern diversification in teleosts. D. rerio kcnj13 mutants develop light and dark stripe regions low in contrast due to pigment cells that lack location-specific shapes and colouration. Regulation of colouration by cell shape transition may point to an important mechanism employed across evolution, where layer-specific and location-specific arrangement of diverse pigment cell types leads to species-specific colouration.

No statistical methods were used to predetermine sample size. The experiments were not randomized. The investigators were aware of allocation during experiments and outcome assessment.

D. rerio and D. aesculapii were maintained as described in Brand et al. (2002). If not newly generated (Table S1), the following lines were used for experiments: D. rerio wild-type Tuebingen (TU), kcnj13^t24ui (Podobnik et al., 2020), kcnj13^td15 (Iwashita et al., 2006), kcnj13^tdxg6 (Irion et al., 2014a), nacre/mitfa^w2 (Lister et al., 1999), pfeffer/csf1ra^tm236b (Odenthal et al., 1996; Parichy et al., 2000b), rose/ednrba^tlf802 (Parichy et al., 2000a), albino/slc45a2^b4 (Dooley et al., 2013b), sparse/kita^b134 (Kelsh et al., 1996), Tg(sox10:mRFP) (Singh et al., 2014), Et(kita:GalTA4,UAS:mCherry)hzm1 (Distel et al., 2009), Tg(sox10:ERT2-Cre);Tg(bactin2:loxP-STOP-loxP-DsRed-express) (Bertrand et al., 2010; Mongera et al., 2013) and D. aesculapii kcnj13^tmp11 (Podobnik et al., 2020). Interspecific hybrids between D. rerio and D. aesculapii were obtained by in vitro fertilizations (Parichy and Johnson, 2001). All species were staged according to the normal table of D. rerio development (Parichy et al., 2009). All animal experiments were performed in accordance with the rules of the State of Baden-Württemberg, Germany, and approved by the Regierungspräsidium Tübingen.

To generate the transgenic rescue lines plasmids with the mitfa promoter sequence from D. rerio (Lister et al., 1999), the coding sequences of kcnj13 from D. rerio or D. aesculapii and the coding sequence of sfGFP was constructed. The construct was subcloned into the Tol2 vector pGEM-T pminiTol2 carrying SV40 elements, a green heart marker cmlc2:Venus and Tol2 restriction sites (Kawakami and Shima, 1999; Yelon et al., 1999). The resulting plasmids were designated as pTol2gh-mitfa-kcnj13^D.rerio-sfGFP (GenBank accession number: OP326275) and pTol2gh-mitfa-kcnj13^D.aesculapii-sfGFP (GenBank accession number: OP326276). Tol2 transgenesis was performed as previously described (Kawakami and Shima, 1999); briefly, a solution (12.5 ng/µl Tol2 mRNA, 50 ng/µl plasmid DNA and 5% Phenol Red) was injected into fertilized eggs of D. rerio kcnj13^t24ui at the one-cell stage. One hundred F0 embryos were selected for marker gene expression at ∼2 days post-fertilization (dpf) and raised to adulthood. Mature F0 founder fish were outcrossed to D. rerio kcnj13^t24ui and F1 larvae positive for marker gene expression were selected to obtain stable transgenic lines. In both cases, lines were identified in which the mutant phenotype was partially rescued. These lines were designated as Tg(mitfa:kcnj13^D.rerio);kcnj13^t24ui and Tg(mitfa:kcnj13^D.aesculapii);kcnj13^t24ui, outcrossed to D. rerio kcnj13^t24ui and selected for marker gene expression in embryos and intact stripe patterns in adults for at least three generations (Table S1).

To generate a D. rerio UAS:Venus line, a plasmid with the coding sequence for the Venus-variant of YFP under the control of the yeast transcription factor GAL4 (six UAS-sites) was constructed (pminiTol2_UAS:Venus, GenBank accession: OP243708); mRNA for the Tol2 transposase was transcribed in vitro from the plasmid pCS2FA-transposase (Kwan et al., 2007) using the mMessageMachine and Poly-A tailing Kits (Invitrogen). TU embryos at the one-cell stage were injected with ∼2-4 nl of injection mix containing 250 ng/µl of in vitro transcribed mRNA and 25 ng/µl of plasmid DNA in PBS with Phenol Red as a tracer dye. The adult F0 fish were crossed to TU and the F1 larvae were screened for expression of the mCherry marker in the heart. From the positive F1 fish a stable line was established by another outcross to TU followed by sibling matings of the F2 fish (Table S1).

For gene knockouts the CRISPR/Cas9 system was applied either as described in Irion et al. (2014b) or according to the guidelines for embryo microinjection of Integrated DNA Technologies (IDT). Briefly, oligonucleotides were cloned into pDR274 to generate the sgRNA vector. sgRNAs were transcribed from the linearized vector using the MEGAscript T7 Transcription Kit (Invitrogen). Alternatively, target-specific crRNAs and universal tracrRNAs were purchased from IDT. Cas9 was expressed as a fusion protein with mCherry in Escherichia coli [BL21(DE)3pLysS] from the plasmid pOPT-Kan_Cas9-mCherry (GenBank accession: OP243709) and purified via double affinity chromatography (His-Tag and Twin-StrepTag) using standard procedures. Before use, the purified protein was dialyzed into PBS containing 300 mM NaCl and 150 mM KCl, aliquoted and stored at −70°C. sgRNAs or crRNA:tracrRNA duplexes were injected as ribonucleoprotein complexes with Cas9 proteins into one-cell-stage embryos. The efficiency of indel generation was tested on eight larvae at 1 dpf by PCR using specific primer pairs and by sequence analysis as described previously (Meeker et al., 2007). The remaining larvae were raised to adulthood. Mature F0 fish carrying indels were outcrossed. Loss-of-function alleles in heterozygous F1 fish were selected to establish homozygous or trans-heterozygous mutant lines (Table S1).

The CRISPR/Cas9-system was used to generate a reporter line for the expression of kcnj13. For the sgRNA template two oligonucleotides (5′-TAGGCCGTCTTTGCTGACCAGG-3′ and 5′-AAACCCTGGTCAGCAAAGACGG-3′) were annealed and cloned into pDR274; the RNA was transcribed in vitro with the MegaScript Kit from Invitrogen. A donor plasmid was constructed containing the KalTA4 variant (Distel et al., 2009) of the GAL4 coding sequence flanked by homology arms and CRISPR target sites (GenBank accession: OP243710). This plasmid (25 ng/µl) was co-injected with Cas9 protein (500 ng/µl) and sgRNA (35 ng/µl) into one-cell-stage embryos from the UAS:Venus line. The resulting F0 fish were backcrossed to UAS:Venus and the F1 larvae were screened for expression of Venus. One founder fish was identified with offspring showing a very strong early signal in the yolk and later also in the pronephros and melanophores, consistent with published expression data (Table S1). To achieve good imaging conditions in this line we generated the loss-of-function allele slc45a2^t22mp, as previously described (Irion et al., 2014b) (Table S1).

Chimeric animals in Figs 2A-D and 4D were generated by transplantations of cells during blastula stage as described in Kane and Kishimoto (2002).

Cre induction was carried out as described in Singh et al. (2014). Labelled clones in Fig. 4E,F were from fish followed over pattern development.

Anaesthesia of postembryonic and adult fish was performed as described previously (Singh et al., 2014). Brightfield images of adult fish in Figs 1A-H and 2A-D were obtained using a Canon 5D Mk II camera. To visualize melanophore protrusions via dispersion of melanosomes using brightfield imaging (Fig. 4G-I), fish were kept in the dark with a final concentration of 100 µM yohimbine (CAS: 65-19-0, Sigma-Aldrich) for 30 min before imaging as described in Hamada et al. (2014). Fish with different pigment patterns vary considerably in contrast, thus requiring different settings for aperture and exposure time, which can result in slightly different colour representations in the pictures. Fluorescence images of postembryonic and adult fish were acquired on a Zeiss LSM 780 NLO confocal (BioOptics Facility, Max Planck Institute for Biology Tübingen, Germany) and a Leica M205 FA stereomicroscope. Repeated imaging of pigment cell clones in metamorphic D. rerio was performed as described in Singh et al. (2014). Maximum intensity projections of confocal scans were uniformly adjusted for brightness and contrast. Images were processed using Adobe Photoshop, Adobe Illustrator CS6 and Fiji (Preibisch et al., 2009).

Images analysis and statistical analyses shown in Fig. S1 were performed in Julia 1.8.5. As necessary, images were scaled down to the same pixel per micron density for both genotypes. Images were cropped or rotated to remove fluorescence from the image that originated from non-melanophore cells, and to approximately position the centre of the stripe in the centre of the image. Cell centres were manually annotated, and regions of interest were then automatically segmented based on the position of cell centres using a watershed method. Each pixel in the region of interest of a was measured for its fluorescence intensity, distance from the cell centre and radial position. Pixels above a baseline intensity threshold were used to compute the radial and distance histograms showing the distance from the cell centre and the number of fluorescent pixels at all possible angles (radial frequency) from the cell centres. The maximum distance and radial frequency for each cell was then used to fit a Gaussian Mixture Model, identifying two distinct clusters of polarized and non-polarized cells. The proportion of cells in each cluster was statistically compared across genotypes using a Chi square test.

We expressed Kcnj13-mCherry with N-terminal His-tags in Sf9-insect cells using a baculovirus/insect cell expression system (Altmannova et al., 2021; Bieniossek et al., 2012). Pink pellets were washed with PBS, stored at −70°C, and later purified at 4°C at all stages. We selected n-Dodecyl-B-D-Maltoside (DDM, Serva Electrophoresis) detergents at ∼2× critical micelle concentration (CMC) and supplied cholesteryl hemisuccinate (CHS, Serva Electrophoresis) lipids for solubilization of the membrane protein. Cell pellets were resuspended in lysis buffer A, treated with a high-pressure homogeniser (Avestin EmulsiFlex-C3) and samples were centrifuged at 100,000 g for 1 h. The supernatant was incubated with Ni-NTA beads for 4 h and applied to a polypropylene column (BioRad) equilibrated in lysis buffer A. The column was washed with buffers B and C, and protein was eluted with buffer D. Fractions were isolated based on pink-marker colouration and concentrated using an Amicon Ultra-15 filter (100 kDa cut-off). The concentrated sample was spun for 1 h on a table-top centrifuge at full speed (21,000 g) and supernatant was applied onto a Superose 6 Increase 5/150 GL column for gel filtration using buffer E. Buffer compositions are provided in Table S2.

Measurements were performed in buffer E (Table S2) using an One^MP mass photometer (Refeyn, Oxford, UK) (Young et al., 2018). Immediately before analysis, the sample was diluted 1:10 with the aforementioned buffer. Molecular mass was determined in the analysis software provided by the manufacturer using a NativeMark (Invitrogen).

The homology model of the tetrameric Kcnj13 channel (Fig. 5E) was built using SWISS-MODEL (Bertoni et al., 2017; Bienert et al., 2017; Studer et al., 2020, 2021; Waterhouse et al., 2018) based on the crystal structure template (2.6-Å resolution) of the potassium channel Kir2.2 from Gallus gallus (PDB ID: 3spg), sharing a sequence similarity of 37% with the target protein Kcnj13 from D. rerio. Similar models with a pTM-based confidence score of ∼60% were generated using AlphaFold-Multimer (Evans et al., 2022 preprint; Jumper et al., 2021).

Reciprocal crosses between species [male D. aesculapii×female D. rerio (pair 1), and male D. rerio×female D. aesculapii (pair 2)] were performed via in vitro fertilization to produce F1 hybrids. Adult parental fish (n=4) and F1 hybrids (n=12; seven hybrids from cross 1, five hybrids from cross 2) were euthanized by exposure to buffered 0.5 g/l MS-222 (Tricaine). Tissues were dissected in ice-cold PBS and collected using TRIzol (Life Technologies). DNA from the parental individuals was isolated from posterior trunk tissue including the fins. RNA was obtained from skin and posterior trunk tissue of F1 hybrids. RNA integrity and quantity were assessed using the Agilent 2100 Bioanalyzer. Metadata is provided in Table S3. Library preparation [DNA/RNA: TruSeq DNA Nano Kit (Illumina); 100ng per sample] and sequencing [NovaSeq 6000 (Illumina), 2×250 bp for DNA and 2×100 bp for RNA] were performed by CeGaT GmbH (Tübingen, Germany). Data are available from PRJEB53585.

All subsequent analyses were based on high-quality clean reads. Quality of the sequencing data was checked using FastQC (version 0.11.9) and adapter sequences were trimmed using fastp (version 0.23.2) (Chen et al., 2018). Genome resequencing reads were aligned to the reference genomes of D. rerio (GRCz11) or D. aesculapii (NCBI ID92583) using BWA-MEM (version 0.7.17-r1188) (Li, 2013 preprint). The aligned SAM files were sorted and converted into BAM files using SAMtools (version 1.11) (Danecek et al., 2021). Then the sorted BAM files were de-realigned and indexed again using Picard (version 2.18.29, https://broadinstitute.github.io/picard/). Transcriptomes were aligned to the reference genomes using STAR aligner (version 2.7.10a) (Dobin et al., 2013). The BAM files directly output by STAR in two-pass mode were deduplicated and indexed by Picard.

To identify species-specific alleles, variant calling was performed according to the best practice pipeline of the Genome Analysis Toolkit (GATK4) (Brouard et al., 2019; McKenna et al., 2010). Specifically, Haplotypecaller was used to detect variants based on genome and transcriptome data. The called variants were joint-genotyped using GentypeGVCFs into a single .vcf file; data from skin and trunk tissue were separately processed. First, SelectVariants was used to filter single nucleotide polymorphisms (SNPs), then the selected SNPs were hard-filtered using Variantfiltration. Specifically, SNPs of ‘QUAL<30.0, QD<2, FS>60, MQ<40, SOR>3, MQRankSum<−12.5 and ReadPosRankSum>−8’ as well as non-biallelic SNPs were filtered out. The remaining SNPs were filtered again using VCFtools (--max-missing 0.8, --maf 0.05). Finally, SNPs shared by genomes and transcriptomes were selected for the subsequent allele-specific expression analysis (ASE) using the intersect function of Bedtools (version 2.30.0) (Quinlan and Hall, 2010).

Read counts for species-specific SNPs were averaged per gene for each hybrid transcriptome using GATK ASEReadCounter (Castel et al., 2015) with default filters enabled. Significant allele-specific expression was defined as ‘fold change’>2 between alleles and adjusted P-values (p-adj)<0.05 from DESeq2 package in R (Love et al., 2014). Finally, the ggplot2 (Wickham, 2016) package in R rendered volcano plots using the data obtained by DESeq2.

We thank Hans-Martin Maischein (now Max Planck Institute for Heart and Lung Research, Bad Nauheim, Germany) and Horst Geiger (Max Planck Institute for Biology, Tübingen, Germany) for help with blastula transplantations; Christian Feldhaus and Aurora Panzera (BioOptics Facility, Max Planck Institute for Biology) for help with confocal microscopy; Veronika Altmannova and Dorota Rousova (Friedrich Miescher Laboratory, Tübingen, Germany) for help with protein purification; and Silke Geiger-Rudolph, Roberta Occhinegro and Reinhard Albrecht for excellent technical assistance (Max Planck Institute for Biology, Tübingen, Germany). The AlphaFold-Multimer model was generated using the BMBF-funded de.NBI Cloud within the German Network for Bioinformatics Infrastructure (de.NBI) (031A532B, 031A533A, 031A533B, 031A534A, 031A535A, 031A537A, 031A537B, 031A537C, 031A537D, 031A538A).

The peer review history is available online at https://journals.biologists.com/dev/lookup/doi/10.1242/dev.201627.reviewer-comments.pdf.