|Home | About | Journals | Submit | Contact Us | Français|
The current gene regulatory network (GRN) for the sea urchin embryo pertains to pregastrular specification functions in the endomesodermal territories. Here we extend gene regulatory network analysis to the adjacent oral and aboral ectoderm territories over the same period. A large fraction of the regulatory genes predicted by the sea urchin genome project and shown in ancillary studies to be expressed in either oral or aboral ectoderm by 24h are included, though universally expressed and pan-ectodermal regulatory genes are in general not. The loci of expression of these genes have been determined by whole mount in situ hybridization. We have carried out a global perturbation analysis in which expression of each gene was interrupted by introduction of morpholino antisense oligonucleotide, and the effects on all other genes were measured quantitatively, both by QPCR and by a new instrumental technology (NanoString Technologies nCounter Analysis System). At its current stage the network model, built in BioTapestry, includes 22 genes encoding transcription factors, 4 genes encoding known signaling ligands, and 3 genes that are yet unknown but are predicted to perform specific roles. Evidence emerged from the analysis pointing to distinctive subcircuit features observed earlier in other parts of the GRN, including a double negative transcriptional regulatory gate, and dynamic state lockdowns by feedback interactions. While much of the regulatory apparatus is downstream of Nodal signaling, as expected from previous observations, there are also cohorts of independently activated oral and aboral ectoderm regulatory genes, and we predict yet unidentified signaling interactions between oral and aboral territories.
Gene regulatory networks (GRNs) are testable, predictive models which can provide comprehensive explanations of why developmental functions occur as they do, in terms of the genomic regulatory code (Davidson, 2006; Oliveri et al., 2008; Smith and Davidson, 2008). Portions of the current sea urchin embryo GRN for pregastrular specification of the endomesoderm are approaching a state of relative maturity, particularly that portion referring to the skeletogenic micromere lineage (Oliveri et al., 2008). A major objective is to extend GRN analysis to the other territories of the embryo. For the pregastrular embryo these are the oral and aboral ectoderm, the boundary region surrounding the oral ectoderm which later becomes the ciliated band, and the neurogenic apical plate (see Fig. 1). Here we present an initial GRN model for specification of the oral and aboral ectoderm and ciliated band territories. The model is based on experimental determination of regulatory gene interactions, as revealed by large scale measurement of the effects of perturbation of regulatory gene expression. The large majority of regulatory genes expressed specifically in either oral or aboral ectoderm territories up to late mesenchyme blastula stage is included in this model, but pan ectodermal genes are not, as our focus is on understanding the specification mechanism by which the various ectodermal territories are distinguished from one another.
There are many valuable returns to be expected as GRN analysis extends to further domains of the sea urchin embryo. The more global the model, the fewer inputs into given domains remain mechanistically unexplained. Or, in other words, the more global the model, the closer we approach a complete solution such that regulatory inputs into any given gene originate in outputs of another gene in the model. In addition, previously unknown inter-territorial signaling interactions may be revealed by GRN analysis that extends across adjacent territories, as we illustrate below.
Specification of the ectodermal territories of this embryo has been far less well explored by experimental embryology than has specification of the endomesodermal territories. Unlike the autonomously specified endomesodermal polarity of the sea urchin egg, the oral-aboral polarity is established only gradually. Lineage labeling experiments showed that by second cleavage some oral-aboral polarity is already evident (Cameron et al., 1990), but this remains easily reversible for some time. Oral vs. aboral regulatory states are not irreversibly committed until late blastula stage. Thus it is only after this stage that the oralizing respecification effects of NiCl2 treatment are no longer seen (Hardin et al., 1992). Long before late blastula, however, specific aboral differentiation genes are transcriptionally activated in the aboral ectoderm, e.g., the spec Ca2+ binding genes (Tomlinson and Klein, 1990), and the cyIIIa cytoskeletal actin gene (Lee et al., 1986). This, together with the lineage labeling data, implied an early asymmetry in the future aboral and oral territories that affects transcriptional activity. Cis-regulatory studies of the cyIIIa gene began in the 1980’s (Hough-Evans et al., 1988; Franks et al., 1990; Kirchhamer and Davidson, 1996; Coffman et al., 1997), and were completed only recently (Brown, 2007). In the early phases of cyIIIa cis-regulatory analysis it became clear that the polarity of expression of the cyIIIa gene is determined by the action of a Zn Finger transcriptional repressor (P3A2) which is maternally encoded (Cutting et al., 1990), but is post-translationally activated only in the oral ectoderm. For what follows, experiments suggesting that the activity of P3A2 is redox sensitive (Coffman and Davidson, 2001) turned out to be prescient.
Three subsequent breakthoughs provided causal linkages between the earliest events in embryogenesis and the initial inputs into the ectodermal GRNs that we derive below. First, Coffman et al. (2004) went on to show that what is polarized in the egg even before cleavage is mitochondrial concentration; that the blastomere inheriting the greatest mitochondrial concentration becomes the polar oral ectoderm founder cell (confirming a redox staining observation of Czihak, 1963); and that forced realignment of the redox gradient results in the predictable respecification of oral-aboral fates (Coffman et al., 2004). Then Duboc et al. (2004) showed that activation of the nodal gene in the oral territory as early as 6th–7th cleavage is required for specification of that territory. They also identified several genes downstream of nodal expression, viz. the bmp2/4 and lefty signaling ligand genes, the bra and gsc transcription factor genes in the oral ectoderm, and the tbx2/3 (Gross et al., 2003) regulatory gene in the aboral ectoderm. Activation of the nodal gene remains the earliest oral ectoderm-specific transcriptional function so far identified. The circle was closed when Nam et al. (2007) showed in a cis-regulatory study of the nodal gene that its initial spatial activation depends on target sites for bzip transcription factors, the activity of which is redox sensitive (confirmed independently by Range et al., 2007). This probably accounts for the activation of nodal only on the future oral side of the embryo. We see here one of the genetic targets by which the redox gradient initiated by mitochondrial concentration asymmetry is transduced into a transcriptional input. The subject of this paper is the genomic regulatory apparatus which governs what happens after the initial oral-aboral transcriptional polarity is set up.
There are specific differences between our experimental approach in solving this GRN and the methods available to us earlier (Davidson et al., 2002a, b; Oliveri and Davidson, 2004). For one thing, we began this analysis with knowledge of the temporal expression profiles of all predicted regulatory genes in the S. purpuratus genome (Howard et al., 2006a, b; Materna et al., 2006; Rizzo et al., 2006; Tu et al., 2006), and whole mount in situ hybridization (WMISH) data on those regulatory genes that are expressed at significant levels up to mid-gastrula stage. This allowed us to begin with a list of most of the possible regulatory players, i.e., most (though not all) regulatory genes expressed in either oral or aboral ectoderm by mid-gastrula. A great deal was left to be determined in regard to the exact spatial patterns of expression of these genes, but at least we knew a priori the identities of the great majority of regulatory genes that would have to be included in a GRN model of oral and aboral ectoderm specification. Secondly, in addition to QPCR, we utilized in this work a new technology for assessing the results of a large fraction of our gene specific perturbation experiments, called the “nCounter Analysis System” (Geiss et al., 2008). In this method the transcripts of interest in a sample are identified and counted automatically by a confocal reader, according to tags hybridized to them bearing diverse fluorescent codes. These differences in approach contributed to a greatly accelerated pace of GRN analysis.
cDNAs of the specifically expressed ectoderm regulatory genes were obtained from RT-PCR or 40 h cDNA library screening, and the sequences were deposited in GenBank (Table 1). Morpholino-substituted antisense oligonucleotides (MASOs) specific to the ectoderm regulatory genes were from Gene Tools (Philomath, OR), the sequences of which are shown in Supplementary Information, Table 1. Of the 32 MASOs used about half had been shown to be functional in earlier work, and the effectiveness of the remainder was confirmed by coinjecting the MASO and mRNA with the MASO target sequence fused in-frame to GFP as described (Hinman et al., 2003). For perturbation analyses, eggs were injected with 300 μM MASO (except Nodal MASO at 50 μM; BMP2/4, chordin, and lefty MASOs at 100 μM). RNA from uninjected control and MASO injected embryos were isolated by RNeasy Micro kit (Qiagen). RNA samples were subjected to nCounter Analysis System or reverse transcribed by iScript cDNA synthesis kit (Bio-Rad) for multiplexed quantitative PCR (QPCR). To monitor the quantitative effects of each perturbation, data were normalized to the amount of ubiquitin mRNA as described (Davidson et al., 2002a, b). QPCR primers used in this study are also listed in Table S1.
Details of CodeSet construction are presented elsewhere (Geiss et al., 2008). In brief, two sequence-specific molecules for each gene of interest were constructed. The capture probe consists of a 35 to 50 base sequence complementary to a particular target mRNA plus a short common sequence coupled to a biotin affinity tag. The reporter probe contains a second 35 to 50 base sequence complementary to the same target mRNA coupled to a color-coded molecular tag. The tag is a single-stranded DNA molecule (in this work, linearized single-stranded M13 DNA) annealed to a series of fluorescently-labeled, complementary RNA segments. The linear order of these differently-colored RNA segments creates a unique code for each gene of interest and provides the detection signal. The expression level for each gene is determined by counting the number of times its corresponding code is counted. The nCounter CodeSet for this study contained probe pairs for 87 test and control genes. Fifty-five probe pairs were specific for S. purpuratus genes, seven were specific for Homo sapiens genes (used as negative controls) and 25 corresponded to various nCounter system controls including a standard curve. Each sample was hybridized in triplicate with approximately 100 ng of total RNA in each reaction. The final concentrations of the hybridization reagents were as follows: 200 pM each capture probe, 40 pM each reporter probe, 5X SSPE (pH 7.5), and 0.1% Tween-20. The final total concentrations of all reporter and capture probes in the multiplexed reactions were 17.5 nM and 3.5 nM, respectively. In addition, in vitro transcribed RNA targets for the nCounter spike-in positive controls were added at final concentrations ranging from 100 fM to 0.1 fM. No target RNA was added for the negative controls. Reagents were mixed and incubated at 65°C in a DNA Engine thermocycler (BioRad) with a heated lid for at least 16 h.
In order to remove unhybridized excess reporter and capture probes before imaging, hybridization reactions were purified by affinity purification. Post-hybridization steps were carried out using the nCounter PrepStation liquid handling robot customized to automate the process end-to-end. Briefly, hybridization reactions were sequentially purified by means of magnetic beads coupled to oligonucleotides complementary to the sequence tags present at both the 3′ and 5′ ends of the hybridized tripartite complex (containing reporter and capture probes bound to the target mRNA molecule). Complexes were initially purified from the 3′ end of the capture probe, washed in low salt buffer, and eluted off the beads at elevated temperature. This step removes excess reporter probes and mRNAs not bound to a capture probe. The eluate is then bound to a second set of beads complimentary to the 5′ end of the reporter molecule in order to remove excess capture probes. After washing and elution, the samples were prepared for binding and imaging.
Purified samples were loaded into the NanoString sample cartridge using the nCoutner PrepStation. The cartridge contains 30 μm deep microfluidic channels. The samples were passed though the channels by hydrostatic pressure and bound to the streptavidin-coated surface via the biotinylated 3′ end of the capture probe. After capture, the surface was washed once, buffer was added to each well, and the bound molecules were stretched and aligned by applying 160 V/cm for 1 min along the fluidic channel. The 5′ ends of the elongated reporters were then immobilized to the surface by addition of a biotinylated oligonucleotide complementary to the 5′ end of each reporter probe. After immobilization, the electrophoresis buffer was removed and replaced with a custom formulation of the anti-photobleaching reagent SlowFade (Invitrogen) for imaging. Slides were imaged on the nCounter Digital Analyzer custom scanner that takes images of each field of view at 60X-magnification in 4 different excitation wavelengths (480, 545, 580 and 622). Images were automatically processed with custom image-processing software. The expression level of a gene is measured by counting the number of times the code for that gene is detected. The details of post-hybridization processing and imaging are described elsewhere (Geiss et al., 2008).
To account for slight differences in hybridization and purification efficiency as well as mRNA content, the raw data for S. purpuratus genes were normalized to the ubiquitin gene transcripts present in all reactions. Each set of samples was normalized to the ubiquitin mRNA levels in their corresponding uninjected control sample. To determine if the counts for each gene were statistically above background, a Student’s T-test against the seven human negatives data was performed. A gene was considered to be above background if the average counts for the S. purpuratus gene was greater than the average counts for the seven H. sapiens negative control genes, and the Student’s T-test P-value was less than 0.05.
Perturbation data were obtained from at least thee independent batches of embryos. In order to compare results from two different measuring methods, data from the nCounter Analysis System were converted to cycles (Ct) assuming the ratio of transcript accumulation per PCR cycle is 1.9. Changes in the prevalence of the transcripts were considered significant if the effect between uninjected control and MASO injected embryos is more than 1.6 Ct (~thee fold). Epistatic interactions were evaluated if the results of two of thee independent experiments showed significant changes. Occasionally MASOs display toxic effects, as indicated when the expression of most genes is depressed significantly upon MASO injection (e.g., the Hes and E2F3 MASOs), and these data are excluded from the analysis. In addition, when the effects of MASOs are observed at early time points but no effects are seen in observations a few hs later, i.e., in early but not late stage embryos, the putative input is excluded (an exception for application of this criterion is nk1 MASO). Data were also ignored if the transcript is too rare for reliable QPCR ratios (>32 cycles) at the time analyzed. The ectoderm GRN model was then formulated from the remaining more robust results.
RNA probes were prepared either from PCR products or cDNA clones. For PCR based probes Sp6 and T7 tailed primers were selected based on their match pattern against an EST database, but in addition the secondary structures, melting temperatures, and possible hairpins were considered. PCR was performed on reverse transcribed cDNA after RNA extraction at certain time points. The PCR products were purified by gel electrophoresis. Gels were stained with SYBR Gold nucleic acid gel stain (Invitrogen Molecular Probes). Bands were visualized, cut, extracted by electrodialysis, and concentrated using a DNA clean up kit (Zymo Research). Resulting probes are of high purity, so that reverse transcription could be done on lower amounts of sample than usual. Digoxigenin (DIG)-labeled probes were purified using an RNeasy Protect Minikit (Qiagen) with on column DNAse digestion. For double WMISH, gsc or lefty were used as a spatial reference for oral expression. Both Spgsc and Splefty were PCR cloned into pGEMT-EZ vector and used as template to prepare dinitrophenol (DNP)-labeled RNA probes as described above.. Double WMISH was performed using a standard method as described by Revillai-Domingo et al. (2006), except that hybridization temperatures were variably adjusted to 59 or 62°° instead of 65°. Probe concentrations were adjusted to 100ng/ml – 350ng/ml for overnight hybridization in order to obtain signal from both DIG and DNP labeled probes. The anti DIG antibody was used for the first antibody incubation. The first staining reaction was carried out using 4.5 μl/ml NBT (N-6876:Sigma-Aldrich), 3.5μl/ml BCIP (Sigma Aldrich). For the second antibody, we incubate with anti-DNP conjugate (Mirrus Gene Transfer). In the second staining reaction NBT was replaced with 4μl/ml INT(Sigma Aldrich) (orange)/4μl/ml BCIP. To increase color contrast the DNP-INT staining step (yellow-orange) was carried out first followed by the DIG-NBT staining. However, the INT/BCIP reaction product is not very stable with respect to the following washing steps, sometimes necessitating adjustment of the hybridization probe concentration.
The GRN model we present here is constructed mainly of regulatory genes and their interactions, plus a few signaling genes and the immediate early regulatory response genes that transduce their signal inputs at the transcriptional level. In the course of the S. purpuratus genome project many regulatory genes not previously studied in sea urchin embryos were observed to be expressed in ectodermal domains. A compilation of almost all regulatory genes specifically expressed in one or another part of the developing ectoderm, including a summary of their domains of expression by 24 h, is given in Table 1. This essentially provides the parts list for the model, though it should be remembered that additional regulatory genes are deployed in the ectodermal domains as time progresses and a network analysis that would extend though the onset of gastrulation at 30 h, and beyond, would include many more genes and their interactive linkages. However, in this work we have only considered genes activated by 24 h, as our first objective was to encompass the processes of ectoderm specification in a GRN that would extend as far into development as does the endomesoderm GRN (24–30 h). Table 1 includes all the genes incorporated in the current model. The model does not include two of the genes of Table 1, e2f3 and hlf, as we could not find any consistent inputs into them from other genes in the model within the period considered. A few additional regulatory genes so far not decisively assigned to either ectodermal domain are not included in Table 1, though they are expressed asymmetrically in the ectoderm by late mesenchyme blastula stage, and these will possibly also need to be included in the next draft of these GRNs. Nonetheless all evidence, genomic as well as earlier experimental evidence (cf. earlier references in Table 1), indicates that the large majority of zygotically expressed regulatory genes that could be involved in the pregastrular specification of oral and aboral ectoderm are linked into the GRN we present in this paper.
The perturbation assays on which the structure of the model is based were organized as matrix analyses, and the model at present includes 29 genes, of which thee are logically implied but are not yet identified. Gene expression was perturbed by MASO injections into the egg, and the effects on transcript levels of every other gene in the system analyzed at different time points by QPCR and/or the nCounter Analysis System (see Materials and methods). Quantitative data for all perturbations for all genes included are presented in Fig. S1, in which replicate perturbation data from >100 individual experiments are presented. Most of the temporal expression profiles for these ectoderm genes were published earlier (Howard et al., 2006a, b; Materna et al., 2006; Tu et al., 2006). In some cases it was necessary to further ascertain whether a gene is expressed in oral or aboral ectoderm by use of double WMISH, which included the oral ectoderm specific gsc or lefty probes for spatial reference. Some representative examples and additional WMISH of less well known ectoderm genes are shown in Fig. 2.
The ectoderm GRN model was based on the perturbation results, taken together with knowledge of the temporal and spatial expression profiles of the regulatory genes that are expressed in oral or aboral ectoderm in the pregastrular embryo. As discussed at length elsewhere, this information, when organized within the framework of the developmental biology of the system, suffices to specify the structure of the GRN model (Davidson, 2006; Oliveri et al., 2008). The dataset is first culled for sources of error, such as generally toxic MASO’s which affect all or almost all genes, or other systematic evidence of unreliable results such as quantitative non-reproducibility (see Materials and methods; all admitted evidence is supported by multiple repeats on different batches of embryos which gave similar results, either by QPCR or the nCounter or both). The surviving data for all runs are shown in Fig. S1. Six regulatory genes that respond to certain of the MASO treatments, and are included in Fig. S1 and Table S1, were for one reason or another not included in this version of the ectoderm GRN model. These are atbf1, arnt, foxk, and thee zinc finger genes. They were excluded because their spatial expression patterns had not been clarified and in most cases MASOs targeting them were not studied. For all the other genes than these we believe we know both the locus of expression and their epistatic relationships, both upstream and downstream.
The data in Fig. S1 plus time and place of expression were then utilized to generate the topology of the GRN model. The perturbation results indicate possible epistatic relations among the genes, but because of the pleiotropic functions of most regulatory genes, a great many of these apparent relations are expected to be indirect. As discussed earlier (Longabaugh et al., 2005; Oliveri et al., 2008; Davidson, 2006), logical criteria can be brought to bear in an effort to judge whether a linkage between interacting regulatory genes is likely to be direct or indirect. For example, a common indicator of indirect effects is non-coincidence of gene expression either in time or space;.if it the case that a regulatory gene directly activates another, they must be expressed in the same cell(s) at about the same times. In some cases we observed genes strongly down-regulated by MASO’s to transcription factors expressed in distinctly different domains. Such observations require the existence of a signaling relationship in which the input gene is upstream of expression of a ligand gene, and the responding gene in the second domain is downstream of the signal transduction apparatus in the receiving cells. The fine scale time courses of gene expression provide an additional aid in organizing the GRN, since genes that are expressed downstream of other genes are activated after several hs delay, as observed in the endomesoderm GRN (Bolouri and Davidson, 2003; Oliveri et al., 2008; Ben-Tabou de-Leon and Davidson, 2008). These and other such “rules” were applied as relevant to each prospective linkage indicated by the data in an effort to determine the likely direct linkages, i.e, to provide specific predictions of direct linkages that can be verified or rejected by cis-regulatory measurements. In the endomesoderm GRN we found that about 90% of similarly predicted direct linkages investigated at the cis-regulatory level are in fact verified, though additional ones are also found (see for example, Yuh et al., 2004; Amore et al., 2006; Lee et al. 2007; Minokawa et al., 2005; Ransick et al., 2006; Smith et al., 2007). To aid the reader in following the train of logic leading to the proposal of each individual linkage in this work, we have constructed Supplemental Information Tables S2 and S3. Here are listed explicitly for every possible interaction the time and place of expression; the linkage path we deduced to be responsible for the perturbation results; and the specific rationale we utilized. Table S2 concerns regulatory genes, and Table S3 signaling factors. As a general parsimonious guideline we impose the criterion that if a linkage could be indirect we assume that it is. The drawback is that this strategy misses some feed forward linkages that have later to be filled in by cis-regulatory analysis or other means, but this does not materially alter the topology of the GRN. For the present dataset we have indicated in Tables S2 and S3 the pathways of indirect linkage as well as the direct ones. In choosing among a set of possible linkages the most likely direct one, we have found that the direct linkages generally produce stronger results when the input is perturbed, while the indirect linkages produce weaker effects.
To compare the two analysis methods used in this study, RNAs from the same batch of embryos were subjected to QPCR and nCounter Analysis System. A sample of the results is reproduced in Fig. 3, and all 106 experiments in which head to head comparisons were included are displayed in Fig. S1. The two methods agree with each other in almost all cases. The only disparities are genes expressed at very low levels at the time analyzed. For example, when using QPCR analysis, it takes more than 30 cycles for nk1 and lhx2.9 to reach a threshold. In the nCounter Analysis System, the counts for nk1 and lhx2.9 are insignificantly different from the background counts (P-value > 0.05). Considering the amount of material needed for the two methods, each QPCR reaction requires cDNA reverse transcribed from ~3 ng of total RNA (the content of one embryo). To measure one gene (four replicates), 12 ng of total RNA are required. For the nCounter Analysis System, 100 ng of total RNA is needed for each hybridization (300 ng for triplicates). In these studies the nCounter Analysis Codeset (the tagged custom made probes identifying each gene) represented fifty-five genes, expression of which was measured simultaneously in the sample. In summary, the nCounter Analysis System is as accurate as QPCR and uses less material.
We believe that the perturbation data and other observations summarized in this paper logically require the GRN structure shown in Fig. 4. On those few linkages where there is prior evidence our conclusions largely agree with this evidence (a single exception is considered below). Because this model is not complete, there remain many unanswered questions, and the several additional asymmetrically expressed regulatory genes not yet included in the model remain to be incorporated. Most importantly, there remains to authenticate this model at the cis-regulatory level by demonstrating the presence and function of the predicted regulatory inputs at the important nodes of the model (Levine and Davidson, 2005). However, it would be impossible to even formulate the next steps of GRN construction absent the provisional GRN shown in Fig. 4. Experience with the endomesoderm GRN shows that the basic structure of the GRN is unlikely to change radically as further linkages and revisions are incorporated, and the model is checked at the cis-regulatory level. Fig. 4 already reveals many interesting, and to some extent surprising, features of oral and aboral ectoderm specification that are unlikely to disappear on further evidence.
The GRN in Fig. 4 allows us to address several issues with greater clarity than would otherwise be possible. These include the exact role of Nodal signaling in ectoderm specification; the additional functions of the initial oral-aboral anisotropy used for nodal activation in the oral ectoderm; and recognition of regulatory circuit design themes which appear in the ectoderm GRNs.
Several unexpected features of the oral ectoderm specification system quickly emerged from the perturbation analysis. The first of these concerned the presence of some genes activated in the prospective oral ectoderm only if nodal is expressed, but others that are activated independently of nodal expression. It had been clearly shown by Lepage and colleagues, as reviewed above, that nodal expression is the initial transcriptional event required for oral ectoderm specification. Indeed, as would be predicted from this, transcription of most of the oral ectoderm regulatory apparatus revealed by this work is strongly depressed when nodal expression is blocked (Fig. S1). However, we discovered that the network of interacting regulatory genes expressed specifically in the oral ectoderm also includes a number of regionally expressed genes that are activated during oral ectoderm specification exactly the same with or without nodal expression. These are the early oral ectoderm and later ciliated band genes hnf6 and otxβ1/2; the lim1 gene, which is expressed along the lower margin of both oral and aboral ectoderm; and the stomodeal gene, foxa. Therefore only a part of the oral ectoderm specification system, though an essential part, is downstream of the Nodal signaling apparatus.
Secondly, it is apparent both from the primary WMISH evidence showing where the various regulatory genes in this analysis are expressed, and from the partitioning of the epistatic relations we observe, that the oral ectoderm actually consists of several spatial regulatory domains which begin to resolve only after the initial specification of the oral territory. A major separate domain is the ciliated band region. From early gastrula stage onwards, genes initially expressed thoughout the oral ectoderm resolve to this trapezoidal band, i.e, the hnf6 (Otim et al., 2004), otxβ1/2 (Yuh et al., 2002), and foxg (Tu et al., 2006) genes (Fig. 1). There remains the main cuboidal epithelium of the oral ectoderm “face” which continues to express the gsc, dri, and hes genes. Expression of foxa and bra genes marks the stomodeal subdomain of the oral face. Finally, there is a separate domain along the interface with the vegetal endomesoderm that extends all across both oral and aboral ectoderm, as indicated by expression of the lim1 gene (Kawasaki et al. 1999), and in the oral part of this, the nk1 (Minokawa et al., 2004), ecr, and foxj1 genes are also expressed. Later, not considered here, are the tightly confined gene expression patterns (e.g., otp, vegf, and fgf) that underlie the positioning of the skeletal rods (Armstrong and McClay, 1994; Di Bernardo et al., 1999; Cavalieri et al., 2003; Duloquin et al., 2007; Röttinger et al., 2008).
The nodal gene participates in a transcriptional “community effect” (Gurdon, 1988; Davidson, 2006). This is an intra-territorial positive feedback relationship such that all the cells of the territory signal to one another. The consequences are that each cell emits more of the same signal, and the territory-specific regulatory state is furthered as well, in that key regulatory genes are also driven by the same signal transduction system. The canonical regulatory circuitry that controls the community effect in the oral ectoderm is that shown in Fig. 5A. Here we see that the Nodal signal causes recipient cells to transcribe the nodal gene though cis-regulatory feedback response elements for the Smad transcription factor, which is activated by Nodal binding to its receptor. This feedback regulation was demonstrated functionally by cis-regulatory site mutation (Nam et al., 2007; Range et al., 2007). In design, the circuitry is identical to the circuitry which produces a Wnt8 community effect in the endomesoderm of the same embryo (Minokawa et al., 2005), and other examples are to be found in GRNs from other developmental systems (Davidson, 2006). In each case the signal ligand gene is activated by its own specific signal transduction system, and all cells of the territory both receive and transcribe the signal. In the case of Nodal the spread of the ligand outside of the oral ectoderm territory is apparently cancelled by the peripheral presence of the diffusible Nodal antagonist Antivin/Lefty (Duboc et al., 2004, 2008). Expression of the lefty gene is also under control of the transcriptional network downstream of nodal expression, as Duboc et al. (2004) showed, and we independently confirm (Fig. 4; Fig. S1). The community effect is an autoregulatory feedback system which drives the cells of the oral ectoderm into a transcriptional lockstep. This is used in turn to drive expression of the initial nodal-dependant regulatory state specifiers of the oral ectoderm. Thus all cells of the future oral ectoderm come to express these genes, following the interpretation of the primordial redox gradient at the bzip sites of the nodalcis -regulatory system (Nam et al., 2007).
In addition to the nodal gene itself (Nam et al., 2007; Range et al., 2007), the two direct targets of the Nodal signal transduction pathway predicted by the analysis are respectively a gene encoding a transcriptional repressor, and a gene encoding a transcriptional activator. The repressor is goosecoid (gsc). This gene has been for some time known to produce an obligate transcriptional repressor which is essential for oral ectoderm specification (Angerer et al., 2001; Amore et al., 2003). Thus we confirm Duboc et al. (2004) that gsc activation depends on nodal expression. But foxg, the activator, is a newly discovered and pivotal early player in oral ectoderm specification originally uncovered in the survey of all S. purpuratusfox genes carried out by Tu et al. (2006). As Fig. 4 shows, we predict that foxg lies high up in the hierarchical oral ectoderm specification GRN.
To summarize, there appear to be four early genetic components producing the initial zygotic transcriptional regulatory state in the specification of the oral ectoderm: the nodal independent hnf6, and otxβ1/2 genes; the nodal dependant gsc repressor, and the nodal dependant foxg activator. The process of setting up a regulatory state is essentially the positive function of activating a novel set of transcriptional regulatory genes, and this can be done either by directly activating them or by repressing a gene which otherwise represses them: the oral ectoderm GRN uses both.
Activation of the skeletogenic micromere specification GRN has been shown to depend on a double negative regulatory gate (Oliveri et al., 2002; 2003; 2008; Revilla-i-Domingo et al., 2007). In that GRN a gene encoding an initial repressor, pmar1, is activated specifically in the micromeres in response to localized anisotropic cues of maternal origin, and its role is to repress zygotic transcription of a second gene encoding a repressor, the hesC gene. The hesC gene is expressed globally, and it represses the upper tier of regulatory state specifiers everywhere except in the micromere lineage, due to pmar1 expression there. This is the mechanism that accounts for the installation of the specific regulatory state of the skeletogenic micromere lineage. In terms of the predicted circuitry, the GRN component downstream of gsc in Fig. 4 has a very similar architecture (abstracted in Fig. 5B). Here gsc plays the role of pmar1: it is an initial regulatory gene activated in the oral ectoderm which encodes a repressor (Angerer et al., 2001). A predicted unknown gene, called repressor A in Figs. 4 and and5B,5B, the target of the Gsc repressor, plays the role of the hesC gene in the skeletogenic micromere GRN. The definitive set of early regulatory genes in the oral ectoderm, including dri, bra,, and hes all are under repressor A control, as evinced by their sharp decrease in expression when gsc mRNA translation is prevented by MASO treatment. In addition, the same gate is used for feedback onto gsc and foxg (Figs. 4, ,5B).5B). The release by this mechanism of this whole set of oral ectoderm target genes from a state of repression would appear from our perturbation data to be what explains their activation in this domain (Fig. S1). It is furthermore consistent with the observation made by Angerer et al. (2001) that ectopic expression of gsc causes ectopic oralization of the ectoderm, even though gsc encodes an obligate repressor. An active search for the repressor A gene is underway in our lab.
Diverse data indicate that the redox anisotropy across the oral–aboral dimension of the early cleavage egg is utilized to promote differential gene expression by more cis-regulatory systems than that of the nodal gene alone. As reviewed above, an initial clue that this anisotropy exists came from study of the cyIIIa cytoskleletal actin gene, a differentiation gene of the aboral ectoderm. This gene is activated on the aboral side specifically because it is repressed on the oral side by the P3A2 factor (Hough-Evans et al., 1988; Franks et al., 1990; Höög et al., 1991), and P3A2 is modified in direct response to the redox gradient (Coffman and Davidson, 2001). The activators of the cyIIIa gene are known (Kirchhamer and Davidson, 1996; Brown, 2007), and they are either ubiquitous or pan-ectodermal; i.e., neither they nor the P3A2 repressor are downstream of the regulatory system in Fig. 4. This shows that the initial cytoplasmic redox input is used at diverse levels of the regulatory hierarchy, extending all the way down to the cyIIIa differentiation gene, a somewhat surprising finding in light of the discovery that the same system is used to activate a key gene at the top of the hierarchy, nodal. The spec1 gene, another differentiation gene of the aboral ectoderm, could be another example, at least in respect to its initial activation. The cis-regulatory system of this gene has been extensively studied, and like cyIIIa, it is driven by widespread activators and is repressed on the oral side of the ectoderm by a specific factor (OER; Yuh et al., 2001 and references therein). It cannot be excluded that OER, like P3A2, is activated on the oral side by the redox system. Eventually OER must come under the control of the oral ectoderm GRN, since treatments which interfere with the operation of the GRN, such as knock-down of dri (Amore et al., 2003), of hnf6 (Otim et al., 2004), or of gsc (Angerer et al., 2001), all cause expression of spec1 to spread around the whole ectoderm. The existence of parallel regulatory pathways indicated by the nodal-independent oral activation of hnf6, and the nodal-independent aboral activation of tbx2/3, suggest that the same mechanism, rooted in the initial redox anisotropy, could be responsible for the early activation of these genes (cf expression times in Table. 1). This hypothesis is summarized in Fig. 5C. We note here that our QPCR data, which display no effect of nodal MASO on tbx2/3 expression, are on the face of it directly contrary to the evidence of Duboc et al. (2004), who show that in gastrula stage embryos of Paracentrotus lividus, blocking expression of nodal eliminates expression of tbx2/3. However, the QPCR observations which show tbx2/3 transcript levels insignificantly affected by nodal MASO treatment (Fig. S1) refer only to the mid and late blastula stages, i.e., the period of tbx2/3 activation. Thus it is not unlikely that the zygotic transcription of the activator that controls tbx2/3 expression is subject to the GRN by gastrula stage studied by Duboc et al. (2004), while the maternally encoded initial form of this activator responds to the redox gradient of the cleavage stage cytoplasm.
The aboral ectoderm, at least so far as current evidence extends, is a simpler territory than is the oral ectoderm. While some of its lineage descendants contribute to the ciliated band according to dextran labeling data (Cameron et al., 1990), the present analysis shows that the ciliated band is mainly to be thought of as a subdivision of the oral ectoderm territory, since ciliated band genes are initially activated in the oral ectoderm but not the aboral ectoderm. Prior to this work only one regulatory gene had been directly implicated in the specification of the aboral ectoderm. This was tbx2/3, which had been shown not only to be expressed specifically in the aboral ectoderm but also required for its specification (Gross et al., 2003). But we now find that there are eight additional regulatory genes, viz., irxA, nk2.2, lhx2.9, dlx (Howard et al., 2006a), hmx (Martinez and Davidson, 1997), msx (Dobias et al., 1997), hox7 (Angerer et al., 1989), all of which encode homeobox regulators, as well as a Zn finger gene, klf7(Z86) (Materna et al., 2006), expressed specifically in aboral ectoderm and linked causally into the GRN (Fig. 4, Table 1). The aboral ectoderm regulators are expressed uniformly thoughout this domain, so far as we can tell. The only subdivision within it is that mentioned above, the lim1 domain immediately adjacent to the vegetal endomesoderm border.
Dynamic lockdown of regulatory state specification by means of feedback circuitry is now an expected and predictable aspect of developmental GRNs (Davidson, 2006). The aboral ectoderm section of Fig. 4 provides an extreme example: the four homeodomain genes, irxA, lhx2.9, dlx, and hox7, are all predicted to be locked together in feedback relations and irxA feeds back on tbx2/3, the first to be expressed, as well. Perhaps this feature underlies the quality of aboral ectoderm specification long ago recognized, viz. that this domain is inflexibly defined. Its role is to produce but a single flattened epithelial cell type. In this it contrasts greatly with the oral ectoderm, which as reviewed above, continues to evolve and produce new regulatory subdivisions well after its initial specification.
A striking observation substantiated in a number of perturbation experiments was the non-autonomy of aboral as well as oral ectoderm specification. It was known earlier that BMP2/4, an oral ectoderm product, is required for aboral ectoderm development (Angerer et al., 2000; Duboc et al., 2004). We confirmed this at the regulatory level, identifying the target genes of the signal transduction system activated by BMP2/4. These we predict to be hox7 and msx, though because of the feedback relations there could be additional direct targets that we have assumed to be indirect effects, even if, as Fig. S1 shows, they are quite strong. However, we predict in addition a second now unknown signal (Signal X) that goes to an overlapping set of target genes in the aboral ectoderm, viz. hox7, msx, and lhx2.9 (dashed blue lines in Fig. 4). This signal is also of oral ectoderm origin but is distinct from the BMP signal, since it is downstream of nk1 expression, which BMP2/4 is not. As if these were not sufficient, other evidence shows that gsc, the major oral ectoderm regulator, requires in turn an input from a signal of aboral ectoderm origin (Signal Z, dashed green lines of Fig. 4). The perturbation evidence suggests that Signal Z feeds back as well on both the bmp4 and nodal genes. This signal is probably controlled by the aboral ectodermal regulatory genes tbx2/3 and hox7. In short, as we discuss below, the image we have is that specification of the aboral and oral ectoderm territories depends reciprocally on signals from one another. On the face of it this is perhaps not in the least surprising, though there is not much prior evidence on this point.
Scientific reports are all more or less accounts of work in progress, and the GRN of Fig. 4 is far from complete. The single most important challenge to be met for this GRN to approach maturity is to add into it the other genes we are aware of which are expressed specifically in oral or aboral ectoderm. In addition to the six ectodermal regulatory genes not included in this model noted above, many others are activated at or after the conclusion of the period encompassed by this GRN, which is only up to 24–30 h, when the perturbation observations in this paper end. The inputs into these genes are not experimentally apparent from perturbation analysis until later times in development. Furthermore, most of the genes already in the GRN eventually, at later times, are expressed specifically only in particular subregions, i.e., the ciliated band, the epithelial oral face, the vegetal border area, or the stomodeal region. A few of the regulatory relationships which account for this just begin to be evident in the GRN of Fig. 4, but until the GRN is extended forward in time (perhaps beyond 40 h) neither the other specifically expressed genes nor the spatial subdivision apparatus can be incorporated into it. Thus the present GRN is essentially an initial specification GRN, though there may be much less difference between initial and later in the aboral ectoderm.
Second, it should be recognized that the GRN does not include ubiquitously expressed or pan-ectodermal genes, though these certainly provide inputs into cis-regulatory modules functioning specifically in either oral or aboral ectoderm. For example CCAAT binding proteins assist in activation of both spec1 and cyIIIa (Yuh et al., 2001; Kirchhammer and Davidson, 1996); soxB1 provides an input into the nodal gene (Range et al., 2007); and there are numerous other examples. But global inputs cannot be the logical drivers of spatial gene expression, and this GRN, like the endomesodermal GRN, is focused on the genomic regulatory logic that accounts for the regional establishment of oral and aboral territories.
Thirdly, in addition to the mysteries noted above regarding yet unidentified signaling interactions predicted in the GRN analysis, a great number of other unresolved issues are now in focus. Most of these are subject to solution at the cis-regulatory level, which is also the direct roadway to general validation of the architecture of Fig. 4. For example, what are in fact the inputs to the nodal -independent oral and aboral ectoderm genes of the GRN? What is the role of sip1, an early ectodermal regulator? This gene is active early in ectoderm specification, then repressed by gsc; but then reactivated by yet undefined inputs. Is there indeed a direct exclusion relation between gsc and tbx2/3 as an extra lockdown of the oral state, analogous to many other such discovered in the sea urchin and other GRNs (Oliveri and Davidson, 2007)? This is suggested by the wipeout of tbx2/3 expression by either nodal mRNA overexpression (MOE), which would cause ectopic gsc expression (Duboc et al., 2004), or by direct gsc MOE (Bradham and McClay, 2006).
But not one of these questions could even be raised were it not for the structure given to the perturbation results by the GRN model, incomplete as it is. Perhaps most importantly, the model provides a suite of testable predictions for the key cis-regulatory modules of the system, e.g., those directing expression of tbx2/3, irxA, hox7, hnf6, dri, gsc, foxg, etc. These predictive targets will vastly increase the efficiency of cis-regulatory analysis, and render it amenable to computational assist.
Supplementary Fig. 1. Perturbation data analysis. Embryos were injected with different MASOs and collected at two different time points after fertilization. The name of MASO and the collection time are shown on top of each panel. RNA samples from different batches of embryos were quantitatively measured by QPCR or the nCounter Analysis System (labeled as NS) for the expression of ectoderm genes (names shown under X-axis). The Y-axis indicates the differences of gene expression (in Ct units) to attain theshold in uninjected control relative to MASO injected embryos. The results regarded as sufficiently robust by the criteria stated in text and hence used to build the ectoderm GRN model are indicated (*). The MASO effects thought to be due to non-specific toxicity, and therefore omitted from the analysis, are labeled with dots (●).
Research was supported by NIH grants HD-37105 and GM-61005 and by the Caltech Beckman Institute.
Publisher's Disclaimer: This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.