A set of simple cell processes is sufficient to model spiral cleavage

Most metazoans start their development via a series of fast cell divisions that partition the zygote into a set of blastomeres. During this cleavage process, a specific spatial cell arrangement, referred to here as the ‘cleavage pattern’, arises in each species. There are several types of cleavage patterns in metazoa (Gilbert and Raunio, 1997). Usually different animal taxa exhibit different cleavage patterns. However, there are several phyla that share a common cleavage pattern. Spiral cleavage is the most abundant cleavage type at the phylum level. It is found in mollusks, annelids and nemerteans. Other lophotrochozoan phyla (platyhelminthes, rotifers, brachiopods, phoronids, gastrotrichs, and bryozoans) also exhibit spiral cleavage in at least some of their species (Hejnol, 2010). The ensemble of phyla with spiralian cleavage has been suggested to form a monophyletic group (Nielsen, 1994; Laumer et al., 2015): the Spiralia. Despite having a very similar cleavage, these phyla have very different adult morphologies.

As in other types of cleavage, typical spiralian cleavage begins with two successive nearly meridional cell divisions that give rise to four large cells (termed A, B, C and D macromeres) that lie in a plane perpendicular to the animal-vegetal (A-V) axis of the egg. The four macromeres then divide usually unequally, along an obliquely equatorial plane, giving rise to four, usually smaller, cells called micromeres. In contrast to the radial cleavage found in many other metazoans, in spiral cleavage each micromere is not placed directly above its sister macromere, but is displaced towards the right (or towards the left, depending on the organism) with respect to its sister macromere. The third cell division, thus, proceeds at an oblique angle relative to the A-V axis. In this process, all micromeres are displaced in the same direction (either all to the right or all to the left with respect to their underlying sister macromeres). If the third cell division produces a micromere to the right, the spiral pattern is classified as ‘dextral’. In this case, the next cell division will proceed to the left, and the ensuing ones will follow in a right-left alternation (the reverse alternation occurs in sinistral spiral cleavage, i.e. if the third division produces a micromere to the left). Owing to this alternation, the four new micromeres appearing between the 4- and 8-cell stages seem to twist clockwise (in the dextral pattern) or counterclockwise (in the sinistral pattern) when viewed from the animal pole (with respect to the four underlying macromeres) (Gilbert and Raunio, 1997; Henry, 2014). Compared with those of other groups, spiralian embryos tend to undergo fewer cell divisions before gastrulation, making it easier to follow the fate of their blastomeres (Gilbert and Raunio, 1997). When cell fates are compared between spiralians, it is often found that the same adult or larval organs in different species arise from the same blastomeres [defined by lineage and relative position in the blastula (Nielsen, 1994; Lyons et al., 2012)]. The first four cell division rounds are synchronous but this synchrony is gradually lost over developmental time and the similarity between groups becomes less obvious (Freeman and Lundelius, 1982; Merkel et al., 2012). In some species with a very derived spiralian cleavage, this synchrony can be lost much earlier (Nielsen, 1994).

There are some studies focusing on the signalling events that, when taking place within a specific spatial blastomere arrangement, lay out the cell fate of each blastomere (Freeman and Lundelius, 1982; Lambert and Nagy, 2003; Kuroda et al., 2009; Grande and Patel, 2008). However, much less is known about how the specific spatial arrangement of blastomeres in spiral cleavage is attained [although there is some work on the early morphogenesis of some invertebrate non-spiralian models, most notably ascidians (Munro et al., 2006)]. A number of developmental processes have been hypothesized to explain spiral cleavage. Some are roughly understood as developmental rules by which the direction of the cell division plane is determined during cleavage (Freeman and Lundelius, 1982), whereas others are processes of mechanical cell interaction that lead to cell displacement during cleavage (Meshcheryakov and Beloussov, 1975; Wandelt and Nagy, 2004; Henley, 2012). Here, we refer to these hypotheses as ‘developmental rules of division plane specification’, or, simply, ‘rules’. We use the term ‘non-directional’ or, simply, ‘cell processes’ for the hypotheses that are not related to the direction of the division plane. Our aim here is to assess which of these previously proposed rules and cell processes, either alone or in combination, are capable of producing the spiral pattern when implemented in a realistic biomechanical model of cleavage. These rules are as follows (see Fig. 1 and section A.5 in Appendix S1 for details).

In many developmental systems, cells tend to divide with their division plane perpendicular to the longest axis of the cell (Minc et al., 2011; Minc and Piel, 2012). It has been proposed that Hertwig's rule could explain part of the spiral cleavage pattern by translating the shape changes of blastomeres into specific directions of the cell division planes (Meshcheryakov, 1978).

Cells can be polarized, which can determine the direction in which cells divide. The direction of this polarization can be affected by factors in the surroundings of the cell or by asymmetries in the oocyte cytosol inherited by the blastomeres (see section A.5.2 in Appendix S1). Therefore, cells tend to divide perpendicularly to the direction of cell polarization. This can happen either because cell growth is biased towards the direction of polarization (Rogulja et al., 2008) (and then, according to Hertwig's rule, cells divide perpendicularly to the direction of polarization) or because polarization directly affects the direction of division (Morin and Bellaïche, 2011). It has been suggested that cell polarization in blastomeres could arise because of a gradient in the distribution of some molecules in the cell cortex (Freeman and Lundelius, 1982; Lu and Johnston, 2013). According to these authors, such an asymmetry would, by promoting a differential attachment of the astral microtubules to the part of the cortex with a higher concentration of those molecules, regulate the tilting of the mitotic apparatus relative to the A-V axis prior to cytokinesis and then bias cell divisions to take place in a specific direction (Freeman and Lundelius, 1982; Lu and Johnston, 2013).

The direction of division in a cell is also known to be affected by which parts of it are in contact with other cells. When this rule applies, cells would tend to divide towards [or away from in some cases (Wang et al., 1997)] the part of the cell contacting other cells. This has been suggested to occur because adhesion in a cell region would modify the underlying cell cortex so that astral microtubules are stabilized in this region, increasing their traction on the mitotic spindle (Hertzler and Wallis, 1992; Goldstein, 1995; Théry and Bornens, 2006). This rule has been proposed to play a role in the cleavage of some spiralian species [e.g. Tubifex worms (Takahasi and Shimizu, 1997)].

In spiralians, any cell division after the third round of division tends to be perpendicular to the previous cell division. This has been suggested as an explanation for the left-right alternation of cell divisions after the 4-cell stage (Guerrier, 1970; Meshcheryakov and Beloussov, 1975; Henley, 2012). This rule has been proposed to be a consequence of the stereotypic duplication and migration of the centrioles (which form a 90° angle between them) between cell divisions, that, in turn, biases the position of the mitotic spindle towards perpendicularity (Théry and Bornens, 2006; Minc and Piel, 2012).

The non-directional cell processes (those not affecting the direction of the division plane but that may affect spiral cleavage) are as follows.

According to Meshcheryakov's and related studies (Meshcheryakov and Beloussov, 1975; Wandelt and Nagy, 2004; Henley, 2012), blastomeres in spiral embryos rotate over themselves immediately after cell division. Rotating over themselves means that all the cell parts radially move with respect to an static rotation axis crossing through the cell centre (see Fig. S7). As a result of this rotation, the cell position does not change but the relative arrangement of its parts does, specially that of its surface. This rotation occurs in the same direction (clockwise or counter-clockwise around the rotation axis, depending on the embryo) in all blastomeres. The aforementioned authors argue that the left/right twist between the 4- and 8-cell stages and the relative placement of blastomeres is produced by this rotation (see section A.5.5 in Appendix S1). This rotation has been observed in many developmental systems (Meshcheryakov and Beloussov, 1975; Danilchik et al., 2006), but its potential role in cleavage remains unclear (Henley, 2012). This cell process does not specify the direction of cell division per se, but it may affect the mechanical interactions between cells after division and lead to changes in the relative positioning of the blastomeres.

Cell adhesion increases the contact surface between cells (the greater the adhesion strength the larger and flatter the contact surface between two cells) (Lecuit and Lenne, 2007; Sandersius and Newman, 2008). In this way, adhesion can lead to cell shape changes that may then affect blastomere shape, which blastomeres are in closer contact, the direction of cell division (e.g. through Hertwig's rule) and, overall, the shape and spatial arrangement of cells in the blastula.

Cell division can give rise to daughter cells of different sizes via several mechanisms. In many embryos, intracellular molecular gradients [e.g. PAR proteins in C. elegans (Cowan and Hyman, 2004)] regulate the relative size of daughter cells (Gilbert and Raunio, 1997; Munro et al., 2006; Salazar-Ciudad et al., 2003). Such gradients promote asymmetric cell division by a differential binding of the spindle microtubules to the cortex that results in asymmetric pulling forces and a displacement of the contractile ring during cytokinesis (Gilbert and Raunio, 1997; Ren and Weisblat, 2006; Lu and Johnston, 2013). Asymmetric cell division may play a role in generating morphological variation among the spiral cleavage patterns because blastomeres of different relative sizes are placed and compacted in different manners due to cell-cell adhesion (Merkel et al., 2012).

In this article we seek to identify a minimal (sufficient) set of rules and cell processes that are capable of generating the spiral pattern. More complex mechanisms could be imagined. For example, recent work in Lymnaea snails (Shibazaki et al., 2004) shows that, during the transition from the 4- to the 8-cell stage, dextral cleavage (which gives rise to adults with dextral shell coiling) exhibits cytoskeletal processes that are not observed in sinistral cleavage. This seems to be a pathological condition and the mutation involved has been identified (Shibazaki et al., 2004; Davison et al., 2016). Nevertheless, if a sinistral spiral cleavage can be produced without such additional cytoskeletal processes, it follows that these cytoskeletal processes are not strictly necessary for spiral cleavage. Instead, these processes would be superimposed on basic (and by default sinistral) processes (Wandelt and Nagy, 2004; Henley, 2012). This basic set of rules and cell processes that are capable of generating a spiral pattern, irrespective of handedness, is what we refer to as a minimal set.

We use a mathematical model, the SpiralMaker, which is simply an incorporation of the above-listed cellular processes into a general modelling framework for animal development previously developed by us (Marín-Riera et al., 2016). Therein each cell is represented as a set of subcellular elements (spheric elastic volumes) in three-dimensional space. The physical interactions between the subcellular elements (hereafter nodes) allows each whole cell to display visco-elastic properties like those observed in real cells (Newman, 2005; Sandersius and Newman, 2008; Marin-Riera et al., 2016). Thus, nodes can adhere to one another (preferentially to nodes from the same cell but also to nodes from other cells) but repel each other if they are too close (reiterating the physical fact that two cells can not occupy exactly the same position in space). These attraction and repulsion forces, together with some noise, lead to node movement and, consequently, to changes in cell shape and spatial location within the blastula. The model also incorporates cell behaviours such as cell division, polarization and cell adhesion.

All our simulations start after the second blastomere division (before this stage, typical spiral cleavage proceeds similar to other, non-spiral, cleavage types). Thus, in all our simulations the initial conditions consist of four blastomeres arranged in a square configuration (Fig. S1). In each blastula simulation, we only specify the initial conditions, and which rules are used and their relative strength. The 3D spatial position and shape of each blastomere over time, and thus the whole blastula pattern, is an output of the model that results from the its dynamics.

Based on the developmental rules and cell processes proposed above, we generate a theoretical developmental morphospace of possible spiral cleavage patterns. Each axis of this morphospace corresponds to one of the rules or cell processes implemented in the model. Along each of those dimensions, the relative contributions of the respective rule or cell processes are quantitatively varied. As we use developmental rules and cell processes, the resulting morphospace is not a geometrical but a generative developmental one: the cleavage patterns resulting from our model are not ordered according to some geometric parameters (e.g. the width/height ratio of the final morphology, as in Raup, 1961). Instead, the cleavage patterns are ordered according to the values of the developmental parameters by which they are generated, irrespective of the morphology of the final patterns they produce. Thus, the distribution of cleavage patterns within morphospace reflects how the different developmental rules and cell processes give rise to different spiralian cleavage patterns.

We explore all combinations of any three rules or cell processes. In addition, all simulations include cell adhesion whose strength is also systematically varied. The distribution of cleavage patterns within the morphospace reflects how the different developmental rules and cell processes give rise to different cleavage patterns. This allows us to evaluate in detail the cleavage pattern-forming capabilities of each rule and cell process, and their combination in a multi-cellular context (see Fig. S5).

As nodes in the model can contain molecules it is possible to implement intracellular molecular gradients. These gradients regulate some of the above explained rules and cell processes (most notably cell polarization and asymmetric division: see sections A.2 and A.5 in Appendix S1 for a more detailed description).

We also compare the spatial arrangement and sizes of the blastoromeres arising from the model with the spatial arrangement and blastomere sizes that we measure in three snail species. We focused on species that were most easily available to us experimentally or for which information about the cleavage pattern has already been published. Those species are the gastropods Crepidula, Planorbella and Lottia. For these three species we obtained 16-cell stage embryos and stained them to measure the volume of each blastomere. In addition, we obtained spiral arrangements from publications for Trochus (one species for each main snail group), the nemertean worm (ribbon worm) Carinoma, and the polychaete worms Nereis and Arenicola (Wilson, 1892; Robert, 1902; Newell, 1948; Freeman and Lundelius, 1982; Maslakova et al., 2004; Goulding, 2009). As individual blastomere volumes were not published, we used, instead, the list of blastomeres that are in physical contact (when available from the literature). These relative contacts between specific blastomeres, which are known to be crucial for the inductive mechanisms acting on later development (Lambert and Nagy, 2003; Grande and Patel, 2008; Kuroda et al., 2009), were then compared between empirical data and simulations. All these species are relatively synchronous in their early cell divisions, being symmetric in some species and asymmetric in others, as we describe later.

There are a number of previous models on cleavage. Some focus on the bio-mechanical properties (rheology) of blastula-like aggregates or make qualitative comparisons between these aggregates and mammalian blastulae (Honda et al., 2008; Sandersius and Newman, 2008). Some others tried, as here, to disentangle the morphogenetic processes responsible for the different types of cleavage that exist in animals. The model that is most similar to ours is that of Kajita et al. (2003) but it applies only to the first two divisions of Caenorhabditis elegans (a non-spiralian) and does not include the entire set of rules and cell processes we include here. Another recent model tries to explain the radial cleavage of the sea urchin (Akiyama et al., 2010). This latter model, however, is only two dimensional and it includes only three of the developmental rules we consider in here.

Fig. 2A shows the only combination of rules and cell processes in which the spiral pattern arises. This combination is the one that includes the cell polarization rule, Sachs' rule, cortical rotation and inter-cellular adhesion (see Fig. 2A). Other combinations of three rules or cell processes are also plotted (see Fig. S2) but they do not produce any cleavage pattern that resembles the spiral one. According to classical embryological descriptions (Gilbert and Raunio, 1997), a cleavage pattern is considered to be spiral only if: (1) the blastomeres are organized in groups of four cells (quartets) along the A-V axis, with the four blastomeres of each quartet forming a square in a plane perpendicular to the A-V axis; (2) the blastomeres closer to the animal pole were in closer contact between one another than the blastomeres closer to the vegetal pole; (3) sister blastomeres were obliquely positioned in respect to the A-V axis; and (4) this oblique positioning was in the same direction for the blastomeres in each quartet (either all to the right or all to the left of their sister blastomere).

The cell polarization rule in our model, if not combined with other rules, leads cell division to be oriented towards the animal pole (this is up and towards the centre of the embryo; Fig. 1A). The results of our simulations show that this rule is strictly required for spiralian cleavage to proceed normally after the third division, although it cannot lead to spiral cleavage on its own. With this rule, the four new micromeres arising in the third cell division (4- to 8-cell stage) are placed close to one another, enabling mechanical interaction between each other. This is important, as we explain below, for further spiral cleavage.

Cell adhesion tends to increase the contact surface between blastomeres. Because of that, each new cell in the model tends to gradually place itself, as it arises, between two other existing cells (in the case of the third division these latter cells are the macromeres). This leads to a relative displacement, or twist, of each micromere with respect to its sister macromere at the resulting 8-cell stage. If no other rules apply, this displacement is random in direction (on average 50% to the left and 50% to the right). As a consequence, each micromere does not necessarily end up correctly positioned between two macromeres as adjacent micromeres may twist in opposite directions, with both of them becoming positioned between the same two macromeres (a feature not found in spiral patterns in nature, Fig. S3A). In fact, at the 4-cell stage of typical spiral cleavage all cell divisions are directed either to the right or to the left. It is then clear that something in addition to the adhesion and polarization rules is required to explain the spiral cleavage pattern.

Our results show that, at the 8-cell stage, cortical rotation, and no other rule or cell process, ensures that all micromeres are displaced in the same direction with respect to their sister macromeres. Owing to cortical rotation, each blastomere at the 8-cell stage rolls over its neighbours in the same sense and places itself correctly between a distinct pair of underlying macromeres. The displacement of the micromeres occurs in all of them at the same time and we found that this requires adhesion between micromeres and macromeres, and between micromeres. In the former, the nodes in the micromeres tend to adhere more strongly to nodes from the macromeres and, as a result, each micromere tends to maximize its contact surface with the macromeres. If cell adhesion in the model is strong enough, this maximization of contact surface can be achieved only if each micromere is placed between two macromeres (as found in the spiral cleavage patterns). Adhesion between micromeres, instead, produces a friction-like effect in the contact areas between rotating micromeres, so that nodes in these contact areas move less than nodes elsewhere in the cell. As a result, the effective rotation of the cell surface introduced by the rotation rule is stronger in those areas than elsewhere in the cell. This causes a coherent displacement of the four micromeres in the same dextral or sinistral direction (Fig. S4). Without strong adhesion the rotation rule simply causes each micromere to rotate slightly around its own axis without any net cell displacement or further morphological changes in the blastula (Fig. S4). The crucial role of cortical rotation in early spiral cleavage is supported by the fact that, as experiments show, when embryos are exposed to actin depolymerization agents (Latrunculin), blastomeres do not rotate, and both dextral and sinistral blastulae are transformed in neutral, radially symmetric blastulae (Shibazaki et al., 2004).

Adhesion and cortical rotation together cannot account for the spatial distribution of blastomeres in the spiral pattern after the 8-cell stage. With these processes, blastomeres always divide towards the animal pole and form cone-shaped blastulae (Fig. S3B). However, our results indicate that the spiral pattern can be produced beyond the 8-cell stage if Sachs' rule is also included (Fig. 2A, encircled). With this rule, the fourth division round occurs at an angle from the A-V axis. This angle is not exactly 90° because, as described, the actual division direction is also affected by the polarization rule.

In the majority of our simulations each combination of rules and cell processes is applied to each embryo during the whole developmental time. To further explore the relative roles of each rule and cell process, we also performed simulations in which some of the rules were interrupted after some stage. This analysis shows that cortical rotation is crucial between the 4- and 8-cell stage to position the micromeres between the macromeres, thus forming an oblique angle in respect to the A-V axis. Thereafter, Sachs' rule leads divisions in successive rounds to be perpendicular to one another. Having the adequate oblique orientation at the 8-cell stage translates into an alternation between right and left oblique cell divisions from that stage onwards (although modified also by the polarization rule). Thus, from the 8-cell stage onwards cortical rotation is not required, as the polarization, Sachs' rules and adhesion are sufficient for the spiral pattern to arise. We also found that the asymmetric cell division did not affect whether spiralian patterns were found, although it was necessary for a more accurate reproduction of the cleavage of some species, as we explain below.

As Fig. 2B and Fig. S2 show, many cleavage patterns from the 8-cell stage onwards (60% of all combinations of rules) lack the tetra-radial symmetry along the A-V axis that is characteristic of many cleavage patterns in animals (spiralians and others). Those asymmetric cleavage patterns are found if Hertwig's rule is applied. According to this rule, the direction of division is parallel to the longest cell axis. As natural blastomeres tend to be spherically shaped, and in our model they are made of a finite number of nodes, small fluctuations around this spherical shape can completely change the direction in which the longest cell axis lies and, thus, the direction of division. The direction of division is then very sensitive to noise (noise is also present in real systems and not only in our model). Thus, when Hertwig's rule is applied, the direction of cell division in the third round of cell division can be different in each of the four macromeres [albeit in general, each micromere appears in the ‘upper’ (animal most) hemisphere of its sister macromere]. As a result, the four new micromeres might not lie at the same relative position along the A-V axis. These relatively small misalignments between macro- and micromeres at the 8-cell stage are amplified during further cell divisions, leading to blastulae with irregular cleavage patterns.

The spiral pattern is not found in the simulations that include the cell-cell contact rule (neither alone nor combined with other rules or cell processes). Our simulations show that with this rule, the direction of cell division tends to point towards the centre of the blastula (Fig. S3C). As no two cells can occupy the same physical space (in the model, cells that are very close repel each other), no extra cells can be placed in the centre of the blastula. As a result, daughter cells are passively displaced towards the periphery of the blastula over developmental time. In the resulting cleavage patterns, micromeres are not densely packed and are arranged in a radial manner that minimizes, or even removes, the contacts between them (Fig. 2C,D and Fig. S1C). This lack of contact between adjacent blastomeres diminishes the effect of cortical rotation because blastomeres cannot roll over their adjacent neighbours (or over their sister blastomeres) without contact. When the cell-cell contact rule is combined with the polarization rule, new micromeres appear towards the animal pole, but they tend to avoid the centre of the blastula (when viewed from the animal pole) so that they form a sort of cavity between them (Fig. 2C). The resulting open configuration prevents the emergence of spiral pattern, in which the four animal-most micromeres are always in contact with one another.

To further explore the accuracy of our model, we systematically varied the parameters in each of the rules and cell processes that we found to lead to spiral patterns and compared the resulting variation in spiral patterns with those of seven different spiralian species. These seven species were chosen because of their experimental availability but happen to have a rather synchronous cleavage. This facilitates our computational analysis for the stages considered in this study. In those simulations, asymmetric cell division was also implemented, so that in each cell division daughter cells would have different sizes according to their position on a gradient along the A-V axis (see Materials and Methods, and section A.5.3.2 in Appendix S1). In this way, we obtained a large set of simulated blastulae that were compared with the blastulae of some spiralian species. Asymmetric cell division was included in here but not in the analysis above because, in the seven species studied, simulated cells differ in size along the A-V axis; however, our simulations above show that asymmetric division per se is not required for achieving the spiral pattern itself.

For the gastropods Crepidula fornicata, Planorbella duryi and Lottia gigantea, we first experimentally obtained 16-cell stage blastulae, stained them and then measured the relative volumes of each blastomere (relative to the embryo volume; see section C.3 in Appendix S1). Within our simulated morphospace, we encountered blastulae closely resembling those of the three species mentioned above. In those simulated blastulae, the volume of each blastomere and its relative positioning is very similar to that of the corresponding blastomere in one of those three species. By corresponding blastomere we mean blastomeres of the same cell lineage in the simulations and in the empirically measured species (Conklin nomenclature, Gilbert and Raunio, 1997) (see Fig. 3A). According to the model, subtle differences in the spiral patterns of these three species arise from differences in the asymmetry of cell division, in adhesion and cortical rotation (see section C.4 in Appendix S1 for the exact quantitative differences). For Planorbella and Lottia, the visual resemblance between the simulated and the real blastulae improved when an eggshell surrounding the whole embryo was introduced (see Fig. 3 and section A.5.7 in Appendix S1). For Crepidula, although the volumes of the real and simulated blastomeres were almost identical, small discrepancies in the blastomere arrangement appear. These discrepancies could be attributed to variations in the mechanical properties (e.g. surface tension) between blastomeres of very disparate size that are not implemented in the model.

The spiral patterns of Trochus niloticus, Carinoma tremaphoros, Nereis diversicolor and Arenicola marina were compared with the ones found in the model, based on the descriptions present in the literature (Wilson, 1892; Robert, 1902; Newell, 1948; Freeman and Lundelius, 1982; Maslakova et al., 2004; Goulding, 2009). From these descriptions, it is possible to determine the overall topology of connections between blastomeres, i.e. which blastomeres (defined by cell lineage) are in physical contact with each other (but not the blastomere volumes) (see Materials and Methods, and section C.2 in Appendix S1). For Arenicola and Nereis, the simulations used initial conditions in which some the cells, the D-blastomere, were bigger than others (as reported previously for such species: Wilson, 1892; Newell, 1948; Freeman and Lundelius, 1982). As Fig. 4B-E shows, we found simulated cleavage patterns resembling those of these four different spiralian species (see section C.4 in Appendix S1 for a detailed description of the parameter values that produce each species blastula). This similarity is at least 90% (90% of the contacts between blastomeres were identical, see section C.2 in Appendix S1). Furthermore, our simulations reveal that the hypervolumes of the parameter space occupied by the different species differ in size and shape (Fig. S5).

Our results indicate that the spatial arrangement of cells characteristic of the spiral cleavage arises mainly due to the combination of an A-V polarization of cell division, Sachs' rule, cortical rotation and adhesion. Cortical rotation is important in the third division, but not later on, to produce a coherent twist or rotation of all micromeres with respect to the macromeres.

Our modelling of realistic bio-mechanics and rules in the context of 3D blastulae has allowed us to discard other rules and cell processes that have been suggested to be responsible for spiral cleavage (Meshcheryakov, 1978; Takahashi and Shimizu, 1997; Aw and Levin, 2009). Nevertheless, these other cell processes may be involved in other non-spiralian cleavage patterns. For example, Cnidarian-like cleavage patterns appear in our model when the direction of cell divisions is mainly determined by Hertwig's rule, and Ctenophora-like patterns arise when the direction of cell division is jointly determined by an A-V gradient and by the cell-contact rule (see Fig. 2B,C). However, whether these processes are involved, or not, in the generation of these and other non-spiralian patterns should be addressed by a specific study in the future (which could be also carried out using SpiralMaker).

Concerning the evolutionary transitions between these spiral and non-spiral patterns, our work provides two new insights. On the one hand, the spiral pattern is accessible from different points of the space of developmental parameters. This may explain why some non-spiralian taxa [e.g. some crustaceans (Nielsen, 1994)] exhibit spiral-like patterns during their early cleavage stages (they have arrived at a similar combination of cell processes via different evolutionary trajectories). On the other hand, the spiral pattern will be lost (as has happened many times in the evolution of Spiralia) whenever one of the cell processes and rules involved does not hold.

Our approach also suggests that up to the 32-cell stage, the formation of the cleavage pattern does not require signalling between cells. A molecular A-V gradient is required and such gradients are known to be common among metazoan oocytes (Raven, 1967; Gilbert and Raunio, 1997; Slack, 2014). This does not imply that signalling is not involved in it, but simply that it is not a logical necessity from what we know about other different cell processes that could also be involved in spiral cleavage. In addition, it is clear that cell signalling may be involved in cell fate determination in many species (Gilbert and Raunio, 1997). In fact, we think that the reason why the spiral cleavage is so conserved is due to this signalling. Invariant cleavage patterns are thought to be adaptive for those metazoans that have fast-developing planktonictrophic (feeding) larva (Salazar-Ciudad, 2010). To avoid sinking away from the plankton, these larvae need to be small but at the same time they need to have functional organs. The tiny organs of these larvae are made of one or a few cells. In these metazoans, fate determination occurs relatively early (compared with larger animals such as vertebrates or arthropods) through short-range signalling between individual cells (Davidson, 1991). As each cell-to-cell signalling event determines whole larval organs, any small change in these signalling events would easily lead to major deleterious consequences on the larval functionality. Because of this, changes in the relative positioning of blastomeres in early development is likely to be highly disadvantageous in spirals but less of a problem in organisms whose body is made of many cell organs (such as many cnidarians and vertebrates; Davidson, 1991; Salazar-Ciudad, 2010). As a consequence, one can expect that the early development and cleavage should be highly conserved (Salazar-Ciudad, 2010).

Our results also show that part of the morphological variation observed in spiralian cleavage can be explained by quantitative variations in few rules and cell processes (i.e. changes in the relative strength of each of such rules). Species-specific patterns, such as those of the seven different spiralian species considered here, arise in response to specific parameter combinations. Thus, in spite of its simplicity, our model seems to be indicative of the developmental processes underlying morphological differences at the level of blastulae between species. However, in some cases the simulated and real blastulae do not look visually identical. This is likely due to the intrinsically provisional or imperfect nature of our model, and of most other models. It could be that cell processes other than those we include in this model explain the 5% or 10% difference between the model and the reality observed. Nevertheless, with this model we can already assess which of the proposed rules and cell processes can, in principle, be involved in spiralian cleavage (and its variation) and which ones cannot.

The distribution of the different species analysed in the parameter space should provide information about which cleavage patterns ought to be easier to produce that other ones by changes in the underlying developmental parameters (Alberch, 1982; Salazar-Ciudad, 2007). We found that some cleavage patterns (e.g. patterns reminiscent of Trochus) can be produced by many different combinations of parameter values (they are found in a large volume of the parameter space), whereas others require a much more restricted combination of parameter values (e.g. 24% of the parameter space leads to patterns reminiscent of Arenicola). We also find that the cleavage patterns of species that are phylogenetically closer (e.g. Arenicola and Nereis) are not closer to each other in the parameter space than to the patterns of other, more phylogenetically distant, species. This suggests that the underlying developmental parameters are relatively easy to change (presumably a small number of mutational changes are required) (Newman, 2011).

The mathematical model used in this study, the SpiralMaker, is a modification of a general model of animal embryonic development, the EmbryoMaker (Marin-Riera et al., 2016), that allows us to implement each of the developmental rules studied in this article. [The EmbryoMaker software implementing the model is available for download (http://www.biocenter.helsinki.fi/salazar/software.html). SpiralMaker is available from the same site.] In the EmbryoMaker, cells are made of subcellular elements (nodes), which can be conceptualized as spherical parts of a cell that occupy a physical volume. Nodes adhere to each other if they touch but repel each other if their centres become to close. A node of a cell adheres with higher affinity to nodes from the same cell than to nodes from different cells (this ensures that cells maintain their physical integrity). Motion is computed in continuous time and 3D space by solving a system of differential equations assuming Langevin over-damped dynamics and some degree of noise (see section A.3 in Appendix S1). As a result of cell division, adhesion, repulsion and noise, the positions of nodes (and thus the position and shape of blastomeres) change over simulation time. This 3D spatial distribution of those nodes represents the morphology of the embryo and within each cell it represents cell shape. Thus, cell shape and position within the embryo, and their variations, are not a free parameter of the model, but result from model biomechanics. In turn, changes in the spatial location and shape of cells configure the overall changes in embryo morphology. Mechanical interactions between nodes, such as the relative cell movement, occur faster than whole-cell processes such as cell division. In the model, this is ensured by choosing long time intervals between cell division rounds.

Because, in most spiralians, cleavage takes place without cell growth, the number of nodes in the SpiralMaker is kept constant during each simulation. Thus, the number of nodes in a cell halves (for symmetric cell division) in each cell division round whereas the number of cells doubles. The SpiralMaker includes only non-epithelial cells, as blastomeres, and only the cell behaviours of cell division, cell polarization and cell adhesion (other types of cells and cell behaviours are included in the EmbryoMaker).

The model includes: (1) an initial condition of the 4-cell stage in which cells have an intracellular molecular gradient along the A-V axis (presumably inherited from the zygote); (2) cells dividing synchronously at regular time intervals; (3) the use of different rules to determine the direction of cell division; (4) the use of different slopes of the intracellular gradient that affects the asymmetry of cell division; (5) the use of different degrees of cortical rotation; and (6) the use of different values of adhesion between nodes. All those things are specific to the SpiralMaker but not to the EmbryoMaker from which it derives.

The developmental rules and cell processes of the SpiralMaker are implemented in the model as simple rules acting on the elements of the cells (see section A.5 in Appendix S1 and Fig. S6). In an initial experiment (see in silico experiment 1 in Appendix S1), we explored whether the combination of these rules and cell processes can generate the spiral pattern or not.

Each simulation used three different rules (plus adhesion, which is also a cell process) and different weighting of the relative importance of each rule and cell process when determining the direction of cell division. Other, high-dimensional combinations of rules and cell processes (e.g. four or five rules and processes by simulation) were not considered. We discarded this possibility because our goal (the emergence of spiral cleavage pattern) was achieved under the simpler combination of only three rules and processes (plus adhesion). For this initial exploration, the number of nodes in the blastula were 1250, cell division was kept approximately symmetric and all simulations were run until the 16-cell stage.

In a more detailed approach (see in silico experiment 2 in section C of Appendix S1) we assessed how variations within the spiral pattern (including the cleavage patterns similar to those of several invertebrate species) can be produced. To achieve this, we took all the combinations of rules that, in experiment 1, were found to produce spiral patterns and performed a more exhaustive sampling of our theoretical developmental morphospace in each rule and cell process parameter (increments of 20% in each parameter from its minimum possible value to its maximum possible value). For this experiment, the number of nodes per simulation was 2500 and simulations were run until the 32-cell stage. In addition, asymmetric cell division was included, and we also considered asymmetry in the size of the blastomeres in initial conditions. We also ran some simulations in which an outer cover (or eggshell) surrounds the blastula (see section 5.7 in Appendix S1).

The output of these simulations, which are realistic 3D representations of 16- and 32-cell stage embryos, were then compared with those of real spiralians. In order to do so, two different quantitative methods were developed (one based on the relative volumes between blastomeres and the other based on the relative contacts between adjacent blastomeres, see section C in Appendix S1). By using these methodologies, simulated cleavage patterns reminiscent of those of several spiralian species (Crepidula, Planorbella, Lottia, Trochus, Carinoma, Nereis and Arenicola) were found within our theoretical developmental morphospace. The data for the first three species (representatives of the main gastropod taxa) were collected from fixed and immunohistochemistry-stained embryos, which were microphotographed with a confocal microscope (Zeiss) and digitally processed to obtain an estimation of the relative volumes between their blastomeres (see section C.3 in Appendix S1 for details about sample collection and handling). Data for the remaining species were collected from bibliographic publication in which clear depictions of both animal and lateral views of the 16-cell stage were available (see section C.2 in Appendix S1).

We thank R. Zimm, I. Salvador-Martinez, P. Hagolani, Angus Davison and two anonymous referees for comments. We also thank the Hybridoma bank at Iowa University for developing the mouse anti-ß-tubulin E7 antibody used in this study and the Marine Biological Laboratory (Woods Hole, MA, USA) (especially Drs Nipam Patel and Lee Niswander) for supporting the collection and fixation of C. fornicata embryos.