The genomic regulatory control of skeletal morphogenesis in the sea urchin

The dynamic anatomical changes that characterize embryogenesis are encoded in the genome. The genomic regulatory control of development can be understood in terms of transcriptional gene regulatory networks (GRNs), which can be defined as dynamic networks of interacting genes that encode transcription factors (i.e. regulatory genes). Such networks are being used to analyze cell specification in diverse organisms (Stathopoulos and Levine, 2005; Ettensohn, 2009; Nikitina et al., 2009; Peter and Davidson, 2009; Davidson, 2010). There has been less progress, however, in using GRNs to explain the complex cell and tissue behaviors that drive changes in embryonic form. Insights in this area will emerge first from model systems in which there is a detailed understanding of both (1) cell specification at the GRN level and (2) morphogenetic mechanisms at the molecular, cellular and tissue levels. Establishing linkages between early cell specification networks and specific morphogenetic processes is crucially important, not just for understanding embryogenesis per se, but also in an evolutionary context, i.e. for understanding the ways in which evolutionary modifications to genetic networks have led to changes in morphological features.

The formation of the embryonic endoskeleton of sea urchins is a powerful experimental model for developing an integrated view of the genomic regulatory control of morphogenesis. The skeleton is the primary determinant of the distinctive, angular shape of the larva and influences its orientation, swimming and feeding (Pennington and Strathmann, 1990; Hart and Strathmann, 1994). The skeleton is a biomineral that consists of calcite and small amounts of occluded proteins. It is secreted by primary mesenchyme cells (PMCs), a specialized population of skeletogenic cells that have a well-defined embryonic lineage. During gastrulation, PMCs undergo a striking sequence of morphogenetic behaviors that includes epithelial-mesenchymal transition, directional cell migration and cell-cell fusion (Wilt and Ettensohn, 2007; Ettensohn, 2009). These cellular behaviors have been analyzed in remarkable detail, largely owing to the optical transparency of the sea urchin embryo and its suitability for in vivo imaging (Gustafson and Wolpert, 1967; Ettensohn and Malinda, 1993; Malinda et al., 1995; Miller et al., 1995; Peterson and McClay, 2003; Hodor and Ettensohn, 2008).

Recently, a GRN that underlies PMC specification has been described (Oliveri et al., 2008; Ettensohn, 2009). This network is initially deployed through the activity of polarized, maternal inputs that activate a small set of early zygotic regulatory genes selectively in the large micromere-PMC lineage. The transcription factors encoded by these genes engage additional layers of regulatory genes; various feedback and feed-forward interactions subsequently stabilize the transcriptional network and drive it forward (Oliveri et al., 2008). Although much information is available concerning the interactions among regulatory genes that are deployed in the network, little is known concerning the downstream circuitry of the network, i.e. the current network model includes few connections to downstream genes that play a direct role in the morphogenetic program of these cells. Some of these downstream morphogenetic effector (‘morphoregulatory’) genes are known, although many important genes in this class remain to be identified.

In previous work, we generated an arrayed cDNA library from PMCs at the mid-gastrula stage, when these cells are engaged in their most prominent morphogenetic behaviors. Initial analysis of this library allowed us to identify several components of the PMC GRN, including delta (Sweet et al., 2002), the transcription factors alx1 and erg (Zhu et al., 2001; Ettensohn et al., 2003), and several biomineralization-related genes (Illies et al., 2002; Cheers and Ettensohn, 2005; Livingston et al., 2006). In the present study, we greatly expand the current PMC GRN model by identifying many additional morphoregulatory genes and by uncovering regulatory inputs into these genes. Our work extends the value of the PMC GRN as a model developmental GRN and establishes a framework for understanding the genomic circuitry that encodes a major anatomical feature.

Strongylocentrotus purpuratus embryos were obtained and cultured at 15°C as described previously (Zhu et al., 2001).

The construction and arraying of the PMC cDNA library has been described (Zhu et al., 2001; Livingston et al., 2006). Briefly, the library was generated from polyA(+) RNA that was isolated from micromeres (presumptive PMCs) that were cultured until sibling control embryos reached the mid-gastrula stage, ∼36 hours post-fertilization. The cDNA library was not normalized or subtracted in any way. Reverse transcription was primed using oligo(dT) and the resulting cDNAs were directionally cloned into pSPORT with an average insert size of 1.5-2 kb. Clones were robotically arrayed in 384-well plates and subjected to Sanger-based sequencing. A total of 51,097 cDNA sequence reads (ESTs) were obtained from randomly selected clones (Genbank entries BG780044-BG789446, DN560823-DN586179 and DN781759-DN810571). In most cases, single 5′-reads were obtained from each clone; in some cases, both 5′ and 3′ sequence reads were obtained.

WMISH was performed as described previously (Zhu et al., 2001).

MOs were injected into fertilized eggs as described by Cheers and Ettensohn (Cheers and Ettensohn, 2004), with the modification that eggs were fertilized in the presence of 0.1% (wt/vol) para-aminobenzoic acid to prevent hardening of the fertilization envelope. MO sequences were: SpAlx1: 5′-TATTGAGTTAAGTCTCGGCACGACA-3′; SpEts1: 5′-GAACAGTGCATAGACGCCATGATTG-3′; SpTbr: 5′-TGTAATTCTTCTCCCATCATGTCTC-3′.

All three MOs were translation-blocking MOs. The SpAlx1 MO has been shown to specifically and effectively block SpAlx1 expression (Ettensohn et al., 2003). The SpEts1 and SpTbr MOs were also described previously (Oliveri et al., 2008). The SpEts1 MO phenocopies the overexpression of a dominant-negative form of Ets1, which blocks PMC specification and ingression (Sharma and Ettensohn, 2010). The SpTbr MO produces a selective effect on skeletogenesis without affecting PMC ingression or migration (Oliveri et al., 2008). MOs were injected at working concentrations of 4 mM (Alx1) or 2 mM (Ets1, Tbr).

QPCR was carried out as described previously (Stamateris et al., 2010). PCR primers that were used in this study are shown in supplementary material Table S4. Technical duplicates were included for each sample and the average of the two C_t values was used to calculate ΔC_t. For temporal QPCR profiles, the absolute numbers of mRNA molecules per embryo at various developmental stages were calculated using Sp-z12 as a standard. The numbers of Sp-z12 transcripts at various stages have been determined by RNA titration (Wang et al., 1995). To determine regulatory inputs into genes, MOs were used to block the expression of ets1, alx1 or tbr and the effects on the expression of downstream skeletogenic genes were analyzed at 28-30 hours post-fertilization by QPCR.

We assessed the expression of each gene in morphant embryos compared with control siblings in three batches of embryos that were derived from independent matings (three biological replicates). We calculated ΔC_t values with respect to Sp-z12 and used Student’s t-test to compare the levels of expression of each gene in control and morphant embryos (n=3). A P value of ≤0.05 (i.e. a 95% confidence value) was taken to be significant and such genes were considered to have an input from the gene that was targeted by the MO. Of the genes that satisfied this criterion, most showed an average level of expression that was <25% of the control level. We also considered an effect on gene expression to be meaningful if, in all three biological replicates, the fold-difference in expression in control embryos relative to morphants was at least threefold and always in the same direction, a criterion applied by Oliveri and co-workers (Oliveri et al., 2008). Applying this criterion caused us to score as regulatory targets a small number of genes (three cases in total) that consistently showed pronounced effects, but which also showed a relatively high level of variability across biological replicates. Lastly, in three cases, genes that showed expression of <33% of the control level in two of the three biological replicates were analyzed independently by WMISH. We considered a gene to have an input if its level of expression was clearly and consistently reduced by MO knockdown, as assessed by WMISH (all three cases).

Sanger-based sequence reads (51,097; average read length >900 nts) from randomly selected PMC cDNAs were mapped to the complete set of GLEAN3 gene predictions (Sea Urchin Genome Sequencing Consortium, 2006). We expanded the GLEAN3 models to include 3′-UTR sequences that were identified through a genome-wide tiling analysis of gene expression (Samanta et al., 2006). Of the 51,097 ESTs, about half (24,238) aligned to GLEAN3 gene models; these ESTs were mapped to 7415 different GLEAN3 genes (supplementary material Table S1). The number of EST matches per gene varied from one to 669 (supplementary material Fig. S1). Approximately one-third of the 7415 GLEAN3 models had a single EST match, approximately one-third had two to four matches, and approximately one-third had five or more matches. The number of EST hits per gene corresponded very roughly to the abundances of the corresponding transcripts at the mid-gastrula stage (36 hours post-fertilization) as measured by QPCR (see below). Manual curation of a randomly selected sample of 100 ESTs that did not match GLEAN3 models showed that ∼50% could be attributed to overly conservative predictions of 3′-UTR sequences, ∼25% represented rRNA sequences or previously unidentified exons, and the remaining 25% represented simple repetitive sequences or other sequences that could not be mapped to the current sea urchin nuclear or mitochondrial genome assemblies.

As one means of assessing the completeness of the collection of PMC-expressed genes, we investigated whether it included all genes that had been shown in previous studies to be expressed selectively by PMCs. Of the 52 such genes that had been identified prior to this study (supplementary material Table S2), only three were not represented in our collection of PMC-expressed genes. In two cases (foxN2/3 and snail), this was attributable to the temporal expression patterns of the genes; in S. purpuratus, foxN2/3 expression is extinguished in the micromere lineage by the gastrula stage (Tu et al., 2006), whereas snail is not expressed until late in gastrulation, and then only at extremely low levels (Oliveri et al., 2008). In the case of fgfr2, closer inspection showed that the collection of PMC ESTs included sequences that corresponded to the 3′UTR of the fgfr2 mRNA, which was incorrectly predicted in the gene model.

Although the catalog of genes that are expressed by PMCs during gastrulation (supplementary material Table S1) is likely to be nearly complete, this list also contains some genes that are not expressed by PMCs but are represented in the catalog because the micromere preparation that was used to generate the original cDNA library was not completely pure and contained small numbers of other cell types (∼5%) (Zhu et al., 2001). For example, we identified several ESTs that corresponded to spec1, a highly abundant, aboral ectoderm-specific transcript. It is significant, however, that many transcripts that are expressed at moderate to high levels at the gastrula stage in various non-skeletogenic lineages were completely absent from the PMC mRNA catalog, or were extremely rare. For example, ESTs that corresponded to the endoderm marker endo16 (GLEAN3_11038/9), the pigment cell marker pks (GLEAN3_02895) and the apical ectoderm marker nk2-1 (GLEAN3_00757) were absent from the PMC EST collection, as expected. Based on these considerations, it seems prudent to conclude that the presence of any specific gene in the PMC gene catalog, particularly if the gene has very few EST matches, is only suggestive of expression in PMCs and confirmation by independent methods is required.

The spatial patterns of expression of 180 genes in the PMC gene catalog were analyzed by whole-mount in situ hybridization (WMISH). These genes were chosen by focusing on relatively abundant transcripts (based on the number of EST matches/gene) and by excluding obvious housekeeping genes or genes that had been previously analyzed by WMISH. Embryonic stages from the fertilized egg to the early pluteus larva were examined. Many of the genes we analyzed by WMISH showed rather general patterns of expression during development; however, we also identified 33 mRNAs that were restricted to the large micromere-PMC lineage or were highly enriched in these cells during gastrulation, when the PMCs are actively engaged in their salient morphogenetic activities (Figs 1, 2). Most of these mRNAs encode proteins that are expected to function in morphogenetic processes and not in cell specification; e.g. cytoskeletal proteins and their regulators, putative cell adhesion and extracellular matrix (ECM) proteins, new biomineralization proteins and novel proteins. We did, however, also identify three new PMC-enriched transcription factors (Fos, Smad1/5/8 and Alx4). The functions of some of the gene products that emerged from the WMISH analysis have been well characterized in other organisms (e.g. Wasp, Arp proteins, Cdc42, etc.) but others are novel proteins or have cellular functions that are less well understood. Information concerning the newly identified PMC-specific and PMC-enriched proteins that emerged from the WMISH analysis is summarized in supplementary material Table S3.

Quantitative polymerase chain reaction (QPCR) was used to determine high-resolution temporal expression profiles of many of the newly identified mRNAs (Fig. 3). Several PMC-specific mRNAs that had been identified in previous studies, but for which detailed expression profiles were not available, were also included in the analysis. Only four of the genes that we examined (casc1, cyp2, smad1/5/8 and jun) were expressed maternally at levels of >1000 transcripts per egg. Zygotic transcription of all of the genes was activated prior to the onset of PMC ingression (24 hours post-fertilization), as shown by an increase in steady-state transcript levels. The timing of zygotic activation, however, appeared to vary considerably. The first genes to be expressed (casc1, cdi, fos, p19, egf-rich/p41, pks2 and stomatin) were activated by 12 hours post-fertilization (early blastula stage) whereas transcripts encoded by late genes (can1 and p16rel1) did not begin to accumulate until shortly before PMC ingression. Maximal levels of expression for most genes were observed during gastrulation (24-48 hours post-fertilization), which includes the developmental stage from which the PMC cDNA library was prepared (mid-gastrula stage; ∼36 hours post-fertilization). More than half of the genes (13/25) showed maximal expression at 42 hours post-fertilization; two exceptions were sm29 and Clectin, which showed maximal transcript levels at 72 hours post-fertilization (the latest stage examined). Calculated on a per-cell basis (32 PMCs per embryo), peak expression levels varied from 30 transcripts per cell (can1) to 600 transcripts per cell (p19).

We used morpholino (MO) knockdowns and QPCR to identify regulatory inputs into many of the newly identified genes. We focused on Ets1, Alx1 and Tbr, three well-characterized transcription factors that provide some of the earliest inputs into the network. For QPCR studies, we restricted our analysis to downstream genes that were expressed exclusively by PMCs, so that potentially subtle changes in transcript levels in PMCs would not be obscured by the presence of the same mRNA in other cell types. We examined the effects of gene knockdowns at 28-30 hours of development (i.e. at the early gastrula stage), when all morphoregulatory genes that we tested were ordinarily expressed at high levels (see Fig. 3).

The effects of Ets1, Alx1 and Tbr knockdowns on the levels of expression of PMC-specific mRNAs are shown in Fig. 4. The regulatory inputs of Ets1, Alx1, and Tbr into their morphoregulatory gene targets are invariably positive inputs. We did, however, detect auto-repression by two of the regulatory genes, alx1 and tbr, as was reported previously (Ettensohn et al., 2003; Oliveri et al., 2008). We also used WMISH to confirm our QPCR results; in every case that we examined, WMISH data supported the assignments of regulatory linkages that were based on QPCR (Fig. 5). We also used WMISH to assess certain inputs into morphoregulatory genes that showed substantial changes in expression in more than one biological replicate but failed to meet the statistical criteria that we established (see Materials and methods). In all three such cases, WMISH analysis supported the assignment of a positive regulatory input (see Fig. 5; effects of Ets1 knockdown on the expression of p58B, sm49 and Clectin).

Figure 6 summarizes the new network connections that were identified in this study. Our findings reveal that: (1) tbr has very few inputs into morphoregulatory genes; (2) ets1 and alx1 provide inputs into the majority (approximately two-thirds) of the morphoregulatory genes, including all of the known biomineralization-related genes that we tested; and (3) the targets of ets1 and alx1 are nearly identical. We detected no essential inputs from ets1, alx1 or tbr into six of the morphoregulatory genes that we tested (casc1, cdi, cyp2, egf-rich/p41, fgfr2 and stomatin). By combining this new network circuitry with previously published work, it is now possible to construct an expanded model of the PMC GRN that bridges upstream regulatory genes to many of the effector genes in the network (Fig. 7).

The formation of the complex and precisely patterned embryonic skeleton is the culmination of a sequence of morphogenetic activities on the part of the biomineral-forming PMCs (Wilt and Ettensohn, 2007). At present, there is a relatively detailed picture of the gene network that underlies early PMC specification; indeed, the micromere-PMC GRN is arguably the most complete developmental GRN in any experimental model (Oliveri et al., 2008; reviewed by Ettensohn, 2009). It is evident that this early specification network impinges on genes that control morphogenesis, as the behaviors of PMCs during gastrulation are dependent upon zygotic transcriptional inputs (Kurokawa et al., 1999; Ettensohn et al., 2003; Wu and McClay, 2007). The connections between the transcriptional GRN and the various cellular activities that underlie skeletogenesis, however, are not well understood. It is precisely this set of regulatory linkages that is most directly related to the assembly and patterning of the larval skeleton.

Our EST analysis and WMISH screen led to the identification of many PMC-specific and PMC-enriched mRNAs that encode cytoskeletal proteins and their regulators, putative adhesion proteins, extracellular matrix proteins and biomineralization proteins (supplementary material Table S3). The developmental functions of many of these proteins have yet to be analyzed in detail, but in other cases their functions are relatively well understood. For example, spicule matrix proteins are a family of secreted proteins that are occluded within the biomineral and influence its growth and physical properties (Wilt and Ettensohn, 2007). A recent proteomic analysis of proteins in the embryonic spicules of S. purpuratus (Mann et al., 2010) has converged on essentially the same suite of C-type lectin domain (CTLD)-containing proteins that has emerged from our analysis of the PMC transcriptome (Illies et al., 2002; Livingston et al., 2006) (this study). It is therefore likely that all of the canonical (CTLD-containing) spicule matrix proteins that are expressed in the embryo have now been identified (supplementary material Fig. S2). Another protein that has an essential function in biomineralization is the novel, PMC-specific transmembrane protein P16 (Cheers and Ettensohn, 2005). P16rel1 and P16rel2, which are similar in structure to P16 and are encoded by closely linked genes, probably play a similar role in biomineralization. A recent functional analysis of P58A and P58B, two related Type I transmembrane proteins that were identified in the present study, has shown that both proteins are required for biomineral formation but not for PMC specification, migration or fusion (Adomako-Ankomah and Ettensohn, 2011). The precise biochemical functions of the P16 and P58 proteins are unknown, although it was recently shown that P16 is phosphorylated and binds to hydroxyapatite (Alvares et al., 2009).

In addition to biomineralization-related proteins, we identified several other secreted or transmembrane proteins that are expressed by PMCs (e.g. P11, EGF-rich/P41, P133, Stomatin and Tsp). P11 and Tsp are both small, secreted proteins; P11 lacks identifiable motifs but Tsp contains a single thrombospondin type I repeat (TSR). In vertebrates, TSR domain-containing proteins, which are typically larger and have more complex domain organizations, are constituents of the ECM and function in cell adhesion, signaling and motility (Adams and Tucker, 2000). Stomatin is a membrane protein; vertebrate stomatins regulate the function of other membrane-associated proteins, including ion channels, by mechanisms that are poorly understood (Lapatsina et al., 2011). P133 is a single-pass transmembrane protein with multiple, interspersed EGF and Laminin G domains, an overall organization that is strikingly similar to that of members of the neurexin family of cell adhesion receptors (Lisé and El-Husseini, 2006). These features suggest that P133 might mediate adhesion between PMCs and substrate molecules in the basal lamina (Hodor et al., 2000), or PMC-PMC adhesion as part of the process of cell fusion.

Our analysis has also identified several PMC-enriched mRNAs that encode cytoskeletal proteins. Almost all of these mRNAs are co-expressed by PMCs and non-skeletogenic mesoderm (NSM) cells, two kinds of mesenchymal cells that exhibit similar kinds of invasive, motile behaviors (see below). Many of these proteins have known biochemical functions and act cooperatively in cells. For example, we identified four components or regulators of the Arp2/3 complex (Arp3 and p21arc/ARPC3, Cdc42 and Wasp) in the set of PMC-enriched mRNAs. Given the well-known role of the Arp2/3 complex in regulating actin dynamics and cell protrusive activity, including the formation of filopodia, it seems likely that this collection of proteins plays a role in regulating the filopodial motility of sea urchin mesenchyme cells (Pollard, 2007; Faix et al., 2009; Mellor, 2010). Cdc42 has been specifically implicated in the formation of filopodia in several cell types (Ahmed et al., 2010). We also identified Lasp1 as a PMC-specific transcript. In mammalian cells, this actin-binding phosphoprotein is enriched in cell protrusions, and Lasp1 has recently been shown to mediate the directional migration of leucocytes by binding to the chemotactic receptor CXCR2 (Raman et al., 2010). Lasp1 is therefore a strong candidate for a motility-related protein in PMCs.

In other experimental systems, emphasis has been placed on the post-translational regulation of cytoskeletal machinery in regulating early morphogenetic processes. For example, studies in Drosophila indicate that spatially controlled, zygotic transcriptional activation of a small number of key proteins that regulate Rho- or G protein-mediated signaling results in the local activation of ubiquitous cytoskeletal machinery, thereby leading to region-specific cell shape changes and cell movements during gastrulation (Dawes-Hoang et al., 2005; Kölsch et al., 2007). Although local, post-translational regulation of the cytoskeleton undoubtedly also contributes to region-specific cell behaviors during sea urchin gastrulation, we also find evidence of a widespread transcriptional upregulation of many cytoskeletal genes in PMCs. In other embryonic and metastatic cell types, invasive cell behavior has been shown to be under transcriptional control. It is of interest that jun and fos have been demonstrated to regulate invasiveness in several cellular contexts (Sherwood et al., 2005; Ozanne et al., 2007) and that both genes are selectively upregulated in PMCs.

Our findings reveal several important features of the regulatory inputs that are provided by alx1, ets1 and tbr into the set of PMC-expressed morphoregulatory genes. (Note, however, that our gene knockdown analysis does not allow us to determine whether these inputs are direct or indirect.) alx1 and ets1 provide positive regulatory inputs into most (approximately two-thirds) of the morphoregulatory genes that we tested. By contrast, tbr provides inputs into a much smaller subset of genes (only lasp1, p16rel1, p16rel2 and sm49), even though tbr is activated early in development; i.e. several hours before PMC ingression and prior to the activation of most of the downstream genes that we analyzed here. The paucity of morphoregulatory genes that receive inputs from tbr is consistent with the finding that tbr also provides few inputs into regulatory genes in the GRN (Oliveri et al., 2008) (SpBase Endomesoderm Gene Regulatory Network). It is significant that tbr has undergone significant, recent evolutionary modifications with respect to its network connectivity (Hinman et al., 2007) and the lack of tbr expression in the skeletogenic centers of adult echinoids and sea stars suggests that this gene was recently co-opted into the micromere-PMC GRN of echinoids (Gao and Davidson, 2008). The limited connectivity of tbr in the PMC GRN probably reflects the recent recruitment of this gene into the network.

One of the most striking findings from our gene knockdown studies is that the targets of alx1 and ets1 are almost identical. The single exception that we identified is p58B, which receives a positive input from ets1 but not from alx1. Even in this case, the data are equivocal; our QPCR measurements suggest that alx1 knockdown might have a modest effect on p58B expression that did not meet the significance criteria that we established. This case aside, we note that all the biomineralization genes that were examined (spicule matrix genes, can1 and members of the p16 and p58 gene families) are co-regulated by ets1 and alx1. Parallel inputs from ets1 and alx1 also regulate vegfr-Ig-10, which encodes a receptor tyrosine kinase with a crucial role in PMC guidance and differentiation (Duloquin et al., 2007).

A variety of mechanisms might explain the parallel connectivity of ets1 and alx1. alx1 knockdown has little, if any, effect on ets1 expression (Ettensohn et al., 2003; Oliveri et al., 2008). By contrast, perturbation of ets1 expression or function, although it does not affect the early phase of alx1 expression, suppresses the later phase (Oliveri et al., 2008; Sharma and Ettensohn, 2010). It seems likely, therefore, that the parallel connectivity of alx1 and ets1 can be attributed partly to the influence of ets1 on alx1 expression. Other evidence indicates, however, that ets1 regulates terminal genes in the network by mechanisms that are independent of alx1 expression. For example, ets1 provides positive inputs into other regulatory genes, including erg and hex, that are required for the expression of differentiation genes but that do not provide regulatory inputs into alx1 (Oliveri et al., 2008) (SpBase Endomesoderm Gene Regulatory Network). Given the evidence that ets1 regulates terminal genes in the network by alx1-independent mechanisms, it is striking that we found no example of a gene that receives an essential input from ets1 but not from alx1, save the one questionable exception noted above. Oliveri and co-workers (Oliveri et al., 2008) proposed that ets1 might cooperate with various regulatory genes in feed-forward loops (A>B, B>C, A>C) that impinge on terminal genes in the network, and experimental evidence for a feed-forward loop that involves ets1 and alx1 has come from analysis of the cis-regulatory control of the cyp1 gene, which receives direct inputs from Dri (a target of alx1) and Ets1 (Amore and Davidson, 2006). If feed-forward regulatory interactions are a common feature of the wiring of the skeletogenic GRN, then one interpretation of our findings is that almost every morphoregulatory gene that receives an input from ets1 is wired with a feed-forward loop that includes alx1.

Although most morphoregulatory genes in the GRN are regulated by ets1 and alx1, six of the genes that we identified (casc1, cdi, cyp2, egf-rich/p41, fgfr2 and stomatin) do not appear to require inputs from any of the three transcription factors that we analyzed. One of the important morphoregulatory genes that is presently unconnected to the GRN is fgfr2, which functions in PMC migration and differentiation (Röttinger et al., 2008). Because the expression of every regulatory gene in the current model of the GRN, with the exception of tel, is influenced, directly or indirectly, by ets1, alx1 or tbr, our findings suggest that these morphoregulatory genes might receive inputs from as yet undiscovered subcircuits within the PMC GRN. These genes might receive inputs from tel, from one or more of the new PMC-enriched transcription factors that were identified in the present study (fos, alx4 and smad1/5/8) or from other regulatory genes in the network that have not yet been subjected to gene knockdown analysis, such as jun.

Analysis of gene expression in cultures of micromeres or dissociated embryonic cells indicates that the initial deployment of the skeletogenic GRN in the micromere-PMC lineage is largely independent of signals from other cell types (Harkey and Whiteley, 1983; Stephens et al., 1989; Page and Benson, 1992). Much evidence, however, demonstrates that, later in development, skeletal growth and patterning is regulated by local, ectoderm-derived cues (Ettensohn, 1990; Kiyomoto and Tsukahara, 1991; Armstrong et al., 1993; Guss and Ettensohn, 1997; Duloquin et al., 2007; Röttinger et al., 2008). Many genes, including many of the downstream effector genes that were identified in the present study, show non-uniform patterns of expression within the PMC syncytium, with high levels of transcripts accumulating at sites of skeletal rod growth (Harkey et al., 1992; Guss et al., 1997). There is evidence from the expression of transgenes within the PMC syncytium that mRNAs remain localized predominantly in the PMC cell body in which they were synthesized (Harkey et al., 1995; Makabe et al., 1995; Wilt et al., 2008), a finding which suggests that local variations in mRNA abundance within the syncytium arise through the local regulation of gene transcription, rather than by an indirect mechanism (e.g. the diffusion and selective trapping of mRNAs). One model, which we favor, is that local, ectoderm-derived signals directly modulate the skeletogenic GRN and control the local availability of gene products (e.g. biomineralization proteins) that control the rate of biomineral deposition. The ectoderm has other important, and probably more subtle, influences on skeletal morphogenesis (for example, on the branching of skeletal rods), which remain mysterious. It will be important in the future to elucidate the various mechanisms by which the ectoderm influences skeletal growth and patterning.

PMCs are a relatively recent evolutionary invention. The appearance of this cell population probably involved at least two steps: (1) the importation of a late larval/adult program of biomineralization into the late (gastrula stage) embryo, and (2) a second heterochronic shift, which shifted the skeletogenic program into the pre-gastrula stage embryo and which was associated with the formation of micromeres and an early-ingressing skeletogenic mesenchyme (Ettensohn, 2009). Our findings shed additional light on the mechanisms that might have accompanied the evolutionary origin of PMCs.

In euechinoid sea urchins (subclass Euechinoidea), a group that includes most modern sea urchins, several regulatory and signaling proteins that are expressed in the micromere-PMC lineage during early development (ets1, erg, foxN2/3, snail, wnt8 and delta) are also deployed in the adjacent NSM territory later in development, during gastrulation. In the present study, we have identified numerous non-regulatory genes (e.g. cytoskeletal proteins and secreted factors) that are expressed selectively in these two territories (Fig. 2). Thus, our current findings provide additional evidence of a striking similarity in the genomic regulatory states of PMCs and NSM cells. These observations support the view that the invention of the micromere-PMC lineage was associated with a heterochronic shift in the deployment of an ancestral GRN that was operative in the late-ingressing skeletogenic mesenchyme of ancestral echinoids. This mode of skeletogenesis is still seen in modern cidaroid urchins, which closely resemble the ancestral stock that gave rise to all extant sea urchins (Wray and McClay, 1989). It seems likely that one vestige of this heterochronic shift that has been retained in modern euechinoid sea urchins is the biphasic temporal expression (i.e. first in PMCs then in NSM cells) of many genes that are co-expressed in these two embryonic territories.

Gene duplication has played a major role in the evolution of the PMC GRN. Many of the genes that are expressed selectively in the micromere-PMC lineage are found in multiple, tandem copies. The stomatin, pmar1 and sm30 loci, for example, consist of five to ten tandem copies of these genes. All msp130 family members are found in small clusters of two to three genes, as are most of the spicule matrix genes (Livingston et al., 2006). The biomineralization gene p16 is adjacent to several related genes, at least one of which (p16rel2) also bears some similarity to spicule matrix genes. In addition, tsp, lasp1, p58 and alx1 are found as pairs of distinct but closely related genes. There is evidence that higher order clustering of these and other PMC-expressed genes is obscured by the incomplete nature of the current S. purpuratus genome assembly. For example, one spicule matrix gene that is expressed in the adult (GLEAN3_13825) (Mann et al., 2008) is closely linked to the msp130-msp130rel1-msp130rel3 gene cluster, indicating that at least some members of these two major gene families are clustered in the genome. These observations point to the active role of gene duplication in expanding the set of terminal morphoregulatory genes in the PMC GRN. Most of the duplicated genes contain introns (the stomatin genes are an exception) and are tandemly arranged, indicating that their duplication has occurred by DNA-mediated processes, such as unequal crossing over, rather than by retrotransposition.

Gene duplications are remarkably common events on an evolutionary time scale, but the usual fate of a duplicate gene pair is the silencing of one copy within a few million years (Hancock, 2005). Several of the gene duplications described above, however, occurred >100 million years ago and the duplicate copies have remained functional (Adomako-Ankomah and Ettensohn, 2011). An adaptive expansion of biomineralization-related genes might have occurred if natural selection favored an increase in their dosage or if a rapid ‘microfunctionalization’ of biomineralization proteins occurred (Hancock, 2005). Extensive duplication of genes that encode secreted, biomineralization proteins has occurred independently in vertebrates (Kawasaki et al., 2009). The duplication of other functional classes of genes might have facilitated the evolution of novel, cell type-specific protein functions or new gene expression patterns. In this regard, we note that the PMC-specific representatives of some protein families, such as the receptor tyrosine kinase family, have atypical structures that might be associated with unusual biochemical properties (Duloquin et al., 2007; Röttinger et al., 2008). Of particular interest is the duplication of transcription factor-encoding genes such as alx1 and pmar1, which might have allowed these genes to acquire new upstream regulators (through changes in cis-regulatory sequences) or new downstream targets (through changes in the biochemical properties of the proteins). Such gene duplication events might have been intimately involved in the heterochronic shifts in the deployment of the skeletogenic GRN that occurred during echinoderm evolution.

The authors thank Angela Bozak, Rachel Schuster and Zhongling Sun (Carnegie Mellon University) for assistance with WMISH analysis; and Melissa Landrum (National Center for Biotechnology Information) for mapping ESTs to the genome assembly.