A lineage-specific transcription profiling of genes downstream of fibroblast growth factor signaling (FGF), was used to find novel players involved in pigment cell formation (Racioppi et al., 2014). Among these, we found a Kelch-like (Klhl) gene family member (Kyoto Hoya gene model KH.L84.23), which exhibited similar expression values of known PCP markers, including Tyr, Tyrp.a and Rab32/38 and that, by reciprocal BLASTs, is similar to vertebrate Kelch-like 21 (Klhl21) and Kelch-like 30 (Klhl30).

To shed light on the evolutionary origins of the Klhl21/30 gene and to gain insights into the poorly studied Klhl family in metazoans, we performed an evolutionary analysis of these proteins in chordates. First, we analyzed their domain organization in PROSITE (de Castro et al., 2006) and Ensembl databases. This revealed the presence of one BTB/POZ domain, one BACK domain and five Kelch repeats (Supplementary Figure S1) in both C. robusta Klhl21/30 and Homo sapiens KLHL21 (Dhanoa et al., 2013). To confirm the orthology of C. robusta Klhl21/30 (KH.L84.23) and study the evolution of Kelch-like genes in chordates, we performed a phylogenetic survey employing a manually curated database (Supplementary Table S1). Our phylogenetic reconstruction included 118 Kelch-like protein sequences from the cephalochordate amphioxus Branchiostoma belcheri, tunicates Oikopleura dioica and Ciona robusta, and human Homo sapiens as representative of vertebrates (Figure 1). Sequences with high degree of molecular divergence have been excluded from phylogeny and listed in the Supplementary Table S2. In total, we identified 25 Klhl subfamilies, conserved throughout chordate evolution from cephalochordates to vertebrates, with a unique exception (Klhl5L) that is tunicate-specific (Figure 1). Genome search and phylogenetic tree defined the complete Kelch-like toolkit of B. belcheri (37), C. robusta (32), O. dioica (14), H. sapiens (42): this confirms the previous survey performed in human (Dhanoa et al., 2013) and represents the first step into understanding Kelch-like evolution in all three chordate subphyla. Though the number of Kelch-like genes in the amphioxus genome has been shaped by various amphioxus-specific duplications (as in the case of Klhl42 duplications), our survey is coherent with relative genomic stasis of this slow-evolving branch (Putnam et al., 2008). In contrast, C. robusta retained 32 Klhl genes and 16 subfamilies, whilst O. dioica maintained 14 genes and 10 subfamilies. This is coherent to the tendency of the fast-evolving Oikopleura dioica to lose large portions of gene families (Berná and Alvarez-Valin, 2014; Albalat and Cañestro, 2016; Martí-Solans et al., 2016; Coppola et al., 2019). Otherwise, we also detected different tunicate-specific duplications of certain Kelch-like genes in either species (Figure 1).

In particular, our phylogenetic tree clustered Ciona KH.L84.23 with one O. dioica protein, two B. belcheri proteins (deriving from a local duplication) and H. sapiens KLHL21 and KLHL30, forming a protein class with robust support that we named Klhl21/30 (orange box; Figure 1). The presence of two related genes in human speaks in favor of a common origin for vertebrate Klhl21 and Klhl30, which could have derived from a local duplication or a whole-genome duplication event in the vertebrate ancestor (Abi-Rached et al., 2002; Dehal and Boore, 2005). In contrast, invertebrate chordates possess a single Klhl21/30 gene, similar to other invertebrates such as the nematode Caenorhabditis elegans and the sea urchin Strongylocentrotus purpuratus (data not shown). Remarkably, Klhl21/30 is one of the few Kelch-like genes not lost by O. dioica and therefore could be functionally relevant for conserved developmental processes. Moreover, we analyzed the Klhl21/30 locus in several available tunicate genomes, finding high degree of local synteny surrounding Klhl21/30 in C. robusta, C. savignyi, Phallusia mammillata, and Halocynthia roretzi (Supplementary Figure S2). Our survey prompted us to describe an ancestral cluster of 11 genes in tunicates, nearly all conserved between C. robusta and P. mammillata. Interestingly, all surveyed species exhibited linkage between Klhl21/30 and Nph4 (nephrocystin-4) genes, which is also observed in human between KLHL21 and NPH4. Conversely, this linkage has been lost for KLHL30 (Supplementary Figure S3). The comparison of Klhl21 and Klhl30 genome environment of species as coelacanth Latimeria chalumnae, spotted gar Lepisosteus oculatus, frog Xenopus tropicalis and human, demonstrated the conservation of flanking genes during gnathostome evolution (Supplementary Figure S3). Moreover, the presence of orthologous genes on both the surveyed chromosomal regions, strongly indicated Klhl21 and Klhl30 as paralogs deriving from one of the events of whole-genome duplication (WGDs) occurred at the root of vertebrate radiation (Dehal and Boore, 2005).

In sum, we have reported here the first evolutionary study dedicated to Kelch-like family in chordates, focused on the seemingly indispensable Klhl21/30 subfamily, which was present in all chordate species examined.

The analysis of the spatio-temporal expression pattern of Klhl21/30, by whole mount in situ hybridization (WISH) in embryos at different developmental stages (Figures 2A–D), showed that during C. robusta embryogenesis, Klhl21/30 is expressed from initial tailbud stage onward in the pigment cell lineage (Figures 2A–D). In particular, at initial tailbud stage it is expressed in two cells (Figures 2A,B), possibly the ones originating the otolith and ocellus pigmented cells of the larva. Later in development, from the middle tailbud stage onward, Klhl21/30 transcript became restricted specifically to one cell (Figures 2C,D). Double WISH using the pigment cell marker, Tyrp1/2a, was performed to better define PCP-specific expression. Detailed and thorough analysis of the expression domain proved that Klhl21/30 is expressed specifically in the a11.193 pair, i.e. the most posterior cells marked by Tyrp1/2a (Haupaix et al., 2014; Racioppi et al., 2014) (white arrow in Figures 2E–F”’) later, starting from middle tailbud stage, the Klhl21/30 expression became restricted to the otolith precursor (white arrow in Figures 2G–G”’), while it disappears from the ocellus pigment cell precursor. Thus, Klhl21/30 represents the earliest-expressed otolith specific gene described so far. Its expression is maintained in the otolith pigment cell precursor while it is downregulated in the ocellus pigment cell precursor, even before expression of the βγ-crystallin, which is expressed specifically in the otolith starting from the larval stage (Shimeld et al., 2005). Klhl21/30 expression in PCPs is consistent with microarray data (Racioppi et al., 2014), in which Klhl21/30 expression was observed to decrease between 8 and 12 hpf. Furthermore, its expression was shown to be dependent on the FGF signaling, which is required for PCPs specification. Coherently, in situ hybridizations performed using Klhl21/30 probe on tailbud embryos electroporated with a construct, Tyrp > dnFGFr, able to drive a dominant-negative form of the FGF receptor in the PCPs (Squarzoni et al., 2011), showed a sharp decrease in the number of tailbud embryos expressing Klhl21/30 in PCPs compared to control (Supplementary Figure S4). Altogether, these data identified Klhl21/30 as the earliest marker of the Ciona otolith and the first Klhl member described in tunicates.

Since Klhl21/30 was identified as a specific marker for PCPs as well as the earliest otolith marker in Ciona robusta, we sought to study its transcriptional regulation to understand the regulatory logics underlying differentiation of these structures (Figure 3). In order to discover the cis-regulatory element involved in the control of Klhl21/30 expression, we took advantage of the genome conservation tracks in the ANISEED genome browser (Brozovic et al., 2018) to select a 949-bp non-coding region upstream Klhl21/30 (“KlA”, −1044 to −95 from the start codon), a region presenting conserved peaks with the sibling species C. savignyi (Figure 3A). This KlA fragment was cloned in a vector containing a GFP reporter gene downstream of a human β-globin minimal promoter (Zeller et al., 2006). Once electroporated the KlA > GFP reporter plasmid into Ciona embryos (Corbo et al., 1997), we detected otolith-specific GFP fluorescence in 70% of larvae at stage 26 (Hotta et al., 2007; Figures 3B,C). By WISH using GFP probe, we registered transcription of KlA > GFP at the middle tailbud stage in one cell of the sensory vesicle, possibly corresponding otolith precursor (Supplementary Figure S5). This indicates that the absence of GFP fluorescence until the larval stage is due to a delay in GFP maturation and accumulation. Our reporter plasmid thus recapitulates the endogenous expression of Klhl21/30 at least in the otolith, confirming the restriction of its expression to the otolith precursor in the late embryogenesis. To dissect the regulatory logics underlying the expression of Klhl21/30, we focused on two smaller, highly conserved regions within the KlA sequence: KlB (385 bp long, −918 to −553) and KlC (441 bp long, −536 to −95). While KlC did not drive GFP expression, 43% of larvae electroporated with KlB > GFP recapitulated the strong GFP signal in otolith, suggesting that KlB fragment retains the minimal regulatory information necessary to drive Klhl21/30 expression (Figure 3C). When we divided KlB roughly into two halves, KlD (199 bp long, −918 to −719) and KlE (166 bp long, −719 to −553), neither was sufficient to drive GFP expression (Figure 3C). However, we found that two smaller fragments centered around KlB, which are KlF (257 bp long, −810 to −553) and KlG (321 bp long, −874 to −553), drove GFP expression in 10 and 22% of larvae, respectively (Figure 3C). These fragments therefore represent the minimal cis-regulatory elements sufficient to recapitulate Klhl21/30 expression. We hypothesized that the KlB region contains transcription factors binding sites (TFBS) crucial for the sustained activation of Klhl21/30 in pigment cell precursors. Therefore, we searched in this fragment putative TFBS using Genomatix software and CIS-BP database with C. intestinalis motifs (Weirauch et al., 2014) and these results have been compared to JASPAR database (Khan et al., 2018; Figure 4). Among several predicted binding sites, we focused our attention on two bHLH-binding motifs. In Ciona 44 bHLH genes have been found (Satou et al., 2003) and, for most of them, expression pattern during development has been described (ANISEED database). We selected Mitf as possible factor binding the bHLH sites in the Klhl21/30 regulatory region (yellow in Figure 4A) because of involvement of Mitf in eye development and in the activation of melanogenic markers Tyr and Tyrp/DCT in vertebrates (Goding, 2000; Curran et al., 2010) and its specific expression in the pigment cell lineage of Ciona (Curran et al., 2010; Abitua et al., 2012) and Halocynthia roretzi (Yajima et al., 2003). We demonstrated that Mitf is co-expressed with Klhl21/30 in both otolith and ocellus pigmented cells precursors at early tailbud stage (Figures 4B–B”’) and that both became restricted to otolith precursor starting from middle tailbud stage (Figures 4C–C”’). Previous experiments showed also that Mitf is downregulated upon perturbation of FGF signaling (Racioppi et al., 2014). To test the relevance of these two putative sites for Mitf (CACGTG), we mutated them individually in the KlB reporter, naming these mutant constructs K_Mitf_mut1 and K_Mitf_mut2 (Figure 4E). When each of these constructs was electroporated, the percentage of larvae showing GFP signal in the otolith was close to zero (2.5 and 3%, respectively, Figure 4E). These data strongly support these binding sites for Mitf as involved in KlB activity, leading to hypothesize that Mitf activates Klhl21/30 in Ciona robusta.

We also identified a well-supported binding site in KlB for Dmrt (TTACAT, violet in Figure 4A), a transcription factor expressed in the early a-line neural plate (Wagner and Levine, 2012). C. savignyi Dmrt mutants present abnormalities in the development of the larval sensory vesicle (Tresser et al., 2010) and the expression of the C. robusta ortholog is under early FGF control (Imai et al., 2006). The KlB fragment harboring a mutated Dmrt site (K_Dmrt_mut) drives GFP expression only in 20% of electroporated larvae, roughly half of the activity of wild-type KlB (Figure 4E). Additionally, two binding sites attributed to the homeodomain transcription factors caught our attention (ATTA, blue in Figure 4A). Among TFs able to bind these sites we focused on Msx because of its early expression in C. robusta pigment cell precursors (Aniello et al., 1999; Russo et al., 2004), and given that it is sharply down-regulated when FGF is blocked (Racioppi et al., 2014). We detected co-expression of Msx and Klhl21/30 in the PCPs at early tailbud stage (Figures 4D–D”’) As development proceeds, the expression of Msx in PCP decrease becoming excluded from most of the a9.49 derivatives in 12 hpf tailbud embryos as already described (Racioppi et al., 2014). This early expression suggests that Msx might be activating early Klhl21/30 transcription, even though Msx factors normally have been described to act as repressors in Ciona (Roure and Darras, 2016). We individually mutated each putative Msx binding site in KlB, resulting in two mutant constructs, which we termed K_Msx_mut1 and K_Msx_mut2. The mutated reporters caused a decrease in the percentage of GFP-expressing larvae of 24 and 19%, respectively, with respect to the wild-type (Figure 4E). Mutating both putative Msx sites simultaneously caused a strong reduction in GFP expression, observed in only 4% of larvae (Figure 4E, K_2xMsx_mut). Taken together, our data suggest that also Dmrt and Msx could be involved in the transcriptional activation of Klhl21/30 with a novel role in Ciona.

To test the potential role of Mitf, Dmrt and Msx as transcriptional activators of Klhl21/30, we used CRISPR/Cas9-mediated mutagenesis to knock out these factors in C. robusta embryos (Stolfi et al., 2014; Gandhi et al., 2017, 2018). To design primers to generate sgRNAs for Mitf, Dmrt, and Msx, we used CRISPOR v4.0 portal (Haeussler et al., 2016), which takes in account the potential single-nucleotide polymorphisms (SNPs) and the specificity of each sgRNA and provides their predicted “efficacy” (predicted efficiency of their ability to induce specific DSBs) by the Fusi/Doench algorithm (Doench et al., 2016). Using a series of candidate sgRNAs, we targeted different parts of the coding regions of the selected genes (Supplementary Table S3), synthesized U6 > sgRNA cassettes through One-Step Overlap PCR (OSO-PCR) and validated them by peakshift analysis (Gandhi et al., 2017; see section “Materials and Methods” for details; Supplementary Table S3).

We co-electroporated each gene-specific U6 > sgRNA cassettes, (targeting either Mitf, Dmrt, or Msx) together with Fog > Cas9, which drives Cas9 expression in the a-line blastomeres through the Fog promoter (Rothbächer et al., 2007), and the KlB > GFP reporter plasmid to verify the effect of the experiments on Klhl21/30 isolated regulatory region. As a control, we co-electroporated Fog > Cas9 and KlB > GFP with an sgRNA cassette targeting the mesoderm-specific transcription factor Mesp.5 (Davidson et al., 2005), which is not expressed in the pigment cell lineage. While the control sample had a proportion of GFP-expressing larvae comparable to electroporation with KlB > GFP alone (43% Figure 5A, compared with Figure 2C), we detected a loss of GFP fluorescence in larvae electroporated with sgRNA cassettes targeting the three putative regulators of Klhl21/30 (Figure 5A). When Mitf was knocked out, KlB > GFP was expressed in only 3% of larvae, a sharp decrease compared to the control condition. Moreover, the knockout of Dmrt and Msx also reduced GFP expression, with only 10 and 4% of larvae showing the KlB > GFP expression, respectively (Figure 5A). Because of the strong effect of Mitf loss-of-function on KlB > GFP activity in C. robusta, and the evolutionary conservation of Mitf function in pigment cell development (Goding, 2000; Levy et al., 2006), we performed a complementary Mitf gain-of-function experiment. We co-electroporated KlB > GFP along with Ebf −2.6 kb > H2B:mCherry and Ebf −2.6 kb/ + 15 STOP > Mitf to overexpress Mitf in Ebf + neural progenitors during development. These resulted in ectopic KlB > GFP expression in Ebf + cells in the central nervous system (Figures 5B–B”), suggesting that Mitf drive the transcription of Klhl21/30 in these cells.

Thus, these data strongly suggest that, downstream the FGF signaling, Klhl21/30 is under the control of Mitf, with Dmrt and Msx acting as important co-activators.

Kelch-like protein sequences from vertebrate Homo sapiens were used as queries in BLASTp and tBLASTn searches in ANISEED (Brozovic et al., 2018), NCBI or Ensembl genome databases of selected species. Orthology was initially assessed by reciprocal best blast hit (RBBH) approach and supported by phylogenetic analyses. Evolutionary reconstructions were performed using ML inferences calculated with PhyML v3.0 and automatic modality of selection of substitution model (Guindon et al., 2010) using protein alignments generated with MUSCLE (Edgar, 2004) and ClustalX (Larkin et al., 2007) programs. Alignment of Supplementary Figure S1 has been carried out employing ClustalX. Full protein sequences were used in our analysis. Accession numbers and names for phylogenetic tree of Figure 1 are listed in Supplementary Table S1, while those excluded are encompassed in Supplementary Table S2. The analysis of synteny conservation (shown in Supplementary Figures S2, S3) was performed by employing ANISEED, Ensembl and Genomicus databases (Louis et al., 2015). To predict putative transcription factor binding sites (TFBS) in the surveyed cis-regulatory regions, we utilized the MatInspector module of the Genomatix Software Suite and CIS-BP employing a Ciona intestinalis (former name for Ciona robusta) DNA-binding-domain classes database (Weirauch et al., 2014). We also used the JASPAR database (Khan et al., 2018) to recognize the potential binding sites.

Adults of Ciona robusta were collected from the Gulf of Naples, or from San Diego, CA, United States, by M-REP. Gametes from many animals were gathered separately for in vitro cross-fertilization followed by dechorionation and electroporation as previously illustrated (Christiaen et al., 2009a; Racioppi et al., 2014). Electroporated plasmid amounts (e.g., 10 μg) were per 700 μl of total volume. Embryos were staged according to the developmental timeline shown in Hotta et al. (2007). To visualize GFP, embryos were fixed in MEM-FA (3.7% methanol-free formaldehyde, 0.1 M MOPS pH 7.4, 0.5 M NaCl, 2 mM MgSO4, 1 mM EGTA) for 30 min and washed several times in PBS-NH₄Cl and in PBS containing 0.05% Triton X-100. Each electroporation of Figures 3, 4 was carried out using 60 μg of plasmid. The statistical significance of electroporations associated to Figure 4 was validated using Fisher exact test. The statistical significance of electroporations of Supplementary Figure S4 was evaluated employing Chi-square test for trend.

Single and double in situ hybridization experiments were performed out essentially as described previously (Christiaen et al., 2009b; Racioppi et al., 2014), using DIG- and FLUO-labeled riboprobes, anti-DIG-POD and anti-FLUO-POD Fab fragments (Roche, Indianapolis, IN), and Tyramide Amplification Signal with Fluorescein (Perkin Elmer, MA). The antisense riboprobes were obtained from plasmids contained in the C. intestinalis gene collection release I: Klhl21/30 (KH2012:KH.L84.23, GC17e22), Tyrp1/2a (KH2012:KH.C8.537, GC31h05), Mitf (KH2012:KH.C10.106, GC28k08), Dmrt (KH2012:KH.S544.3, GC02f18), and Msx (KH2012:KH.C2.957, GC42h24) (Satou et al., 2002).

The cis-regulatory elements upstream Klhl21/30 were PCR-amplified from genomic DNA and their localization on upstream sequence is shown in Supplementary Figure S6. Insertion of the products into pSP72 vector containing GFP (Zeller et al., 2006) was carried out using TOPO-TA Cloning kit (Invitrogen). The QuickChange Site-Directed Mutagenesis Kit (Agilent) was employed to generate the mutations inside the putative binding sites (Mitf, Dmrt, Msx) identified in the sequence of the KlB element (Figure 4). All the oligos used for cloning the putative regulatory regions and for mutational experiments are listed in Supplementary Table S4.

Cas9 and sgRNA expression vectors were constructed or used as previously described (Stolfi et al., 2014; Gandhi et al., 2017). The predictive algorithm used for designing in vivo-transcribed sgRNAs has been the Fusi/Doench (Doench et al., 2016), available on CRISPOR portal (Haeussler et al., 2016). The target sequences have been selected not too close to the translational start, to have the higher impact on the function of protein of interest (Gandhi et al., 2017). DNA oligos used to generate sgRNAs are listed in Supplementary Table S3. One-step overlap PCR (OSO-PCR) was employed for the fast synthesis of a U6 > sgRNA cassettes through a single PCR reaction, performed using Platinum Pfx Polymerase (Invitrogen). The products were cloned into pCESA plasmid using Gibson Assembly Cloning Kit (NEB). Before the electroporation, the products were purified utilizing AMPure XP (Agencourt).

Mitf, Msx, and Dmrt sgRNAs were validated by PCR amplification and Sanger sequencing of the targeted region as previously described (Gandhi et al., 2018). Pooled larvae were lysed for 30 min in 180 μl buffer + 5 μl of Proteinase K from QIAamp DNA Micro Kit (Qiagen) and eluted in 20 μl of water. Approximately 200 ng/μl of genomic DNA extracted from hatched larvae was used for “touchdown” PCRs with Platinum Pfx Polymerase (Invitrogen), as described in Gandhi et al. (2017). The genomic oligos (“peakshift oligos”) employed for “touchdown” PCRs were selected at 150–500 bp away from the primer (Supplementary Table S3), to ensure a proper peakshift. Electroporations were performed as single biological replicates. Electroporation mix recipes can be found in the Supplementary Figure S7. The statistical significance of electroporations associated to Figure 5 was validated using Fisher exact test. Mitf overexpression construct (Ebf-2.6kb/ + 15 STOP > Mitf) was designed as described in Supplementary Figure S7: STOP indicates the presence of stop codon between the ATG of Ebf and the NotI site where Mitf was cloned (no amino acid sequence of Ebf is included in Mitf protein). Images were captured using Confocal ZEISS LSM 700 or ZEISS Apotome.2 compound microscopes.