pmar1/phb homeobox genes and the evolution of the double-negative gate for endomesoderm specification in echinoderms

The rewiring of gene regulatory networks (GRNs) is essential for morphological evolution. However, not all modifications of a GRN alter the associated morphology. In some cases, different upstream developmental pathways result in a conserved developmental output, a phenomenon known as ‘developmental system drift’ (True and Haag, 2001; Kalinka and Tomancak, 2012), which is supported by the studies focusing on the ‘hourglass model’ (e.g. Irie and Kuratani, 2014; Hu et al., 2017). However, it is still a mystery why early developmental processes are more variable than later processes. Because morphological evolution is tightly linked with the rewiring of GRNs, a better understanding of the flexibility of GRNs would contribute to a deeper understanding of morphological evolution.

One of the factors that may contribute to early developmental diversification is lineage-specific genes, i.e. genes found in the particular lineage. The regulation of early developmental process by lineage-specific genes has been frequently observed; examples of these genes include bicoid in Drosophila (Frohnhöfer and Nüsslein-Volhard, 1986), SPILE in spiralians (Paps et al., 2015; Morino et al., 2017), Siamois in amphibians (Lemaire et al., 1995), dharma in zebrafish (Yamanaka et al., 1998) and CRX-related genes in mammals (Töhönen et al., 2015; Maeso et al., 2016). However, it is not fully understood how these lineage-specific genes acquired their functions in early embryogenesis, except in the case of the dipteran bicoid gene (e.g. Stauber et al., 1999; Kotkamp et al., 2010; Liu et al., 2018). Therefore, we focused on the sea urchin (echinoid) pmar1/micro1 genes, which are key upstream factors involved in endomesoderm specification during embryogenesis, but have been identified in only one of the two echinoid lineages.

The class Echinoidea is classified into two subclasses: Cidaroida and Euechinoida (see Fig. 8A). These are estimated to have diverged 268.8 million years ago (Thompson et al., 2015). Most developmental research in this group has been performed using species of Euechinoida, such as Strongylocentrotus purpuratus, Lytechinus variegatus, Paracentrotus lividus and Hemicentrotus pulcherrimus (McClay, 2011). There are slight differences in early developmental processes between cidaroids and euechinoids; e.g. whereas mesodermal skeletogenic cells ingress before gastrulation in euechinoids, cidaroid skeletogenic mesenchyme cells delaminate after gastrulation (for further details, see Yamazaki et al., 2014). Despite these apparent differences, skeletogenic cells differentiate from cells in the vegetal region in cidaroids, resulting in pluteus larvae with a similar morphology. pmar1/micro1 (hereafter, pmar1) genes have been isolated from a variety of euechinoids, including the sand dollar and heart urchin (Di Bernardo et al., 1995; Kitamura et al., 2002; Oliveri et al., 2002; Ettensohn et al., 2007; Yamazaki et al., 2010; Yamazaki and Minokawa, 2015). However, no pmar1 gene ortholog has been identified in cidaroids (Yamazaki et al., 2012; Erkenbrack and Davidson, 2015; Dylus et al., 2016) or in other echinoderm species, such as sea cucumber (holothuroid) (McCauley et al., 2012; Thompson et al., 2017), brittle star (ophiuroid) (Dylus et al., 2016) or starfish (asteroid) (McCauley et al., 2010). Dylus et al. (2016) reported the pplx gene as a pmar1-related gene in brittle star, which we will discuss later.

Euechinoids possess multiple copies of pmar1 genes, which are tandemly arrayed in the genome. Each of these genes encodes a transcription factor with a paired-type homeodomain that functions as a transcriptional repressor with engrailed homology region 1-like (eh1-like) motifs. Although unusual expression of a P. lividus ortholog (hbox12) in the ectodermal region has been reported (Di Bernardo et al., 1995), all of the other pmar1 genes are transiently expressed in the micromere-skeletogenic cell lineage at the vegetal pole from the 16-cell stage to the mid-blastula stage as far as examined (Kitamura et al., 2002; Oliveri et al., 2002; Yamazaki et al., 2010; Yamazaki and Minokawa, 2015); i.e. pmar1 is one of the earliest zygotically expressed genes. In embryos injected with pmar1 mRNA, almost all cells develop into skeletogenic cells. A comprehensive GRN has been established for the skeletogenic cell lineage of euechinoids; pmar1 is the most upstream zygotic factor in the network and is directly activated by maternal β-catenin. Pmar1 promotes the specification of skeletogenic cells by repressing a hairy gene, hesC, which represses downstream regulatory genes involved in skeletogenesis, such as alx1, tbr, ets1 and delta (Revilla-i-Domingo et al., 2007; see reviews by Oliveri et al., 2008; Minokawa, 2017; Shashikant et al., 2018; see Fig. 8B, euechinoid). However, recent studies on two cidaroids, Prionocidaris baculosa and Eucidaris tribuloides, indicated that HesC does not repress skeletogenic regulatory gene orthologs during the early phase of endomesoderm specification; in the hesC knockdown cidaroid embryos, expression of alx1, tbr and ets1 was not affected at the mid-blastula stage (Yamazaki et al., 2014), although HesC seems to exhibit a regulatory function related to alx1 expression at the relatively later stage (Erkenbrack and Davidson, 2015), suggesting the absence of the typical double-negative gate of Pmar1 and HesC in cidaroids. Considering the recently proposed hypothesis that the larval skeletogenic cells arose in the common ancestor of eleutherozoan echinoderms (all echinoderms except crinoids) (Erkenbrack and Thompson, 2019), the Pmar1-HesC double-negative gate was likely established after the acquisition of a larval skeleton; i.e. this gate likely evolved independently of one of the novel morphologies acquired during echinoderm evolution.

The purpose of this study is to reveal the evolutionary history of the establishment of the double-negative gate of Pmar1 and HesC in echinoderms, which importantly occurred without changing the expression pattern of key endomesodermal developmental genes, such as alx1, ets1 and delta. In the course of the comparative analysis of the GRNs among echinoderm species, we unexpectedly identified cidaroid pmar1 gene orthologs through temporal RNA-sequencing (RNA-seq) analysis. This prompted us to further examine the pmar1-related genes of other echinoderms. Based on the results of expression and functional analyses of these genes, we discuss the evolution of the endomesoderm gene network in echinoderms.

As endomesoderm regulatory genes, alx1 and ets1 are expressed in skeletogenic cells in cidaroids (Yamazaki et al., 2014; Erkenbrack and Davidson, 2015). However, the upstream regulatory mechanism has not yet been revealed. To screen candidate genes responsible for regulating the onset of skeletogenic regulatory gene expression in the cidaroid, we performed RNA-seq analysis using embryos of the cidaroid P. baculosa at the two-cell (2 hours postfertilization; h), 16-cell (4 h), ∼64-cell (6 h), ∼240-cell (10 h) and ∼500-cell (14 h) stages. In these embryos, the zygotic expression of Pb-alx1 was first observed at 10 h. Based on changes in the obtained FPKM (fragments per kilobase of transcript per million mapped reads) values, 43 candidate transcription factor genes that were activated before or simultaneously with alx1 were selected (Table S1). The detailed criteria for selecting the candidate genes are described in the Materials and Methods section. As shown in Fig. S1, we examined the spatial expression patterns of the candidate genes and found some genes showing mesoderm-specific expression, such as kruppel-like1 (krl-like1) and kruppel-like2 (krl-like2), the euechinoid ortholog of which is not required for skeletogenic cell specification (Yamazaki et al., 2008). In addition, we identified a sequence that showed remarkable similarity to S. purpuratus pmar1c. This was unexpected because pmar1 was thought to have emerged in the common ancestor of euechinoids (Erkenbrack and Davidson, 2015; Thompson et al., 2017). Its transient and relatively low expression may have hidden its existence in previous transcriptome data. We identified a similar sequence in another cidaroid, Eucidaris tribuloides, via a BLAST search using the P. baculosa sequence as a query against the E. tribuloides genome 1.0 sequence at EchinoBase (www.echinobase.org) (Kudtarkar and Cameron, 2017). The sequence was not found in the E. tribuloides transcriptome data obtained from EchinoBase.

Given the existence of the pmar1 gene in cidaroids, before moving on to the functional analyses of cidaroid pmar1, we examined the molecular evolutionary history of pmar1. In the euechinoid S. purpuratus, the pmar1-related phb1 gene was identified as a paired-class homeobox gene by Howard-Ashby et al. (2006), but a detailed analysis has not been performed. Dylus et al. (2016) reported that phb1-related sequences exist in cidaroids, other echinoderms and acorn worms, and demonstrated that phb1-like genes and pmar1 are closely related but are recognized as distinct classes of paired homeobox genes according to phylogenetic analysis. Based on a BLAST search of the assembled transcriptome and genomic sequences (details are described in the Materials and Methods section), we identified 1, 6, 2 and 1 pmar1/phb1-like sequences in sea cucumbers, brittle stars, starfishes and feather stars (crinoids), respectively.

To evaluate the relationships between these obtained sequences and euechinoid Pmar1, we performed phylogenetic analysis using deduced homeodomain (HD) sequences with several classes of other paired-type genes according to the method of Dylus et al. (2016) (Fig. 1A). Our results showed a monophyletic clade, including the cidaroid pmar1-like sequences and euechinoid pmar1 sequences with a high support value (93%). Accordingly, we designated the genes obtained from P. baculosa and E. tribuloides as Pb-pmar1 and Et-pmar1, respectively. In addition, the monophyly of the clade including these sequences, phb1, and pmar1 was supported with significant values, suggesting that these genes are paralogs. However, our analyses did not resolve the relationships between these genes, which is often the case when performing phylogenetic analyses with only 60 amino acids of a homeodomain. In particular, the long branches of pmar1 genes made the tree less resolvable. Some sequences from brittle stars form a clade with sea urchin pmar1 genes, although this clade is not supported by sufficient values, possibly owing to artificial long branch attraction. Indeed, pmar1/phb1-related genes show extensive gene duplication in the brittle star lineage, and some of the brittle star genes present an accelerated substitution rate, as reflected by their long branches (such as Ak-phbC and Afi-pplx). We designated the identified genes from nonurchin echinoderms as follows: for sea cucumbers, phb; for the brittle star Amphipholis kochii, phbA to phbF; for starfishes, phbA and phbB; and for the feather star Oxycomanthus japonicus, phb.

The euechinoid Pmar1 proteins commonly contain a HD in the N terminus and two engrailed homology region 1-like (eh1-like) motifs in the C terminus (Fig. 1B). The eh1 motif is a repression motif that interacts with the co-repressor groucho (Copley, 2005). Our previous study demonstrated that these eh1-like motifs are responsible for the repressive function of Pmar1 (Micro1) (Yamazaki et al., 2009). The eh1-like motifs were found only in sea urchin Pmar1 and Phb1, including those of cidaroids, but not in any of the Phb sequences of nonurchin echinoderms (Fig. 1B,C). The cidaroid Pmar1 contains only one eh1-like motif, which is highly conserved compared with euechinoid Pmar1 (Fig. 1B,C).

The paired-type HDs are classified into three subclasses according to the 50th amino acid (glutamine, Q; lysine, K; and serine, S), which is involved in binding sequence preference (Wilson et al., 1993). HDs with a Q50 are shared by euechinoid Pmar1 and Phb1, whereas pplx from another brittle star encodes a HD with an irregular 50th amino acid, histidine (Dylus et al., 2016). Among the Pmar1/Phb sequences identified in this study, one Phb from the brittle star A. kochii (Ak-PhbC) includes a K50, whereas the others share a Q50.

In addition, similar to euechinoid pmar1 (micro1) genes (Nishimura et al., 2004; Ettensohn et al., 2007; Cavalieri et al., 2017), the cidaroid pmar1 gene is extensively duplicated in the genome (Fig. 1D). In the genomic sequence of E. tribuloides obtained from EchinoBase, we found two long scaffolds containing sequences similar to the pmar1 gene (scaffolds 7904 and 10027). Both scaffolds include four copies of pmar1-related sequences, whose orientations are different, and one short pmar1-like sequence without a start codon is present at the corresponding position in scaffold 7904. In summary, we identified multicopy pmar1 genes in cidaroids and pmar1/phb1-related phb genes from nonurchin echinoderms.

The euechinoid pmar1 (micro1) genes are commonly expressed transiently in the micromeres in the 16-cell stage and the descendant skeletogenic cells at the vegetal pole (e.g. Yamazaki et al., 2010; Yamazaki and Minokawa, 2015). To examine whether the expression patterns of P. baculosa pmar1 are similar to those of euechinoid pmar1, we performed expression analysis through quantitative PCR (qPCR) and whole-mount in situ hybridization. Consistent with the RNA-seq analysis (Table S1), qPCR demonstrated the transient activation of Pb-pmar1 during early stages; the expression level of Pb-pmar1 reached a peak during the ∼64-cell (6 h) and ∼120-cell stages (8 h) (Fig. 2A). Whole-mount in situ hybridization showed the vegetal expression of Pb-pmar1. Almost no embryos (1/11) at the 16-cell stage (4 h) showed the whole-mount in situ hybridization signal (Fig. 2B), although a low level of the transcript was detected by qPCR. At the 32-cell stage (5 h), a subset of embryos (11/24) showed a signal in the smaller blastomeres (arrowheads in Fig. 2C) located in the vicinity of the vegetal pole (Yamazaki et al., 2012). The pmar1-expressing cells may differentiate into larval skeletogenic cells because previous lineage-tracing experiments in another cidaroid, E. tribuloides, demonstrated that smaller vegetal blastomeres develop into skeletogenic cells (Wray and McClay, 1988). The signal continued to be detected at 8 h but not at 14 h (Fig. 2D,E). Thereafter, the transcripts remained undetectable up to 36 h by qPCR (Fig. 2A). This expression pattern is very similar to that of euechinoid pmar1 genes.

We also examined the expression patterns of phb genes in the embryos of the starfish Patiria pectinifera (Fig. 2F-O), the sea cucumber Apostichopus japonicus (Fig. 2P, Fig. S2) and the brittle star A. kochii (Fig. 2Q-S, Fig. S2). The expression of the starfish P. pectinifera phbA and phbB genes (Ppe-phbA and Ppe-phbB) was first detected at the 200- to 400-cell stage (6 h) (Fig. 2G,L), and both genes showed expression at the vegetal pole of the hatched blastula (Fig. 2H,I,M,N). At the mid-gastrula stage (24 h), phbA expression was detected in the region encircling the blastopore (Fig. 2J), which seems to be endoderm lineage, but phbB expression was no longer detected (Fig. 2O). The sea cucumber and brittle star phb genes showed similar expression patterns: their expression was first detected at the cleavage stage or early blastula stage, and was subsequently maintained in either the mesoderm or endoderm lineage of cells (Fig. 2P-S, Fig. S2). Thus, our analyses showed that all phb genes examined are expressed in the endomesoderm region at the vegetal pole. Some of the phb genes (Ppe-phbA, Aj-phb and Ak-phbA/C) were detected in the presumptive endoderm region during relatively later stage, which is clearly distinct from the expression patterns of euechinoid pmar1 genes.

The result of the above expression analysis suggests that the function of cidaroid Pmar1 is similar to that of euechinoid Pmar1. To examine whether Pmar1 also controls skeletogenic cell specification in cidaroid embryos, we performed overexpression analysis using P. baculosa embryos (Fig. 3). The phenotype of pmar1 mRNA-injected cidaroid embryos was not identical to that observed in euechinoids (i.e. fate conversion to the skeletogenic cell phenotype in almost all cells), although excess mesoderm cell differentiation was observed. When the control embryos developed into elongated swimming blastulae, Pmar1-overexpressing embryos showed a rather spherical morphology (16 h; Fig. 3A,E). During gastrulation, a broader area of the vegetal side invaginated in overexpressing embryos (Fig. 3B,F). To examine the effects of skeletogenic and nonskeletogenic mesenchyme cell formation, we counted the mesenchyme cell number at the late gastrula stage (40 h) (Fig. 3C,D,G-J). The number of presumptive skeletogenic cells expressing the skeletogenic cell marker P4 did not increase significantly in embryos overexpressing Pmar1 (Fig. 3I). On the other hand, the number of mesenchyme cells showed a significant increase in Pmar1-overexpressing embryos (Fig. 3J).

The observed phenotypic difference suggests that the regulatory function of Pmar1 differs between cidaroids and euechinoids. Thus, we examined the effect of cidaroid Pmar1 on skeletogenic gene orthologs by assessing the expression of hesC, alx1, tbr, ets1 and delta, as well as the putative endoderm regulatory gene foxA (Erkenbrack et al., 2018) by whole-mount in situ hybridization in the blastula stage (Fig. 3K-N). Pb-pmar1 mRNA-injected embryos showed global activation of hesC (Fig. 3Ka,La,Ma,Na). This effect on hesC expression was opposite to that observed in euechinoids, in which overexpression of Pmar1 suppresses the expression of hesC (Revilla-i-Domingo et al., 2007). On the other hand, the effects on other genes were similar to those found in euechinoids. The expression of tbr, ets1 and delta was expanded throughout the entire embryo (Fig. 3Kc-e,Lc-e,Mc-e,Nc-e). In contrast, foxA expression disappeared in Pmar1-overexpressing embryos (Fig. 3Kf,Lf,Mf,Nf). The expression of alx1 was also moderately expanded in Pmar1-overexpressing embryos at the earlier stage (Fig. 3Kb,Lb), but ectopic expression was not observed at the later stage (Fig. 3Mb,Nb), which is consistent with the lack of an increase in skeletogenic cell numbers in Pmar1-overexpressing embryos at the gastrula stage. It should be noted that the expression pattern of alx1 was different from those of the other genes in the cidaroid embryos. During the blastula stage, the expression of tbr, ets1 and delta was detected uniformly in the vegetal region of normal P. baculosa embryos (Fig. S3B-D,F-H), whereas patchy expression of alx1 was frequently observed (Fig. S3A,E). This suggests that an additional mechanism, possibly related to a molecular or mechanical bias in the earlier stage, exists for alx1 regulation in this species. In summary, cidaroid Pmar1 promotes the activation of endomesodermal/skeletogenic regulatory gene orthologs (alx1, tbr, ets1 and delta) but not by repressing hesC.

Because the GRN downstream of Pmar1 is likely to differ between cidaroids and euechinoids, we asked whether any biochemical features of Pmar1 have changed during sea urchin evolution. To address this issue, we examined whether cidaroid Pmar1 can perform similar functions in euechinoid embryos. We injected the cidaroid P. baculosa pmar1 mRNA into the euechinoid H. pulcherrimus and observed the resultant phenotype (Fig. 4). The phenotype of embryos overexpressing Pb-Pmar1 was identical to that of embryos overexpressing the euechinoid Pmar1. The Pmar1-overexpressing embryos showed global expression of alx1 and ets1 at the hatched blastula stage (Fig. 4F-H), whereas these genes were expressed specifically in the skeletogenic cell region of control embryos (Fig. 4A-C). Until the gastrula stage, almost all cells of the Pmar1-overexpressing embryos developed into mesodermal mesenchyme cells that expressed the skeletogenic cell marker P4 (Fig. 4D,E,I,J). These observations suggest that cidaroid Pmar1 exhibits euechinoid Pmar1-like activity in euechinoid embryos, probably through the repression of hesC, and may function as a repressor.

To estimate the function of the ancestral genes of pmar1 and phb, we further examined the functions of two phb genes (phbA and phbB) in the starfish P. pectinifera (Fig. 5). We performed knockdown and overexpression analyses using morpholino antisense-oligos (MOs) and synthesized mRNAs, respectively. In embryos injected with phbA and/or phbB MOs, gastrulation and mesenchyme formation were inhibited (Fig. 5E-H) compared with these processes in control embryos (Fig. 5A-D,M). At the gastrula stage during normal development, the archenteron is subdivided into two regions: the endoderm region, which shows alkaline phosphatase (AP) activity; and the AP-negative mesodermal region (Kuraishi and Osanai, 1994). All mesenchyme cells express the antigen of the MC5 monoclonal antibody (Hamanaka et al., 2011). In the knockdown embryos, AP activity was significantly reduced (Fig. 5C,G), and the total number of mesenchyme cells recognized by the MC5 antibody decreased (Fig. S4). Double-knockdown caused more-severe effects (see the detailed observations of archenteron and mesenchyme cell formation in Fig. S4), implying that the two phb genes function redundantly. In contrast, the embryos overexpressing PhbA and PhbB formed an enlarged AP-positive region and subsequently developed into an exogastrula (Fig. 5I-L). These observations suggest that, similar to echinoid Pmar1, the starfish Phb proteins are required for the formation of vegetal tissues.

Cidaroid Pmar1 leads to the activation of endomesoderm regulatory genes (e.g. alx1 and ets1) and hesC, as mentioned above. To evaluate the regulatory function of the starfish Phb proteins, we also analyzed the expression of endomesoderm regulatory gene orthologs in PhbA/B-perturbed P. pectinifera embryos. We examined the expression of hesC, ets1/2 (ets1 ortholog), tbr, delta and foxA in experimental embryos at the hatched blastula stage (12 h) (Fig. 5M-O). In the Phb-knockdown embryos, the expression of hesC, ets1/2, delta and foxA was significantly reduced at the vegetal pole (Fig. 5Na,b,d,e), whereas tbr expression was not affected (Fig. 5Nc). Conversely, in Phb-overexpressing embryos, hesC and delta expression was expanded throughout the embryos (Fig. 5Oa,d), which is similar to that in Pmar1-overexpressing cidaroid embryos (see Fig. 3). However, ets1/2 and tbr expression did not appear to be affected (Fig. 5Ob,c), and expansion of foxA expression was observed in the Phb-overexpressing starfish embryos (Fig. 5Oe), suggesting that the regulatory functions of starfish PhbA/B for foxA, ets1/2 and tbr are distinct from that of cidaroid Pmar1. Nonetheless, these observations suggest that, similar to cidaroid Pmar1, starfish PhbA/B leads to the activation of hesC, ets1/2 and delta, although an additional factor(s) is needed for ets1/2 expression. Based on these results, we suggest that starfish PhbA/B exhibit a regulatory function similar to that of cidaroid Pmar1.

To determine whether the biochemical activity of starfish Phb proteins is comparable with that of echinoid Pmar1, we examined the activity of starfish Phb when expressed in euechinoid embryos. We injected the P. pectinifera phbA and phbB mRNAs into the eggs of the euechinoid H. pulcherrimus (Fig. S5). When the control embryos developed into mesenchyme blastulae, the embryos overexpressing starfish PhbA showed a moderate increase in skeletogenic cells (Fig. S5A,B). In contrast, no obvious effects were observed in the embryos expressing starfish PhbB (Fig. S5C). These observations suggest that at least one starfish Phb protein can exhibit some degree of Pmar1-like activity in euechinoid embryos.

Because starfish Phbs have no typical eh1-like motifs (see Fig. 1B), we asked whether starfish Phb proteins function as repressors. We overexpressed the mRNAs encoding the two types of proteins: PhbA/B fused to the Drosophila Engrailed repression domain (EnR) or to the VP16 activation domain (VP16AD) (Fig. 6). The Phb proteins fused with the EnR domain caused phenotypes similar to those caused by the wild-type proteins (Fig. 6C,D). In contrast, overexpression of PhbA- or PhbB-VP16AD retarded the development of endomesodermal tissues (Fig. 6E-H). These results suggest that starfish PhbA and PhbB function as repressors similar to euechinoid Pmar1.

Our data indicate that, in contrast to euechinoid Pmar1, Pmar1/Phb leads to the activation of hesC in both cidaroid and starfish embryos, although this promotion probably occurs indirectly because both Pmar1 and Phb may function as repressors. Erkenbrack and Davidson (2015) demonstrated that a Delta signal is required for hesC expression in the vegetal embryo of another cidaroid, E. tribuloides. Here, we confirmed the same result in the cidaroid P. baculosa and the starfish P. pectinifera. To estimate the function of Delta signaling, we treated the embryos with N-[N-(3,5-difluorophenacetyl)-L-alanyl]-S-phenylglycine t-butyl ester (DAPT), which inhibits Notch signaling, according to Erkenbrack et al. (2018) (Fig. 7). In the DAPT-treated P. baculosa embryos, the vegetal expression of hesC disappeared, while tbr, ets1 and foxA expression was not affected (Fig. 7F,H-J). Because the experimental embryos showed a slightly increased signal intensity of alx1 (compare Fig. 7G with Fig. 7B), we counted the number of cells expressing alx1 (Fig. 7K). The average number increased moderately but significantly; the control and experimental embryos showed alx1 expression in 12.2±1.3 cells (n=13) and 14.5±2.1 cells (n=13) on average, respectively. We also found that the number of skeletogenic cells expressing P4 increased moderately in Delta signaling-deficient embryos (Fig. S6A-E), suggesting that P. baculosa HesC represses skeletogenic fate in the cells around the vegetal pole, as noted by Erkenbrack and Davidson (2015).

A significant decrease in hesC mRNA at the vegetal pole was also found in the DAPT-treated starfish embryos (Fig. 7L,O), whereas the expression of ets1/2 and foxA was not altered (Fig. 7M,N,P,Q). We also confirmed the results of Hinman and Davidson (2007), showing that depletion of Delta function resulted in the activation of ets1/2 in the starfish gastrula (compare Fig. S6I and Fig. S6N). In the DAPT-treated embryos, cell conversion into globular mesenchyme cells was observed in the whole upper region of the archenteron in the mid-gastrula stage (Fig. S6G,L). This implies that Delta-Notch signaling affects the genes responsible for the epithelial-mesenchymal transition of the mesoderm region in starfish embryos but may be regulated independently of HesC because HesC-knockdown in starfish embryos had no effect on mesoderm differentiation, as further discussed. In summary, the above data suggest that cidaroid and starfish hesC genes are regulated by Delta-Notch signaling. This system appears to be the ancestral mode of hesC regulation in eleutherozoans.

The sea urchin Pmar1-HesC double-negative gate for endomesoderm specification has been considered to have been newly acquired during echinoid evolution. However, it is still unknown how the system was established because pmar1-related genes and their upstream regulators have not been examined in noneuechinoid echinoderms. Through experiments using various echinoderms, we provide a hypothetical evolutionary scenario for the pmar1 gene (Fig. 8A) and models of the upstream GRN for endomesoderm in starfish and cidaroids comparable with that of the euechinoid S. purpuratus (Fig. 8B). The Pmar1-HesC regulatory system provides a remarkable opportunity to understand the diversification of the early developmental GRN.

Our screening of upstream regulators of endomesoderm development in noneuechinoid echinoderms revealed that cidaroids have pmar1 genes and that other echinoderms possess pmar1/phb1-related phb genes (Fig. 1). Our phylogenetic analyses supported the hypothesis that duplication of an ancient phb1 gene (referred to as phb genes in this study) led to the emergence of the pmar1 gene, which was reported by Dylus et al. (2016) (Fig. 8A). We note that the eh1-like motif is shared by Phb1 and Pmar1 in echinoids. Thus, we prefer the evolutionary history that pmar1 and phb1 were derived from an ancestral gene in which the eh1-like motif evolved in the common ancestor of echinoids, although the relationships between pmar1, phb1 and other phb genes have not been revealed by phylogenetic analysis. After the emergence of the pmar1 gene, tandem duplications of the pmar1 gene occurred in the genome, and the substitution rate of Pmar1 was accelerated, a phenomenon known as asymmetric evolution (Holland et al., 2017). This asymmetric evolution made our phylogenetic analyses less resolvable. Thus, we cannot exclude the alternative phylogenetic history in which the common ancestors of echinoderms possessed both phb and pmar1, and gene loss occurred in multiple lineages, although this scenario is less parsimonious.

Although pmar1 shows an accelerated substitution rate, its basic biochemical nature and developmental function have not changed. Pmar1 remains a transcriptional repressor, and it may have obtained an additional eh1-like motif in the echinoid lineage. In addition to the asymmetric evolution of pmar1, another important evolutionary change is the positioning of hesC as a downstream target gene, as mentioned below. However, our analysis revealed that there was a time lag between the extensive duplication of phb/pmar1 and the acquisition of hesC as a target gene. The asymmetric gene duplication occurred before the cidaroid-euechinoid divergence, but HesC regulation was placed under the control of Pmar1 in the common ancestor of euechinoids.

Our models of the endomesoderm GRN in cidaroids and starfish provide valuable information for understanding the evolutionary scenario of the upstream GRN (Fig. 8B). Our results support the idea that in the common ancestor of starfish and sea urchins (eleutherozoan), the phb gene exhibited a regulatory function in endomesoderm specification by regulating ets1 (ets1/2) and delta (Fig. 8A). Because Phb likely functions as a transcriptional repressor, starfish ets1 and delta may be regulated by an unknown repressor that is repressed by Phb proteins (X in Fig. 8B). Because HesC also acts as a repressor of alx1 expression in cidaroids, this interaction of HesC with alx1 dates back to the common ancestor of echinoids, even though hesC is activated indirectly, and not repressed, by Pmar1. Thus, the presumptive GRN in cidaroids seems to represent an intermediate state in the evolution of the Pmar1-HesC double-negative gate; i.e. HesC is just beginning to regulate alx1 to specify mesodermal skeletogenic cells in this GRN. After euechinoids diverged from the cidaroid lineages, hesC regulation was placed under the control of Pmar1. Our cross-species analysis suggests that the positioning of hesC downstream of Pmar1 occurred through the modification of the cis-regulatory sequence of hesC and not through the alteration of the coding sequence of Pmar1.

These results bring us to the idea that HesC was one of the genes whose developmental role was most drastically changed during echinoderm evolution. McCauley et al. (2010) demonstrated that, in embryos of the starfish P. miniata, HesC knockdown caused no obvious effects on the expression of ets1/2 or tbr or endomesoderm formation. Similarly, we observed that the embryos of the starfish P. pectinifera injected with two distinct hesC-specific MOs showed no defects in the vegetal tissues until the late gastrula stage (Fig. S7), suggesting that HesC is not essential for early endomesoderm specification. On the other hand, in cidaroid embryos, HesC appears to repress alx1 expression to some extent and skeletogenic cell fate (Erkenbrack and Davidson, 2015; this study). It is only in the euechinoid lineage that multiple other endomesodermal regulatory genes, such as ets1 and delta, are regulated by HesC (Fig. 8B). The next issue to address to understand the changes in the GRN will be how the multiple endomesodermal regulatory genes were placed under the control of HesC, which may have occurred through the modification of cis-regulatory motifs of target genes or those in the coding sequences of HesC. It should be noted that this event is expected to have occurred in a coordinate manner with the addition of hesC regulation under Pmar1 control.

Regarding the evolutionary change in the Pmar1-target gene repertoire, we consider the possibility that the target may not have been simply transferred from an unknown factor (X) to HesC. The present study shows that Phb/Pmar1 repress an unknown repressor other than HesC in starfish and cidaroid embryos (Fig. 8B). On the other hand, our previous study using two distantly related euechinoid species suggested the existence of an additional repressor downstream of euechinoid Pmar1. In euechinoid embryos injected with pmar1 (micro1) mRNA, the fate conversion of almost all cells into putative skeletogenic cells is commonly observed (Oliveri et al., 2002; Nishimura et al., 2004; Yamazaki et al., 2009), which suggests that Pmar1 represses the global repressor(s). In contrast, hesC MO-embryos show only a moderate expansion of skeletogenic cell region, although the transient global activation of some micromere/skeletogenic regulatory genes (delta and ets1) is observed at the earlier stage (Yamazaki and Minokawa, 2016). Furthermore, the Pmar1-overexpressing embryos show almost no expression of a nonskeletogenic regulatory gene gcm at the blastula stage, whereas hesC MO-embryos of two euechinoid species show expanded expression of gcm (Yamazaki and Minokawa, 2016), i.e. the regulatory states in these embryos are clearly distinguishable. This difference implies that an additional repressor of skeletogenic regulatory genes is present to repress skeletogenic cell fate in the animal region of euechinoids, which is supported by structure-function correlation analysis of the Pmar1 (Micro1) protein (Yamazaki et al., 2009). Accordingly, we predict that there is an unknown repressor (X) shared by eleutherozoans downstream of Phb/Pmar1; i.e. an additional target of Pmar1, hesC, had been added in the euechinoid lineage. Our data suggest that the gene encoding the unknown repressor in the cidaroid is not a member of the hairy gene family because we found no hairy genes showing expression patterns similar to euechinoid hesC in the cidaroid P. baculosa (i.e. nonvegetal ectodermal expression during the blastula stage). Two hairy genes (hesA and hesD) showed zygotic expression during the early stages of P. baculosa (Fig. S1); strong expression of hesA was detected in whole embryos, whereas hesD showed only faint expression at the blastula stage. To understand how the upstream GRN has been rewired, it is crucial to identify the unknown repressor X in these animals. Further study of the sea urchin phb1 genes may also be informative to determine how the ancestral function of phb was modified after the emergence of pmar1.

Our results illuminate the evolutionary history of the echinoderm GRN, in which the upstream GRN recruited a new component (i.e. hesC) without changing the developmental outcome. The stepwise rewiring of transcription networks via the addition of target genes is one of the general evolutionary pathways observed in the transcription network of yeast. Li and Johnson (2010) proposed that transition through intermediate states would not decrease fitness, which is the key to understanding the evolutionary process of GRN rewiring. In the case of GRN modification in euechinoids, the following two changes must have occurred: (1) recruitment of new repressor targets under the control of HesC; and (2) recruitment of hesC as a repressor target of Pmar1, irrespective of the order. The future questions that we need to address using the experimental system of echinoderms are as follows: what is the intermediate state that enabled GRN rewiring without causing catastrophe, and what sort of molecular evolution occurred during this stepwise process?

The collection of adult P. baculosa and H. pulcherrimus and the handling of gametes and embryos were performed according to a previously described method (Yamazaki et al., 2010; Hibino et al., 2019). The collection and handing of gametes of P. pectinifera and A. kochii were performed according to Koga et al. (2010). The adults of A. japonicus were collected around the Misaki Marine Biological Station, University of Tokyo. The gametes were obtained according to the method of Kikuchi et al. (2015).

To survey the candidate upstream genes in the endomesoderm GRN of the cidaroid P. baculosa, we performed screening based on temporal expression patterns. We first performed RNA-seq using samples from embryos at five developmental stages (2 h, 4 h, 6 h, 10 h, and 14 h). The criteria for the selection of candidate genes were as follows: (1) zygotic activation (more than a fivefold increase in the FPKM value compared with the value of maternal expression) earlier than or simultaneously with expression of alx1, (2) encoding a DNA-binding domain according to the HMMER search (Johnson et al., 2010), and (3) an FPKM value greater than 5 at the onset of Pb-alx1 activation (10 h). To identify these obtained sequences, we checked the top BLAST hit sequence in the sea urchin S. purpuratus gene database (named Sp genes) of EchinoBase (www.echinobase.org/Echinobase/). The amplified sequence of P. baculosa pmar1 has been deposited in DDBJ (accession number LC483152), and the sequences of the other genes are shown in the supplementary Materials and Methods. We further identified pmar1, phb1 and phb genes from the other echinoids and echinoderms. See supplementary Materials and Methods for more details.

The genes used for the outgroup were selected according to the previous analysis (Dylus et al., 2016). See supplementary Materials and Methods for more details.

Whole-mount in situ hybridization was performed as described previously (Koga et al., 2010; Yamazaki et al., 2010; Morino et al., 2012). See Supplementary Materials and Methods and Table S2 for more details.

The overexpression and translational perturbation analyses were performed using synthesized mRNAs and MOs, respectively. See Supplementary Materials and Methods for more details.

The qPCR analysis was performed as described previously (Yamazaki et al., 2012) using CFX Connect Real-Time PCR Detection System (Bio-Rad). The EF1alpha gene was used for the internal reference standard according to Koga et al. (2016). The primer sequences used are follows: Pb-pmar1-qF, 5′-CGATATCGACGTGCGAGAAA-3′; Pb-pmar1-qR, 5′-TGAAACCAGACCTGTATTCTC-3′; Pb-EF1a-qF, 5′-GCGTGAGCGAGGTATCACAAT-3′; and Pb-EF1a-qR, 5′-ACAATCAGCACCGCACAATC-3′.

The differentiation of skeletogenic cells, mesenchyme cells and endoderm was evaluated using fixed embryos of P. baculosa, H. pulcherrimus and P. pectinifera. See Supplementary Materials and Methods for more details.

We thank Masashi Noguchi, Setsuo Kiyomoto, Norihiko Deguchi and Toshimitsu Fukuhata for the collection of adult P. baculosa. We also thank the staff of Noto Marine Laboratory of Kanazawa University: Nobuo Suzuki, Toshio Sekiguchi, Hiroyasu Kamei and Shozo Ogiso for the culturing system of adult P. baculosa and laboratory equipment; Misa Yamaguchi and Sae Sugino for the support in a part of embryological experiments using P. baculosa; staff of Tateyama Marine Laboratory of Ochanomizu University and Research Center for Marine Biology of Tohoku University for collecting P. pectinifera and H. pulcherrimus; Gen Hamanaka for providing valuable technical advices for microinjection in P. pectinifera; Naoki Irie for his help in identifying genes from A. japonicus and O. japonicus; and Masato Kiyomoto and Hiroyuki Kaneko for providing antibodies. We thank two anonymous reviewers for their helpful comments and advice.