A central challenge of developmental and evolutionary biology is to understand the transformation of genetic information into morphology. Elucidating the connections between genes and anatomy will require model morphogenetic processes that are amenable to detailed analysis of cell/tissue behaviors and to systems-level approaches to gene regulation. The formation of the calcified endoskeleton of the sea urchin embryo is a valuable experimental system for developing such an integrated view of the genomic regulatory control of morphogenesis. A transcriptional gene regulatory network (GRN) that underlies the specification of skeletogenic cells (primary mesenchyme cells, or PMCs) has recently been elucidated. In this study, we carried out a genome-wide analysis of mRNAs encoded by effector genes in the network and uncovered transcriptional inputs into many of these genes. We used RNA-seq to identify >400 transcripts differentially expressed by PMCs during gastrulation, when these cells undergo a striking sequence of behaviors that drives skeletal morphogenesis. Our analysis expanded by almost an order of magnitude the number of known (and candidate) downstream effectors that directly mediate skeletal morphogenesis. We carried out genome-wide analysis of (1) functional targets of Ets1 and Alx1, two pivotal, early transcription factors in the PMC GRN, and (2) functional targets of MAPK signaling, a pathway that plays an essential role in PMC specification. These studies identified transcriptional inputs into >200 PMC effector genes. Our work establishes a framework for understanding the genomic regulatory control of a major morphogenetic process and has important implications for reconstructing the evolution of biomineralization in metazoans.

The progressive changes in form that characterize embryogenesis are encoded in the genome. The properties of cells that drive these changes in form, like other specialized cellular properties, arise as a consequence of differential gene expression. Programs of differential gene expression can be viewed as dynamic networks of regulatory genes (genes that encode transcription factors, or TFs), and the cis-regulatory DNA elements to which TFs bind. Such gene regulatory networks (GRNs) are proving to be powerful tools for analyzing cell specification and the evolution of development (Stathopoulos and Levine, 2005; Davidson, 2010; Peter and Davidson, 2011; Van Nostrand and Kim, 2011; Wunderlich and DePace, 2011). A current limitation of this conceptual approach to development, however, is that we have a poor understanding of the connections between transcriptional networks and the morphogenetic processes that build tissues and organs. A marriage of regulatory network biology and morphogenesis will require experimental models that are amenable both to systems-level approaches and to detailed analysis of morphogenetic mechanisms. Integrating transcriptional networks and morphogenesis will also be crucial in an evolutionary context, i.e. for understanding how evolutionary modifications to genetic programs have supported changes in animal anatomy (Ettensohn, 2013).

The endoskeleton of the sea urchin embryo provides an opportunity to elucidate the genetic circuitry that underlies the formation of a major anatomical feature. The skeleton is a biomineral composed of calcium carbonate (in the form of calcite) and small amounts of occluded proteins. It is secreted by primary mesenchyme cells (PMCs), a population of cells derived from the large micromeres (LMs) of the cleavage stage embryo. During gastrulation, PMCs undergo a sequence of morphogenetic behaviors that includes ingression (epithelial-mesenchymal transition), directional migration, cell-cell fusion and biomineral formation (Wilt and Ettensohn, 2007; Ettensohn, 2013). These cellular behaviors have been analyzed in detail in living embryos (Gustafson and Wolpert, 1967; Malinda et al., 1995; Miller et al., 1995; Guss and Ettensohn, 1997; Peterson and McClay, 2003; Hodor and Ettensohn, 2008; Adomako-Ankomah and Ettensohn, 2013). The skeleton has several important functions; it influences the shape, orientation, swimming and feeding of the larva (Pennington and Strathmann, 1990; Hart and Strathmann, 1994), and its growth during larval development is responsive to environmental cues (Adams et al., 2011).

A GRN that underlies skeletogenic specification is activated in the LM-PMC lineage by localized maternal factors (Emily-Fenouil et al., 1998; Wikramanayake et al., 1998; Logan et al., 1999; Weitzel et al., 2004; Ettensohn, 2006). These maternal inputs function cell-autonomously to drive the zygotic expression of a small number of lineage-specific TFs, including Ets1 (Kurokawa et al., 1999) and Alx1 (Ettensohn et al., 2003). Early TFs in the GRN engage additional layers of regulatory genes, and various feedback and feed-forward interactions stabilize the network and drive it forward (Oliveri et al., 2008). Although considerable information is available concerning interactions among regulatory genes in the network, we have a very limited understanding of the downstream effector genes that execute skeletal morphogenesis and their transcriptional control.

In previous studies, we used an in situ hybridization screen to identify candidate effector genes in the PMC GRN and analyzed the developmental functions and regulatory control of several of these genes (Illies et al., 2002; Cheers and Ettensohn, 2005; Livingston et al., 2006; Adomako-Ankomah and Ettensohn, 2011; Rafiq et al., 2012; Adomako-Ankomah and Ettensohn, 2013). Here, we expand this analysis to a genome-wide level by carrying out an RNA-seq-based analysis of effector genes in the PMC GRN and by identifying key regulatory inputs (direct or indirect) into many of these genes. We increase by approximately an order of magnitude the number of known PMC-enriched transcripts. The great majority of these encode effector proteins, many with known or predicted functions, whereas others encode newly identified, PMC-specific TFs. We find that Ets1 and/or Alx1 provide essential regulatory inputs into >50% of the genes differentially expressed by PMCs at the early gastrula stage, pointing to the pivotal role of these TFs in controlling the cell-specific identity of PMCs. Genome-wide mRNA profiling of embryos treated with the MEK inhibitor U0126, which blocks PMC specification by inhibiting the phosphorylation of Ets1 (Fernandez-Serra et al., 2004; Röttinger et al., 2004), reveals that the PMC GRN is a major target of MAPK signaling during early embryogenesis and shows that Ets1 and Alx1 are key mediators of MAPK inputs into the GRN. Overall, this work greatly expands our understanding of the architecture and regulation of the PMC GRN and extends the utility of this experimental system as a model for developing an integrated view of the genomic regulatory control of morphogenesis.

RNA-seq analysis of mRNAs differentially expressed by PMCs at the onset of gastrulation

Our previous work focused on a subset of highly abundant, PMC-enriched transcripts (Rafiq et al., 2012). To obtain a more global perspective, we used RNA-seq to compare the abundance of transcripts in PMCs and a non-PMC fraction at the early mesenchyme blastula stage 24 hours post-fertilization (hpf). At this stage of development, PMCs are the only cells that have ingressed into the blastocoel. Thus, we enriched PMCs by isolating basal lamina bags from embryos at this stage (Harkey and Whiteley, 1980). Most effector genes in the PMC GRN are expressed at 24 hpf (Rafiq et al., 2012).

We compared the expression of ∼21,000 distinct S. purpuratus transcripts (Tu et al., 2012) in PMCs and the non-PMC fraction. Most mRNAs were expressed at similar levels (Fig. 1; R2=0.91, P<2×10-16). Cuffdiff analysis identified 420 transcripts with expression levels that differed significantly in the PMC and non-PMC samples (supplementary material Tables S1 and S2). All but five of these mRNAs were more abundant in PMCs than in the non-PMC sample. We refer to the genes that encode this collection of 420 mRNAs as the differentially expressed (DE) gene set. A summary of information concerning the 420 DE genes is presented in supplementary material Table S1 and quantitative expression values for all S. purpuratus transcripts from the RNA-seq analysis are provided in supplementary material Table S2. Overall, RNA-seq-based gene expression profiling increased by approximately an order of magnitude the number of known PMC-enriched mRNAs and therefore provided a far more complete picture of the output of this transcriptional network than was previously available.

Fig. 1.

Linear scatter plot of FPKM values derived from RNA-seq analysis of genes expressed at 24 hpf in PMCs and a non-PMC fraction isolated from S. purpuratus embryos. FPKM values shown are the means of two biological replicates and range from <1 to 5564 in purified PMCs and <1 to 3524 in the non-PMC fraction (‘other cells’), for the ∼21,000 transcripts detected. R2=0.91, P<2x10-16.

Fig. 1.

Linear scatter plot of FPKM values derived from RNA-seq analysis of genes expressed at 24 hpf in PMCs and a non-PMC fraction isolated from S. purpuratus embryos. FPKM values shown are the means of two biological replicates and range from <1 to 5564 in purified PMCs and <1 to 3524 in the non-PMC fraction (‘other cells’), for the ∼21,000 transcripts detected. R2=0.91, P<2x10-16.

To assess the completeness of our analysis (i.e. the false-negative rate), we examined a set of 36 mRNAs that were previously reported to be restricted to PMCs at this developmental stage, based on a WMISH screen and literature survey [see table S2 in the supplementary material of Rafiq et al. (Rafiq et al., 2012)]. Of these transcripts, only 4/36 (11%) were not found in the collection of DE genes. The four mRNAs that were not identified (ctd, p19, sm37 and stomatin) all yielded FPKM values in PMCs that were higher than in the non-PMC cell fraction, but the data failed to meet the significance criteria of the Cuffdiff analysis. This sampling indicates that, although the collection of DE genes is not exhaustive, it is likely to have captured the great majority of transcripts that are differentially expressed by PMCs at the mesenchyme blastula stage.

Because the significance threshold of the Cuffdiff analysis was relatively stringent (estimated false discovery rate=0.05), it seems likely that the DE set contains few false positives. We took two approaches to further assess the frequency of false positives in the DE gene set. First, we examined 50 genes chosen at random from those S. purpuratus genes annotated with Gene Ontology (GO) terms associated with metabolism, DNA replication, protein translation and other likely housekeeping functions. None of these genes was found to be differentially expressed in our analysis. Second, we used WMISH to analyze the expression patterns of 41 DE genes (these were selected because examination of the predicted gene products suggested a possible role in skeletal morphogenesis, as discussed below). WMISH analysis confirmed that 25 of these mRNAs were enriched in PMCs at the mesenchyme blastula stage (and, in most cases, at later developmental stages) (Fig. 2). The remaining WMISH probes, most of which were directed against low-abundance transcripts, showed uniformly low levels of staining. Based on analysis of >200 PMC-enriched mRNAs (Rafiq et al., 2012 and this study), we found that the threshold for WMISH detection was an FPKM value of 5-10 in whole embryo samples, which corresponded to 4-7 transcripts/PMC (assuming 32 PMCs/embryo). In general, RNA-seq data agreed well with WMISH analysis, i.e. transcripts that were (1) relatively abundant as indicated by a high FPKM value and (2) highly enriched in PMCs as reflected by a high log2-fold difference (supplementary material Table S1), yielded robust WMISH patterns that were restricted to PMCs at the mesenchyme blastula stage.

Fig. 2.

WMISH analysis of mRNAs differentially expressed in PMCs. Lateral views are shown for all stages and gastrula stage embryos are also viewed along the animal-vegetal axis (far right panels; ventral side down). The gene names shown are those assigned in SpBase. For genes that have not been assigned descriptive names, SPU identifiers (also from SpBase) are shown. Also included is a previously unannotated gene, located on Scaffold17:88148-92454 (v.3.1), which encodes a member of a family of PMC-specific Ig/TM proteins. BL, blastula; MB, mesenchyme blastula; EG, early gastrula; MG, mid-gastrula; LG, late gastrula; G-vv, mid- to late gastrula-vegetal view.

Fig. 2.

WMISH analysis of mRNAs differentially expressed in PMCs. Lateral views are shown for all stages and gastrula stage embryos are also viewed along the animal-vegetal axis (far right panels; ventral side down). The gene names shown are those assigned in SpBase. For genes that have not been assigned descriptive names, SPU identifiers (also from SpBase) are shown. Also included is a previously unannotated gene, located on Scaffold17:88148-92454 (v.3.1), which encodes a member of a family of PMC-specific Ig/TM proteins. BL, blastula; MB, mesenchyme blastula; EG, early gastrula; MG, mid-gastrula; LG, late gastrula; G-vv, mid- to late gastrula-vegetal view.

Characterization of DE genes

Most of the 420 DE genes were expressed at relatively low levels at the mesenchyme blastula stage (<10 transcripts/PMC, assuming 16 PMCs/embryo at this stage; supplementary material Fig. S1). Only 39 (∼10%) of the PMC-enriched mRNAs were expressed at levels greater than 25 transcripts/PMC. Not surprisingly, many of these abundant transcripts encoded biomineralization proteins, including three members of the Msp130 family (Msp130, Msp130r1 and Msp130r2), six spicule matrix proteins [Sm20 (Clect_14), Sm29, Sm30, Sm49 (C-lectin/PMC1), Sm32/50 and C-lectin], and other biomineralization proteins such as P16, P16rel2 (Hypp_2998), P58A (FcgbpL) and carbonic anhydrase (Cara7LA) (Livingston et al., 2006; Rafiq et al., 2012).

The temporal expression profiles of the 420 DE genes were extracted from transcriptome profiling data of Tu et al. (Tu et al., 2012) and analyzed by hierarchical clustering (Fig. 3). This analysis revealed subsets of DE genes with coordinated temporal expression profiles, including (from top to bottom in Fig. 3): (1) genes with high levels of maternal expression; (2) a small set of genes that showed a sharp spike in expression during late cleavage (10 hpf); (3) genes that showed maximal expression during late gastrulation and post-gastrula stages (48-72 hours); and (4) a major class of genes (almost half) that were expressed at very low levels early in development, peaked in expression during the late blastula-gastrula stages (18-40 hpf), and then declined in expression (genes in the lower half of Fig. 3).

Fig. 3.

Hierarchical clustering of temporal expression patterns of genes differentially expressed in PMCs. The temporal expression profiles of 420 genes differentially expressed in PMCs were obtained from the RNA-seq data of Tu et al. (Tu et al., 2012), available on SpBase. Each gene is represented by a single row and each time point by a single column. The color scale ranges from deep red (2.5-fold higher than mean expression) to deep blue (2.5-fold lower than mean expression). White indicates the mean expression value. For reference, 24, 48 and 72 hpf correspond to the mesenchyme blastula, late gastrula and late prism stages, respectively.

Fig. 3.

Hierarchical clustering of temporal expression patterns of genes differentially expressed in PMCs. The temporal expression profiles of 420 genes differentially expressed in PMCs were obtained from the RNA-seq data of Tu et al. (Tu et al., 2012), available on SpBase. Each gene is represented by a single row and each time point by a single column. The color scale ranges from deep red (2.5-fold higher than mean expression) to deep blue (2.5-fold lower than mean expression). White indicates the mean expression value. For reference, 24, 48 and 72 hpf correspond to the mesenchyme blastula, late gastrula and late prism stages, respectively.

To gain insight into the possible roles of DE genes, we first examined the functional assignments of these genes as annotated in SpBase (Cameron et al., 2009). Forty-six percent of the DE genes have been assigned to functional categories based on hand annotations (Sodergren et al., 2006) or primary GO terms derived by blast2go (Tu et al., 2012). Fig. 4 shows the distribution of these functional classes. Consistent with the skeletogenic function of PMCs, one of the largest classes of DE genes was the biomineralization set. It seems likely that many of the DE genes without functional annotations (as well as genes currently annotated as ‘novel’) also encode biomineralization-related proteins. In addition, some genes in the ‘calcium toolkit’, ‘kinase’ and ‘metalloprotease’ classes play important roles in biomineral formation. As an independent means of assessing the subset of DE genes with functions in biomineralization, we examined a set of ∼200 proteins identified in a recent proteomic analysis of partially purified embryonic spicules (Mann et al., 2010) and identified 62 gene products that were common to the two sets (supplementary material Table S1). Adhesion-related proteins constituted another sizable functional category, and many of these proteins are likely to be involved in PMC migration and/or fusion (see below).

Fig. 4.

Distribution of DE genes by functional class. The distribution is based on the primary functional assignments of DE genes in their public annotations (SpBase). Functional assignments are based on hand annotation (Sodergren et al., 2006) and, where lacking, on primary GO terms derived by blast2go (Tu et al., 2012). Out of the 420 DE genes, 194 have been assigned to functional categories. Novel genes, biomineralization genes and adhesion genes constitute almost half of this set. The y-axis indicates the number of genes in each functional class.

Fig. 4.

Distribution of DE genes by functional class. The distribution is based on the primary functional assignments of DE genes in their public annotations (SpBase). Functional assignments are based on hand annotation (Sodergren et al., 2006) and, where lacking, on primary GO terms derived by blast2go (Tu et al., 2012). Out of the 420 DE genes, 194 have been assigned to functional categories. Novel genes, biomineralization genes and adhesion genes constitute almost half of this set. The y-axis indicates the number of genes in each functional class.

We examined the complete set of 368 annotated TFs in the S. purpuratus genome and found that more than half were detectable in basal lamina bag-purified PMCs at levels >1 transcript/cell. Several TFs in this set, however, including gcm, foxa, scl and gataC, are restricted to cell types other than PMCs at the mesenchyme blastula stage (Ransick and Davidson, 2006; Croce and McClay, 2010; Flynn et al., 2011) and their identification in our analysis reflected the low level of contamination of the bag preparations with cell types other than PMCs. When we restricted our analysis to the set of 420 DE genes, we identified only 11 TFs (Table 1). One of these (evi1) was expressed at lower levels in PMCs than in other cells; the other ten TF mRNAs (alx1, alx4, cebpa, foxB, foxO, mef2, mitf, nk7, smad2/3 and tbr) were enriched in PMCs relative to the non-PMC population to varying degrees, ranging from 3.5-fold (foxO) to 15-fold (nk7). WMISH data are available for seven of these genes [alx1 (Ettensohn et al., 2003), alx4 (Rafiq et al., 2012), foxb (Minokawa et al., 2004), foxO (Tu et al., 2006), smad2/3 (Poustka et al., 2007) tbr (Fuchikami et al., 2002) and nk7 (this study, Fig. 2)] and in all cases confirms that expression is enriched in PMCs at the mesenchyme blastula stage.

Table 1.

Transcription factors identified in this analysis

Transcription factors identified in this analysis
Transcription factors identified in this analysis

Hand curation of the set of DE genes revealed many new candidate effectors of skeletal morphogenesis, some of which are highlighted below.

Biomineralization proteins

Transport/channel proteins

The deposition of CaCO3 by PMCs is associated with the uptake of Ca2+ and HCO -3 ions from the blastocoel (Stumpp et al., 2012). We identified five solute carrier proteins (members of the Slc4, Slc5, Slc10, Slc24 and Slc26 families) that might mediate these transport functions.

Secreted metalloproteases

Metalloprotease inhibitors reversibly block spiculogenesis by PMCs in vivo and in vitro (Roe et al., 1989; Ingersoll and Wilt, 1998). The DE genes include a suite of four matrix metalloprotease genes, arranged in tandem on a single chromosome, that are likely to encode the relevant enzymes. Fig. 2 shows WMISH data for two of these genes (mmpl2 and mmpl6).

Fam20C

This secreted kinase was recently shown to phosphorylate extracellular biomineralization proteins in vertebrates (Ishikawa et al., 2012; Tagliabracci et al., 2012).

Otopetrin

Otopetrin (otop2L) is a transmembrane (TM) protein essential for the development of otoliths/otoconia, which are extracellular calcium carbonate-containing crystals that mediate vestibular mechanosensory function in vertebrates (Hurle et al., 2003; Hughes et al., 2004; Söllner et al., 2004).

Adhesion/migration proteins

Nephronectin

During their migration, PMCs interact with extracellular matrix (ECM) fibers that contains the sea urchin ortholog of vertebrate Frem2 (Hodor et al., 2000). Frem2 and related proteins have been implicated in epithelium-mesenchyme adhesion and mutations in these genes underlie Fraser’s syndrome (Smyth and Scambler, 2005). Frem proteins are required for the proper incorporation of nephronectin, an integrin ligand required for kidney morphogenesis, into the ECM (Kiyozumi et al., 2012). The sea urchin ortholog of nephronectin (npnt) is expressed selectively by PMCs during gastrulation and might play an important role in PMC-substrate interactions.

Adhesion receptors

PMCs selectively express several type I TM proteins with variable numbers of extracellular Ig, Egf, Lrr and Fn3 repeats, an organization which suggests that these proteins might function as adhesion receptors. Examples include Lrr/Igr_10, Fn3/Egff_1 and Fn3f_9.

Aquaporin-9

One abundant, PMC-specific transcript (aqp9) encodes a member of the aquaporin family of TM, water channel proteins, which have recently been implicated in regulating the protrusive activity and migration of cancer cells (Verkman, 2011).

Cell-cell fusion proteins

The dynamics of PMC fusion have been analyzed in vivo (Hodor and Ettensohn, 1998), but molecules that mediate fusion have not been identified. In Drosophila, TM proteins with multiple extracellular Ig domains (Sns, Rst and Duf) are required for myoblast interactions prior to fusion (Abmayr and Pavlath, 2012). We have identified four PMC-specific, type I TM proteins with multiple extracellular Ig repeats that are the closest relatives of Sns/Rst/Duf in the sea urchin genome and strong candidates for regulators of PMC fusion. WMISH data for three of these genes (kirre1L, SPU_026000 and Scaffold17:88148-92454) are shown in Fig. 2.

Transcriptional inputs into DE genes

To identify regulatory inputs into the 420 DE genes, we used RNA-seq to analyze changes in gene expression following knockdown of Ets1 or Alx1 (supplementary material Tables S3 and S4). Ets1 and Alx1 are pivotal early TFs in the PMC GRN (Kurokawa et al., 1999; Ettensohn et al., 2003). RNA-seq was used to profile gene expression in controls and morphants at 28-30 hpf (early gastrula). We chose this stage because the severity of the morphant phenotypes could be unambiguously scored (see Materials and Methods) and because an earlier, more limited, analysis of gene expression changes in Ets1 and Alx1 morphants was carried out at this stage (Rafiq et al., 2012). Most effector genes in the PMC GRN are robustly expressed at 28-30 hpf (Rafiq et al., 2012) (Fig. 3).

To assess the reliability of our RNA-seq analysis, we compared QPCR data from a previous study that examined the effects of Ets1 and Alx1 knockdowns on the expression of ∼20 effector genes (Rafiq et al., 2012) with Nanostring and RNA-seq-based expression data obtained in the present study for the same set of genes at the same developmental stage. This analysis showed that effects of knockdowns on gene expression were highly reproducible across these experiments, which were carried out using embryos derived from three different male-female matings and which used three different methods of transcript quantification (supplementary material Fig. S2).

RNA-seq-based gene expression profiling showed that 223/420 DE genes (53%) were significantly affected by knockdown of Ets1 and/or Alx1 (Fig. 5; supplementary material Table S1). This demonstrated the pivotal role of these TFs in controlling the cell-specific identity of PMCs. Of the DE genes with inputs from Ets1 or Alx1, most (144/223, or ∼2/3) were downregulated in both classes of morphants.

Fig. 5.

Venn diagram showing overlapping distributions of genes affected by Ets1 knockdown, Alx1 knockdown, or U0126 treatment among the 420 genes differentially expressed by PMCs. More than half (223/420, 53%) of DE genes are affected by knockdown of Ets1 and/or Alx1; the great majority of these inputs (∼90%) are positive. Of these 223 DE genes, ∼2/3 (144/223) are affected in both classes of morphants. 101 DE genes are sensitive to U0126, a number that includes more than half of all U0126 targets genome-wide. Most of the U0126-sensitive DE genes have inputs (direct or indirect) from Ets1 and/or Alx1.

Fig. 5.

Venn diagram showing overlapping distributions of genes affected by Ets1 knockdown, Alx1 knockdown, or U0126 treatment among the 420 genes differentially expressed by PMCs. More than half (223/420, 53%) of DE genes are affected by knockdown of Ets1 and/or Alx1; the great majority of these inputs (∼90%) are positive. Of these 223 DE genes, ∼2/3 (144/223) are affected in both classes of morphants. 101 DE genes are sensitive to U0126, a number that includes more than half of all U0126 targets genome-wide. Most of the U0126-sensitive DE genes have inputs (direct or indirect) from Ets1 and/or Alx1.

We compared the temporal expression profiles of two cohorts of DE genes: (1) those that were affected both by Ets1 and Alx1 knockdowns, and (2) those that were not regulated by either TF (Fig. 6). Hierarchical clustering revealed that the Ets1/Alx1-regulated gene set contained few genes that exhibited high levels of maternal transcripts and, more strikingly, the majority of these genes had a strong spike in expression between the late blastula and mid-gastrula stages (18-30 hpf) (Fig. 6A). By contrast, DE genes that were not regulated by Ets1 or Alx1 showed a much broader distribution of expression patterns, with peak expression levels distributed relatively evenly across all developmental stages, and many genes showed high levels of maternal expression (Fig. 6B).

Fig. 6.

Distinct temporal gene expression profiles of Ets1/Alx1 co-regulated targets and non-target genes in the DE set. Hierarchical clustering of the temporal expression profiles of (A) 143 DE genes that are sensitive both to Ets1 and Alx1 knockdown and (B) 198 DE genes that are not regulated by either Ets1 or Alx1 (see Fig. 3 legend for details). The Ets1/Alx1-regulated gene set contains few genes that exhibit high levels of maternal transcripts and most genes show maximal expression between the late blastula and mid-gastrula stages (18-30 hpf). By contrast, DE genes that are not affected by Ets1 or Alx1 knockdowns show a much broader distribution of temporal expression patterns, including many cases of high maternal expression.

Fig. 6.

Distinct temporal gene expression profiles of Ets1/Alx1 co-regulated targets and non-target genes in the DE set. Hierarchical clustering of the temporal expression profiles of (A) 143 DE genes that are sensitive both to Ets1 and Alx1 knockdown and (B) 198 DE genes that are not regulated by either Ets1 or Alx1 (see Fig. 3 legend for details). The Ets1/Alx1-regulated gene set contains few genes that exhibit high levels of maternal transcripts and most genes show maximal expression between the late blastula and mid-gastrula stages (18-30 hpf). By contrast, DE genes that are not affected by Ets1 or Alx1 knockdowns show a much broader distribution of temporal expression patterns, including many cases of high maternal expression.

In parallel studies, we carried out RNA-seq transcriptional profiling of 28- to 30-hour embryos that had been treated from the 2-cell stage with U0126, a MEK inhibitor that blocks PMC specification by inhibiting the phosphorylation of Ets1 (Fernandez-Serra et al., 2004; Röttinger et al., 2004). Genome-wide, we identified 180 transcripts that showed significant changes of expression in U0126-treated embryos (supplementary material Table S5). Remarkably, the majority of these mRNAs (101/180, 56%) were also DE genes, suggesting that the PMC GRN is the principal target of MAPK signaling during early embryogenesis (Fig. 5). All 101 of the DE mRNAs significantly affected by U0126 treatment were downregulated in the presence of the inhibitor. Of these, the large majority (83%) were also regulated by Ets1 and/or Alx1 (67% were affected in both classes of morphants), pointing to these two TFs as key mediators of MAPK inputs into the GRN (Fig. 5).

Non-DE genes

The expression of many (>1500) genes not in the DE set was also significantly affected by Ets1/Alx1 knockdowns and/or U0126 treatment (see supplementary material Tables S3-S5). These might include early targets of Ets1/MAPK in the non-skeletogenic mesoderm (NSM) (Fernandez-Serra et al., 2004; Röttinger et al., 2004), but most of the gene expression changes are probably indirect and reflect the additive effects of complex tissue interactions. As an initial step in analyzing these targets, we focused primarily on the suite of all regulatory genes, as these are well annotated and WMISH data are available for almost all regulatory genes expressed at detectable levels during embryogenesis.

Oliveri et al. (Oliveri et al., 2008) showed that one function of alx1 in the LM progeny is to repress gcm, a regulatory gene ordinarily expressed by adjacent, presumptive pigment cells. Our RNA-seq analysis confirmed an increase in gcm expression in Alx1 morphants. We found that Alx1 knockdown resulted in a significant upregulation of 23 other regulatory genes. WMISH data are available for eight of these genes, and a surprisingly large fraction (6/8) are expressed selectively in the NSM during normal development. Four genes (scl, lmo2t, rxr/Z177 and sna) are expressed by blastocoelar cells (a population of presumptive immunoctyes), whereas six1/2 and soxE are expressed by pigment cells and coelomic pouch cells (probably prospective germ cells), respectively. In addition, two of the genes for which WMISH data are unavailable, irf4 and nfil3, have vertebrate orthologs that play important roles in immune cell development, suggesting that these mRNAs might also be restricted to the blastocoelar cell lineage during normal development. We also examined the set of ∼100 non-DE effector genes that are upregulated both in Alx1 and Ets1 morphants and identified several proteins that are predicted to function in immune system development or physiology, including two Toll receptors (Sp-TlrL_9 and Sp-Tlr072) and a leukocyte receptor cluster member (Sp-Leng9L). Although further analysis of the expression patterns of these regulatory and effector genes in control embryos and morphants will be required, these findings support the view that a key function of alx1 is to repress multiple, alternative mesodermal regulatory states, including the blastocoelar cell fate, in the LM progeny.

Surprisingly, Alx1 and Ets1 morphants exhibited a significant downregulation of hox7, a regulatory gene expressed in the aboral ectoderm, as well as spec2c and spec2ce1-3, two aboral ectoderm differentiation markers (supplementary material Tables S3 and S4). Two other aboral ectoderm regulatory genes, dlx and msx, were downregulated in Ets1 morphants. Our data therefore point to a previously unsuspected interaction between LM progeny and the aboral ectoderm that occurs before the early gastrula stage (i.e. the stage at which we analyzed gene expression).

A complex sequence of PMC behaviors underlies the morphogenesis of the embryonic skeleton (Wilt and Ettensohn, 2007; Ettensohn, 2013). These behaviors require zygotic transcriptional inputs (Kurokawa et al., 1999; Ettensohn et al., 2003; Wu and McClay, 2007). Our work has provided the most complete picture to date of the effector genes that direct skeletogenesis and has revealed important features of the transcriptional control of these genes.

The identification of morphogenetic effector genes

The morphogenetic functions of some PMC effector genes are well understood. The spicule matrix proteins are a family of 15-20 closely related proteins occluded within the biomineral that influence its growth and physical properties, probably by regulating the conversion of amorphous calcium carbonate to the crystalline state (Wilt and Ettensohn, 2007; Gong et al., 2012; Rafiq et al., 2012). Non-fibrillar collagens produced by PMCs serve as an essential substrate for the cells (Wessel et al., 1991). Several PMC-specific, type I TM proteins, including P16, P58A and P58B, play essential roles in biomineral deposition (Cheers and Ettensohn, 2005; Adomako-Ankomah and Ettensohn, 2011). The precise biochemical functions of the P16 and P58 proteins are unknown, although P16 is phosphorylated and binds to hydroxyapatite (Alvares et al., 2009). A PMC-specific, GPI-anchored carbonic anhydrase is likely to be involved in biomineral remodeling (Livingston et al., 2006). All these proteins are associated with biomineralization, the major specialized function of the PMCs. Only one effector, VEGFR-Ig10, has been shown to mediate other aspects of the morphogenetic program of PMCs. This signaling receptor plays an essential role in PMC guidance; in addition, local VEGF signals acting through VEGFR-Ig10 control regional patterns of skeletal growth (Duloquin et al., 2007; Adomako-Ankomah and Ettensohn, 2013).

Our RNA-seq-based analysis has increased by approximately an order of magnitude the number of known PMC-enriched mRNAs and therefore provides a more complete picture of the output of this transcriptional network than was previously available. We identified a large number of putative effectors of skeletal morphogenesis that are candidates for further functional studies. In some cases (e.g. Fam20C and Otopetrin, see below), functions can be inferred from information concerning the vertebrate counterparts of these genes. Our work has also revealed specific proteins that are likely to account for pharmacological evidence that metalloproteases (Roe et al., 1989; Ingersoll and Wilt, 1998), calcium channels (Hwang and Lennarz, 1993) and ion transporters (Yasumasu et al., 1985; Mitsunaga et al., 1986; Fujino et al., 1987; Stumpp et al., 2012) are essential effectors of skeletogenesis.

Regulatory inputs into effector genes

Our findings have revealed important features of the regulatory inputs into the set of 420 effector genes. We focused on two key TFs in the PMC specification network: Ets1 and Alx1. These TFs provide regulatory inputs near the top of the regulatory network and are essential for PMC specification. Knockdown of Ets1 or Alx1 causes LM descendants to take on alternative mesodermal fates (Kurokawa et al., 1999; Ettensohn et al., 2003; Ettensohn et al., 2007; Oliveri et al., 2008). Zygotic expression of both TFs is restricted to the LM lineage early in development; alx1 transcription is activated selectively in LMs in the first interphase after these cells are born.

We found that of the 420 DE genes, more than half (223/420, 53%) received essential inputs from Ets1 and/or Alx1 (Fig. 5), the great majority of which (∼90%) were positive. When only the most abundant mRNAs are considered, this value increased to 74% (i.e. 74/100 DE transcripts with the highest FPKM values in purified PMCs). We also noted that of the DE genes annotated with the GO terms ‘biomineralization’ or ‘metalloprotease’, 84% (32/38) were subject to regulatory inputs from one or both of these TFs. These findings demonstrate the central role of Ets1 and Alx1 in controlling the cell-specific identity of PMCs. At the same time, our analysis identified 197 DE genes that were not significantly affected by Ets1 or Alx1 knockdown. This number is likely to be inflated by the stringency of the Cuffdiff analysis; for example, many of these mRNAs showed modest changes in expression in morphants (e.g. 50-75% reduction in mRNA level) that were scored as non-significant. More importantly, we can assume that the MO knockdowns were incomplete. With these caveats in mind, we identified ∼150 DE mRNAs, many of which were very abundant, that showed changes in expression of <50% in both Ets1 and Alx1 morphants relative to controls. These findings indicate that Ets1/Alx1-independent circuits also make contributions to the specialized molecular properties of PMCs.

One of the most striking findings from this and previous work (Rafiq et al., 2012) is that many effector genes are regulated positively by both Ets1 and Alx1. Of the 223 DE genes with inputs from Ets1 and/or Alx1, ∼2/3 (144/223) were affected in both classes of morphants (Fig. 5). Several mechanisms might underlie this apparent co-regulation. First, Ets1 might regulate effector genes indirectly through its effect on Alx1 expression. Perturbation of Ets1 function does not affect the early phase of alx1 expression, but suppresses the later phase (Oliveri et al., 2008; Sharma and Ettensohn, 2010). In our study, Ets1 knockdown reduced alx1 expression by 80%, whereas Alx1 knockdown had no effect on ets1 expression. Moreover, of the 170 DE genes that were regulated by Ets1, 85% showed significant changes in expression following Alx1 knockdown. Thus, most of the effects of Ets1 knockdown might be explained through the effect of Ets1 on Alx1 expression. It was reported previously that forced expression of Alx1 is unable to rescue the effects of Ets1 knockdown (Oliveri et al., 2008), which appears inconsistent with this model, but subsequent studies have shown that the effects of Alx1 are highly dosage dependent (Ettensohn et al., 2007; Damle and Davidson, 2011). Second, Ets1 and Alx1 might regulate effector genes in concert, via a feed-forward mechanism (e.g. Ets1>Alx1, Alx1>Effector X, Ets1>Effector X). Experimental evidence in support of this model has come from analysis of the cis-regulatory control of cyp1, which receives direct inputs from Dri (a target of alx1) and Ets1 (Amore and Davidson, 2006) and sm50, which receives direct inputs from Ets1 (Yajima et al., 2010). Third, Ets1 and Alx1 might regulate the expression of a common intermediary TF that provides essential inputs into many effector genes. If this is the case, then knockdowns of Ets1 and Alx1 would be expected to produce similar effects. Candidates include regulatory genes in the DE gene set that are downregulated both in Ets1 and Alx1 morphants (alx4, cebpa, foxB and nk7).

MAPK signaling and the PMC GRN

The MAPK pathway plays a crucial role in PMC specification. Previous studies documented a transient, localized activation of ERK in the LM lineage shortly before ingression and showed that U0126, a selective MEK inhibitor, blocks PMC ingression and the expression of several terminal differentiation genes (Fernandez-Serra et al., 2004; Röttinger et al., 2004). The ERK/MAPK pathway is not required for the maternally driven activation of the network or the initial expression of early regulatory genes such as alx1, but becomes active at the blastula stage, when it functions to maintain the expression of alx1 and possibly other regulatory genes (Fernandez-Serra et al., 2004; Sharma and Ettensohn, 2010). The activation of ERK in the LM-PMC lineage does not require signals from other cell populations (Fernandez-Serra et al., 2004; Röttinger et al., 2004). Significantly, Röttinger et al. (Röttinger et al., 2004) showed that Ets1 (which contains consensus MAPK phosphorylation and ERK docking sites) is a direct target of MEK/ERK signaling.

Our RNA-seq analysis showed that a surprisingly small fraction of the transcriptome is dependent upon MAPK signaling during early embryogenesis. We identified only 180 transcripts that exhibited significant changes in expression at 28-30 hpf in response to MEK inhibition. Strikingly, more than half of these transcripts (101/180) were contained in the DE collection. Our data are consistent with immunostaining studies indicating that ERK is selectively activated in the LM lineage and support the view that the PMC GRN is the principal target of MAPK signaling during early development. Later in gastrulation, p-ERK is also detected at the tip of the archenteron, where it plays a role in the specification of non-skeletogenic mesoderm (Fernandez-Serra et al., 2004; Röttinger et al., 2004). Genes that are sensitive to U0126 but not differentially expressed in PMCs (79/180 genes) might be direct targets of MAPK signaling in NSM cells or indirect targets in tissues that are dependent upon PMCs for their normal development. Of the 101 DE genes sensitive to U0126, most (74/101, 73%) were also found to be regulated by Ets1, strongly supporting the view that Ets1 is the key mediator of MAPK inputs into the PMC GRN. We also identified a small number of DE transcripts (17) that were sensitive to U0126 but not to knockdown of either Ets or Alx1; these included tbr and mitf, mRNAs that encode PMC-restricted TFs. The mechanism by which MAPK/ERK signaling regulates the expression of these genes is unknown, although it is possible that the maternal pool of Ets1 protein (Yajima et al., 2010), which is not affected by MO knockdown and might be activated by MAPK, could be responsible. Our analysis also showed that >50% of the genes within the DE set that were sensitive to Ets1 knockdown were insensitive to U0126 (Fig. 5), suggesting that ERK-mediated phosphorylation is required for only a subset of the regulatory functions of Ets1. Lastly, our studies define a discrete, signal-dependent submodule of the larger genetic circuitry that controls PMC identity, represented by the subset (∼1/4) of DE genes sensitive to MEK inhibition.

The evolution of biomineralization

The further elucidation of the genetic network that underlies skeletogenic specification and morphogenesis in echinoderms has important implications for reconstructing the evolution of biomineralization in metazoans. The fossil record documents a widespread and relatively synchronous emergence of biomineralization in many metazoan lineages during the Cambrian period (Knoll, 2003; Murdock and Donohue, 2011). It is widely accepted that biomineralized structures, in the strictest sense, appeared independently in these lineages. For example, the first true mineralized vertebrate skeletons are thought to have appeared in ostracoderms, a group of stem gnathostomes, as a dermal skeleton, independently of the echinoderm skeleton (Donoghue and Sansom, 2002; Murdock and Donohue, 2011). An important unanswered question, however, concerns the extent to which this occurred by exploiting a common ‘toolkit’, i.e. a set of ancestral biochemical and developmental pathways that was independently co-opted for biomineral formation in diverse animal taxa (Westbroek and Marin, 1998; Jackson et al., 2007; Murdock and Donohue, 2011).

Our findings reveal new and surprising connections between genes that control biomineralization in modern echinoderms and vertebrates, despite the difference in biomineral content and micro-architecture in these taxa (Bottjer et al., 2006). We found that PMCs selectively express the single sea urchin member of the otopetrin family. Otopetrin 1 is required for the formation of calcite otoliths/otoconia in vertebrates (Hurle et al., 2003; Hughes et al., 2004; Söllner et al., 2004). The precise biochemical function of this 12-pass TM protein is unknown, but it might play a role in regulating cytosolic Ca2+ levels in response to extracellular signals (Kim et al., 2010; Kim et al., 2011). We also identified in the DE gene set the S. purpuratus ortholog of Fam20C, an extracellular kinase that phosphorylates extracellular biomineralization proteins in vertebrates (Ishikawa et al., 2012; Tagliabracci et al., 2012), as a protein differentially expressed by PMCs. Other classes of proteins with conserved functions in biomineralization in echinoderms and vertebrates include collagens, matrix metalloproteases and carbonic anhydrases (Livingston et al., 2006; Krane and Inada, 2008; Wuthier and Lipscomb, 2011). Alx1 family members play conserved roles as upstream transcriptional regulators of skeletogenesis in both taxa (Ettensohn et al., 2003). Our studies therefore reveal an extensive, common biomineralization toolkit that was likely to be present in the ancestral deuterostome and might have been exploited in diverse animal lineages.

Adult animals and embryo culture

Adult Strongylocentrotus purpuratus were obtained from Patrick Leahy (California Institute of Technology, USA). Spawning was induced by intracoelomic injection of 0.5 M KCl and embryos were cultured in artificial seawater (ASW) at 15°C in a temperature-controlled incubator.

PMC isolation

PMCs were isolated from early mesenchyme blastula stage embryos at 24 hours post-fertilization (hpf), as previously described (Harkey and Whiteley, 1980). Briefly, embryos were washed three times in calcium- and magnesium-free ASW (CMFSW), twice in 1 M glycine, and resuspended in bag isolation medium (per liter: 400 ml 1 M dextrose, 400 ml CMFSW, 200 ml distilled water). Embryos were dissociated by gentle pipetting. Basal lamina bags containing PMCs were collected using a sucrose step gradient. A ‘non-PMC’ (or ‘other cell’) fraction was collected from the same batch of embryos, also as described (Harkey and Whiteley, 1980). The same dissociation procedure was used except that embryos were washed only once in 1 M glycine to minimize rupturing of basal lamina bags. After resuspension in bag isolation medium, the sample was centrifuged at 650 g for 10 minutes, and the supernatant containing the non-PMC fraction was collected. The purity of isolated PMCs was >95% as assessed by immunostaining with monoclonal antibody 6a9 (Ettensohn and McClay, 1988). For analysis of transcript levels by RNA-seq, the PMC and non-PMC samples were isolated from two embryo cultures, derived from separate matings, which served as biological replicates.

Morpholino (MO) injections

MOs (Gene Tools) were injected into fertilized eggs as previously described (Cheers and Ettensohn, 2004), with the modification that eggs were fertilized in the presence of 0.1% (w/v) para-aminobenzoic acid to prevent hardening of the fertilization envelope. MO sequences (5′-3′) were: SpAlx1, TATTGAGTTAAGTCTCGGCACGACA; SpEts1, GAACAGTGCATAGACGCCATGATTG; control, CCTCTTACCTCAGTTACAATTTATA.

The Ets1 and Alx1 MOs, both of which are translation-blocking, have been shown to be specific and effective (Ettensohn et al., 2003; Oliveri et al., 2008; Rafiq et al., 2012). MOs were injected at an initial concentration of 2 mM (Ets1) or 4 mM (Alx1). Injection solutions also contained 20% (v/v) glycerol and 0.16% (w/v) Texas Red dextran. The control MO was injected at the same concentration as the corresponding translation-blocking MO. For comparisons of transcript levels in controls and Ets1/Alx1 morphants by RNA-seq, 500 embryos were pooled at 28-30 hpf for each sample. Because there was some embryo-to-embryo variability in morphant phenotypes, we hand selected embryos that lacked PMCs at the start of invagination. A single embryo culture was used for the analysis of Alx1 knockdown and a separate culture for the analysis of Ets1 knockdown.

RNA-seq

Total RNA was extracted using the NucleoSpin RNA II Kit (Clontech) and precipitated with ethanol. For comparisons of transcript levels in PMCs and other cells, RT-PCR was performed using the RETROScript Kit (Clontech) and primers for several PMC-specific transcripts (fc2, p133 and can1) and one housekeeping gene (z12), in order to confirm the expected difference in gene expression between the PMC and non-PMC samples. RNA samples were provided to the USC Epigenome Center and their quality was assessed using a BioAnalyzer. 600 ng-1 μg of total RNA was used for the construction of each Illumina HiSeq library. Sequencing was carried out with an Illumina HiSeq2000 machine (50 cycles, paired-end reads) with four to five indexed libraries in each lane. Approximately 40 million reads were obtained per sample. All data were analyzed using the open-source Tuxedo Suite (Langmead et al., 2009; Trapnell et al., 2012) with default parameters. TopHat (2.0.8b) was used to map sequence reads to the S. purpuratus transcriptome (Tu et al., 2012). The relative abundance of transcripts, represented by their FPKM (fragments per kilobase of transcript per million mapped reads) values, was estimated using Cufflinks (2.0.1). Cuffdiff (part of the Cufflinks package) was used to identify significant differences in the abundance of transcripts between samples (false discovery rate=0.05), and CummeRbund (2.0.0) (R/Bioconductor) was used for scatterplot analysis. All sequences were deposited in the NCBI Sequence Read Archive (SRA accession numbers SRP033427 and SRP031836).

Whole-mount in situ hybridization (WMISH)

Embryos were fixed for 1 hour at room temperature in 4% paraformaldehyde in ASW and stored at 4°C in 100% methanol. WMISH was carried out as previously described (Lepage et al., 1992; Duloquin et al., 2007).

Nanostring analysis

Quantitative analysis of transcript levels was carried out with a Nanostring nCounter system and a codeset corresponding to ∼90 genes in the PMC GRN, as previously described (Adomako-Ankomah and Ettensohn, 2013).

Analysis of temporal expression profiles

Temporal expression profiles of genes were obtained from the transcriptome data of Tu et al. (Tu et al., 2012) and were based on raw FPKM values at 0, 10, 18, 24, 30, 40, 48, 56, 64 and 72 hpf. Temporal expression patterns were analyzed by hierarchical clustering using Euclidean distances (MATLAB, MathWorks).

U0126 treatment

Embryos were treated with 7 μM U0126 (in DMSO) continuously from the 2-cell stage and sibling control embryos were treated with DMSO alone. At 28-30 hpf, total RNA was collected for RNA-seq analysis. Several control and UO126-treated embryos were immunostained with monoclonal antibody 6a9 to confirm that PMC specification was blocked by U0126 treatment, as reported previously (Fernandez-Serra et al., 2004; Röttinger et al., 2004).

We thank Noah Most and Matthew Early for their contributions to the WMISH analysis.

Author contributions

K.R., T.S. and C.A.E. designed the experiments. K.R. and T.S. conducted the experiments and prepared the figures. C.A.E. wrote the manuscript.

Funding

This work was supported by a National Science Foundation grant [IOS-1021805] to C.A.E.

Abmayr
S. M.
,
Pavlath
G. K.
(
2012
).
Myoblast fusion: lessons from flies and mice
.
Development
139
,
641
656
.
Adams
D. K.
,
Sewell
M. A.
,
Angerer
R. C.
,
Angerer
L. M.
(
2011
).
Rapid adaptation to food availability by a dopamine-mediated morphogenetic response
.
Nat. Commun.
2
,
592
.
Adomako-Ankomah
A.
,
Ettensohn
C. A.
(
2011
).
P58-A and P58-B: novel proteins that mediate skeletogenesis in the sea urchin embryo
.
Dev. Biol.
353
,
81
93
.
Adomako-Ankomah
A.
,
Ettensohn
C. A.
(
2013
).
Growth factor-mediated mesodermal cell guidance and skeletogenesis in the sea urchin embryo
.
Development
140
,
4214
4225
.
Alvares
K.
,
Dixit
S. N.
,
Lux
E.
,
Veis
A.
(
2009
).
Echinoderm phosphorylated matrix proteins UTMP16 and UTMP19 have different functions in sea urchin tooth mineralization
.
J. Biol. Chem.
284
,
26149
26160
.
Amore
G.
,
Davidson
E. H.
(
2006
).
cis-Regulatory control of cyclophilin, a member of the ETS-DRI skeletogenic gene battery in the sea urchin embryo
.
Dev. Biol.
293
,
555
564
.
Bottjer
D. J.
,
Davidson
E. H.
,
Peterson
K. J.
,
Cameron
R. A.
(
2006
).
Paleogenomics of echinoderms
.
Science
314
,
956
960
.
Cameron
R. A.
,
Samanta
M.
,
Yuan
A.
,
He
D.
,
Davidson
E.
(
2009
).
SpBase: the sea urchin genome database and web site
.
Nucleic Acids Res.
37
,
D750
D754
.
Cheers
M. S.
,
Ettensohn
C. A.
(
2004
).
Rapid microinjection of fertilized eggs
.
Methods Cell Biol.
74
,
287
310
.
Cheers
M. S.
,
Ettensohn
C. A.
(
2005
).
P16 is an essential regulator of skeletogenesis in the sea urchin embryo
.
Dev. Biol.
283
,
384
396
.
Croce
J. C.
,
McClay
D. R.
(
2010
).
Dynamics of Delta/Notch signaling on endomesoderm segregation in the sea urchin embryo
.
Development
137
,
83
91
.
Damle
S.
,
Davidson
E. H.
(
2011
).
Precise cis-regulatory control of spatial and temporal expression of the alx-1 gene in the skeletogenic lineage of s. purpuratus
.
Dev. Biol.
357
,
505
517
.
Davidson
E. H.
(
2010
).
Emerging properties of animal gene regulatory networks
.
Nature
468
,
911
920
.
Donoghue
P. C.
,
Sansom
I. J.
(
2002
).
Origin and early evolution of vertebrate skeletonization
.
Microsc. Res. Tech.
59
,
352
372
.
Duloquin
L.
,
Lhomond
G.
,
Gache
C.
(
2007
).
Localized VEGF signaling from ectoderm to mesenchyme cells controls morphogenesis of the sea urchin embryo skeleton
.
Development
134
,
2293
2302
.
Emily-Fenouil
F.
,
Ghiglione
C.
,
Lhomond
G.
,
Lepage
T.
,
Gache
C.
(
1998
).
GSK3beta/shaggy mediates patterning along the animal-vegetal axis of the sea urchin embryo
.
Development
125
,
2489
2498
.
Ettensohn
C. A.
(
2006
).
The emergence of pattern in embryogenesis: regulation of beta-catenin localization during early sea urchin development
.
Sci. STKE
2006
,
pe48
.
Ettensohn
C. A.
(
2013
).
Encoding anatomy: developmental gene regulatory networks and morphogenesis
.
Genesis
51
,
383
409
.
Ettensohn
C. A.
,
McClay
D. R.
(
1988
).
Cell lineage conversion in the sea urchin embryo
.
Dev. Biol.
125
,
396
409
.
Ettensohn
C. A.
,
Illies
M. R.
,
Oliveri
P.
,
De Jong
D. L.
(
2003
).
Alx1, a member of the Cart1/Alx3/Alx4 subfamily of Paired-class homeodomain proteins, is an essential component of the gene network controlling skeletogenic fate specification in the sea urchin embryo
.
Development
130
,
2917
2928
.
Ettensohn
C. A.
,
Kitazawa
C.
,
Cheers
M. S.
,
Leonard
J. D.
,
Sharma
T.
(
2007
).
Gene regulatory networks and developmental plasticity in the early sea urchin embryo: alternative deployment of the skeletogenic gene regulatory network
.
Development
134
,
3077
3087
.
Fernandez-Serra
M.
,
Consales
C.
,
Livigni
A.
,
Arnone
M. I.
(
2004
).
Role of the ERK-mediated signaling pathway in mesenchyme formation and differentiation in the sea urchin embryo
.
Dev. Biol.
268
,
384
402
.
Flynn
C. J.
,
Sharma
T.
,
Ruffins
S. W.
,
Guerra
S. L.
,
Crowley
J. C.
,
Ettensohn
C. A.
(
2011
).
High-resolution, three-dimensional mapping of gene expression using GeneExpressMap (GEM)
.
Dev. Biol.
357
,
532
540
.
Fuchikami
T.
,
Mitsunaga-Nakatsubo
K.
,
Amemiya
S.
,
Hosomi
T.
,
Watanabe
T.
,
Kurokawa
D.
,
Kataoka
M.
,
Harada
Y.
,
Satoh
N.
,
Kusunoki
S.
, et al. 
. (
2002
).
T-brain homologue (HpTb) is involved in the archenteron induction signals of micromere descendant cells in the sea urchin embryo
.
Development
129
,
5205
5216
.
Fujino
Y.
,
Mitsunaga
K.
,
Yasumasu
I.
(
1987
).
Inhibitory effects of omeprazole, a specific inhibitor of H+, K+-ATPase, on spicule formation in sea urchin embryos and in cultured micromere-derived cells
.
Dev. Growth Differ.
29
,
591
597
.
Gong
Y. U.
,
Killian
C. E.
,
Olson
I. C.
,
Appathurai
N. P.
,
Amasino
A. L.
,
Martin
M. C.
,
Holt
L. J.
,
Wilt
F. H.
,
Gilbert
P. U.
(
2012
).
Phase transitions in biogenic amorphous calcium carbonate
.
Proc. Natl. Acad. Sci. USA
109
,
6088
6093
.
Guss
K. A.
,
Ettensohn
C. A.
(
1997
).
Skeletal morphogenesis in the sea urchin embryo: regulation of primary mesenchyme gene expression and skeletal rod growth by ectoderm-derived cues
.
Development
124
,
1899
1908
.
Gustafson
T.
,
Wolpert
L.
(
1967
).
Cellular movement and contact in sea urchin morphogenesis
.
Biol. Rev. Camb. Philos. Soc.
42
,
442
498
.
Harkey
M. A.
,
Whiteley
A. H.
(
1980
).
Isolation, culture, and differentiation of echinoid primary mesenchyme cells
.
Rouxs Arch. Dev. Biol.
189
,
111
122
.
Hart
M. W.
,
Strathmann
R. R.
(
1994
).
Functional consequences of phenotypic plasticity in echinoid larvae
.
Biol. Bull.
186
,
291
299
.
Hodor
P. G.
,
Ettensohn
C. A.
(
1998
).
The dynamics and regulation of mesenchymal cell fusion in the sea urchin embryo
.
Dev. Biol.
199
,
111
124
.
Hodor
P. G.
,
Ettensohn
C. A.
(
2008
).
Mesenchymal cell fusion in the sea urchin embryo
.
Methods Mol. Biol.
475
,
315
334
.
Hodor
P. G.
,
Illies
M. R.
,
Broadley
S.
,
Ettensohn
C. A.
(
2000
).
Cell-substrate interactions during sea urchin gastrulation: migrating primary mesenchyme cells interact with and align extracellular matrix fibers that contain ECM3, a molecule with NG2-like and multiple calcium-binding domains
.
Dev. Biol.
222
,
181
194
.
Hughes
I.
,
Blasiole
B.
,
Huss
D.
,
Warchol
M. E.
,
Rath
N. P.
,
Hurle
B.
,
Ignatova
E.
,
Dickman
J. D.
,
Thalmann
R.
,
Levenson
R.
, et al. 
. (
2004
).
Otopetrin 1 is required for otolith formation in the zebrafish Danio rerio
.
Dev. Biol.
276
,
391
402
.
Hurle
B.
,
Ignatova
E.
,
Massironi
S. M.
,
Mashimo
T.
,
Rios
X.
,
Thalmann
I.
,
Thalmann
R.
,
Ornitz
D. M.
(
2003
).
Non-syndromic vestibular disorder with otoconial agenesis in tilted/mergulhador mice caused by mutations in otopetrin 1
.
Hum. Mol. Genet.
12
,
777
789
.
Hwang
S. P.
,
Lennarz
W. J.
(
1993
).
Studies on the cellular pathway involved in assembly of the embryonic sea urchin spicule
.
Exp. Cell Res.
205
,
383
387
.
Illies
M. R.
,
Peeler
M. T.
,
Dechtiaruk
A. M.
,
Ettensohn
C. A.
(
2002
).
Identification and developmental expression of new biomineralization proteins in the sea urchin Strongylocentrotus purpuratus
.
Dev. Genes Evol.
212
,
419
431
.
Ingersoll
E. P.
,
Wilt
F. H.
(
1998
).
Matrix metalloproteinase inhibitors disrupt spicule formation by primary mesenchyme cells in the sea urchin embryo
.
Dev. Biol.
196
,
95
106
.
Ishikawa
H. O.
,
Xu
A.
,
Ogura
E.
,
Manning
G.
,
Irvine
K. D.
(
2012
).
The Raine syndrome protein FAM20C is a Golgi kinase that phosphorylates bio-mineralization proteins
.
PLoS ONE
7
,
e42988
.
Jackson
D. J.
,
Macis
L.
,
Reitner
J.
,
Degnan
B. M.
,
Wörheide
G.
(
2007
).
Sponge paleogenomics reveals an ancient role for carbonic anhydrase in skeletogenesis
.
Science
316
,
1893
1895
.
Kim
E.
,
Hyrc
K. L.
,
Speck
J.
,
Lundberg
Y. W.
,
Salles
F. T.
,
Kachar
B.
,
Goldberg
M. P.
,
Warchol
M. E.
,
Ornitz
D. M.
(
2010
).
Regulation of cellular calcium in vestibular supporting cells by otopetrin 1
.
J. Neurophysiol.
104
,
3439
3450
.
Kim
E.
,
Hyrc
K. L.
,
Speck
J.
,
Salles
F. T.
,
Lundberg
Y. W.
,
Goldberg
M. P.
,
Kachar
B.
,
Warchol
M. E.
,
Ornitz
D. M.
(
2011
).
Missense mutations in Otopetrin 1 affect subcellular localization and inhibition of purinergic signaling in vestibular supporting cells
.
Mol. Cell. Neurosci.
46
,
655
661
.
Kiyozumi
D.
,
Takeichi
M.
,
Nakano
I.
,
Sato
Y.
,
Fukuda
T.
,
Sekiguchi
K.
(
2012
).
Basement membrane assembly of the integrin α8β1 ligand nephronectin requires Fraser syndrome-associated proteins
.
J. Cell Biol.
197
,
677
689
.
Knoll
A. H.
(
2003
).
Biomineralization and evolutionary history
.
Reviews in Mineralogy and Geochemistry
54
,
329
356
.
Krane
S. M.
,
Inada
M.
(
2008
).
Matrix metalloproteinases and bone
.
Bone
43
,
7
18
.
Kurokawa
D.
,
Kitajima
T.
,
Mitsunaga-Nakatsubo
K.
,
Amemiya
S.
,
Shimada
H.
,
Akasaka
K.
(
1999
).
HpEts, an ets-related transcription factor implicated in primary mesenchyme cell differentiation in the sea urchin embryo
.
Mech. Dev.
80
,
41
52
.
Langmead
B.
,
Trapnell
C.
,
Pop
M.
,
Salzberg
S. L.
(
2009
).
Ultrafast and memory-efficient alignment of short DNA sequences to the human genome
.
Genome Biol.
10
,
R25
.
Lepage
T.
,
Ghiglione
C.
,
Gache
C.
(
1992
).
Spatial and temporal expression pattern during sea urchin embryogenesis of a gene coding for a protease homologous to the human protein BMP-1 and to the product of the Drosophila dorsal-ventral patterning gene tolloid
.
Development
114
,
147
163
.
Livingston
B. T.
,
Killian
C. E.
,
Wilt
F.
,
Cameron
A.
,
Landrum
M. J.
,
Ermolaeva
O.
,
Sapojnikov
V.
,
Maglott
D. R.
,
Buchanan
A. M.
,
Ettensohn
C. A.
(
2006
).
A genome-wide analysis of biomineralization-related proteins in the sea urchin Strongylocentrotus purpuratus
.
Dev. Biol.
300
,
335
348
.
Logan
C. Y.
,
Miller
J. R.
,
Ferkowicz
M. J.
,
McClay
D. R.
(
1999
).
Nuclear beta-catenin is required to specify vegetal cell fates in the sea urchin embryo
.
Development
126
,
345
357
.
Malinda
K. M.
,
Fisher
G. W.
,
Ettensohn
C. A.
(
1995
).
Four-dimensional microscopic analysis of the filopodial behavior of primary mesenchyme cells during gastrulation in the sea urchin embryo
.
Dev. Biol.
172
,
552
566
.
Mann
K.
,
Wilt
F. H.
,
Poustka
A. J.
(
2010
).
Proteomic analysis of sea urchin (Strongylocentrotus purpuratus) spicule matrix
.
Proteome Sci.
8
,
33
.
Materna
S. C.
,
Nam
J.
,
Davidson
E. H.
(
2010
).
High accuracy, high-resolution prevalence measurement for the majority of locally expressed regulatory genes in early sea urchin development
.
Gene Expr. Patterns
10
,
177
184
.
Miller
J.
,
Fraser
S. E.
,
McClay
D.
(
1995
).
Dynamics of thin filopodia during sea urchin gastrulation
.
Development
121
,
2501
2511
.
Minokawa
T.
,
Rast
J. P.
,
Arenas-Mena
C.
,
Franco
C. B.
,
Davidson
E. H.
(
2004
).
Expression patterns of four different regulatory genes that function during sea urchin development
.
Gene Expr. Patterns
4
,
449
456
.
Mitsunaga
K.
,
Fujino
Y.
,
Yasumasu
I.
(
1986
).
Change in the activity of Cl-,HCO3(-)-ATPase in microsome fraction during early development of the sea urchin, Hemicentrotus pulcherrimus
.
J. Biochem.
100
,
1607
1615
.
Murdock
D. J. E.
,
Donoghue
P. C.
(
2011
).
Evolutionary origins of animal skeletal biomineralization
.
Cells Tissues Organs
194
,
98
102
.
Oliveri
P.
,
Tu
Q.
,
Davidson
E. H.
(
2008
).
Global regulatory logic for specification of an embryonic cell lineage
.
Proc. Natl. Acad. Sci. USA
105
,
5955
5962
.
Pennington
J. T.
,
Strathmann
R. R.
(
1990
).
Consequences of the calcite skeletons of planktonic echinoderm larvae for orientation, swimming, and shape
.
Biol. Bull.
179
,
121
133
.
Peter
I. S.
,
Davidson
E. H.
(
2011
).
Evolution of gene regulatory networks controlling body plan development
.
Cell
144
,
970
985
.
Peterson
R. E.
,
McClay
D. R.
(
2003
).
Primary mesenchyme cell patterning during the early stages following ingression
.
Dev. Biol.
254
,
68
78
.
Poustka
A. J.
,
Kühn
A.
,
Groth
D.
,
Weise
V.
,
Yaguchi
S.
,
Burke
R. D.
,
Herwig
R.
,
Lehrach
H.
,
Panopoulou
G.
(
2007
).
A global view of gene expression in lithium and zinc treated sea urchin embryos: new components of gene regulatory networks
.
Genome Biol.
8
,
R85
.
Rafiq
K.
,
Cheers
M. S.
,
Ettensohn
C. A.
(
2012
).
The genomic regulatory control of skeletal morphogenesis in the sea urchin
.
Development
139
,
579
590
.
Ransick
A.
,
Davidson
E. H.
(
2006
).
cis-regulatory processing of Notch signaling input to the sea urchin glial cells missing gene during mesoderm specification
.
Dev. Biol.
297
,
587
602
.
Roe
J. L.
,
Park
H. R.
,
Strittmatter
W. J.
,
Lennarz
W. J.
(
1989
).
Inhibitors of metalloendoproteases block spiculogenesis in sea urchin primary mesenchyme cells
.
Exp. Cell Res.
181
,
542
550
.
Röttinger
E.
,
Besnardeau
L.
,
Lepage
T.
(
2004
).
A Raf/MEK/ERK signaling pathway is required for development of the sea urchin embryo micromere lineage through phosphorylation of the transcription factor Ets
.
Development
131
,
1075
1087
.
Sharma
T.
,
Ettensohn
C. A.
(
2010
).
Activation of the skeletogenic gene regulatory network in the early sea urchin embryo
.
Development
137
,
1149
1157
.
Smyth
I.
,
Scambler
P.
(
2005
).
The genetics of Fraser syndrome and the blebs mouse mutants
.
Hum. Mol. Genet.
14
Suppl. 2
,
R269
R274
.
Sodergren
E.
,
Weinstock
G. M.
,
Davidson
E. H.
,
Cameron
R. A.
,
Gibbs
R. A.
,
Angerer
R. C.
,
Angerer
L. M.
,
Arnone
M. I.
,
Burgess
D. R.
,
Burke
R. D.
, et al. 
;
Sea Urchin Genome Sequencing Consortium
(
2006
).
The genome of the sea urchin Strongylocentrotus purpuratus
.
Science
314
,
941
952
.
Söllner
C.
,
Schwarz
H.
,
Geisler
R.
,
Nicolson
T.
(
2004
).
Mutated otopetrin 1 affects the genesis of otoliths and the localization of Starmaker in zebrafish
.
Dev. Genes Evol.
214
,
582
590
.
Stathopoulos
A.
,
Levine
M.
(
2005
).
Genomic regulatory networks and animal development
.
Dev. Cell
9
,
449
462
.
Stumpp
M.
,
Hu
M. Y.
,
Melzner
F.
,
Gutowska
M. A.
,
Dorey
N.
,
Himmerkus
N.
,
Holtmann
W. C.
,
Dupont
S. T.
,
Thorndyke
M. C.
,
Bleich
M.
(
2012
).
Acidified seawater impacts sea urchin larvae pH regulatory systems relevant for calcification
.
Proc. Natl. Acad. Sci. USA
109
,
18192
18197
.
Tagliabracci
V. S.
,
Engel
J. L.
,
Wen
J.
,
Wiley
S. E.
,
Worby
C. A.
,
Kinch
L. N.
,
Xiao
J.
,
Grishin
N. V.
,
Dixon
J. E.
(
2012
).
Secreted kinase phosphorylates extracellular proteins that regulate biomineralization
.
Science
336
,
1150
1153
.
Trapnell
C.
,
Roberts
A.
,
Goff
L.
,
Pertea
G.
,
Kim
D.
,
Kelley
D. R.
,
Pimentel
H.
,
Salzberg
S. L.
,
Rinn
J. L.
,
Pachter
L.
(
2012
).
Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks
.
Nat. Protoc.
7
,
562
578
.
Tu
Q.
,
Brown
C. T.
,
Davidson
E. H.
,
Oliveri
P.
(
2006
).
Sea urchin Forkhead gene family: phylogeny and embryonic expression
.
Dev. Biol.
300
,
49
62
.
Tu
Q.
,
Cameron
R. A.
,
Worley
K. C.
,
Gibbs
R. A.
,
Davidson
E. H.
(
2012
).
Gene structure in the sea urchin Strongylocentrotus purpuratus based on transcriptome analysis
.
Genome Res.
22
,
2079
2087
.
Van Nostrand
E. L.
,
Kim
S. K.
(
2011
).
Seeing elegance in gene regulatory networks of the worm
.
Curr. Opin. Genet. Dev.
21
,
776
786
.
Verkman
A. S.
(
2011
).
Aquaporins at a glance
.
J. Cell Sci.
124
,
2107
2112
.
Weitzel
H. E.
,
Illies
M. R.
,
Byrum
C. A.
,
Xu
R.
,
Wikramanayake
A. H.
,
Ettensohn
C. A.
(
2004
).
Differential stability of beta-catenin along the animal-vegetal axis of the sea urchin embryo mediated by dishevelled
.
Development
131
,
2947
2956
.
Wessel
G. M.
,
Etkin
M.
,
Benson
S.
(
1991
).
Primary mesenchyme cells of the sea urchin embryo require an autonomously produced, nonfibrillar collagen for spiculogenesis
.
Dev. Biol.
148
,
261
272
.
Westbroek
P.
,
Marin
F.
(
1998
).
A marriage of bone and nacre
.
Nature
392
,
861
862
.
Wikramanayake
A. H.
,
Huang
L.
,
Klein
W. H.
(
1998
).
beta-Catenin is essential for patterning the maternally specified animal-vegetal axis in the sea urchin embryo
.
Proc. Natl. Acad. Sci. USA
95
,
9343
9348
.
Wilt
F. H.
,
Ettensohn
C. A.
(
2007
).
The morphogenesis and biomineralization of the sea urchin larval skeleton
. In
Handbook of Biomineralization: Biological Aspects and Structure Formation
(ed.
Bauerlein
E.
), pp.
183
210
.
Weinheim
:
Wiley-VCH Press
.
Wu
S. Y.
,
McClay
D. R.
(
2007
).
The Snail repressor is required for PMC ingression in the sea urchin embryo
.
Development
134
,
1061
1070
.
Wunderlich
Z.
,
DePace
A. H.
(
2011
).
Modeling transcriptional networks in Drosophila development at multiple scales
.
Curr. Opin. Genet. Dev.
21
,
711
718
.
Wuthier
R. E.
,
Lipscomb
G. F.
(
2011
).
Matrix vesicles: structure, composition, formation and function in calcification
.
Front. Biosci. (Landmark Ed.).
16
,
2812
2902
.
Yajima
M.
,
Umeda
R.
,
Fuchikami
T.
,
Kataoka
M.
,
Sakamoto
N.
,
Yamamoto
T.
,
Akasaka
K.
(
2010
).
Implication of HpEts in gene regulatory networks responsible for specification of sea urchin skeletogenic primary mesenchyme cells
.
Zoolog. Sci.
27
,
638
646
.
Yasumasu
I.
,
Mitsunaga
K.
,
Fujino
Y.
(
1985
).
Mechanism for electrosilent Ca2+ transport to cause calcification of spicules in sea urchin embryos
.
Exp. Cell Res.
159
,
80
90
.

Competing interests

The authors declare no competing financial interests.

Supplementary information