Genome-wide analysis of the skeletogenic gene regulatory network of sea urchins Free

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.

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; R²=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.

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 log₂-fold difference (supplementary material Table S1), yielded robust WMISH patterns that were restricted to PMCs at the mesenchyme blastula stage.

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).

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).

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.

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

The deposition of CaCO₃ by PMCs is associated with the uptake of Ca²⁺ and HCO ^-₃ 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.

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).

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

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).

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.

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.

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).

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.

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.

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).

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).

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 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.

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).

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 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 Ca²⁺ 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 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.

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.

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.

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).

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).

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).

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).

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.