Abstract STUDY QUESTION What is the prevalence, reproducibility and biological significance of transcriptomic differences between sister blastomeres of the mouse 2-cell embryo? SUMMARY ANSWER Sister 2-cell stage blastomeres are distinguishable from each other by mRNA analysis, attesting to the fact that differentiation starts mostly early in the mouse embryo; however, the interblastomere differences are poorly reproducible and invoke the combinatorial effects of known and new mechanisms of blastomere diversification. WHAT IS KNOWN ALREADY Transcriptomic datasets for single blastomeres in mice have been available for years but have never been systematically analysed together, although such an analysis may shed light onto some unclarified topics of early mammalian development. Two unknowns that remain are at which stage embryonic blastomeres start to diversify from each other and what is the molecular origin of that difference. At the earliest postzygotic stage, the 2-cell stage, opinions differ regarding the answer to these questions; one group claims that the first zygotic division yields two equal blastomeres capable of forming a full organism (totipotency) and another group claims evidence for interblastomere differences reminiscent of the prepatterning found in embryos of lower taxa. Regarding the molecular origin of interblastomere differences, there are four prevalent models which invoke (1) oocyte anisotropy, (2) sperm entry point, (3) partition errors of the transcript pool and (4) asynchronous embryonic genome activation in the two blastomeres. STUDY DESIGN, SIZE, DURATION Seven transcriptomic studies published between 2011 and 2017 were eligible for retrospective analysis, since both blastomeres of the mouse 2-cell embryo had been analysed individually regarding the original pair associations and since the datasets were made available in public repositories. Five of these studies, encompassing a total of 43 pairs of sister blastomeres, were selected for further analyses based on high interblastomere correlations of mRNA levels. A double cut-off was used to select mRNAs that had robust interblastomere differences both within and between embryos (hits). The hits of each study were compared and contrasted with the hits of the other studies using Venn diagrams. The hits shared by at least four of five studies were analysed further by bioinformatics. PARTICIPANTS/MATERIALS, SETTING, METHODS PubMed was systematically examined for mRNA expression profiles of single 2-cell stage blastomeres in addition to publicly available microarray datasets (GEO, ArrayExpress). Based on the original normalizations, data from seven studies were screened for pairwise sample correlation at the gene level (Spearman), and the top five datasets with the highest correlation were subjected to hierarchical cluster analysis. Interblastomere differences of gene expression were expressed as a ratio of the higher to the lower mRNA level for each pair of blastomeres. A double cut-off was used to make the call of interblastomere difference, accepting genes with mRNA ratios above 2 when observed in at least 50% of the pairs, and discarding the other genes. The proportion of interblastomere differences common to at least four of the five datasets was calculated. Finally, the corresponding gene, pathway and enrichment analyses were performed utilizing PANTHER and GORILLA platforms. MAIN RESULTS AND THE ROLE OF CHANCE An average of 17% of genes within the datasets are differently expressed between sister blastomeres, a proportion which falls to 1% when considering the differences that are common to at least four of the five studies. Housekeeping mRNAs were not included in the 17% and 1% gene lists, suggesting that the interblastomere differences do not occur simply by chance. The 1% of shared interblastomere differences comprise 100 genes, of which 35 are consistent with at least one of the four prevalent models of sister blastomere diversification. Bioinformatics analysis of the remaining 65 genes that are not consistent with the four models suggests that at least one more mechanism is at play, potentially related to the endomembrane system. Although there are many dimensions to the issue of reproducibility (biological, experimental, analytical), we consider that the sister blastomeres are poised to escape high interblastomere correlations of mRNA levels, because at least five sources of diversity superimpose on each other, accounting for at least 25 = 32 different states. As a result, interblastomere mRNA differences of a given 2-cell embryo are necessarily difficult to reproduce in another 2-cell embryo. LARGE SCALE DATA Data were as provided by the original studies (GSE21688, GSE22182, GSE27396, GSE45719, GSE57249, E-MTAB-3321, GSE94050). LIMITATIONS, REASONS FOR CAUTION The original studies present similarities (e.g. fertilization in vivo after ovarian stimulation) as well as differences (e.g. mouse strains, method and timing of blastomere separation). We identified robust mRNA differences between the sister blastomeres, but these differences are underestimated because our double cut-off method works with thresholds and affords more protection against false positives than false negatives. Regarding the false negatives, transcriptome analysis may have captured only part of the interblastomere differences due to: (1) the 2-fold cut-off not being sensitive enough to detect the remaining part of the interblastomere differences, (2) the detection limit of the transcriptomic methods not being sufficient, or (3) interblastomere differences being oblivious to transcriptomic identification because transcriptional changes are oscillatory or because differences are mediated non-transcriptionally or post-transcriptionally. Regarding the false positives, it seems unlikely that a difference was found just by chance for the same group of transcripts due to the same technical error, given that different laboratories produced the data. WIDER IMPLICATIONS OF THE FINDINGS It is clear that the sister blastomeres are distinguishable from each other by mRNA analysis even at the 2-cell stage; however, efforts to identify large stable patterns may be in vain. This elicits thoughts about the wisdom of adding new transcriptomic datasets to the ones that already exist; if all transcriptomic datasets produced so far show a reproducibility of 1%, then any future study would probably face the same issue again. Possibly, a solid identification of the ‘large stable pattern that should be there but was not found’ requires an even larger dataset than the sum of the seven datasets considered here. Conversely, small stable patterns may be easier to identify, but their biological relevance is less obvious. Alternatively, interblastomere differences may not be mediated by nucleic acids but by other cellular components. STUDY FUNDING/COMPETING INTEREST(S) This study was supported by the Deutsche Forschungsgemeinschaft (grant DFG BO 2540-4-3 to M.B. and grant NO 413/3-3 to V.N.). The authors declare that they have no competing financial interests. reproducibility, mouse, 2-cell stage embryo, blastomere, transcriptome, interblastomere differences Introduction A basic issue in mammalian development concerns the status of the 2-cell stage embryo. In contrast to what is known about the early development of many non-mammalian organisms, mechanisms of early mammalian cell differentiation remain controversial regarding the stage at which embryonic blastomeres are demonstrably different from each other and the molecular players involved. The 2-cell stage blastomeres are traditionally viewed as totipotent, i.e. having the individual ability, if separated from the embryo, to form a full organism in a permissive environment. At the same time, being two cells, they also create a formal opportunity of reciprocal diversity. There is not necessarily a conflict between being different and being totipotent, as long as the differences do not restrict the potentiality. The totipotency of a blastomere in mice, for example, is revealed by the number of epiblast cells formed by the blastomere progeny being >4 (Morris et al., 2012). An objection may be put forward that these progenies started out equally, but became different under the influence of the environment, whereby a direct molecular examination of the starting blastomeres is desirable. However, as we are going to show, the published molecular data on 2-cell stage blastomeres in mice present some inconsistencies. This retrospective analysis is devoted to the issue of interblastomere similarities and differences, in gene transcription terms. One line of evidence supports the proposal that there are interblastomere differences in the mammalian 2-cell embryo (Piotrowska et al., 2001; Gardner, 2007) reminiscent of the patterning found in embryos of lower taxa, for example, Caenhorabditis, Drosophila, ascidians, Xenopus or Danio rerio (reviewed in Piotrowska and Zernicka-Goetz, 2002). By contrast, another line of evidence supports the proposal that the 2-cell stage blastomeres in mammals are equal; when kept together in the same embryo, their progenies contribute about half of the blastocyst’s cells (Alarcon and Marikawa, 2005; Waksmundzka et al., 2006), and when the two blastomeres are separated from each other, they are capable of full development as monozygotic twins. However, it became apparent by more careful inspection of the records of monozygotic gestations that the twin embryos formed a single newborn more often than they formed two (Willadsen, 1979; Tsunoda and McLaren, 1983; Allen and Pashen, 1984; Togashi et al., 1987; Matsumoto et al., 1989; Wang et al., 1997; Sotomaru et al., 1998; Mitalipov et al., 2002; Tagawa et al., 2008). Continuing this series, the latest large-scale study of monozygotic mouse embryos (n = 1210) revealed functional differences of cell lineage formation and post-blastocyst development in two-thirds of the twin pairs (Casser et al., 2017). Given the prevalence of these observations, it should be feasible to find gene expression correlates of the recurring functional differences between sister blastomeres. To address this, we focus here on the comparison of transcriptome composition between the sister 2-cell blastomeres in the mouse. The choice of analysing the transcriptome is made because the mRNAs carry the information contained in the gene into the cytoplasm of cells; thus, different populations of mRNAs read from the genome are likely to make a different spectrum of proteins, which are the molecules ultimately responsible for functional traits and, thereby, functional differences. We here ask if interblastomere differences of mRNA composition at the 2-cell stage are frequent and reproducible, and we elaborate on which origin and mechanism as well as biological significance they could have. These questions are instrumental to define when postzygotic differentiation starts, not only as a fundamental step in embryogenesis, but also as metrics to evaluate cellular totipotency. We would expect that the frequency of mRNA differences between sister blastomeres compares well to that of the functional differences revealed by monozygotic twin experiments, which are conspicuous (Casser et al., 2017). We would also expect that the interblastomere differences at the 2-cell stage are more reproducible compared to later stages, since the former do not include exogenous effects accumulated over a longer period of time, for example, effects of the embryo culture environment. Regarding the origin of interblastomere differences, the scientific literature suggests four prevalent origins and, thereby, mechanisms: unequal partition of molecular anisotropies already existing in the mammalian oocyte (Held et al., 2012), preferential distribution of sperm-associated material in one blastomere (Piotrowska and Zernicka-Goetz, 2002), unequal partition of low-abundance mRNAs present in the zygote (Shi et al., 2015) or asynchronous embryonic genome activation (EGA) between the sister blastomeres (Biase et al., 2014). For convenience, we refer to them as models 1, 2, 3 and 4, respectively. We do not know how much of the reality can be explained by these models, i.e. whether interblastomere differences also originate in other ways besides the four mechanisms listed above. This analysis of the sister blastomeres’ transcriptomes may bring clarity into the subject. Although transcriptomic datasets for single blastomeres have been available in mice since 2011, they have not been analysed systematically together, due to various reasons. Firstly, at the beginning, there were simply not enough studies in which both blastomeres of the same embryo had been preserved and analysed side by side. Secondly, a simple metrics appropriate to discern the interblastomere difference was not available. Thirdly and finally, statistical methods that intuitively seem to be suited for analysing series of data pairs, such as the t-test, may not be applied to the blastomeres of a 2-cell embryo: it is not known which blastomere is the ‘number 1’ and which blastomere is the ‘number 2’, making it impossible to build the two groups that are required for the t-test. This problem was recognized years ago (VerMilyea et al., 2011) and still persists. An extension of the t-test into a one-way analysis of variance (ANOVA) would reveal if the total variation originates more from within the 2-cell embryos than from between the 2-cell embryos (Biase et al., 2014). However, when elaborate statistical methods are applied to multiple datasets taken together, the risk of violating the methods’ assumptions increases. A minimal assumption approach has recently been designed to analyse the blastomere transcriptomes in the particular case of the sister blastomeres at the 2-cell stage (Tang et al., 2011; Shi et al., 2015; Zheng et al., 2016). This approach uses the ratio of the gene expression signals (interblastomere mRNA ratio) as simple metrics. In this study, we analysed the single blastomere transcriptomes of seven published datasets to assess the prevalence and reproducibility of transcriptomic differences between sister blastomeres of the mouse 2-cell embryo. With this information in hand, we may be able to uncover patterns of gene expression and to elaborate on their biological significance. Using the interblastomere mRNA ratio as metrics and a double cut-off to make the call, we determined that 17% of genes (total N = 10 050) within datasets are differently expressed between sister blastomeres; however, the proportion of differently expressed genes fell to 1% (n = 100 genes) when crossing the datasets. Thus, the sister blastomeres are broadly distinguishable from each other by mRNA expression, but with little reproducibility. Each blastomere pair seems to be different in its own way. Of the shared 1% interblastomere differences, one-third was assigned, while two-thirds could not be assigned to the four models that have been proposed to explain how interblastomere differences arise. Clearly, it is necessary to invoke additional biological mechanisms than the ones featured in the literature so far to explain the interblastomere variability observed. Regardless of their exact nature, when multiple, for example 5, sources of natural variability superimpose on each other and each of them has a binary outcome (equal vs. unequal blastomeres, coded as ‘1’ and ‘0’), they then form a five-digit binary code, which has 25 = 32 possible states. Thus, the interblastomere differences of a given 2-cell embryo are necessarily difficult to reproduce in another embryo. We propose that while interblastomere mRNA differences clearly exist at the 2-cell stage, the search for underlying mechanisms has not been exhausted with four models; and even if we knew precisely all the mechanisms that are at play, efforts to identify large stable patterns would be in vain, due to the combinatorial effects of known and new mechanisms of blastomere diversification. Despite the fact that large stable transcriptomic patterns of interblastomere difference are elusive, some of the genes that recur more often among the interblastomere differences (Eomes, Tead1, Phlda2, Cops3) play roles in cell lineage specification or proliferation, consistent with the functional observations made in monozygotic twin blastocysts (Casser et al., 2017). Materials and Methods Identification of potentially eligible datasets We systematically mined the PubMed database for RNA expression profiling of single 2-cell stage blastomeres. The following keywords and their combinations were used: mouse, embryo, 2-cell stage, blastomere and transcriptome. In addition, to ensure no relevant studies were missed, we searched publicly available microarray datasets through August 2017 in two public repositories: NCBI Gene Expression Omnibus (GEO) (http://www.ncbi.nlm.nih.gov/geo/) and ArrayExpress database of the European Molecular Biology Laboratory–European Bioinformatics Institute (http://www.ebi.ac.uk/arrayexpress/). The following information was extracted from each study identified: genetic background and number of embryos, method of embryo production after natural ovulation or superovulation, method of embryo bisection and repository accession number (Table I). We inspected the studies manually to ensure that the blastomere pair associations had been respected. Table I Synopsis of the studies considered in this retrospective analysis. Publication Roberts et al. VerMilyea et al. Tang et al. Deng et al. Biase et al. Goolam et al. Casser et al. Repository GSE21688 GSE27396 GSE22182 GSE45719 GSE57249 E-MTAB-3321 GSE94050 Mouse strain CF1 x CF1 B6D2F1 x B6D2F1 MF1 x MF1 CAST/EiJ (CAST) x C57BL/6J (C57) C57BL/6 × C57BL/6 (C57BL/6xCBA) x (C57BL/6xCBA) B6C3F1 x CD1 Stimulated ovarian cycle × ✓ ? ✓ ✓ ✓ ✓ In-vivo fertilization ✓ ✓ ✓ ✓ ✓ ✓ ✓ In-vitro culture in KSOM(aa) ✓ ✓ ✓ Oxygen in vivo ✓ ✓ ✓ ✓ 20% ✓ ✓ ✓ Method of blastomere separation: mechanic ✓ ✓ chemical ✓ ✓ combined ✓ ✓ ✓ Time of blastomere separation About half way through the second cell cycle Within 10 min of the first mitotic division E 1.5 Early, middle, late 2-cell stage 48 h post-hCG late 2-cell (48 h after hCG) 2 h after the first cleavage Blastomere pairs analysed 6 22 4 15 10 8 10 Transcript analysis Agilent-014 868 Whole Mouse Genome Microarray 4×44 K Operon OpArray Mus musculus 38k v4.0.3 SOLiD system Illumina HiSeq 2000 HiSeq 2000 or HiSeq 2500 Illumina HiSeq 2500 Affymetrix GeneChip® Mouse Transcriptome Array 1.0 Detected genes 20498 36008 21118 19754 21310 24129 26597 PMID 21076082 21468028 21731673 24408435 25096407 27015307 28811525 Publication Roberts et al. VerMilyea et al. Tang et al. Deng et al. Biase et al. Goolam et al. Casser et al. Repository GSE21688 GSE27396 GSE22182 GSE45719 GSE57249 E-MTAB-3321 GSE94050 Mouse strain CF1 x CF1 B6D2F1 x B6D2F1 MF1 x MF1 CAST/EiJ (CAST) x C57BL/6J (C57) C57BL/6 × C57BL/6 (C57BL/6xCBA) x (C57BL/6xCBA) B6C3F1 x CD1 Stimulated ovarian cycle × ✓ ? ✓ ✓ ✓ ✓ In-vivo fertilization ✓ ✓ ✓ ✓ ✓ ✓ ✓ In-vitro culture in KSOM(aa) ✓ ✓ ✓ Oxygen in vivo ✓ ✓ ✓ ✓ 20% ✓ ✓ ✓ Method of blastomere separation: mechanic ✓ ✓ chemical ✓ ✓ combined ✓ ✓ ✓ Time of blastomere separation About half way through the second cell cycle Within 10 min of the first mitotic division E 1.5 Early, middle, late 2-cell stage 48 h post-hCG late 2-cell (48 h after hCG) 2 h after the first cleavage Blastomere pairs analysed 6 22 4 15 10 8 10 Transcript analysis Agilent-014 868 Whole Mouse Genome Microarray 4×44 K Operon OpArray Mus musculus 38k v4.0.3 SOLiD system Illumina HiSeq 2000 HiSeq 2000 or HiSeq 2500 Illumina HiSeq 2500 Affymetrix GeneChip® Mouse Transcriptome Array 1.0 Detected genes 20498 36008 21118 19754 21310 24129 26597 PMID 21076082 21468028 21731673 24408435 25096407 27015307 28811525 Table I Synopsis of the studies considered in this retrospective analysis. Publication Roberts et al. VerMilyea et al. Tang et al. Deng et al. Biase et al. Goolam et al. Casser et al. Repository GSE21688 GSE27396 GSE22182 GSE45719 GSE57249 E-MTAB-3321 GSE94050 Mouse strain CF1 x CF1 B6D2F1 x B6D2F1 MF1 x MF1 CAST/EiJ (CAST) x C57BL/6J (C57) C57BL/6 × C57BL/6 (C57BL/6xCBA) x (C57BL/6xCBA) B6C3F1 x CD1 Stimulated ovarian cycle × ✓ ? ✓ ✓ ✓ ✓ In-vivo fertilization ✓ ✓ ✓ ✓ ✓ ✓ ✓ In-vitro culture in KSOM(aa) ✓ ✓ ✓ Oxygen in vivo ✓ ✓ ✓ ✓ 20% ✓ ✓ ✓ Method of blastomere separation: mechanic ✓ ✓ chemical ✓ ✓ combined ✓ ✓ ✓ Time of blastomere separation About half way through the second cell cycle Within 10 min of the first mitotic division E 1.5 Early, middle, late 2-cell stage 48 h post-hCG late 2-cell (48 h after hCG) 2 h after the first cleavage Blastomere pairs analysed 6 22 4 15 10 8 10 Transcript analysis Agilent-014 868 Whole Mouse Genome Microarray 4×44 K Operon OpArray Mus musculus 38k v4.0.3 SOLiD system Illumina HiSeq 2000 HiSeq 2000 or HiSeq 2500 Illumina HiSeq 2500 Affymetrix GeneChip® Mouse Transcriptome Array 1.0 Detected genes 20498 36008 21118 19754 21310 24129 26597 PMID 21076082 21468028 21731673 24408435 25096407 27015307 28811525 Publication Roberts et al. VerMilyea et al. Tang et al. Deng et al. Biase et al. Goolam et al. Casser et al. Repository GSE21688 GSE27396 GSE22182 GSE45719 GSE57249 E-MTAB-3321 GSE94050 Mouse strain CF1 x CF1 B6D2F1 x B6D2F1 MF1 x MF1 CAST/EiJ (CAST) x C57BL/6J (C57) C57BL/6 × C57BL/6 (C57BL/6xCBA) x (C57BL/6xCBA) B6C3F1 x CD1 Stimulated ovarian cycle × ✓ ? ✓ ✓ ✓ ✓ In-vivo fertilization ✓ ✓ ✓ ✓ ✓ ✓ ✓ In-vitro culture in KSOM(aa) ✓ ✓ ✓ Oxygen in vivo ✓ ✓ ✓ ✓ 20% ✓ ✓ ✓ Method of blastomere separation: mechanic ✓ ✓ chemical ✓ ✓ combined ✓ ✓ ✓ Time of blastomere separation About half way through the second cell cycle Within 10 min of the first mitotic division E 1.5 Early, middle, late 2-cell stage 48 h post-hCG late 2-cell (48 h after hCG) 2 h after the first cleavage Blastomere pairs analysed 6 22 4 15 10 8 10 Transcript analysis Agilent-014 868 Whole Mouse Genome Microarray 4×44 K Operon OpArray Mus musculus 38k v4.0.3 SOLiD system Illumina HiSeq 2000 HiSeq 2000 or HiSeq 2500 Illumina HiSeq 2500 Affymetrix GeneChip® Mouse Transcriptome Array 1.0 Detected genes 20498 36008 21118 19754 21310 24129 26597 PMID 21076082 21468028 21731673 24408435 25096407 27015307 28811525 Data citations The datasets analysed in this study (GSE21688, GSE22182, GSE27396, GSE45719, GSE57249, E-MTAB-3321 and GSE94050) are as provided by the original studies. Selection of datasets eligible for further analysis We adopted the normalization performed in the original studies as retrieved from the repository. In the case of microarray data, the values were retrieved from the repository as normalized already (Roberts et al., 2011; Casser et al., 2017). In the case of RNA-Seq, some datasets were retrieved from the repository as ‘reads per kilobase per million mapped reads’ (RPKM)-normalized (Deng et al., 2014), while other datasets were retrieved as counts (Tang et al., 2011; Goolam et al., 2016). Therefore, we chose to apply the RPKM normalization to all RNA-Seq datasets considered in this study. Genes that presented zero counts in all blastomeres were eliminated. Of the remaining genes, some had zero counts in some blastomeres, but not in all. All counts were added to 1 to avoid divisions by zero later in the workflow. These counts were normalized based on the library size (1 million) and transcript length (in kilobases) and used as RPKM without applying threshold filters. We subjected these normalized data to pairwise sample correlation at the gene level (Spearman) and hierarchical cluster analysis (Ward) using the statistical programme JMP® 13.1.0 to confirm, respectively, the expected value of r very close to 1 and the matching of the sister blastomeres. From this point on, analysis was identical for all datasets. Intra-study analysis of interblastomere differences of gene expression We were interested in genes that are consistently differently expressed between the sister blastomeres. Since it is not possible to know a priori which blastomere is ‘first’ and which blastomere is ‘second’, the call of differential gene expression between sister blastomeres was made by dividing the higher value by the lower value across each pair of blastomeres, as described previously (Tang et al., 2011; Shi et al., 2015; Zheng et al., 2016). Briefly, one blastomere will show Ehigher and the other will show Elower (where E represents mRNA expression level) for any mRNA detected, and the ratio Ehigher/Elower represents the adimensional parameter on which to base our analysis. If the sister blastomeres are molecularly equal, then the statistical distribution of Ehigher/Elower will be centred on 1, whereas it will depart from 1 if the sister blastomeres are unequal. Our definition of ‘consistent’ was that 50% or more of the blastomere pairs would feature a ratio >2. Inter-study analysis of interblastomere differences of gene expression We used a Venn diagram tool, InteractiVenn (Heberle et al., 2015), to first identify the transcripts that are common to the five studies, and then compare, within this common domain, the five studies for the mRNAs that present interblastomere differences >2 in at least 50% of the pairs. Comparison of genes featuring interblastomere differences with the genes expressed in oocytes and sperm It has been proposed that when the number of mRNA molecules in the zygote is lower than 180 (36 RPKM; 1 RPKM = 5 molecules), the probability of mRNA molecules being equally partitioned between the daughter blastomeres falls below 80% (Shi et al., 2015). We identified these least-abundant mRNAs’ unfertilized metaphase II oocytes (Huang et al., 2017) and determined the 93rd percentile as the cut-off. The RNA-Seq datasets for metaphase II mouse oocytes were retrieved in Fastq format from dataset GSE 81216 (Huang et al., 2017) and converted to RPKM by Exiqon (https://xplorerna.exiqon.com/). We then analysed the transcript abundances of early, mid and late 2-cell stage blastomeres (Deng et al., 2014) to see which mRNAs fall below the cut-off in blastomeres. We also inspected the blastomere datasets to identify which mRNAs were already present in spermatozoa. The RNA-Seq datasets for mouse spermatozoa were obtained from dataset E-MTAB-5210 (SpermBase; Schuster et al., 2016). Gene ontology and gene set enrichment analyses of interblastomere differences conserved across studies The list of common genes featuring, for all studies, interblastomere differences >2 in at least 50% of the pairs was loaded into PANTHER (Mi et al., 2017), for gene ontology categorization and pathway analysis, and into GORILLA (Eden et al., 2009), for gene set enrichment analysis (GSEA), using this list as the target and the list of common detected transcripts (N = 10 050) as the background. Results Identification of transcriptomic studies featuring individual 2-cell stage blastomeres regarding the original pair associations The goal of this study was to assess the prevalence, reproducibility and biological significance of transcriptomic differences, if any, between sister blastomeres at the mouse 2-cell stage using publicly available datasets. Accordingly, our analysis can only be as good as the original studies in it. Suitable datasets were limited to studies where both blastomeres (1) were isolated with a non-damaging intervention, and (2) were subjected to single-cell transcriptomics, whether by microarray or RNA-Seq, respecting the blastomere pair associations. We identified seven studies (Roberts et al., 2011; Tang et al., 2011; VerMilyea et al., 2011; Biase et al., 2014; Deng et al., 2014; Goolam et al., 2016; Casser et al., 2017) as candidates for inclusion in this retrospective analysis. An eighth study was identified but was not eligible because the gene expression data had not been made available in a public repository (Yang et al., 2017). While there are differences between the seven studies, primarily the genetic background of the embryos and their timing and the transcriptomic platforms used, the strong unifying element lies in the high-throughput single cell analysis of the 2-cell stage blastomeres while preserving the knowledge of the original pair associations (75 pairs). An overview of the studies included here is shown in Table I. Selection of eligible transcriptomic studies for detailed analysis The seven datasets were screened for high correlation of gene expression between the sister blastomeres, since the blastomeres are of the same cell type, derived from the same cell. Even if molecular differences are expected to exist naturally between sister blastomeres (given the immense number of molecules in a cell, the chance of them becoming equally partitioned upon zygotic division is inconceivably small; Elsasser, 1984), substantial departures from a 1:1 ratio should be the exception, resulting in a linear correlation very close to 1. Accordingly, we calculated the Spearman correlation coefficients between the transcriptomes of sister blastomeres, using the normalized transcript levels provided by the studies (Fig. 1A). The studies with high interblastomere correlations were five of the seven, encompassing a total of 43 pairs of sister blastomeres (accession numbers: GSE21688 (Roberts et al., 2011), GSE22182 (Tang et al., 2011), GSE45719 (Deng et al., 2014), E-MTAB-3321 (Goolam et al., 2016), GSE94050 (Casser et al., 2017); Fig. 1). Of note, in one of these five studies (Deng et al., 2014), the sister blastomeres were collected not only respecting the pair associations, but also at different time points during the second cell cycle (early, mid and late). In this particular study, the statistical distribution of the correlation values (GSE45719 in Fig. 1A) was narrow in spite of the fact that the embryos were fertilized at different times and varied among each other regarding the kinetics of mRNA synthesis and degradation (Piko and Clegg, 1982) and the onset of EGA (Nothias et al., 1995). The narrow statistical distribution suggests that either the variability between embryos is unproblematic or, as seems equally possible, calculating the Spearman correlation coefficients between the transcriptomes of sister blastomeres exerted a normalizing effect (see next section). Figure 1 View largeDownload slide Screening of the transcriptomic datasets considered for this retrospective study. (A) Box plot representations of the distributions of the linear correlation coefficient (Spearman) between the transcriptomes of sister blastomeres. Among the seven eligible studies, five studies had high coefficients and two studies had lower coefficients. (B) Clustering of the blastomeres’ transcriptomes from the five studies selected from the seven of this retrospective analysis, using a constellation diagram output. In two studies (datasets GSE22182 (Tang et al., 2011) and E-MTAB-3321 (Goolam et al., 2016)), all 12 sister blastomeres were assigned to the expected pairs, while the blastomeres in the other three studies (datasets GSE21688 (Roberts et al., 2011), GSE45719 (Deng et al., 2014) and GSE94050 (Casser et al., 2017)) were matched with the correct counterpart in 15 of 31 pairs. Circles represent the correct matching of the sister blastomeres. (C1) A first cut-off was applied to interblastomere ratios to retain ratios >2. (C2) A second cut-off was applied to the ratios >2 to retain those genes that occurred in at least 50% of blastomere pairs. Figure 1 View largeDownload slide Screening of the transcriptomic datasets considered for this retrospective study. (A) Box plot representations of the distributions of the linear correlation coefficient (Spearman) between the transcriptomes of sister blastomeres. Among the seven eligible studies, five studies had high coefficients and two studies had lower coefficients. (B) Clustering of the blastomeres’ transcriptomes from the five studies selected from the seven of this retrospective analysis, using a constellation diagram output. In two studies (datasets GSE22182 (Tang et al., 2011) and E-MTAB-3321 (Goolam et al., 2016)), all 12 sister blastomeres were assigned to the expected pairs, while the blastomeres in the other three studies (datasets GSE21688 (Roberts et al., 2011), GSE45719 (Deng et al., 2014) and GSE94050 (Casser et al., 2017)) were matched with the correct counterpart in 15 of 31 pairs. Circles represent the correct matching of the sister blastomeres. (C1) A first cut-off was applied to interblastomere ratios to retain ratios >2. (C2) A second cut-off was applied to the ratios >2 to retain those genes that occurred in at least 50% of blastomere pairs. Secondly, we wanted to determine whether the original pair associations would be recognized again a posteriori. The five datasets with high interblastomere correlation were analysed individually in an unsupervised, hierarchical cluster analysis of the blastomeres’ transcriptomes. It appears that most, but not all, pairs were recognized (see circles in Fig. 1B). If human error (e.g. incorrect sample labelling) is discounted, then not being able to match the sister blastomeres can be interpreted in two opposite ways: (1) it may indicate that the molecular analysis failed to recognize the biological bond of the sister blastomeres due to, for example, technical confounders (RNA amplification bias) or (2) it may just confirm the high quality of the analysis in the presence of a factor that marks one blastomere but not the other, such as the sperm-entry point according to model 2. It should be noted that the example of model 2 does not exclude other models, and that model 2 also requires the absence of a spatial redistribution and overall random localization of sperm-associated materials. In other words, it can be argued both ways: datasets GSE22182 and E-MTAB-3321 give the correct representation because they recognize all pairs, or datasets GSE21688, GSE45719 and GSE94050 give the correct representation because they detect the asymmetry introduced by the sperm entry point. There is also a third possibility, namely, a mixture of equal and unequal pairs due to the presence of functional subpopulations among the 2-cell embryos, as indeed has already been proposed (Casser et al., 2017); the more embryos which are examined, the higher the probability of sampling from both groups. For these reasons, there is no justification for discarding data, and we chose to retain the five datasets for detailed analyses. Within-study identification of genes that are differently expressed between sister blastomeres A convenient approach to analyse the differences in gene expression between sister blastomeres has been described (Tang et al., 2011; Shi et al., 2015; Zheng et al., 2016), in which the higher value of gene expression across each pair is divided by the lower value (ratio Ehigher/Elower) (see Materials and Methods). Calculating the ratio addresses the biological problem of the differences of post-fertilization age between embryos. Embryos can differ from each other in the time elapsed from fertilization, which means that the progression of cellular processes, for example, mRNA synthesis and degradation can vary from embryo to embryo. However, when we apply the interblastomere ratio, time is the same for the numerator and the denominator, which leads de facto to a normalization for time. We consider that the ratio also addresses the technical problem of the differences of the analytical platform among studies. The distribution of the ratios Ehigher/Elower (Fig. 1C1) was adopted here and analysed further. We directed our interest toward the ratios lying convincingly far from 1, looking for interblastomere differences that show up consistently in many pairs. We defined an arbitrary but robust double cut-off criterion to make the call of ‘different’ gene expression between sister blastomeres, namely: >2-fold intra-pair ratio (cut-off 1) observed in ≥50% of the blastomere pairs (cut-off 2). A 2-fold ratio is used as the threshold because it has been shown that mRNAs with ratios up to 1.5 were more frequent in technical replicates than in biological replicates (Shi et al., 2015), so that, after adding a security margin, we suggest that differences of more than 2-fold are no longer capturing noise. A 50% cut-off for the blastomere pairs was chosen because each study contains multiple pairs, and we reasoned that the interblastomere difference should be observed in most of them to qualify as conspicuous. Even though choosing cut-off values can introduce false negatives (e.g. genuine, and yet discarded, 1.5-fold differences), it puts all five datasets on the same footing. The double cut-off criterion also appears to be effective at removing unwanted information; when we examined genes that are expected to show no interblastomere differences, for example, housekeeping genes (Mamo et al., 2007), they were removed by the cut-off (see next paragraph). Overall, our results show that an average of 17 ± 14% of the mRNAs detected satisfy the double cut-off criterion (Fig. 1C2). This does not mean that the substance is essentially the same (17% « 100%), but it means that robust differences unlikely to be false positives (>2-fold in at least 50% of the pairs examined) occur in 17% of the mRNAs detected. Inter-study identification of genes that are differently expressed between sister blastomeres We used a Venn diagram tool to identify genes that are consistently differently expressed between sister blastomeres across studies (consisting of 43 blastomere pairs in total). Because the transcriptomic depths of the datasets are different, it follows that there is a bias in the size of the intersections between the differently expressed genes of the five studies. To address this, the total gene lists of each dataset were uploaded to the InteractiVenn (Heberle et al., 2015) application (http://www.interactivenn.net/), yielding Venn diagrams in Edwards layout (Fig. 2). The Venn diagram intersection of the five transcriptomes resulted in a set of 10 050 mRNAs (Fig. 2A). We used this common set of genes to interrogate each study for the mRNAs that present >2-fold intra-pair difference observed in ≥50% of the blastomere pairs (Fig. 2B). This analysis returned five genes that are shared by all five studies (1700010D01Rik, Fstl5, Pnrc2, Zfp472 and Tead1; Supplementary Table 1). The ratios of these five genes across the 43 pairs of blastomeres (Fig. 2C) were clearly distinct from the ratios of housekeeping genes (Fig. 2D; Supplementary Fig. 1). Venn diagram analysis also returned 95 genes that are shared by four of five studies, for a total of 100 genes (Supplementary Table 1). We chose not to consider genes that were shared by less than four studies, even though this would return more genes, including false positives. The 100 genes were a sufficiently large group for a gene set enrichment analysis (GSEA). Figure 2 View largeDownload slide Analysis of the eligible transcriptomic datasets. (A) Overlap of the transcriptomes of the five studies considered in this retrospective analysis, identifying a shared domain of 10 050 genes. (B) Conservation of the differently expressed mRNAs (ratio >2 in at least 50% of blastomere pairs) in the shared domain of 10 050 genes. (C) According to the double cut-off, five genes were differently expressed in all five studies. (D) Using the same double cut-off, housekeeping genes were classified as not differently expressed. Grey horizontal lines in (C) and (D) represent the overall mean across the five datasets. Figure 2 View largeDownload slide Analysis of the eligible transcriptomic datasets. (A) Overlap of the transcriptomes of the five studies considered in this retrospective analysis, identifying a shared domain of 10 050 genes. (B) Conservation of the differently expressed mRNAs (ratio >2 in at least 50% of blastomere pairs) in the shared domain of 10 050 genes. (C) According to the double cut-off, five genes were differently expressed in all five studies. (D) Using the same double cut-off, housekeeping genes were classified as not differently expressed. Grey horizontal lines in (C) and (D) represent the overall mean across the five datasets. As we mentioned before, the structure of the data is such that the sister blastomere pairs are recognized in some embryos, while they go unrecognized in others (Fig. 1B). There could be reasons to discard either the matched or the unmatched pairs, which we refrained from doing because we considered that this mix reflects biological reality. However, if we hypothesize that the interblastomere differences had been affected by human or technical errors, such as incorrect labelling of the samples or bias of mRNA amplification, then discarding part of the data would have been appropriate. We, therefore, repeated the Venn analysis on a subset of the blastomere pairs. We excluded the 16 pairs in which a blastomere could not be matched with its sister blastomere according to our previous hierarchical cluster analysis of the transcriptomes (Fig. 1B). The remaining 27 pairs, distributed among all five studies, were subjected to Venn analysis, returning four of the five genes (1700010D01Rik, Fstl5, Pnrc2 and Tead1). This means that our conclusions are stable and that the core of reproducible interblastomere differences is very small. We put the emphasis on ‘reproducible’; we do not state that the interblastomere differences are few, but that the robust and reproducible ones are few. A functional annotation and a GSEA were performed using PANTHER (Thomas et al., 2003; Mi et al., 2017) and GORILLA (Eden et al., 2009), respectively, to shed light on the functional roles of the 100 genes that were shared by at least four studies, including the five genes that were shared by all studies (referred to as ‘meta-genes’). In terms of functional annotation, these genes are involved in 13 biological processes and 13 signalling pathways. The integrin signalling is prominent among the 100 genes, followed by EGF, FGF, TGF and WNT signalling. In terms of GSEA, there appear to be few if any themes that interconnect the 100 genes according to GORILLA: significant enrichment is found in GO:0002115 (‘store-operated calcium entry’) and GO:0097084 (‘vascular smooth muscle cell development’). This paucity could be genuine or due to the small enrichment signal that is found when gene sets are small. Understanding interblastomere mRNA differences invokes the combinatorial effects of known and new mechanisms of blastomere diversification Although consistent interblastomere differences of gene expression are not an infallible indicator of importance, genes featuring consistent interblastomere differences do undoubtedly attract attention. We compared and contrasted the 100 genes that were shared by at least four studies with the four models that try to explain the genesis of interblastomere differences at the 2-cell stage (see Introduction and Fig. 3B). Briefly, these differences could arise from unequal partition of molecular anisotropies already existing in the unfertilized oocyte (model 1), from sperm-associated material being retained in one blastomere (model 2), from unequal partition of low-abundance mRNAs present in the zygote (model 3) or from asynchronous EGA between the sister blastomeres (model 4). Figure 3 View largeDownload slide Distribution of ‘cellular component’ GO terms among the 100 genes with interblastomere mRNA differences shared between studies. (A) The 100 genes with shared interblastomere differences were broken down into genes with (n = 35) versus without (n = 65) assignment to (B) the four prevalent models that explain the origin of interblastomere differences at the 2-cell stage. (C) The 35 assigned and the 65 unassigned genes were analysed in the GO ‘cellular component’. Shown are the proportions (%) of GO terms associated with the main cellular structures (ER and Golgi, GO:0005783, GO:0005794; Nucleus, GO:0005634; Organelle, GO:0043227, GO:0043228; Plasma membrane, GO:0005886; Cytoskeleton, GO:0005856). The high number of unassigned genes indicates the necessity to invoke additional mechanisms of interblastomere diversification, potentially related to the endomembrane systems (ER and Golgi). Figure 3 View largeDownload slide Distribution of ‘cellular component’ GO terms among the 100 genes with interblastomere mRNA differences shared between studies. (A) The 100 genes with shared interblastomere differences were broken down into genes with (n = 35) versus without (n = 65) assignment to (B) the four prevalent models that explain the origin of interblastomere differences at the 2-cell stage. (C) The 35 assigned and the 65 unassigned genes were analysed in the GO ‘cellular component’. Shown are the proportions (%) of GO terms associated with the main cellular structures (ER and Golgi, GO:0005783, GO:0005794; Nucleus, GO:0005634; Organelle, GO:0043227, GO:0043228; Plasma membrane, GO:0005886; Cytoskeleton, GO:0005856). The high number of unassigned genes indicates the necessity to invoke additional mechanisms of interblastomere diversification, potentially related to the endomembrane systems (ER and Golgi). Of the 100 gene transcripts, nine were ascribed to model 1, since they were detected in metaphase II oocytes (Huang et al., 2017), but not in mature sperm (SpermBase; Schuster et al., 2016) and since they declined during the 2-cell stage (Deng et al., 2014) consistent with maternal product degradation: Pnrc2, Zfp472, Ppp2r5a, Eomes, Sirt5, Il17rb, Sdc4, Gfra3 and St8sia6 (Supplementary Table 2). Of the 100 gene transcripts, four were ascribed to model 2, since they were detected in mature sperm (Schuster et al., 2016), but not in oocytes (Huang et al., 2017): Gm5941, Rnf32, Nkapl and Stk17b (Supplementary Table 2). In model 3, mRNAs that are least abundant in the zygote are more prone to random partition errors and, thereby, unequal levels in daughter blastomeres (Shi et al., 2015). Since every zygote has a different age post-fertilization, we used the more homogeneous unfertilized MII oocytes (Huang et al., 2017) in lieu of zygotes to define a low threshold (see Materials and methods). We found that 13 of the 100 gene transcripts fall below this threshold and, thus, are consistent with model 3: Pnrc2, Eomes, 1810030O07Rik, Rras, Tcf19, Parp14, 1700029F12Rik, Pmvk, Ttyh1, Cxcl16, Osr2, Zfp566 and Wrap53 (Supplementary Table 2). In model 4, interblastomere differences arise as a result of asynchronous EGA producing unequal transcript accumulation in the sister blastomeres. We identified 50 mRNAs that increase during the 2-cell stage (Deng et al., 2014): 1810030O07Rik, Rras, Tcf19, Parp14, 1700029F12Rik, Ttyh1, Cxcl16, Zfp566, Wrap53, gm5941, Nkapl, Stk17b, 1700010D01Rik, Cxcl11, Lgals1, Tead1, Hist1h4i, Zswim6, Tfpi, Zfp369, Far1, Bcs1l, Slc25a17, Akr1b8, Bet1l, Hmgn3, Lin52, Spink2, Rwdd2a, Wdr25, H2-K1, Spata24, Ramp2, Tbc1d22a, Stim1, Hmg20b, Hes1, Zfp81, Phlda2, Serinc5, Pls3, Nmb, Ids, Dusp6, Rabggtb, Actr6, Cdkl5, Lage3, Cops3 and Mrps18c; of these 50 mRNAs, 17 are also sensitive to embryo treatment with alpha-amanitin (Zeng and Schultz, 2005), consistent with model 4 (Supplementary Table 2). The genes listed above are too many for individual analysis, but we note that four of them, Eomes, Tead1, Phlda2 and Cops3, are relevant for cell lineage specification or proliferation (see Discussion). Adding up the numbers of genes that are consistent with the models returns a total of 9 + 4 + 13 + 17 = 43, which is lower than 100. Naturally, we expect that the models above are not mutually exclusive, i.e. the same gene may be featured in more than one model. This proved to be the case with 1700010D01Rik, for example, one of the five ‘meta-genes’ that are consistently differently expressed between sister blastomeres of all studies. 1700010D01Rik could be associated with spermatogenesis (Hong et al., 2005; Ellis et al., 2011), consistent with model 2; however, it is also alpha-amanitin-sensitive in mouse zygotes (Zeng and Schultz, 2005), consistent with model 4. Due to such overlaps, the number of 43 genes is reduced to 35 (Fig. 3A), which makes the gap with the 100 genes even wider. The remaining 65 genes could not be assigned to any of the four models, suggesting that at least one additional mechanism of sister blastomere diversification must be in play. We reasoned that a GO analysis of the 65 unassigned genes could point us at the additional mechanism/model. While not reaching statistical significance, GO analysis in the cellular component (GO CC) revealed a three-times higher proportion of terms related to the ‘endoplasmic reticulum’ (GO:0005783) and ‘Golgi apparatus’ (GO:0005794) in the 65 unassigned genes compared to the 35 assigned genes (Fig. 3C). This is consistent with the previous GO analysis conducted in the ‘biological process’ (GP BP) on the 100 shared genes, which returned a term related to Calcium ion channels that are found on the endoplasmic reticulum and on the plasma membrane (Prakriya and Lewis, 2015). Of note, the genes featured in the GO BP analysis of the 100 genes coincide with those featured in the GO CC analysis of the 65 genes. In summary, it is necessary to invoke additional biological mechanisms than the ones featured in the four models to understand the interblastomere variability observed, and the endomembrane system stands as a candidate. Discussion The sister blastomeres of the mammalian 2-cell stage embryo are traditionally considered to be totipotent, but there is mounting experimental evidence that they are not equally able to form monozygotic twins. This is a tricky situation in developmental biology, because it gives contradictory answers to the basic question: ‘At which stage do embryonic blastomeres start to differentiate from each other?’, with implications for the reproducibility of scientific results (Munafo et al., 2017). To investigate this ambiguity, we conducted a retrospective analysis of single-cell datasets, assessing the prevalence, reproducibility and biological significance of transcriptomic differences between sister blastomeres of the mouse 2-cell embryo. As we are going to discuss, cell diversification has clearly occurred at the 2-cell stage, but its pattern is necessarily difficult to reproduce at the mRNA level, because there are multiple sources of interblastomere variability acting in various combinations. The evidence of functional differences between sister 2-cell stage blastomeres has been out there for a long time (Willadsen, 1979; Tsunoda and McLaren, 1983; Allen and Pashen, 1984; Togashi et al., 1987; Matsumoto et al., 1989; Wang et al., 1997; Sotomaru et al., 1998; Mitalipov et al., 2002; Tagawa et al., 2008) and recently it grew even stronger owing to a study that featured more than one thousand monozygotic twins under one and the same set of laboratory conditions (Casser et al., 2017). However, the molecular foundations of these functional differences are elusive, and our present analysis of the transcriptomes of sister blastomeres did not find a large univocal pattern of interblastomere differences based on existing models. Our search was complicated by several factors which could account for at least some variation in the datasets. The mouse strains featured in the original studies were multiple (CF1 x CF1; B6D2F1 x B6D2F1; MF1 x MF1; CAST/EiJ (CAST) x C57Bl/6J; C57Bl/6J x C57Bl/6J; (C57Bl/6J x CBA) x (C57Bl/6J x CBA); B6C3F1 x CD1). Additionally, the embryos were sampled at different time points (the time points were defined more accurately only in VerMilyea et al., 2011, and in Deng et al., 2014). It is tempting to propose that the original studies were perhaps not homogeneous enough, and that a search for a blastomere pattern would be better served by using one and the same biological system. However, the more you standardize, the less you can generalize, which is of value in biology. We also acknowledge that part of the difference may lie below the detection limit of the transcriptomic methods, and that our method of analysis is more likely to produce false negatives than false positives (for example, differences of 1.8-fold or differences observed in 40% of the pairs are discarded even if perhaps they are ‘true’). However, these criteria are justified, because, for example, differences of up to 1.5-fold are often due to noise (Shi et al., 2015) and not accepting a difference unless it occurs in most of the blastomere pairs is a way to protect against false positives. Therefore, our analysis is likely to provide an underestimation of the mRNA differences between the sister blastomeres. In addition, there could be differences that are oblivious to transcriptomic identification because they are mediated post-transcriptionally or non-transcriptionally. Even though the aforementioned confounders are present, we note that interblastomere differences are still observed even in the presence of a superior experimental design (Goolam et al., 2016), in which a spike control was used, whereas the behaviour of the housekeeping genes was very similar between the sister blastomeres. In our opinion, these considerations point at biological, not technical, origins for the interblastomere differences. Following on from these considerations, three major questions are: (1) how could the interblastomere differences arise biologically, (2) why were they difficult to reproduce, and (3) what functional implications could they have (biological significance)? Regarding the origin of interblastomere mRNA differences at the 2-cell stage, a most simple explanation that would ascribe them to cell cycle asynchrony between sister blastomeres seems unlikely, since key cell cycle genes, such as Atr, Chk1 and Wee1 (Artus and Cohen-Tannoudji, 2008), were not featured among the 100 genes of our analysis. This is consistent with the report which found that DNA methylation dynamics between the blastomeres within the same 2-cell embryo is highly synchronized (Guo et al., 2017). While the scatter of global gene expression widened from the early-to-mid and from the mid-to-late 2-cell stage, the intra-stage correlation remained very high (Piras et al., 2014), supporting the existence of a highly synchronized cell cycle in the blastomeres within the same embryo. Aside from cell cycle asynchrony, we have grouped the other possibilities into four models for simplicity and we have compared and contrasted the 100 hits of our analysis with these models. Out of 100 genes, 35 are consistent with at least one of the models (the remaining 65 genes are consistent with none). Model 4 (EGA) is the one that received the most support from the current set of data, as well as support from another recent study (Biase et al., 2014), followed by model 3 (partition error of low-copy-number mRNAs), model 1 (maternal products) and model 2 (sperm factors) (Supplementary Table 2). It should be noted that EGA occurs at the 2-cell stage in mice but later in other mammalian species, which limits the general validity of our analysis, and it remains to be explained why one blastomere has an edge over the other when it comes to the EGA. A recent hypothesis envisions that differential histone or DNA methylation marks are imposed on the maternal and paternal chromosomes during gametogenesis, and these marks segregate at the first mitosis of the zygote, giving rise to two epigenetically distinguishable sister blastomeres that support transcription differently (Takaoka and Hamada, 2014). This would be earlier than the stage when the sister blastomeres can become different due to allelic segregation (Ito et al., 1988). Additionally, there could be differences in the polyadenylation and differential degradation of maternal mRNAs leading to a faster or slower transition into the embryonic programme of gene expression (Svoboda et al., 2015; Sousa Martins et al., 2016; Sha et al., 2017). The length of the polyA tail can influence the efficiency of reverse transcription during the initial steps of the transcriptomic analysis and, thereby, the perceived (as opposed to the actual) transcript abundance. It cannot be deduced from the studies whether the different transcriptomic signal is due to the number of mRNA molecules, which varies substantially during preimplantation development and features a minimum at the 2-cell stage (Fan et al., 2015) or to their polyA tail lengths. The sole way to solve this conundrum would be to apply PAL sequencing analysis (Subtelny et al., 2014), which was not performed in any of the studies considered here. In contrast to model 4, model 2 is the one that has the least support from the current set of data. In spite of its simplicity, the possibility of a higher concentration of former sperm-associated mRNAs in one blastomere seems unrealistic; a mouse sperm head has a length and width of 7.9 and 3.2 μm, respectively (Cummins and Woodall, 1985) and, thus, sperm content is diluted more than 2000 times in the mouse ooplasm. The question remains whether the sperm entry point may leave a mark that is not decoded at the 2-cell stage, but later. The fact that 65 of the 100 genes could not be assigned to the four models suggests that the repertoire of possible mechanisms is not exhausted with four. As a fifth model, the comparison of the 35 assigned genes with the 65 unassigned genes suggests a role for the partition of membranous systems between the sister blastomeres. These systems could be those of the endoplasmic reticulum and Golgi apparatus, consistent with older investigations (Maro et al., 1985). Of note, the endoplasmic reticulum was featured in a genetic mouse model of unequal blastomere partition, in which the sister blastomeres differ in size by 10% at the 2-cell stage (Gao et al., 2018). In turn, the endoplasmic reticulum mediates Calcium ion homoeostasis (Kim et al., 2014) including Calcium oscillations that may regulate further progression into embryogenesis. Regarding the modest reproducibility of interblastomere mRNA differences at the 2-cell stage, we consider that the aforementioned multiple origins (at least four plus one) also provide a key to address this aspect. If the various natural origins superimpose on each other and occur in any number of combinations, then the interblastomere differences of a given 2-cell embryo are necessarily difficult to reproduce in another 2-cell embryo. In a very simplified scenario, should each of the multiple origins unfold independently of the others and give rise to binary states (equal vs. unequal blastomeres), this would result in 25 = 32 combined states, and each state would, thus, have a mere 3% (100/32) probability of being reproduced exactly in another 2-cell embryo. The idea of a combinatorial model for blastomere diversification is not new. It was proposed in Caenhorabditis, for example, that at least three events can control the specification of cells at the 8-cell stage, where each event appears to be an independent regulatory input that controls a binary decision. Therefore, these three regulatory inputs act in various combinations and each unique combination appears to specify a unique cell identity (Moskowitz et al., 1994). Finally, the modest reproducibility of interblastomere mRNA differences at the 2-cell stage could also depend on transcriptional oscillations (Kiessling et al., 2010; Tang et al., 2011), which are a very poorly understood phenomenon. One of the 100 genes, Hes1, encodes a canonical effector of Notch signalling and features oscillating expression with a period of about two hours in various cell types including neural progenitor cells (Baek et al., 2006) and embryonic stem cells (Kobayashi and Kageyama, 2011). It is here a matter of speculation, whether transcriptional oscillations of Hes1 and other genes take place in 2-cell stage blastomeres, or not. Should transcriptional oscillations exist, and should they influence the expression of target genes in blastomeres, this could decrease the reproducibility of interblastomere differences even further. As to the functional implications of interblastomere mRNA differences at the 2-cell stage, the exact consequences of interblastomere differences are difficult to predict, but certain effects seem likely, even with the lack of consensus about the possible causes. Four among the 100 genes, for example, are relevant for cell lineage specification or proliferation (Eomes, Tead1, Phlda2 and Cops3). Eomes is required for mouse trophoblast development and mesoderm formation (Russ et al., 2000). In Danio rerio, Eomes is required for endoderm induction and is non-uniformly localized in the 1–2 cell stage embryo (Bjornson et al., 2005). Tead1 is expressed in both mouse trophectoderm and inner cell mass (Nishioka et al., 2008), where it interacts with Yap (Sawada et al., 2008), suggesting common roles in Hippo signalling. Phlda2 is expressed in mouse trophectoderm and extraembryonic ectoderm, where it modulates placental growth (Frank et al., 2002; Salas et al., 2004). Cops3 is crucial for the maintenance of cell proliferation in the mouse epiblast (Yan et al., 2003), which was the main contributor to the phenotypic discordance of monozygotic twin blastocysts (Casser et al., 2017). It is therefore tempting to envision that the interblastomere differences (whatever the cause) might influence the fate of daughter cells at the blastocyst stage (Casser et al., 2017). Less is known about the other four genes that, together with Tead1, have interblastomere differences of mRNA expression that are shared by all five studies: 1700010D01Rik, Fstl5, Pnrc2 and Zfp472. We already noted the behaviour of 1700010D01Rik, which is associated with spermatogenesis (Hong et al., 2005; Ellis et al., 2011) and is also alpha-amanitin-sensitive in mouse zygotes (Zeng and Schultz, 2005). Fstl5 may encode a calcium ion-binding protein which has been detected in human oocytes (Virant-Klun et al., 2016). Pnrc2 is involved in regulation of transcription, with evidence for nonsense-mediated mRNA decay (Cho et al., 2009). Zfp472 (Krim-1) codes for a KRAB box containing zinc finger protein that interacts with c-Myc (Hennemann et al., 2003). In conclusion, our analysis provides a compass for those who are interested in the earliest phase of embryo development and how the cells start to become different from each other. Even if mRNA differences between the sister blastomeres are underestimated by our double cut-off method, our results leave little doubt that the sister blastomeres are already distinguishable from each other at the 2-cell stage; however, efforts to identify large stable patterns based on mRNA may be in vain, due to combinatorial contributions to diversity. This elicits thoughts about the wisdom of adding new transcriptomic datasets to the ones already available, with the aim of exposing systematic patterns in sister 2-cell blastomeres; if all transcriptomic datasets to date show a reproducibility of 1%, which is comparable to the 3% predicted on the basis of simple combinatorial considerations, then any new study would probably face the same problem. Possibly, a solid identification of the ‘large stable pattern that should be there but was not found’ requires an even larger dataset than the sum of the seven datasets considered here. Conversely, small stable patterns may be easier to identify, but their biological relevance is less obvious. Alternatively, interblastomere differences may be mediated not only by mRNAs, i.e. nucleic acids, but also by other cellular components, for example, proteins, lipids or other types of molecules. It is especially desirable to learn whether the interblastomere mRNA differences project on the protein level. This would have to be tested in the rather cumbersome way of immunocytochemistry, since proteomics is currently not sensitive enough regarding the minimum amount of cellular material present in a single germ cell (Virant-Klun et al., 2016). Although there are multiple dimensions to the issue of interblastomere diversity (biological, experimental, analytical), we propose that the sister blastomeres are poised to escape reproducibility of their molecular make-up, and that the mRNA dimension may be insufficient to study how the blastomeres start to become different from each other. This has implications for the validity of using the molecular make-up of one blastomere to infer what made the companion blastomere developmentally successful or unsuccessful (Held et al., 2012), placing constraints on the usefulness of embryo biopsy-based strategies to study what regulates developmental competence. Supplementary data Supplementary data are available at Molecular Human Reproduction online. Acknowledgements We thank Andreas Huge who gave us precious advice on the analysis of transcriptomic data. K. John McLaughlin, Ursula Eichenlaub-Ritter, Hans-Werner Denker and Georg Fuellen commented on earlier versions of the manuscript. Amy Pavlak and Philip Saunders checked the manuscript for grammar. Last but not least, we thank the Max Planck Institute for Molecular Biomedicine and its Director, Prof. Hans R. Schöler, for infrastructural support. Authors’ roles M.B., V.N. and S.S. conceived the study. M.B., E.C. and S.I. analysed the data. M.B. wrote the initial draft. M.B. and E.C. wrote the manuscript in its final form with feedback from all authors. All authors approved the final version. Funding This study was supported by the Deutsche Forschungsgemeinschaft (grant DFG BO 2540-4-3 to M.B. and grant NO 413/3-3 to V.N.). Conflict of interest The authors declare that they are not aware of any conflict of interest that may affect the study. References Alarcon VB , Marikawa Y . Unbiased contribution of the first two blastomeres to mouse blastocyst development . Mol Reprod Dev 2005 ; 72 : 354 – 361 . Allen WR , Pashen RL . Production of monozygotic (identical) horse twins by embryo micromanipulation . J Reprod Fertil 1984 ; 71 : 607 – 613 . Artus J , Cohen-Tannoudji M . Cell cycle regulation during early mouse embryogenesis . Mol Cell Endocrinol 2008 ; 282 : 78 – 86 . Baek JH , Hatakeyama J , Sakamoto S , Ohtsuka T , Kageyama R . Persistent and high levels of Hes1 expression regulate boundary formation in the developing central nervous system . Development 2006 ; 133 : 2467 – 2476 . Biase FH , Cao X , Zhong S . Cell fate inclination within 2-cell and 4-cell mouse embryos revealed by single-cell RNA sequencing . Genome Res 2014 ; 24 : 1787 – 1796 . Bjornson CR , Griffin KJ , Farr GH III , Terashima A , Himeda C , Kikuchi Y , Kimelman D . Eomesodermin is a localized maternal determinant required for endoderm induction in zebrafish . Dev Cell 2005 ; 9 : 523 – 533 . Casser E , Israel S , Witten A , Schulte K , Schlatt S , Nordhoff V , Boiani M . Totipotency segregates between the sister blastomeres of two-cell stage mouse embryos . Sci Rep 2017 ; 7 : 8299 . Cho H , Kim KM , Kim YK . Human proline-rich nuclear receptor coregulatory protein 2 mediates an interaction between mRNA surveillance machinery and decapping complex . Mol Cell 2009 ; 33 : 75 – 86 . Cummins JM , Woodall PF . On mammalian sperm dimensions . J Reprod Fertil 1985 ; 75 : 153 – 175 . Deng Q , Ramskold D , Reinius B , Sandberg R . Single-cell RNA-seq reveals dynamic, random monoallelic gene expression in mammalian cells . Science 2014 ; 343 : 193 – 196 . Eden E , Navon R , Steinfeld I , Lipson D , Yakhini Z . GOrilla: a tool for discovery and visualization of enriched GO terms in ranked gene lists . BMC Bioinformatics 2009 ; 10 : 48 . Ellis PJ , Yu Y , Zhang S . Transcriptional dynamics of the sex chromosomes and the search for offspring sex-specific antigens in sperm . Reproduction 2011 ; 142 : 609 – 619 . Elsasser WM . Outline of a theory of cellular heterogeneity . Proc Natl Acad Sci USA 1984 ; 81 : 5126 – 5129 . Fan X , Zhang X , Wu X , Guo H , Hu Y , Tang F , Huang Y . Single-cell RNA-seq transcriptome analysis of linear and circular RNAs in mouse preimplantation embryos . Genome Biol 2015 ; 16 : 148 . Frank D , Fortino W , Clark L , Musalo R , Wang W , Saxena A , Li CM , Reik W , Ludwig T , Tycko B . Placental overgrowth in mice lacking the imprinted gene Ipl . Proc Natl Acad Sci U S A 2002 ; 99 : 7490 – 7495 . Gao Z , Zhang X , Yu X , Qin D , Xiao Y , Yu Y , Xiang Y , Nie X , Lu X , Liu W et al. . Zbed3 participates in the subcortical maternal complex and regulates the distribution of organelles . J Mol Cell Biol 2018 ; 10 : 74 – 88 . Gardner RL . The axis of polarity of the mouse blastocyst is specified before blastulation and independently of the zona pellucida . Hum Reprod 2007 ; 22 : 798 – 806 . Goolam M , Scialdone A , Graham SJL , Macaulay IC , Jedrusik A , Hupalowska A , Voet T , Marioni JC , Zernicka-Goetz M . Heterogeneity in Oct4 and Sox2 targets biases cell fate in 4-cell mouse embryos . Cell 2016 ; 165 : 61 – 74 . Guo F , Li L , Li J , Wu X , Hu B , Zhu P , Wen L , Tang F . Single-cell multi-omics sequencing of mouse early embryos and embryonic stem cells . Cell Res 2017 ; 27 : 967 – 988 . Heberle H , Meirelles GV , da Silva FR , Telles GP , Minghim R . InteractiVenn: a web-based tool for the analysis of sets through Venn diagrams . BMC Bioinformatics 2015 ; 16 : 169 . Held E , Salilew-Wondim D , Linke M , Zechner U , Rings F , Tesfaye D , Schellander K , Hoelker M . Transcriptome fingerprint of bovine 2-cell stage blastomeres is directly correlated with the individual developmental competence of the corresponding sister blastomere . Biol Reprod 2012 ; 87 : 154 . Hennemann H , Vassen L , Geisen C , Eilers M , Moroy T . Identification of a novel Kruppel-associated box domain protein, Krim-1, that interacts with c-Myc and inhibits its oncogenic activity . J Biol Chem 2003 ; 278 : 28799 – 28811 . Hong S , Choi I , Woo JM , Oh J , Kim T , Choi E , Kim TW , Jung YK , Kim DH , Sun CH et al. . Identification and integrative analysis of 28 novel genes specifically expressed and developmentally regulated in murine spermatogenic cells . J Biol Chem 2005 ; 280 : 7685 – 7693 . Huang Y , Kim JK , Do DV , Lee C , Penfold CA , Zylicz JJ , Marioni JC , Hackett JA , Surani MA . Stella modulates transcriptional and endogenous retrovirus programs during maternal-to-zygotic transition . Elife 2017 ; 6 . doi:10.7554/eLife.22345 . Ito K , McGhee JD , Schultz GA . Paternal DNA strands segregate to both trophectoderm and inner cell mass of the developing mouse embryo . Genes Dev 1988 ; 2 : 929 – 936 . Kiessling AA , Bletsa R , Desmarais B , Mara C , Kallianidis K , Loutradis D . Genome-wide microarray evidence that 8-cell human blastomeres over-express cell cycle drivers and under-express checkpoints . J Assist Reprod Genet 2010 ; 27 : 265 – 276 . Kim B , Zhang X , Kan R , Cohen R , Mukai C , Travis AJ , Coonrod SA . The role of MATER in endoplasmic reticulum distribution and calcium homeostasis in mouse oocytes . Dev Biol 2014 ; 386 : 331 – 339 . Kobayashi T , Kageyama R . Hes1 oscillations contribute to heterogeneous differentiation responses in embryonic stem cells . Genes (Basel) 2011 ; 2 : 219 – 228 . Mamo S , Gal AB , Bodo S , Dinnyes A . Quantitative evaluation and selection of reference genes in mouse oocytes and embryos cultured in vivo and in vitro . BMC Dev Biol 2007 ; 7 : 14 . Maro B , Johnson MH , Pickering SJ , Louvard D . Changes in the distribution of membranous organelles during mouse early development . J Embryol Exp Morphol 1985 ; 90 : 287 – 309 . Matsumoto K , Miyake M , Utsumi K , Iritani A . Production of identical twins by separating two-cell rat embryos . Gamete Res 1989 ; 22 : 257 – 263 . Mi H , Huang X , Muruganujan A , Tang H , Mills C , Kang D , Thomas PD . PANTHER version 11: expanded annotation data from Gene Ontology and Reactome pathways, and data analysis tool enhancements . Nucleic Acids Res 2017 ; 45 : D183 – D189 . Mitalipov SM , Yeoman RR , Kuo HC , Wolf DP . Monozygotic twinning in rhesus monkeys by manipulation of in vitro-derived embryos . Biol Reprod 2002 ; 66 : 1449 – 1455 . Morris SA , Guo Y , Zernicka-Goetz M . Developmental plasticity is bound by pluripotency and the Fgf and Wnt signaling pathways . Cell Rep 2012 ; 2 : 756 – 765 . Moskowitz IP , Gendreau SB , Rothman JH . Combinatorial specification of blastomere identity by glp-1-dependent cellular interactions in the nematode Caenorhabditis elegans . Development 1994 ; 120 : 3325 – 3338 . Munafo MR , Nosek BA , Bishop DVM , Button KS , Chambers CD , Percie du Sert N , Simonsohn U , Wagenmakers E-J , Ware JJ , Ioannidis JPA . A manifesto for reproducible science . Nat Hum Behav 2017 ; 1 : 0021 EP . Nishioka N , Yamamoto S , Kiyonari H , Sato H , Sawada A , Ota M , Nakao K , Sasaki H . Tead4 is required for specification of trophectoderm in pre-implantation mouse embryos . Mech Dev 2008 ; 125 : 270 – 283 . Nothias JY , Majumder S , Kaneko KJ , DePamphilis ML . Regulation of gene expression at the beginning of mammalian development . J Biol Chem 1995 ; 270 : 22077 – 22080 . Piko L , Clegg KB . Quantitative changes in total RNA, total poly(A), and ribosomes in early mouse embryos . Dev Biol 1982 ; 89 : 362 – 378 . Piotrowska K , Wianny F , Pedersen RA , Zernicka-Goetz M . Blastomeres arising from the first cleavage division have distinguishable fates in normal mouse development . Development 2001 ; 128 : 3739 – 3748 . Piotrowska K , Zernicka-Goetz M . Early patterning of the mouse embryo – contributions of sperm and egg . Development 2002 ; 129 : 5803 – 5813 . Piras V , Tomita M , Selvarajoo K . Transcriptome-wide variability in single embryonic development cells . Sci Rep 2014 ; 4 : 7137 . Prakriya M , Lewis RS . Store-operated calcium channels . Physiol Rev 2015 ; 95 : 1383 – 1436 . Roberts RM , Katayama M , Magnuson SR , Falduto MT , Torres KE . Transcript profiling of individual twin blastomeres derived by splitting two-cell stage murine embryos . Biol Reprod 2011 ; 84 : 487 – 494 . Russ AP , Wattler S , Colledge WH , Aparicio SA , Carlton MB , Pearce JJ , Barton SC , Surani MA , Ryan K , Nehls MC et al. . Eomesodermin is required for mouse trophoblast development and mesoderm formation . Nature 2000 ; 404 : 95 – 99 . Salas M , John R , Saxena A , Barton S , Frank D , Fitzpatrick G , Higgins MJ , Tycko B . Placental growth retardation due to loss of imprinting of Phlda2 . Mech Dev 2004 ; 121 : 1199 – 1210 . Sawada A , Kiyonari H , Ukita K , Nishioka N , Imuta Y , Sasaki H . Redundant roles of Tead1 and Tead2 in notochord development and the regulation of cell proliferation and survival . Mol Cell Biol 2008 ; 28 : 3177 – 3189 . Schuster A , Tang C , Xie Y , Ortogero N , Yuan S , Yan W . SpermBase: a database for sperm-borne RNA contents . Biol Reprod 2016 ; 95 : 99 . Sha QQ , Dai XX , Dang Y , Tang F , Liu J , Zhang YL , Fan HY . A MAPK cascade couples maternal mRNA translation and degradation to meiotic cell cycle progression in mouse oocytes . Development 2017 ; 144 : 452 – 463 . Shi J , Chen Q , Li X , Zheng X , Zhang Y , Qiao J , Tang F , Tao Y , Zhou Q , Duan E . Dynamic transcriptional symmetry-breaking in pre-implantation mammalian embryo development revealed by single-cell RNA-seq . Development 2015 ; 142 : 3468 – 3477 . Sotomaru Y , Kato Y , Tsunoda Y . Production of monozygotic twins after freezing and thawing of bisected mouse embryos . Cryobiology 1998 ; 37 : 139 – 145 . Sousa Martins JP , Liu X , Oke A , Arora R , Franciosi F , Viville S , Laird DJ , Fung JC , Conti M . DAZL and CPEB1 regulate mRNA translation synergistically during oocyte maturation . J Cell Sci 2016 ; 129 : 1271 – 1282 . Subtelny AO , Eichhorn SW , Chen GR , Sive H , Bartel DP . Poly(A)-tail profiling reveals an embryonic switch in translational control . Nature 2014 ; 508 : 66 – 71 . Svoboda P , Franke V , Schultz RM . Sculpting the transcriptome during the oocyte-to-embryo transition in mouse . Curr Top Dev Biol 2015 ; 113 : 305 – 349 . Tagawa M , Matoba S , Narita M , Saito N , Nagai T , Imai K . Production of monozygotic twin calves using the blastomere separation technique and Well of the Well culture system . Theriogenology 2008 ; 69 : 574 – 582 . Takaoka K , Hamada H . Origin of cellular asymmetries in the pre-implantation mouse embryo: a hypothesis . Philos Trans R Soc Lond B 2014 ; 369 :20130536. doi:10.1098/rstb.2013.0536 Tang F , Barbacioru C , Nordman E , Bao S , Lee C , Wang X , Tuch BB , Heard E , Lao K , Surani MA . Deterministic and stochastic allele specific gene expression in single mouse blastomeres . PLoS ONE 2011 ; 6 : e21208 . Thomas PD , Campbell MJ , Kejariwal A , Mi H , Karlak B , Daverman R , Diemer K , Muruganujan A , Narechania A . PANTHER: a library of protein families and subfamilies indexed by function . Genome Res 2003 ; 13 : 2129 – 2141 . Togashi M , Suzuki H , Miyai T , Okamoto MT . Production of monozygotic twins by splitting of 2-cell stage embryos in mice . Jpn J Anim Reprod 1987 ; 33 : 51 – 57 . Tsunoda Y , McLaren A . Effect of various procedures on the viability of mouse embryos containing half the normal number of blastomeres . J Reprod Fertil 1983 ; 69 : 315 – 322 . VerMilyea MD , Maneck M , Yoshida N , Blochberger I , Suzuki E , Suzuki T , Spang R , Klein CA , Perry AC . Transcriptome asymmetry within mouse zygotes but not between early embryonic sister blastomeres . EMBO J 2011 ; 30 : 1841 – 1851 . Virant-Klun I , Leicht S , Hughes C , Krijgsveld J . Identification of maturation-specific proteins by single-cell proteomics of human oocytes . Mol Cell Proteomics 2016 ; 15 : 2616 – 2627 . Waksmundzka M , Wisniewska A , Maleszewski M . Allocation of cells in mouse blastocyst is not determined by the order of cleavage of the first two blastomeres . Biol Reprod 2006 ; 75 : 582 – 587 . Wang M , Kato Y , Tsunoda Y . Effects of several factors on the monozygotic twin production in the mouse . J Reprod Dev 1997 ; 43 : 91 – 95 . Willadsen SM . A method for culture of micromanipulated sheep embryos and its use to produce monozygotic twins . Nature 1979 ; 277 : 298 – 300 . Yan J , Walz K , Nakamura H , Carattini-Rivera S , Zhao Q , Vogel H , Wei N , Justice MJ , Bradley A , Lupski JR . COP9 signalosome subunit 3 is essential for maintenance of cell proliferation in the mouse embryonic epiblast . Mol Cell Biol 2003 ; 23 : 6798 – 6808 . Yang L , Ma Z , Cao C , Zhang Y , Wu X , Lee R , Hu B , Wen L , Ge H , Huang Y et al. . MR-seq: measuring a single cell’s transcriptome repeatedly by RNA-seq . Sci Bull 2017 ; 62 : 391 – 398 . Zeng F , Schultz RM . RNA transcript profiling during zygotic gene activation in the preimplantation mouse embryo . Dev Biol 2005 ; 283 : 40 – 57 . Zheng Z , Li H , Zhang Q , Yang L , Qi H . Unequal distribution of 16S mtrRNA at the 2-cell stage regulates cell lineage allocations in mouse embryos . Reproduction 2016 ; 151 : 351 – 367 . © The Author 2018. Published by Oxford University Press on behalf of the European Society of Human Reproduction and Embryology. All rights reserved. For Permissions, please email: email@example.com This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/about_us/legal/notices)
Molecular Human Reproduction – Oxford University Press
Published: May 9, 2018
It’s your single place to instantly
discover and read the research
that matters to you.
Enjoy affordable access to
over 18 million articles from more than
15,000 peer-reviewed journals.
All for just $49/month
Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly
Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place.
Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals.
Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more.
All the latest content is available, no embargo periods.
“Hi guys, I cannot tell you how much I love this resource. Incredible. I really believe you've hit the nail on the head with this site in regards to solving the research-purchase issue.”Daniel C.
“Whoa! It’s like Spotify but for academic articles.”@Phil_Robichaud
“I must say, @deepdyve is a fabulous solution to the independent researcher's problem of #access to #information.”@deepthiw
“My last article couldn't be possible without the platform @deepdyve that makes journal papers cheaper.”@JoseServera