Dynamic transcriptional symmetry-breaking in pre-implantation mammalian embryo development revealed by single-cell RNA-seq

Mammalian pre-implantation embryos provide a unique opportunity to study how a single cell diverges from a totipotent state into different fates within limited rounds of cell division (Zernicka-Goetz et al., 2009). Morphologically, distinct cell lineage patterning first emerges during the fourth cleavage division to generate a 16-cell embryo. At this stage, some blastomeres divide symmetrically, contributing two daughter cells to the outside region of the embryo, whereas others divide asymmetrically and contribute one daughter cell to the outside and another to the inside region (Rossant and Tam, 2009). This difference in cell allocation has long been considered to be the first emergence of asymmetry in early embryo development; the ‘inside cells’ will contribute to the inner cell mass (ICM), whereas the ‘outside cells’ will contribute to the trophectoderm (TE) (Rossant and Tam, 2009; Takaoka and Hamada, 2012; Zernicka-Goetz et al., 2009).

However, contrary to this ‘late-deciding’ model, recent emerging evidence from studies on morphologically indistinguishable blastomeres before the 16-cell stage has revealed distinctive molecular markers in blastomeres during the four- to eight-cell stages that are indicative of future lineage segregations, leading to the idea that the late-onset event of ‘inside’ and ‘outside’ cell divergence is not a random choice but rather an anticipated result rooted in the embryonic division history (Torres-Padilla et al., 2007). As early as the four-cell embryo stage, different levels of CARM1-mediated histone 3 arginine methylation can be found in blastomeres, with a higher level of methylation directed to the ICM lineages (Torres-Padilla et al., 2007). At the four- to eight-cell stage, two distinct patterns of OCT4 (POU5F1) kinetics are observed between blastomeres, reflecting the different levels of OCT4 binding affinity to DNA, with the slower OCT4 kinetics (higher binding ability to DNA) tending to give rise to the ICM (Plachta et al., 2011). Further evidence obtained using the state-of-the-art Rainbow lineage-tracing system has shown that four-cell blastomeres do indeed display developmental biases in contributing to either the ICM or TE (Tabansky et al., 2013). These advancements have promoted the concept that the molecular differences in blastomeres emerge much earlier (as early as the four-cell stage) than the morphological clue.

Despite these observations, the open question remains regarding when the very first molecular asymmetry emerges, and how such initial asymmetry develops with the potential to direct lineage bifurcation. The emerging technology of single-cell RNA-seq has provided the opportunity to study the transcriptome-wide dynamics of the transcriptional symmetry-breaking process during pre-implantation mammalian development. In the present study, combining single-blastomere transcriptome analysis from human and mouse pre-implantation embryos (Deng et al., 2014; Yan et al., 2013), we show that random segregation at the first cleavage division, known as ‘partitioning errors’, provides an important source of initial blastomere-to-blastomere heterogeneity. During the two- to 16-cell stages, zygotic transcriptional activation generates different molecular feedbacks that minimize or enhance the initial blastomere-to-blastomere biases, resulting in different gene clusters with either a ‘monostable’ (ubiquitous expression between blastomeres) or ‘bistable’ (strongly asymmetric expression between blastomeres) pattern. For those genes with a bistable pattern, the relative ratio of opposing lineage specifiers in each blastomere creates an inclined ‘lineage strength’, providing a transcriptional clue to guide future cell fates ahead of morphological distinction. We believe that these analyses represent a major step in understanding the very first mammalian embryonic symmetry-breaking process, as a result of both chance distribution and dynamic transcriptional regulation before morphological distinction.

When considering the source of cell-to-cell heterogeneity, there are two major theoretical models. One is partitioning error, which occurs when, at cell division, cellular substances are unevenly distributed into the two daughter cells (Huh and Paulsson, 2011a). The other is ‘stochastic transcription noise’, which is derived from random transcriptional starts within each cell cycle (Cai et al., 2006; Elowitz et al., 2002; Raj et al., 2006). These two models are usually mixed-up in a real cell system, both contributing to cell-to-cell heterogeneity, and are therefore very difficult to analyze separately (Huh and Paulsson, 2011b). However, unlike normal mitotic cell division during which gene transcription is always active, mammalian zygotic transcription is mostly silent before the first embryo cleavage (Li et al., 2013; Tadros and Lipshitz, 2009). Therefore, the first mammalian embryo cleavage division provides a unique system in which partitioning error comprises the dominant source of biases between two-cell blastomeres, and its contribution to initiation of blastomere-to-blastomere biases can be calculated mathematically and tested by single-blastomere sequencing data.

By considering mammalian zygotes as a homogenous sphere and the first embryonic cleavage division as a binomial distribution system, we calculated the theoretical chance (Eqn 1 in Materials and Methods) of a biased distribution (with 10%, 15% and 20% differences between two daughter blastomeres) for a putative substance with regard to their initial copy numbers, and showed that the chance of an uneven distribution depends on their initial quantity (Fig. 1A): the substances present in small quantities are more likely to be differentially distributed into the two daughter cells, whereas those present in larger quantities have a greater chance of being equally divided (Fig. 1A). To test this binomial model in a real mammalian embryo cleavage division, we analyzed single-blastomere RNA-seq data from mouse and human two-cell blastomeres (Tang et al., 2011; Yan et al., 2013). As shown in a representative two-cell mouse embryo, the overall comparative transcriptome profile between two-cell blastomeres showed that highly expressed genes are usually similarly expressed, whereas differential expression is more often observed for genes with lower expression (Fig. 1B), which is similar to the mathematical prediction. To test more stringently the accuracy of this binomial model in a two-cell embryo, especially the extent of bias distribution with regard to the copy number of a gene transcript, we next estimated the transcript numbers in each two-cell mouse blastomere by combining qPCR and single-cell transcriptome reads (Materials and Methods). This estimate showed that every five transcripts approximately equaled one RPKM (reads per kilobase transcriptome per million reads) value in a two-cell blastomere, which was similar to what we previously reported (Tang et al., 2009). Then, by using the estimated transcript number for each gene, the blastomere-to-blastomere heterogeneity at the transcript level could be tested against the theoretical curve under binominal partition (Materials and Methods). The statistics of transcript partitioning at the two-cell embryo stage showed that the RNA copy numbers detected in two-cell blastomeres were approximately aligned with the binomial prediction (Fig. 1C,D; Fig. S1A,B), especially when the copy numbers were low (<100) (Fig. 1D; Fig. S1B). This analysis demonstrated that during the first embryo cleavage division, small biases initially arise following a binomial distribution pattern, with the low-quantity transcripts being more likely to be asymmetrically distributed.

Once the initial bias is generated in two-cell blastomeres under binomial partitioning, where will it end? It has been previously shown that transcription noise is an important source of cell-to-cell heterogeneity (Raj et al., 2006). Therefore, it is our hypothesis that the increasing level of zygotic transcriptional activation after the middle-to-late two-cell-embryo stage should increase the level of blastomere-to-blastomere bias. To test this hypothesis, we utilized two methods to analyze quantitatively the stage-specific blastomere-to-blastomere heterogeneity using single-blastomere RNA-seq datasets of two-cell mouse embryos at the early, middle and late stages (Deng et al., 2014) (Fig. 2A).

First, we calculated the expression noise between RNA-seq technical repeats of one oocyte/blastomere and those from two-cell blastomeres based on a previously published formula (Swain et al., 2002) (Eqn 2 in Materials and Methods). This calculation clearly showed that the extent of blastomere-to-blastomere expression bias from the early two-cell embryo was higher than that of the single-oocyte/blastomere technical repeats (Fig. 2B). It can be assumed that this difference represents the contribution of partitioning error because early two-cell blastomeres have only just been formed by cleavage division and zygotic transcription is almost silent at this stage (Lee et al., 2014) so the contribution of transcriptional noise is minimal. The early-to-late embryo progression showed an ever increasing expression bias, which presumably represents the contribution of transcriptional noise (Fig. 2B) because during this period there are no cleavage divisions, therefore transcriptional noise becomes the only source of increased asymmetry. These analyses dissect the relative contribution of partitioning errors and transcriptional noise that contribute to transcriptional asymmetry between two-cell blastomeres in a stage-dependent manner, supporting the hypothesis illustrated in Fig. 2A.

Our second method for representing blastomere-to-blastomere expression asymmetry is as follows: for a gene detected in a two-cell embryo, one blastomere will show E_higher and the other will show E_lower (where E represents gene expression level), and the value of E_higher/E_lower will represent the extent of asymmetry for this gene between the two blastomeres. When the two examined blastomeres are more ‘identical’, the calculated E_higher/E_lower values for more genes should be close to 1, whereas if the two blastomeres are more ‘distinct’, more genes will show increased E_higher/E_lower values. In this context, the numbers of genes with 1<E_higher/E_lower<1.5 and E_higher/E_lower>2 were calculated for each two-cell embryo (early, middle and late stages) and for the technical repeats (Fig. 2C). The results clearly showed that the number of genes with 1<E_higher/E_lower<1.5 was highest in the technical repeats and continually decreased during the early-to-late embryo progression, whereas this trend was reversed for the number of genes with E_higher/E_lower>2 (Fig. 2C). These analyses confirmed the conclusion drawn from the data in Fig. 2B, both showing that the first embryo cleavage division generates blastomere-to-blastomere bias as a result of partitioning error, and that during the early-to-late two-cell embryo progression, the blastomere-to-blastomere bias continually increases owing to zygotic transcriptional activation.

To analyze further the expression asymmetry in four-cell and eight-cell embryos, we performed a pairwise comparison for all blastomeres within an embryo to reveal the expression asymmetry between each pair of blastomeres. As shown in Fig. 3A,B and Fig. S2, the extent of blastomere-to-blastomere asymmetry showed an overall increase from the two-cell stage to the eight-cell stage in both mouse and human embryos, which could be due to continuous transcriptional activities and multiplied partitioning error (Eqns 3-5 in Materials and Methods; Fig. 3A) during embryo development.

In addition to the observed overall increase in gene expression asymmetry along embryonic development (two-cell to eight-cell stage) (Fig. 3B; Fig. S2), it is notable that during continuous embryonic cleavage, the available amount of embryonic sample for each RNA-seq decreases and this will cause increased technical noise, especially for genes expressed at a low level (Streets et al., 2014), which might partially account for the observed asymmetry increase. Because highly expressed genes are less affected by both technical noise and partition errors, we further selected a subgroup of genes from each mouse embryonic stage with both a strong expression bias (top 20%) and high expression level (top 20%), which we believe will more faithfully represent the level of transcription and minimize the potential influence of technical bias (which affects to a greater degree genes with lower expression). For these selected genes (highly expression level with strong bias) (Table S1), we analyzed the trends of blastomere-to-blastomere asymmetry during the two- to 16-cell stages, and found that the extent of asymmetry at a certain stage proceeds with two distinct patterns: for some genes, the extent of asymmetry tends to be minimized between blastomeres, whereas other genes become increasingly asymmetrically distributed (Fig. 4A-C). Notably, at the eight-cell stage, only very few genes with high asymmetry showed decreased asymmetry onward (Fig. 4C), suggesting that, compared with two- to four-cell embryos, the transcriptional bias between eight-cell blastomeres has become increasingly irreversible and deterministic. Interestingly, for those genes with strong bias distribution at the 16-cell stage, gene ontology analysis revealed that most of these genes are cataloged in, for example, ‘positive regulation of cell differentiation’, ‘positive regulation of developmental process’ (Fig. 4D), implicating their potential in directing different cell fates. In addition to the overall analysis, the actual asymmetry trends of three typical genes were examined during the two-cell to 16-cell stages: a housekeeping gene (Tubb2c, which encodes the tubulin beta-2C chain; also known as Tubb4b), a transcription factor for self-renewal (Pou5f1, also known as Oct4), and a lineage specifier with a defined role starting at the four-cell stage (Carm1, the histone arginine methylase). Each of these genes exhibits a different pattern of gene expression asymmetry, which may be described as overall-flat (Tubb2c), wave-like (Pou5f1) or ever-increasing (Carm1) (Fig. 4E), suggesting dynamic transcriptional regulation for each gene.

The observed distinct asymmetry dynamics for different genes are, in general, supposed to be regulated by series of molecular chain reactions (MacArthur et al., 2012) with either negative or positive transcriptional circuits. This could be simplified as two distinct models involving a negative feedback-controlled ‘monostable pattern’ (Fig. 4F) or a positive feedback-controlled ‘bistable pattern’ (Fig. 4G) (Materials and Methods; Fig. S3), with the latter (such as the pattern of Carm1) being a potential driving force for transcriptional symmetry-breaking in embryo development.

In a real biological system, the divergence of lineage is usually decided by more than one factor, and recent evidence from stem cell research has revealed that lineage specifiers with counteracting effects work in balance to maintain pluripotency, whereas a tilted ratio leads to lineage differentiation (Montserrat et al., 2013; Shu et al., 2013). Such a scenario similarly exists in the context of early mouse embryos, in which biased lineage-driving forces are crucial in guiding future cell fates (Bruce and Zernicka-Goetz, 2010). Here, we proposed that the relative ratio of a pair of opposing lineage specifiers is influenced by both partitioning error at cleavage division, and transcriptional regulation within the cell cycle. By supposing that two lineage specifiers have an initial ratio that maintains the cell in an undifferentiating state (Fig. 5A), we first calculated the probability of tilting the initial ratio (to an extent of 10%) of opposing lineage specifiers based on partitioning error (Eqn 6 in Materials and Methods), and showed that the odds of breaking the initial ratio are inversely correlated with the level of either one or both substances before cleavage (Fig. 5B), suggesting that the copy number of lineage specifiers deposited in the previous cycle can influence their distribution pattern (relative ratio with other lineage specifiers) in the daughter cells, thus contributing to the bifurcation of future cell fate.

Transcriptional regulation is another important contributor that can tilt the ratio of lineage specifiers. For example, Cdx2 and Carm1 are a known pair of lineage specifiers, which guide different cell fates in mouse early embryo. The expression of Cdx2 at the two- to four-cell stages is very low (mostly inherited from maternal storage), but its expression increases drastically in the eight-cell embryo (Fig. 5C); thus, de novo transcription at this stage becomes the dominant factor affecting the ratio of Cdx2 to Carm1. Whereas for Carm1, transcription of which increases only mildly at the eight-cell stage (Fig. 5C), the Carm1 to Cdx2 ratio is affected by both transcripts inherited from previous divisions and those newly transcribed in the current cell cycle.

As Carm1 and Cdx2 have been reported to direct distinct cell fates (Jedrusik et al., 2008; Torres-Padilla et al., 2007), it is our hypothesis that their relative ratio within an eight-cell blastomere provides clues for the future lineage specification. We therefore analyzed the relative expression levels of both Carm1 and Cdx2 in each blastomere from eight-cell mouse and human embryos, and showed that different eight-cell blastomeres indeed have a ‘wax and wane’ phenomenon regarding the relative ratio of Carm1 and Cdx2. As shown in Fig. 5D-G, most blastomeres are either dominated by Carm1 or Cdx2, suggesting that a predisposed lineage strength has been formed in these blastomeres. It is notable that some blastomeres showed a similar expression level of both Carm1 and Cdx2, thus representing an undecided, bi-potent state as suggested previously (Bruce and Zernicka-Goetz, 2010). According to our analysis, when these bi-potent blastomeres proceed into next cell cycle, the relative ratio of Carm1 and Cdx2 could be further tilted as a result of both partitioning error and transcriptional regulation, eventually creating a situation in which one side wins out on the other.

Moreover, we have further provided additional lists of genes showing the ‘wax and wane’ phenomenon with regard to Carm1 and Cdx2 in eight-cell blastomeres (Table S2). These genes could be potential lineage specifiers and their function warrants further investigations.

Finally, based on our single-blastomere transcriptome data analysis, a scenario of early transcriptional symmetry-breaking and lineage specification is summarized in Fig. 6. In this model, a static view before the final lineage destination will give an impression of an undecided lineage combat (Fig. 6), showing a heterogeneous level of opposing transcription factors (i.e. a ‘stochastic appearance’) that could be either inherited from a previous division or produced by de novo transcription, i.e. also having a ‘traceable origin’.

The mechanisms of early mammalian cell fate determination remain debatable as it has not yet been established whether the first bifurcation of cell fate in mammalian embryo emerges randomly at the morula stage (Dietrich and Hiiragi, 2007; Wennekamp and Hiiragi, 2012) or if the molecular clues of differentiation emerge before morphological distinction (Plachta et al., 2011; Tabansky et al., 2013; Torres-Padilla et al., 2007; Zernicka-Goetz, 2013). In the present study, our single-blastomere transcriptome analysis of early embryo development have showed that the earliest blastomere-to-blastomere biases emerge with the first cleavage division, owing to random segregation, but subsequent zygotic transcriptional activation triggers transcriptional regulation that fine-tunes these small biases in a more defined manner, minimizing or amplifying the initial biases with negative or positive feedback mechanisms. We believe that only those transcriptional biases that finally develop into a bistable pattern (strongly asymmetric between blastomeres), bear the potential to guide lineage fates. Moreover, our analysis supports a scenario in which opposing lineage specifiers within an early blastomere constantly compete with each other based on their relative ratio, pushing the blastomere onto a predisposed, yet flexible, lineage track before morphological distinction. These analyses revealed mammalian early embryo symmetry-breaking as a continuous process rather than a sudden emergence, with the driving force involving both chance separation and transcriptional circuits.

As shown in our analysis, at the first cleavage division, small biases will inevitably arise in a binomial distribution pattern, with lower-quantity substances bearing a greater chance of being asymmetrically distributed (Fig. 1A-D). Such quantity-dependent asymmetric distribution in two-cell blastomeres, in our opinion, is an important step in resolving the dilemma of how two-cell-stage blastomeres can be both ‘identical’ (in totipotency) and ‘different’ (in guiding future lineages) (Zernicka-Goetz, 2006). In a real biological system such as a zygote, the substance found in greater quantity is usually required for maintaining basic cell properties, whereas those present in small amounts are usually elements with fine-tuning regulatory functions (transcriptional factor, non-coding RNAs, etc.). Therefore, these regulatory factors tend to be unequally distributed in two-cell blastomeres as a result of partitioning errors, which may trigger downstream events that influence future cell fate. Moreover, apart from the biased transcript distribution, it could be repeatedly observed that, among the two-cell blastomeres (which are generated at the same time), one cell always divides before the other (Fig. S4), supporting the idea that at the late two-cell stage, the two-cell blastomeres already bear slight biological differences regarding their intrinsic properties.

Recently, it is notable that two groups also used single-blastomere transcriptome analysis to study the early symmetry-breaking process. One report showed that as early as the two-cell embryo stage, some genes have already become strongly asymmetrically expressed in different blastomeres and it has been suggested that these genes may be involved in directing future lineages (Biase et al., 2014). The other report further showed that the blastomeres have become strikingly variable at eight-cell stage owing to the amplification effect of transcriptional activities (Piras et al., 2014). These conclusions are in general in accordance to our analysis, as we showed an overall increase in blastomere-to-blastomere biases driven by both zygotic transcriptional activation and multiplied partitioning errors during ongoing embryonic cleavage (Fig. 2B,C and Fig. 3A,B). However, when we performed more stringent analysis by filtering out genes with low expression, which in large ruled out the influence of technical biases, we found a more dynamic framework of transcriptional regulation during the two- to 16-cell stages, showing two distinct patterns of progression: for some genes, the extent of blastomere-to-blastomere asymmetry tends to be minimized, displaying a monostable pattern [in general, the single stable equilibrium in gene expression dynamics is a result of negative feedback regulation (Alon, 2007)], whereas for other genes, the extent of asymmetry becomes increasingly larger, displaying a bistable pattern, as a result of positive feedback regulation as previously pointed out (Smits et al., 2006).

These analyses revealed that although strong biases of transcripts could be observed in early blastomeres as early as the two- to four- cell stages, it is not a guarantee that these early observed biases will be kept at the later embryo stages; alternatively, only those crucial molecules (mRNA, protein or non-coding RNA) with the ability to trigger self-reinforcement or gene regulatory feedbacks (Davidson, 2010) bear the ability to drive initial small biases to produce tangible biological differences.

Despite the observed dynamic transcriptional symmetry-breaking during early embryo development, we would like to emphasize that even if a lineage specifier is asymmetrically distributed between early blastomeres, they could only ‘guide’, but not ‘decide’ the lineage track until other more definitive clues such as inner-outer position, cell polarity, cell-cell contacts (Lorthongpanich et al., 2012) emerge. However, once a blastomere has acquired dominant lineage specifiers (for example, Carm1) at the four-cell stage, the chance of keeping the dominant position in its daughter generation will be higher than those blastomeres with low level expression. This is, in principle, consistent with the biological observations that these Carm1-dominant blastomeres showed a biased cell fate as previously published (Tabansky et al., 2013; Torres-Padilla et al., 2007). Also, previous experiments have shown that if the transcript of a lineage specifier is exogenously injected in a subset of blastomeres, it will guide the lineage direction in a more definite manner (Cockburn et al., 2013; Jedrusik et al., 2008; Torres-Padilla et al., 2007), most probably because the injection of exogenous transcripts secured its dominant expression level, which overwhelms the fluctuations caused by other factors (such as partitioning error and transcription regulation) in its daughter cells.

Regarding the lineage competition by opposing lineage specifiers in the early embryo, it remains unclear at what time the competition first begins and by what mechanisms the winner steadily establishes its cell lineage. In the present study, we revealed an interesting phenomenon that the relative ratio of Carm1 and Cdx2 (which guide ICM or TE specifications, respectively) in each eight-cell blastomere is different, suggesting that inclined lineage strengths might have been formed in these eight-cell blastomeres. However, the lineage track in these blastomeres have not been fully decided, because in the next round of division, the ratio of Carm1 to Cdx2 could be reversed as a result of both random segregation and transcriptional regulation. Nonetheless, such a ratio-dependent lineage battle will continue, until one side finally becomes strong enough to consolidate its dominant position and define the lineage fate (Fig. 6) by introducing other more definitive clues, such as inner-outer position, cell polarity, etc. Similar lineage competition and wax and wane of opposing lineage specifiers (Fig. 6) may also exist in the second wave of cell specification of the ICM into the epiblast and primitive endoderm (Ohnishi et al., 2014), which would involve both a traceable cell origin from their lineage ancestors (Morris et al., 2010) and the plasticity to change cell fates before the final differentiation (Yamanaka et al., 2010).

In summary, our single-blastomere transcriptome analysis revealed a dynamic transcriptional symmetry-breaking process occurring far earlier than morphological distinction, contributed to by both random segregations at each embryo cleavage division and transcriptional regulatory feedbacks. In the future, more detailed verification of our analysis with regard to lineage specification will depend on our improved ability to examine precisely the quantity of molecules (proteins, RNAs) in an amplification-free manner, such as by single-molecule fluorescent in situ hybridization detection of selected genes at single-cell resolution (Itzkovitz et al., 2012) with time-lapse observation, as well as functional examination of candidate genes regarding their roles in lineage specification.

We analyzed single-blastomere RNA-seq datasets from pre-implantation human embryos (Yan et al., 2013; data available in Gene Expression Omnibus under accession number GSE36552). Single-blastomere RNA-seq data from mouse pre-implantation embryos have been previously published (Tang et al., 2011; Deng et al., 2014; data available in Gene Expression Omnibus under accession numbers GSE22182 and GSE45719, respectively). Gene expression was calculated as RPKM (reads per kilobase transcriptome per million reads; Audic and Claverie, 1997; Mortazavi et al., 2008).

Eight- to ten-week-old C57BL/6 female mice were stimulated to superovulate by standard methods (Wan et al., 2013) and then mated with DBA2 male mice. The animal-use protocols in the present study were approved by the Animal Research Committee of the Institute of Zoology, Chinese Academy of Sciences. Two-cell embryos were collected and cultured in vitro in M2 medium under mineral oil, and their progression to four-cell stage was video recorded by time-lapse microscopy (UltraVIEW VoX 3D Live-Cell Imaging System). The time lag between the three-cell and four-cell stages of each embryo was recorded and analyzed.

Genes with both high noise (top 20%) and relative high expression level (top 20%) at each embryonic stage were selected, and the pattern of these genes' noise during the progression from two-cell to 16-cell stages was clustered by Genesis (Sturn et al., 2002), using a hierarchical clustering method, with ‘Complete linkage clustering’ agglomeration parameter in the software.

GO enrichment was analyzed using DAVID (http://david.abcc.ncifcrf.gov/). A hypergeometric test was performed using the default parameters to adjust the P value.

The real count of Hprt in one MII mouse oocyte is ∼4829.2, according to a previous report (Steuerwald et al., 2000), and the relative gene expression of Hprt between an MII oocyte and a two-cell-stage blastomere was obtained by real-time PCR. The timing of oocyte and two-cell embryo retrieval was the same as we previously described (Tang et al., 2011). The calculation showed that 1 RPKM that corresponded to 4.465 counts for one two-cell blastomere.

The effect of negative feedback regulation is keeping the expression level around one stable state (monostable), whereas an auto-activating positive-feedback loop is able to exhibit bistable states, which may have three fixed points, including two local stable states (one high and one low) and one unstable state.

We simulated the monostable and bistable process based on the following assumptions:

1. The whole process is divided into two parts: determinate dynamics behavior during the cell cycle and mock process during the cell division period (Fig. S3). Fluctuations from gene expression are neglected, i.e. the process of gene expression is considered to be determinate as is its ordinary differential equation (ODE) dynamics.

2. Duration (or period) of cell cycle T is long enough that the system state can reach the corresponding local stable state in each cell cycle. In mathematical terms, that is γ, Γ≪T.

3. Production of gene expression, i.e. protein and mRNA, are both considered.

4. The total volumes of the embryo are unchanged in the first six cell divisions (from one cell to 32 cells). Only this cell period is considered.

The parameter values we selected in simulation were: K=30 h⁻¹, k₀=0.1 h⁻¹, ϒ=30 h⁻¹, Γ=1 h⁻¹, k_H=10 (Vilar et al., 2002; Zheng et al., 2011).

Array B_i contains the percentages between reference gene b and reference gene a of each blastomere in order number i embryo. In the present manuscript, genes a and b can be Carm1(Cdx2) and Cdx2(Carm1). The genes with correlation coefficient R between array A_i and B_i >0.7 in every embryo were selected, deemed as highly correlated with the expression pattern of reference gene b.

Statistical analysis was conducted using GraphPad Prism 6.0 software. One-way ANOVA was used for statistical analysis. For all statistical analyses, a value of P<0.05 was considered statistically significant.

We thank Drs Lei Li and Haibin Wang at the Institute of Zoology, Chinese Academy of Sciences for helpful discussions; and Mr Ming Ge at the Institute of Zoology, Chinese Academy of Sciences and Ms Yue Wang at the National Institute of Biological Sciences, Beijing, for their help in mouse breeding and embryo manipulation.