The phylogenomics and comparative functional genomics of avian species were investigated in the Bird 10,000 Genomes (B10K) project because of the important evolutionary position of birds and their value as a research model. However, the systematic profiling of transcriptional changes prior to oviposition has not been investigated in avian species because of the practical difficulties in obtaining pre-oviposited eggs. In this study, a total of 137 pre-oviposited embryos were collected from hen ovaries and oviducts and subjected to RNA-sequencing analyses. Two waves of chicken zygotic genome activation (ZGA) were observed. Functionally distinct developmental programs involving Notch, MAPK, Wnt and TGFβ signaling were separately detected during cleavage and area pellucida formation. Furthermore, the early stages of chicken development were compared with the human and mouse counterparts, highlighting chicken-specific signaling pathways and gradually analogous gene expression via ZGA. These findings provide a genome-wide understanding of avian embryogenesis and comparisons among amniotes.

INTRODUCTION

Avian species have been used as valuable research models for virology, cancer, immunology, developmental biology and biotechnology (Streit et al., 2000; Skromne and Stern, 2001; Stern, 2005; Han, 2009; Torlopp et al., 2014). Owing to their important position in the comparative genomics of vertebrate species and wide relevance across many research fields, the Bird 10,000 Genomes (B10K) project was recently initiated. Through this project, a phylogenetic hierarchy of avian species involving comparative genomics regarding flight and functional adaptations was generated (Jarvis et al., 2014; Zhang et al., 2014a, 2015). Although historically significant, no genome-wide study on the initial and crucial events of avian embryogenesis prior to oviposition has yet been reported, mainly because of the practical difficulties in accessing pre-oviposited embryos.

Based on the criteria of Eyal-Giladi and Kochav (EGK) (Eyal-Giladi and Kochav, 1976; Kochav et al., 1980), intrauterine embryonic development before oviposition consists of ten stages, which include separate periods of cleavage and area pellucida formation (Sheng, 2014). After the first cleavage at EGK.I, representative processes, such as rapid asymmetric cellularization and an increased number of cell layers, take place until EGK.VI. Embryonic polarity, by which tilting and gravity forces lead to the reorientation of the yolk and formation of the posterior side, and the expression of germ layer markers in a salt-and-pepper manner are initiated from EGK.VII to EGK.X, during area pellucida formation (Kochav and Eyal-Giladi, 1971; Mak et al., 2015). However, from more than a morphological perspective, investigating the temporal regulation of genes controlling fundamental cellular programming and signaling pathways is essential in order to understand the early embryonic development of each species.

Transcriptional regulation of early developmental processes, including zygotic genome activation (ZGA), functional networks and cellular programming, has already been characterized using next-generation sequencing technology in representative model organisms, such as mouse, cow, frog and zebrafish, as well as in humans (Aanes et al., 2011; Cantone and Fisher, 2013; Tan et al., 2013; Xue et al., 2013; Yan et al., 2013; Graf et al., 2014; Lee et al., 2014), but not in avian species. Recently, molecular studies on the cellularization process, pluripotent state, and germ cell specification in avian species have been enabled using a non-invasive method for acquiring intrauterine chick embryos and zebra finch embryos at stages EGK.VI-VII (Lee et al., 2013; Mak et al., 2015; Nagai et al., 2015). Using this method, we performed RNA sequencing (RNA-seq) in oocytes, zygotes and intrauterine embryos of chickens to investigate comprehensively the transcriptional changes involved in early embryogenesis. To our knowledge, this is the first transcriptomic analysis approach investigating avian early embryonic development. Through this study, we sought to fill the gap in the comparative genomics of model organisms by identifying avian developmental mechanisms.

RESULTS

Early chicken embryo acquisition and RNA-seq data processing

Based on the EGK morphological staging system (Eyal-Giladi and Kochav, 1976; Kochav et al., 1980), we sampled embryos from the shell gland at the following stages: EGK.I (first cleavage), EGK.III (ZGA), EGK.VI (last stage of the cleavage period), EGK.VIII (progression of area pellucida formation) and EGK.X (oviposition) (Fig. 1A). Additionally, we sampled oocytes from the ovary and zygotes from the magnum to assess changes upon fertilization (Fig. 1A). Intrauterine eggs were retrieved from oviposition time-checked hens by non-invasive collection methods (Fig. 1A; see Materials and Methods). In total, 137 early embryos from oocyte to EGK.X-stage blastoderms (15 oocytes from ovaries, 9 zygotes from the magnum, and 17 EGK.I, 19 EGK.III, 20 EGK.VI, 27 EGK.VIII and 30 EGK.X embryos from the shell glands of hens) were collected according to the EGK morphological criteria (Fig. S1). Between three and ten embryos of each stage were randomly pooled, and the total RNA content per embryo was assessed in each pooled sample (Table S1). A single embryo contained 1457 ng total RNA, on average, through intrauterine development. Based on the total RNA, RNA-seq was performed using a total of 21 pooled samples (three biological replications per stage from oocyte to EGK.X) on the Illumina NextSeq500 system. As a result, 58,931,937 paired-end reads were available after adapter sequence trimming, and the poor-quality reads were filtered out (Table S2) from the originally generated 1279 million reads (386 Gb). Finally, 35,291,134 paired-end reads were mapped to the Ensembl chicken reference genome (galGal4; Table S3), and those reads were mapped to 22,010 genes, after filtering out non-mapped genes across all samples.

Fig. 1.

RNA-seq analysis of early chicken embryos. (A) Schematic of early development and embryo acquisition in chicken. Intrauterine embryos from EGK.I, showing first cleavage, to EGK.VIII, in which area pellucida formation and cell layer reduction were in progress from the posterior side (light gray), were collected from the shell gland. EGK.X embryos in which the area pellucida was fully formed (light gray) and the area opaca present (outlined in dark gray) were obtained after oviposition. Germinal vesicle (represented by the round circle in the center) oocytes in the ovary and zygotes without any cleavage in the magnum were also retrieved according to egg-laying time. h, staying duration for each part of oviduct and hours after fertilization for each stage of embryos. (B) Multidimensional scaling (MDS) plot, based on gene expression of the whole transcriptome in pre-oviposited chicken embryos. Each point represents an RNA-seq sample, and the dashed ovals represent clustered samples. (C) Number of significantly up- and downregulated genes between consecutive stages (FDR-adjusted P<0.05).

Fig. 1.

RNA-seq analysis of early chicken embryos. (A) Schematic of early development and embryo acquisition in chicken. Intrauterine embryos from EGK.I, showing first cleavage, to EGK.VIII, in which area pellucida formation and cell layer reduction were in progress from the posterior side (light gray), were collected from the shell gland. EGK.X embryos in which the area pellucida was fully formed (light gray) and the area opaca present (outlined in dark gray) were obtained after oviposition. Germinal vesicle (represented by the round circle in the center) oocytes in the ovary and zygotes without any cleavage in the magnum were also retrieved according to egg-laying time. h, staying duration for each part of oviduct and hours after fertilization for each stage of embryos. (B) Multidimensional scaling (MDS) plot, based on gene expression of the whole transcriptome in pre-oviposited chicken embryos. Each point represents an RNA-seq sample, and the dashed ovals represent clustered samples. (C) Number of significantly up- and downregulated genes between consecutive stages (FDR-adjusted P<0.05).

Transcriptomic dynamics representing two waves of ZGA in early chicken development

Based on quantified gene expression, the relationships among the samples and stages were investigated using multidimensional scaling methods (Fig. 1B). From a stage standpoint, most intrauterine developmental stages were clearly distinguishable, but three stages (zygote, EGK.I and EGK.III) were not. Statistical tests were performed to identify differentially expressed genes (DEGs) between adjacent stages according to developmental stage (i.e. oocyte-zygote, zygote-EGK.I, EGK.I-III, EGK.III-VI, EGK.VI-VIII and EGK.VIII-X). As a result, pairwise comparisons of gene expression between consecutive stages revealed no dynamic changes among these developmental stages (EGK.I and EGK.III) at the 5% significance level with false discovery rate (FDR) multiple testing adjustment, whereas large numbers of DEGs were observed in other consecutive stages (Fig. 1C).

Among the six consecutive stage comparisons, the greatest changes in gene expression levels (4819 upregulated and 4129 downregulated DEGs) were observed between the oocyte and zygote stages (Fig. 1C). In most cases, equivalent numbers of up- and downregulated genes were observed; however, the significant DEGs between EGK.III and EGK.VI were skewed towards upregulation (1287 upregulated versus 13 downregulated). These results could indicate that there are two waves of ZGA, from oocyte to zygote and from EGK.III to EGK.VI, designated as the first (first ZGA) and second (second ZGA) waves of ZGA, based on previous reports in chicken and on minor and major ZGA in mammals (Lee et al., 2013; Xue et al., 2013; Nagai et al., 2015).

Developmental networks and signaling pathways associated with cell division after fertilization in chicken

To identify the functional characteristics of the DEGs between consecutive stages (i.e. oocyte-zygote, zygote-EGK.I, EGK.I-III, EGK.III-VI, EGK.VI-VIII and EGK.VIII-X), enrichment tests were performed for both significantly up- and downregulated genes based on the gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases (FDR-adjusted P<0.05). From these analyses, we expected to identify transition-specific and pattern-enriched functions using these gene lists, respectively. As a result, diverse transition-specific genes (Tables S4-S7) were identified (enrichment test P-value<0.05). Because only four transition points, oocyte-zygote, EGK.III-VI, EGK.VI-VIII and EGK.VIII-X, were observed among the seven developmental stages, only their corresponding DEG lists were analyzed further. The large number of DEGs among the four transitions (Fig. 1B) produced a broad spectrum of functional terms in our enrichment analysis.

Next, functional term clustering analysis was performed using significantly detected GO and KEGG terms derived from DAVID analysis, considering N:M relationships among the combination of various biological terms and genes. Based on the significant terms and their associated gene information (see details in Materials and Methods), correlations were calculated to estimate co-occurrence and co-relationships between biological terms (an FDR-adjusted P-value<0.01 derived from correlation tests was considered significant). The identified relationships were visualized as networks to demonstrate the co-occurring terms.

In the results, many biological terms were changed between the oocyte and zygote stages (Fig. 2A), which is the first ZGA after fertilization, exhibiting the most dynamic changes during early development. Of the observed terms, phosphate metabolic process (GO:0006796), histone acetylation (GO:0016573), cell cycle (gga04110), DNA replication (gga03030), MAPK signaling pathway (gga04010), regulation of small GTPase-mediated signal transduction (GO:0051056), Wnt signaling pathway (gga04310) and Notch signaling pathway (gga04330) were upregulated. Transmembrane transport (GO:0055085), apoptosis (GO:0006915) and ribosome (gga03010) were downregulated. These results indicate that the biological networks induced by the first ZGA might be involved in the cleavage period and remain constant until EGK.VI.

Fig. 2.

Functional network of significantly expressed genes between oocyte and zygote stages during the cleavage period. (A) Functionally clustered network using significantly detected biological terms, such as biological processes and pathways, from the GO and KEGG databases between the oocyte and zygote stages. Each node represents a biological term, and an edge represents a significant correlation between two terms (FDR-adjusted P<0.01). (B) The expression patterns of the gene sets involved in cell cycle, DNA replication and MAPK signaling after fertilization. (C) Heatmap showing the expressions of the signaling pathways involved in small GTPase, WNT and Notch signaling during the early cleavage period.

Fig. 2.

Functional network of significantly expressed genes between oocyte and zygote stages during the cleavage period. (A) Functionally clustered network using significantly detected biological terms, such as biological processes and pathways, from the GO and KEGG databases between the oocyte and zygote stages. Each node represents a biological term, and an edge represents a significant correlation between two terms (FDR-adjusted P<0.01). (B) The expression patterns of the gene sets involved in cell cycle, DNA replication and MAPK signaling after fertilization. (C) Heatmap showing the expressions of the signaling pathways involved in small GTPase, WNT and Notch signaling during the early cleavage period.

We investigated the upregulated functional terms and signaling pathways in depth. Many cyclins and cyclin-dependent kinases and a representative S-phase marker, PCNA, were expressed and grouped into cell cycle and DNA replication functions (Fig. 2B; Table S8). Genes involved in MAPK signaling, such as RAS, RAF and several MAP kinases, were also activated (Fig. 2B; Table S8). Additionally, the genes encoding GTPase cycle regulators of small GTPase-mediated signal transduction, ligands and receptors of the Wnt signaling pathway, downstream effectors of the Wnt/planar cell polarity (PCP) pathway, and Notch signaling pathway-related genes were observed (Fig. 2C; Table S9). We examined the spatial expression of Wnt and Notch signaling-related genes including WNT4, WNT6, NOTCH1 and NOTCH2 during cleavage stages by in situ hybridization (Fig. 3). As a result, we showed that WNT4, WNT6 and NOTCH1 were expressed from EGK.I to EGK.VI, and NOTCH2 was expressed in EGK.I and EGK.III embryos (Fig. 3A,B). Furthermore, asymmetric expression of these genes was observed in EGK.III embryos (Fig. 3C). These biological terms seem to be related to the rapid asymmetric cellularization that occurs during the cleavage period.

Fig. 3.

In situ hybridization of WNT and Notch signaling-related genes during cleavage stages. (A,B) Dorsal views of whole-mount in situ hybridization of WNT4, WNT6 (A), NOTCH1 and NOTCH2 (B) from EGK.I to EGK.VI. (C) Longitudinal section views of asymmetric expression of WNT4, WNT6, NOTCH1 and NOTCH2 in EGK.III. Scale bars: 500 µm.

Fig. 3.

In situ hybridization of WNT and Notch signaling-related genes during cleavage stages. (A,B) Dorsal views of whole-mount in situ hybridization of WNT4, WNT6 (A), NOTCH1 and NOTCH2 (B) from EGK.I to EGK.VI. (C) Longitudinal section views of asymmetric expression of WNT4, WNT6, NOTCH1 and NOTCH2 in EGK.III. Scale bars: 500 µm.

The definite functional signatures of morphological changes from the EGK.VI stage

In contrast to the preceding stages, distinct programs were observed from EGK.VI to EGK.X. The functional networks are characterized and shown in Fig. 4A-C. By the second ZGA, between EGK.III and EGK.VI, the upregulated signaling pathways, such as the Wnt signaling pathway (gga04310) and TGFβ signaling pathway (gga04350), were clearly observed (Fig. 4A). The downregulated functional networks, which were involved in cleavage such as phosphate metabolic process, histone acetylation, MAPK signaling pathway, regulation of small GTPase-mediated signal transduction, and Notch signaling pathway, initially appeared from EGK.VIII to EGK.X (Fig. 4A-C).

Fig. 4.

Functional network for the area pellucida formation period after the EGK.VI stage. (A-C) Functionally clustered networks using significantly detected biological terms between stages EGK.III and EGK.VI (A), EGK.VI and EGK.VIII (B), and EGK.VIII and EGK.X (C). Each node represents a biological term, and an edge represents a significant correlation between two terms (FDR-adjusted P<0.01). (D) Expression of the gene sets associated with WNT and TGFβ signaling. (E) Expression of the terms involved in development of the three germ layers and in ECM-receptor interaction during the area pellucida formation period.

Fig. 4.

Functional network for the area pellucida formation period after the EGK.VI stage. (A-C) Functionally clustered networks using significantly detected biological terms between stages EGK.III and EGK.VI (A), EGK.VI and EGK.VIII (B), and EGK.VIII and EGK.X (C). Each node represents a biological term, and an edge represents a significant correlation between two terms (FDR-adjusted P<0.01). (D) Expression of the gene sets associated with WNT and TGFβ signaling. (E) Expression of the terms involved in development of the three germ layers and in ECM-receptor interaction during the area pellucida formation period.

Notably, regulators of Wnt and TGFβ signaling, including various ligands, inhibitors, receptors and downstream mediators, were found (Fig. 4D; Table S10). Specifically, Wnt signaling-related genes were grouped into three clusters based on expression levels (WNT7A, WNT7B, WNT9B and WNT3 in one cluster, FZD7 in the second, and WNT11, WNT3A, WNT5A, WNT9A, DKK1, WNT8C, WNT8B and others in the third), but we were unable to determine the functional categorization of the ligands in each cluster, because canonical (WNT3, WNT3A, WNT7B, WNT8C and WNT8B) and non-canonical (WNT5A, WNT11 and WNT7A) ligands were both present. These signaling pathways were increased consistently until EGK.X, indicating that the candidate signaling pathways responsible for onset of area pellucida formation are activated after EGK.VI (Fig. 4A-C). Furthermore, several biological processes potentially related to lineage specification during intrauterine development were observed, such as neuron development (GO:0048666), heart development (GO:0007507), lung development (GO:0030324) and ECM–receptor interactions (gga04512), which seemed to increase continuously from EGK.VI to EGK.X (Fig. 4E; Table S11). These regulators, which were distinct from those associated with the cleavage stages, might control the area pellucida formation period.

Multi-species comparisons among vertebrate embryos, including pre-oviposited chicken, pre-implantation human and mouse, and early zebrafish embryos

First, to identify the distinct developmental programs involved in cleavage and area pellucida formation in chicken embryos, a co-expression network analysis based on correlations was performed using protein-coding genes known to be important in early embryonic development (Hamatani et al., 2004; Wang et al., 2004; Yanai et al., 2011). The expression of genes, including representative signaling genes encoding ligands, receptors, DNA-binding factors, and kinases of the Notch, MAPK, Wnt and TGFβ signaling pathways along with significantly expressed well-known transcription factors in early chicken development, was analyzed (Fig. 5A; Table S12). In the co-expression analysis, pairwise correlation tests were employed, and a 1% significance level after FDR multiple-testing adjustments was considered to indicate a co-expressed relationship between two genes. The network showed distinctions between the cleavage and area pellucida formation stages. For example, genes encoding all Notch signaling pathway, MAPK signaling genes except MAP2K2, five Wnt ligands and receptors such as WNT4, WNT6, FZD1, RYK and CCNYL1, and ETS2 were enriched only in the subnetwork of cleavage, indicating a significant function in rapid cellularization. In contrast, pluripotency markers, such as NANOG, POUV, MYC and SOX2, and the other five Wnt ligands and receptors (WNT5B, WNT8C, DKK1, FZD7 and FZD10) were positively enriched in the subnetwork of area pellucida formation. Considering the functional characteristics of the genes contained in this subnetwork, it is possible that they have key roles in area pellucida formation and the maternal-to-zygotic transition (MZT). These genes are thus candidates for the clearly distinct genetic networks involved in cleavage and area pellucida formation period during early chicken embryogenesis.

Fig. 5.

The collective network of distinct developmental programs and comparison of transcriptomic profiles between pre-oviposited chicken embryos and early mammalian embryos. (A) The mRNA co-expression network, composed of diverse transcriptomes, includes representative signaling pathway-related genes and well-known transcription factors. The node size represents centrality, determined based on statistically significant relationships (FDR-adjusted P<0.01). (B) Heat map obtained from a comparison between early chicken and early human and mouse embryos. The relative similarity measure was used, based on gene expression in the transcriptome. The color from red to blue shows the relative similarity between stages. The x-axis represents the early developmental stages of human or mouse, and the y-axis represents the profiling of the early developmental stages of chicken. The dotted and solid boxes, respectively, indicate the first and second transitions in chicken versus the other two species.

Fig. 5.

The collective network of distinct developmental programs and comparison of transcriptomic profiles between pre-oviposited chicken embryos and early mammalian embryos. (A) The mRNA co-expression network, composed of diverse transcriptomes, includes representative signaling pathway-related genes and well-known transcription factors. The node size represents centrality, determined based on statistically significant relationships (FDR-adjusted P<0.01). (B) Heat map obtained from a comparison between early chicken and early human and mouse embryos. The relative similarity measure was used, based on gene expression in the transcriptome. The color from red to blue shows the relative similarity between stages. The x-axis represents the early developmental stages of human or mouse, and the y-axis represents the profiling of the early developmental stages of chicken. The dotted and solid boxes, respectively, indicate the first and second transitions in chicken versus the other two species.

Because our transcriptome analysis showed the involvement of several signaling pathways such as Notch, MAPK, Wnt and TGFβ signaling in two defined periods of chicken early development, we presumed that gene expression in other well-known species would be different because of the difference in ontogeny between aves and mammals. Additionally, we considered that the counterpart developmental stages in other amniotes needed to be compared with those in the chicken embryo. To address these issues, we performed a comparative analysis of gene expression in pre-implantation embryos in humans and mice (GSE44183 in the GEO database) (Xue et al., 2013). In those species, gene expression from the oocyte to the morula stages (oocyte, pronucleus, zygote, 2-cell, 4-cell, 8-cell and morula) was assessed as pre-implantation developmental stages. For the multi-species comparison, 1:1:1 orthologous genes were defined using the Ensembl database. Next, a reads per kilobase per million mapped reads (RPKM) normalization method was used to assess different gene lengths among the species. Finally, a relative similarity measure was used to estimate similarity among the various developmental stages (Fig. 5B; see details in Materials and Methods). Relative similarities tended to increase gradually with developmental progression in the three species. In the comparison between human and mouse, similar patterns of relative similarity with embryonic development were observed (Fig. S2). In addition, minor and major ZGAs are processed in human (zygote and 8-cell) and mouse (zygote and 2-cell) (Lee et al., 2014), which is consistent with our relative similarity analysis results. Furthermore, the three observed transitions in chicken species were closely correlated with our observations in the transcriptomic analysis. In addition to mammals, we compared the transcriptomes of early embryos in zebrafish [2-cell, 64-cell, 3.3 hours post-fertilization (hpf), 6 hpf and 9 hpf] with those of chickens. After identifying 1:1 orthologs between chicken and zebrafish, the relative similarities were examined (Fig. S3). As with mammals, the values for zebrafish embryos also increased continuously as development progressed and increased sharply with ZGA.

In addition to the relative similarity analysis among species, we compared orthologous gene expression related to the cell signaling networks indicated in Fig. 5A during representative stages: oocyte, 1-cell zygote or pronucleus, post-major ZGA (EGK.VI in chicken, 8-cell stage in human and 4-cell stage in mouse) and post-MZT (EGK.X in chicken and morula in human and mouse) (Fig. 6; Table S13). From these analyses, Notch signaling genes, such as DLL4, NOTCH1, NOTCH2 and RBPJ, were strongly expressed at all stages in chicken, but not in human or mouse (Fig. 6A). Wnt ligands, including WNT6 during cleavage and WNT5B during area pellucida formation, were also significantly expressed in chicken compared with mammals (Fig. 6B). MAP kinases, MAPK1, MAP2K3 and MAPK14, which exhibited low expression during early cleavage in mammals, were predominantly expressed during the cleavage period (Fig. 6C). In the TGFβ family, TGFB3 expression was specific to the area pellucida formation period, and its expression was lower at later stages in mammals (Fig. 6D). These results indicate that some genes related to signaling, such as Notch, MAPK, Wnt and TGFβ, are differentially expressed between chicken and mammals during the early developmental program.

Fig. 6.

The comparison of orthologous genes during representative stages among chicken, human and mouse. (A-D) Gene expression related to Notch (A), Wnt (B), MAPK (C) and TGFβ (D) cell signaling. Representative stages including oocyte, 1-cell zygote or pronucleus, post-major ZGA (EGK.VI in chicken, 8-cell in human and 4-cell in mouse), and post-MZT (EGK.X in chicken and morula in human and mouse) were compared. Normalized RPKM values of genes in each replicate are represented by colored dots and mean values by shaded lines for analogous stages in chicken (red), human (green) and mouse (blue).

Fig. 6.

The comparison of orthologous genes during representative stages among chicken, human and mouse. (A-D) Gene expression related to Notch (A), Wnt (B), MAPK (C) and TGFβ (D) cell signaling. Representative stages including oocyte, 1-cell zygote or pronucleus, post-major ZGA (EGK.VI in chicken, 8-cell in human and 4-cell in mouse), and post-MZT (EGK.X in chicken and morula in human and mouse) were compared. Normalized RPKM values of genes in each replicate are represented by colored dots and mean values by shaded lines for analogous stages in chicken (red), human (green) and mouse (blue).

All of this evidence, including the significantly detected genes, and associated functional terms, during the zygote and EGK.VI stages and zygote and EGK.VIII stages, is illustrated in Fig. 7. Early chicken, human and mouse embryos exhibit parallel developmental processes such as cell cycle, transcription, and translation, according to the results of a previous study (Xue et al., 2013) and the transcriptomic analysis in this study, despite differences in ontogenetic physiology and expression of signaling pathway-related genes between aves and mammals.

Fig. 7.

Summary and complementation of early chicken embryo developmental processes with those of human and mouse. Parallel stages of development for chicken, human and mouse were mapped based on the analysis of functional terms and developmental processes including ZGA. Signaling pathways such as Notch, MAPK, Wnt and TGFβ signaling were shown to be avian specific.

Fig. 7.

Summary and complementation of early chicken embryo developmental processes with those of human and mouse. Parallel stages of development for chicken, human and mouse were mapped based on the analysis of functional terms and developmental processes including ZGA. Signaling pathways such as Notch, MAPK, Wnt and TGFβ signaling were shown to be avian specific.

DISCUSSION

From the moment of fertilization, precise regulation of gene expression and dynamic courses occurs during early embryogenesis in various organisms. However, most of the early developmental programs in avian species have not been determined. In this study, RNA-seq helped identify the steps occurring in pre-ovipositional embryonic development in the chicken. To examine in a systematic manner the genetic regulation involved in early chicken embryogenesis compared with other species, transcriptomic data were analyzed from various perspectives.

Of the many observations from our transcriptomic analysis, the most interesting was the very large number of actively expressed transcripts between the oocyte and zygote stages (Fig. 1B). This phenomenon occurs between the oocyte and zygote stages in other mammals, the so-called minor ZGA. In previous research, transcriptional activity was measured in terms of bromouridine triphosphate incorporation, and greater incorporation was seen in the male pronucleus from the 1-cell embryo stage (Aoki et al., 1997). Additionally, soon after fertilization in the mouse, epigenetic modification of the paternal pronucleus changes vigorously (Aoshima et al., 2015). From this evidence, we specified this transcriptional event in the chicken as the first ZGA, which remains to be elucidated in depth. Surprisingly, no transcriptional differences were observed among the zygote, EGK.I and EGK.III stages (Fig. 1B). Indeed, these results suggest that transcriptional silencing occurs during early chicken embryonic development, similar to the cleavage stages (from the 1- to 4-cell stage) of human embryos (Braude et al., 1988; Xue et al., 2013). Even though central cells positive for phosphorylation of the RNA polymerase II C-terminal domain (p-PolII) were apparent during the late EGK.II to early EGK.III stage according to a previous report (Nagai et al., 2015), significant expression was not found at this point, when cells are still rapidly dividing, which was confirmed at the transcriptome level in this study. In contrast, from the EGK.III to the EGK.VI stage, 1287 transcripts with biological functions involved in transcription were upregulated, and only 13 transcripts were downregulated; this difference indicated the presence of another ZGA. This may be accompanied by lengthening of the cell cycle, as observed in Xenopus (Newport and Kirschner, 1982). Based on these observations, we suggest that the first and second ZGA events at the transcript level occurred between the oocyte and zygote stages and between the EGK.III and EGK.VI stages, respectively.

In the first and second ZGA, genes that seemed to be correlated functionally with the cleavage and area pellucida formation periods, respectively, were induced (Fig. 5A). From the functional network analysis results (Fig. 2A; Fig. 4A-C), we found relevance between the early developmental courses and transcriptomic signatures. First, the genes activated during the first ZGA were largely associated with the cleavage period. During these stages, we found several functional processes enriched after fertilization. Cell cycle, DNA replication and MAPK signaling appear to be involved in mitosis in early embryos. In sea urchin embryos after fertilization, the MAPK pathway, via Cdk2 activation, regulates DNA replication and entry into mitosis (Zhang et al., 2005; Kisielewska et al., 2009). The cleavage stages during intrauterine development in the chicken appear with rapid cellularization from a single cell to approximately 8000 cells within 10 h of the first cleavage event. Mitotic modulators from the first ZGA might regulate sustainable and rapid cell division in the early chicken embryo.

Notch signaling plays important roles in cell–cell interactions, mediated by Delta ligands and Notch receptors. As an example, in Caenorhabditis elegans, Notch signaling is involved in early cell fate specification for asymmetric differentiated cell functions (Priess, 2005). Similarly, chicken NOTCH1 and NOTCH2 were expressed asymmetrically in cleavage stage embryos (Fig. 3C). Further investigation of the Notch signaling pathway during the first ZGA is needed to clarify whether it actually shows cellular diversity in early chicken embryos. Moreover, the spatiotemporal expression of marker genes, which are related to the three germ layers during later stages (Fig. 4E; Table S11), could be investigated to elucidate early fate mapping during the cleavage stages in chicken. Wnt signaling and small GTPase signaling might interact to regulate cell polarity and cleavage patterns. During embryogenesis, maternal Wnt signaling is necessary for axis formation via the canonical pathway in Xenopus (Tao et al., 2005). A recent finding suggested that the non-canonical Wnt/STOP pathway of maternal Wnt signaling plays a role in the progression of cell cleavage during early Xenopus embryogenesis (Huang et al., 2015). Additionally, the small GTPase RhoA acts as a regulator of spindle formation and cytokinesis during early embryonic development in C. elegans and Sus scrofa (Tse et al., 2012; Zhang et al., 2014b). Furthermore, the Wnt/PCP pathway, via RhoA interaction with the actin cytoskeleton, regulates the orientation of cell division (Castanon et al., 2013). During early cleavage events in chicken embryos, formation of the cleavage furrow in the embryonic central region is asymmetrical, occurring in a radial manner (Lee et al., 2013). Moreover, the increase in layer number during the cleavage period appears to be caused by cell division orientation, to some extent (Nagai et al., 2015). In this regard, maternal Wnt ligands in the chicken, such as WNT4 and WNT6, which differ from maternal WNT11 in Xenopus (Tao et al., 2005), might regulate cleavage formation and orientation during early development. During cleavage in EGK.III, WNT4 and WNT6 showed asymmetric expression, suggesting the diversification of embryonic cells in this period (Fig. 3C). However, WNT4 and WNT6 need to be further elucidated in terms of their interaction with their receptor, because of their characteristics as species-specific maternal Wnt ligands. Consequently, the representative processes during the cleavage period, such as rapid asymmetric cellularization and increase in layer number, could be regulated by genes expressed during the first ZGA.

Otherwise, the maternal Wnt ligands were downregulated and replaced with WNT5B, WNT8C and the antagonist DKK-1, which could regulate axis formation, during the second ZGA (Fig. 5A). In addition, TGFβ signaling might also cooperate with Wnt signaling in axis formation and lineage specification. Even though the interaction between Wnt and TGFβ/activin signaling is known to regulate the initiation of primitive streak formation during post-ovipositional development (Skromne and Stern, 2001), the synergism between these pathways in pre-ovipositional development has yet to be elucidated, primarily for the more highly expressed Wnt ligands including WNT11, WNT3A, WNT5A and WNT9A, in terms of both canonical and non-canonical pathways (Fig. 4D). Along with these signals, we also observed, by functional enrichment analysis, simultaneous expression of the three germ layers in neuronal development (ectoderm), heart development (mesoderm) and lung development (endoderm) after the EGK.VI stage. Moreover, the ECM–receptor interaction, which might play an important role in tissue and organ morphogenesis, was identified during the area pellucida formation stage (Fig. 4A-C). Four genes related to the ectoderm and endoderm, SOX3, SOX17, GATA4 and GATA6, were expressed during stages EGK.VII to EGK.VIII in the zebra finch (Mak et al., 2015). However, the simultaneous expression associated with the three germ layers in lineage specification during intrauterine development in the chicken is difficult to demonstrate and requires further study. Collectively, embryonic polarity and gene expression related to the hypoblast and epiblast appear to be activated by the second ZGA during area pellucida formation from the EGK.VI stage.

Finally, we compared gene expression patterns in chickens with other amniotes, human and mouse, during early embryogenesis. The ontogenetic physiology of early development in aves and mammals, such as cell division patterns and lineage segregation, seems to be different (Sheng, 2014). Even though some signaling genes were found to be restricted to chicken development, unexpectedly, all stages of the chicken showed gradually increasing similarity in gene expression to the mammalian stages after the minor and major ZGA, which occur in the pronucleus and 8-cell stage in humans and the 2-cell stage in mice, respectively (Lee et al., 2014). Furthermore, similar patterns were observed in chicken and zebrafish regarding ZGA (Fig. S3). Additionally, the relative similarities were greater after the major ZGA than after the minor ZGA (Fig. 5B), indicating that the earliest genes of the first ZGA or minor ZGA differed evolutionarily, as shown in the comparison of species in a previous study (Heyn et al., 2014). This observation suggests that the global transcriptome and basic cellular functions, including cell cycle, small GTPase signaling, transcription and translation, in avian and mammalian embryos increasingly overlap and resemble each other over the course of development, despite the divergence in their timing of activation. In human and mouse embryos, fundamental programs such as cell cycle before the major ZGA (4-cell stage in humans and 1-cell stage in mice), transcription after the major ZGA (8-cell stage in humans and 2-cell stage in mice), and translation (morula in humans and 8-cell stage and morula in mice) were identified in a previous study (Xue et al., 2013). In this study, we found cell cycle before EGK.VI, transcription after EGK.VI, and translation after EGK.VIII to be enriched based on functional annotation (Fig. 2A, Fig. 3A-C). Considering all of the available evidence, in this study, we summarized the chicken pre-ovipositional developmental stages and correlated the chicken stages with those of early mammalian early development, although differential gene expression related to signaling pathways existed between chicken and mammals (Fig. 7). In conclusion, we report here the first temporal transcriptomic analysis during early chicken development. Our results provide an explanation of the developmental signaling features governed by two waves of ZGAs during early developmental events in avian species, at the transcriptomic level. From these results, we foresee important studies in early avian development, including fundamental transcriptome and systematic signaling pathways, from an evo-devo perspective. Furthermore, our results will contribute to a greater understanding of early avian embryogenesis and bridge gaps in comparative studies with other amniotes in various developmental respects.

MATERIALS AND METHODS

Experimental animals and animal care

The care and experimental use of chickens were approved by the Institute of Laboratory Animal Resources, Seoul National University (SNU-150827-1). Chickens were maintained according to a standard management program at the University Animal Farm, Seoul National University, Korea. The procedures for animal management, reproduction and embryo manipulation adhered to the standard operating protocols of our laboratory.

Chicken early embryo preparation

The egg-laying times of the white leghorn (WL) hens were recorded, and intrauterine eggs from EGK stages I-VIII were harvested (Fig. 1A) using an abdominal massage technique (Lee et al., 2013). Briefly, the abdomen was pushed gently until the shell gland was exposed; the surface of the shell gland expanded when an egg was present for egg shell formation. After expansion of the shell gland surface, massaging was used to move the egg gently towards the cloaca until the intrauterine egg was released. EGK stage X blastoderms were collected from WL hens after oviposition. To collect oocytes and zygotes, WL hens were sacrificed and the follicles were collected. Oocyte and zygote were collected from one hen simultaneously. Because of the little transcriptomic difference between pre- and post-ovulatory oocyte in the previous study (Elis et al., 2008), and the infeasible simultaneous acquisition of post-ovulatory oocyte and zygote from one hen, we isolated the only pre-ovulatory large F1 oocyte. Only the zygote embryos not showing cleavage and located in the magnum were collected within 1 h after fertilization according to the recorded egg-laying times. All embryos were classified according to the morphological criteria (Fig. S1). Shortly after collection, the embryos were separated from the egg using sterilized paper, and the shell membrane and albumen were detached from the yolk. A piece of square filter paper (Whatman) with a hole in the center was placed over the germinal disc. After cutting around the paper containing the embryo, it was gently turned over and transferred to saline to remove the yolk and vitelline membrane to allow embryo collection.

RNA isolation and RNA-seq library preparation

Total RNA was isolated from early embryos using TRIzol reagent (Invitrogen). The quality of the extracted total RNA was determined using the DropSense96 system (Trinean), Ribogreen (Invitrogen) and the Agilent 2100 Bioanalyzer (Agilent Technologies) (Table S1). Total RNA was used for construction of cDNA libraries using the TruSeq Stranded Total RNA Sample Preparation kit (Illumina). The resulting libraries were subjected to transcriptome analysis using the Illumina Nextseq 500 platform to produce paired 150 bp reads.

RNA-seq data pre-processing

Based on the RNA libraries generated, paired-end sequencing (150 bp read length and ∼150-200 bp insert size) was performed using the NextSeq500 platform (Illumina). To generate clean reads, Trimmomatic (ver. 0.32) (Bolger et al., 2014) was used. Per-base sequence qualities were checked using FastQC (ver. 0.11.2) (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/) and filtered fastq files.

Alignment and quantification of mapped reads

The galGal4 Ensembl genome was used as a reference, and filtered reads were aligned onto this reference genome using Tophat2 (ver. 2.0.12) (Trapnell et al., 2009). SAMtools (ver. 0.1.18.0) (Li et al., 2009) was used to convert SAM to BAM files. To count the mapped reads, HTSeq-count (Anders et al., 2015) was used with the merged gene annotation file (.GTF), with mRNAs derived from Ensembl.

Detection of differentially expressed genes

As some recent statistical research suggests, it is appropriate to consider RNA-seq gene expression data as a negative binomial distribution, rather than a normal distribution, in statistical analyses (Robinson and Oshlack, 2010). Based on this, a negative binomial-based generalized linear model (GLM) was used to detect DEGs between stages. In total, six statistical tests, oocyte versus zygote, zygote versus EGK.I, EGK.I versus EGK.III, EGK.III versus EGK.VI, EGK.VI versus EGK.VIII, and EGK.VIII versus EGK.X, were performed using the edgeR package implemented in R (Robinson et al., 2010). Based on edgeR package, one-way analysis of deviance model was firstly fitted for estimating dispersion on the samples, as follows:
formula
(1)
where the Stage represents multiple group containing seven developmental stages across the 21 RNA-seq samples. Because we are only interested in differences between two adjacent stages of development, a total of six two-group comparisons (Oocyte-Zygote, Zygote-EGK.I, EGK.I-EGK.III, EGK.III-EGK.VI, EGK.VI-EGK.VIII and EGK.VIII-EGK.X) were performed using contrast matrix on the model described by Eqn. 1. Under the null hypothesis, H0:Before stage=After stage, a likelihood ratio test was performed with FDR-adjusted P<0.05 considered significant (Benjamini and Hochberg, 1995).

Functional analysis of differentially expressed genes

Biological processes from the GO and KEGG databases were used for functional annotation of the genes shown to be differentially expressed between adjacent stages. An enrichment test (hypergeometric test) was performed for each term using DAVID (Huang et al., 2007); here, P<0.05 was considered to indicate a significantly enriched term. Significantly detected GO and KEGG terms were visualized as a correlation-based network. To visualize the network, a distance-matrix is required, thus generating a gene association matrix (if a biological term associates a gene it is designated 1; if not, it is designated 0). Then, correlation coefficients can be calculated based on the generated gene-membership vectors by considering N:M relationships between biological terms and genes (estimated correlation is always positive value because the vector contains only 0 or 1). A correlation value close to 1 between two biological terms means that the two terms co-occur and a value closer to 0 indicates that the two terms are highly independent of one another. From these correlation-based distance matrixes, a network analysis was performed using the ‘qgraph’ package in R with ‘spring’ layout for placement based on centrality.

In situ hybridization

To make hybridization probes, the total RNA from each blastodermal stage was reverse transcribed, and the cDNA was amplified using chicken WNT4, WNT6, NOTCH1 and NOTCH2 primers (Table S14). The PCR products of the correct size were cloned with the pGEM-T Easy Vector System (Promega). After sequence verification, the recombinant plasmids containing the genes of interest were amplified with T7 (5′-TGTAATACGACTCACTATAGGG-3′) and SP6 (5′-CTATTTAGGTGACACTATAGAAT-3′) specific primers to prepare the templates for labeling with hybridization probes (Rengaraj et al., 2008). Digoxigenin (DIG)-labeled sense and antisense hybridization probes of each gene were transcribed in vitro using the DIG RNA Labeling Kit (Roche Diagnostics). Whole-mount in situ hybridization was performed following the standard published protocol for chickens (Stern, 1998). In addition, intrauterine embryos were embedded in paraffin and sectioned at 10 μm on an HM 355S automatic microtome (Thermo Fisher Scientific). After deparaffinization, rehydration and antigen retrieval, each slide was mounted with Vectashield Antifade Mounting Medium with DAPI (Vector Laboratories). The whole-mount embryos and sections were observed and imaged under a Ti-U fluorescence microscope (Nikon).

Correlation-based co-expressed network analysis of coding genes

In the co-expression analysis, a total of 42 genes were used including representative genes encoding proteins and transcription factors that are known to be important in early embryonic developmental of other model organisms (Table S12). First, mRNAs included in four representative signaling pathways (Notch, MAPK, Wnt and TGFβ signaling pathways) were used. Second, the differentially expressed significant transcription factors (ETS1, ETS2, c-myb, NANOG, POUV, SOX2 and MYC) were included. Based on the gene expression matrix of these 42 genes, pairwise correlation tests were performed and correlation coefficients were calculated for all possible combinations of those genes. Based on the correlation matrix (42×42), a network plot was generated using the qgraph package in R with the spring layout for considering clustering patterns of co-expressed genes. In this plot, only significant relationships (FDR-adjusted P<0.01 from the pairwise correlation test) were visualized as edges to identify relationships more conservatively. Finally, node size was determined based on the centrality of each gene (number of significant connections at FDR-adjusted P<0.01).

Species comparison of early embryonic gene expression patterns in human, mouse, zebrafish and chicken

In various analyses, early embryonic transcriptomic changes in chickens were identified. Although early embryonic transcriptomic changes have been investigated in human, mouse and zebrafish, no such reports have been conducted in chicken. Thus, one of our aims was to map the developmental stages among different species based on the gene expression patterns identified in this study. For the species comparison, early embryo RNA-seq data of other species were obtained from public transcriptome databases (Gene Expression Omnibus, accession number GSE44183, for human and mouse; European Nucleotide Archive, accession number ERP001280, for zebrafish). In this series, human (oocyte, pronucleus, zygote, 2-cell, 4-cell, 8-cell and morula), mouse (oocyte, pronucleus, 2-cell, 4-cell, 8-cell and morula) and zebrafish (2-cell, 64-cell, 3.3 hpf, 6 hpf and 9 hpf) samples were obtained. The raw sequenced files (.fastq) were processed in the same manner as chicken RNA-seq data in order to measure gene expression levels. Quality control, alignment and quantification methods were the same, but reference genome and gene-annotation files were different. For the reference genomes, genome builds of the hg19, mm9 and zv10 were used for human, mouse and zebrafish, respectively. For the between species comparison, RPKM normalization was performed to account for the different gene lengths in each species. In addition, 1:1:1 orthologous genes among human, mouse and chicken were identified using Ensembl Biomart tool to match different gene annotations of three species. As a result, three types of gene expression matrices (including 9810 genes) were generated for human, mouse and chicken, respectively. In the case of zebrafish, which is evolutionarily distant, only 1:1 orthologous genes between zebrafish and chicken were considered independently. A total of 6221 Ensembl genes were identified as 1:1 orthologous genes in two species. Based on these matrixes, we attempted to calculate relative similarities using the following formula:
formula
where rg(X) and rg(Y) are rank values. The score of relative similarity represents extended Spearman's correlation coefficients ([−1 to 1] to [0 to 1] scale). ‘0’ and ‘1’ values, respectively, represent the lowest and highest similarities given pairwise similarities. This statistic was used to identify transition points of similarities in another species given the similarity of all transcripts of other species based on the correlation in different developmental stages of different species. By linear transforming (stretching) the correlations of the whole gene expression levels to a range of values from 0 to 1 given the general similarities of the other species mapping the developmental process, this statistic makes it possible to identify transition points. If there are transition points, as in other species, this statistic will show a notable gap.

Acknowledgements

We thank the members of the Jae Yong Han's lab for suggestions and inspirational discussions.

Footnotes

Author contributions

Conceptualization: Y.S.H., J.Y.H.; Methodology: Y.S.H., M.S.; Software: M.S.; Validation: Y.S.H., B.R.L., H.J.L., Y.H.P., J.Y.H.; Formal analysis: Y.S.H., M.S.; Investigation: Y.S.H., M.S., S.K.K., H.J.C., H.K.; Resources: Y.S.H., H.K., J.Y.H.; Data curation: Y.S.H., M.S.; Writing - original draft: Y.S.H., M.S., H.K., J.Y.H.; Writing - review & editing: Y.S.H., M.S., B.R.L., H.J.L., Y.H.P., H.C.L., J.Y., H.K., J.Y.H.; Visualization: Y.S.H., M.S.; Supervision: H.K., J.Y.H.; Project administration: J.Y.H.; Funding acquisition: J.Y.H.

Funding

This work was supported by the Creative Research Initiatives Program (Center for Avian Germ Cell Modulation and Cloning) from a National Research Foundation of Korea (NRF) grant, funded by the Korean government (MSIP; no. 2015R1A3A2033826).

Data availability

Chicken RNA sequencing data were deposited in Gene Expression Omnibus (GSE86592). Human and mouse early embryo RNA-seq data were taken from Gene Expression Omnibus (GSE44183). Zebrafish RNA-seq data were taken from the European Nucleotide Archive (ERP001280).

References

Aanes
,
H.
,
Winata
,
C. L.
,
Lin
,
C. H.
,
Chen
,
J. P.
,
Srinivasan
,
K. G.
,
Lee
,
S. G. P.
,
Lim
,
A. Y. M.
,
Hajan
,
H. S.
,
Collas
,
P.
,
Bourque
,
G.
, et al. 
(
2011
).
Zebrafish mRNA sequencing deciphers novelties in transcriptome dynamics during maternal to zygotic transition
.
Genome Res.
21
,
1328
-
1338
.
Anders
,
S.
,
Pyl
,
P. T.
and
Huber
,
W.
(
2015
).
HTSeq – a Python framework to work with high-throughput sequencing data
.
Bioinformatics
,
31
,
166
-
169
.
Aoki
,
F.
,
Worrad
,
D. M.
and
Schultz
,
R. M.
(
1997
).
Regulation of transcriptional activity during the first and second cell cycles in the preimplantation mouse embryo
.
Dev. Biol.
181
,
296
-
307
.
Aoshima
,
K.
,
Inoue
,
E.
,
Sawa
,
H.
and
Okada
,
Y.
(
2015
).
Paternal H3K4 methylation is required for minor zygotic gene activation and early mouse embryonic development
.
EMBO Rep.
16
,
803
-
812
.
Benjamini
,
Y.
and
Hochberg
,
Y.
(
1995
).
Controlling the false discovery rate: a practical and powerful approach to multiple testing
.
J. R. Stat. Soc. Ser. B
,
57
,
289
-
300
.
Bolger
,
A. M.
,
Lohse
,
M.
and
Usadel
,
B.
(
2014
).
Trimmomatic: a flexible trimmer for Illumina sequence data
.
Bioinformatics
30
,
2114
-
2120
.
Braude
,
P.
,
Bolton
,
V.
and
Moore
,
S.
(
1988
).
Human gene expression first occurs between the four- and eight-cell stages of preimplantation development
.
Nature
332
,
459
-
461
.
Cantone
,
I.
and
Fisher
,
A. G.
(
2013
).
Epigenetic programming and reprogramming during development
.
Nat. Struct. Mol. Biol.
20
,
282
-
289
.
Castanon
,
I.
,
Abrami
,
L.
,
Holtzer
,
L.
,
Heisenberg
,
C. P.
,
van der Goot
,
F. G.
and
González-Gaitán
,
M.
(
2013
).
Anthrax toxin receptor 2a controls mitotic spindle positioning
.
Nat. Cell Biol.
15
,
28
-
39
.
Elis
,
S.
,
Batellier
,
F.
,
Couty
,
I.
,
Balzergue
,
S.
,
Martin-Magniette
,
M.-L.
,
Monget
,
P.
,
Blesbois
,
E.
and
Govoroun
,
M. S.
(
2008
).
Search for the genes involved in oocyte maturation and early embryo development in the hen
.
BMC Genomics
9
,
110
.
Eyal-Giladi
,
H.
and
Kochav
,
S.
(
1976
).
From cleavage to primitive streak formation: a complementary normal table and a new look at the first stages of the development of the chick. I. General morphology
.
Dev. Biol.
49
,
321
-
337
.
Graf
,
A.
,
Krebs
,
S.
,
Zakhartchenko
,
V.
,
Schwalb
,
B.
,
Blum
,
H.
and
Wolf
,
E.
(
2014
).
Fine mapping of genome activation in bovine embryos by RNA sequencing
.
Proc. Natl. Acad. Sci. USA
111
,
4139
-
4144
.
Hamatani
,
T.
,
Carter
,
M. G.
,
Sharov
,
A. A.
and
Ko
,
M. S. H.
(
2004
).
Dynamics of global gene expression changes during mouse preimplantation development
.
Dev. Cell
6
,
117
-
131
.
Han
,
J. Y.
(
2009
).
Germ cells and transgenesis in chickens
.
Comp. Immunol. Microbiol. Infect. Dis.
32
,
61
-
80
.
Heyn
,
P.
,
Kircher
,
M.
,
Dahl
,
A.
,
Kelso
,
J.
,
Tomancak
,
P.
,
Kalinka
,
A. T.
and
Neugebauer
,
K. M.
(
2014
).
The earliest transcribed zygotic genes are short, newly evolved, and different across species
.
Cell Rep.
6
,
285
-
292
.
Huang
,
D. W.
,
Sherman
,
B. T.
,
Tan
,
Q.
,
Kir
,
J.
,
Liu
,
D.
,
Bryant
,
D.
,
Guo
,
Y.
,
Stephens
,
R.
,
Baseler
,
M. W.
and
Lane
,
H. C.
(
2007
).
DAVID Bioinformatics Resources: expanded annotation database and novel algorithms to better extract biology from large gene lists
.
Nucleic Acids Res.
35
Suppl. 2
,
W169
-
W175
.
Huang
,
Y.-L.
,
Anvarian
,
Z.
,
Döderlein
,
G.
,
Acebron
,
S. P.
and
Niehrs
,
C.
(
2015
).
Maternal Wnt/STOP signaling promotes cell division during early Xenopus embryogenesis
.
Proc. Natl. Acad. Sci. USA
112
,
5732
-
5737
.
Jarvis
,
E. D.
,
Mirarab
,
S.
,
Aberer
,
A. J.
,
Li
,
B.
,
Houde
,
P.
,
Li
,
C.
,
Ho
,
S. Y.
,
Faircloth
,
B. C.
,
Nabholz
,
B.
,
Howard
,
J. T.
, et al. 
(
2014
).
Whole-genome analyses resolve early branches in the tree of life of modern birds
.
Science
346
,
1320
-
1331
.
Kisielewska
,
J.
,
Philipova
,
R.
,
Huang
,
J.-Y.
and
Whitaker
,
M.
(
2009
).
MAP kinase dependent cyclinE/cdk2 activity promotes DNA replication in early sea urchin embryos
.
Dev. Biol.
334
,
383
-
394
.
Kochav
,
S.
and
Eyal-Giladi
,
H.
(
1971
).
Bilateral symmetry in chick embryo determination by gravity
.
Science
171
,
1027
-
1029
.
Kochav
,
S.
,
Ginsburg
,
M.
and
Eyal-Giladi
,
H.
(
1980
).
From cleavage to primitive streak formation: a complementary normal table and a new look at the first stages of the development of the chick. II. Microscopic anatomy and cell population dynamics
.
Dev. Biol.
79
,
296
-
308
.
Lee
,
H. C.
,
Choi
,
H. J.
,
Park
,
T. S.
,
Lee
,
S. I.
,
Kim
,
Y. M.
,
Rengaraj
,
D.
,
Nagai
,
H.
,
Sheng
,
G.
,
Lim
,
J. M.
and
Han
,
J. Y.
(
2013
).
Cleavage events and sperm dynamics in chick intrauterine embryos
.
PLoS ONE
8
,
e80631
.
Lee
,
M. T.
,
Bonneau
,
A. R.
and
Giraldez
,
A. J.
(
2014
).
Zygotic genome activation during the maternal-to-zygotic transition
.
Annu. Rev. Cell Dev. Biol.
30
,
581
-
613
.
Li
,
H.
,
Handsaker
,
B.
,
Wysoker
,
A.
,
Fennell
,
T.
,
Ruan
,
J.
,
Homer
,
N.
,
Marth
,
G.
,
Abecasis
,
G.
and
Durbin
,
R.
(
2009
).
The sequence alignment/map format and SAMtools
.
Bioinformatics
25
,
2078
-
2079
.
Mak
,
S.-S.
,
Alev
,
C.
,
Nagai
,
H.
,
Wrabel
,
A.
,
Matsuoka
,
Y.
,
Honda
,
A.
,
Sheng
,
G.
and
Ladher
,
R. K.
(
2015
).
Characterization of the finch embryo supports evolutionary conservation of the naive stage of development in amniotes
.
eLife
4
,
e07178
.
Nagai
,
H.
,
Sezaki
,
M.
,
Kakiguchi
,
K.
,
Nakaya
,
Y.
,
Lee
,
H. C.
,
Ladher
,
R.
,
Sasanami
,
T.
,
Han
,
J. Y.
,
Yonemura
,
S.
and
Sheng
,
G.
(
2015
).
Cellular analysis of cleavage-stage chick embryos reveals hidden conservation in vertebrate early development
.
Development
142
,
1279
-
1286
.
Newport
,
J.
and
Kirschner
,
M.
(
1982
).
A major developmental transition in early Xenopus embryos: II. Control of the onset of transcription
.
Cell
30
,
687
-
696
.
Priess
,
J. R.
(
2005
).
Notch signaling in the C. elegans embryo
.
WormBook
,
1
-
16
.
Rengaraj
,
D.
,
Kim
,
D. K.
,
Zheng
,
Y. H.
,
Lee
,
S. I.
,
Kim
,
H.
and
Han
,
J. Y.
(
2008
).
Testis-specific novel transcripts in chicken: in situ localization and expression pattern profiling during sexual development
.
Biol. Reprod.
79
,
413
-
420
.
Robinson
,
M. D.
and
Oshlack
,
A.
(
2010
).
A scaling normalization method for differential expression analysis of RNA-seq data
.
Genome Biol.
11
,
1
.
Robinson
,
M. D.
,
McCarthy
,
D. J.
and
Smyth
,
G. K.
(
2010
).
edgeR: a Bioconductor package for differential expression analysis of digital gene expression data
.
Bioinformatics
26
,
139
-
140
.
Sheng
,
G.
(
2014
).
Day-1 chick development
.
Dev. Dyn.
243
,
357
-
367
.
Skromne
,
I.
and
Stern
,
C. D.
(
2001
).
Interactions between Wnt and Vg1 signalling pathways initiate primitive streak formation in the chick embryo
.
Development
128
,
2915
-
2927
.
Stern
,
C. D.
(
1998
).
Detection of multiple gene products simultaneously by in situ hybridization and immunohistochemistry in whole mounts of avian embryos
.
Curr. Top. Dev. Biol.
36
,
223
-
243
.
Stern
,
C. D.
(
2005
).
The chick; a great model system becomes even greater
.
Dev. Cell
8
,
9
-
17
.
Streit
,
A.
,
Berliner
,
A. J.
,
Papanayotou
,
C.
,
Sirulnik
,
A.
and
Stern
,
C. D.
(
2000
).
Initiation of neural induction by FGF signalling before gastrulation
.
Nature
406
,
74
-
78
.
Tan
,
M. H.
,
Au
,
K. F.
,
Yablonovitch
,
A. L.
,
Wills
,
A. E.
,
Chuang
,
J.
,
Baker
,
J. C.
,
Wong
,
W. H.
and
Li
,
J. B.
(
2013
).
RNA sequencing reveals a diverse and dynamic repertoire of the Xenopus tropicalis transcriptome over development
.
Genome Res.
23
,
201
-
216
.
Tao
,
Q.
,
Yokota
,
C.
,
Puck
,
H.
,
Kofron
,
M.
,
Birsoy
,
B.
,
Yan
,
D.
,
Asashima
,
M.
,
Wylie
,
C. C.
,
Lin
,
X.
and
Heasman
,
J.
(
2005
).
Maternal wnt11 activates the canonical wnt signaling pathway required for axis formation in Xenopus embryos
.
Cell
120
,
857
-
871
.
Torlopp
,
A.
,
Khan
,
M. A. F.
,
Oliveira
,
N. M. M.
,
Lekk
,
I.
,
Soto-Jiménez
,
L. M.
,
Sosinsky
,
A.
and
Stern
,
C. D.
(
2014
).
The transcription factor Pitx2 positions the embryonic axis and regulates twinning
.
eLife
3
,
e03743
.
Trapnell
,
C.
,
Pachter
,
L.
and
Salzberg
,
S. L.
(
2009
).
TopHat: discovering splice junctions with RNA-Seq
.
Bioinformatics
25
,
1105
-
1111
.
Tse
,
Y. C.
,
Werner
,
M.
,
Longhini
,
K. M.
,
Labbe
,
J.-C.
,
Goldstein
,
B.
and
Glotzer
,
M.
(
2012
).
RhoA activation during polarization and cytokinesis of the early Caenorhabditis elegans embryo is differentially dependent on NOP-1 and CYK-4
.
Mol. Biol. Cell
23
,
4020
-
4031
.
Wang
,
Q. T.
,
Piotrowska
,
K.
,
Ciemerych
,
M. A.
,
Milenkovic
,
L.
,
Scott
,
M. P.
,
Davis
,
R. W.
and
Zernicka-Goetz
,
M.
(
2004
).
A genome-wide study of gene activity reveals developmental signaling pathways in the preimplantation mouse embryo
.
Dev. Cell
6
,
133
-
144
.
Xue
,
Z.
,
Huang
,
K.
,
Cai
,
C.
,
Cai
,
L.
,
Jiang
,
C.-Y.
,
Feng
,
Y.
,
Liu
,
Z.
,
Zeng
,
Q.
,
Cheng
,
L.
,
Sun
,
Y. E.
, et al. 
(
2013
).
Genetic programs in human and mouse early embryos revealed by single-cell RNA sequencing
.
Nature
500
,
593
-
597
.
Yan
,
L.
,
Yang
,
M.
,
Guo
,
H.
,
Yang
,
L.
,
Wu
,
J.
,
Li
,
R.
,
Liu
,
P.
,
Lian
,
Y.
,
Zheng
,
X.
,
Yan
,
J.
, et al. 
(
2013
).
Single-cell RNA-Seq profiling of human preimplantation embryos and embryonic stem cells
.
Nat. Struct. Mol. Biol.
20
,
1131
-
1139
.
Yanai
,
I.
,
Peshkin
,
L.
,
Jorgensen
,
P.
and
Kirschner
,
M. W.
(
2011
).
Mapping gene expression in two Xenopus species: evolutionary constraints and developmental flexibility
.
Dev. Cell
20
,
483
-
496
.
Zhang
,
W. L.
,
Huitorel
,
P.
,
Glass
,
R.
,
Fernandez-Serra
,
M.
,
Arnone
,
M. I.
,
Chiri
,
S.
,
Picard
,
A.
and
Ciapa
,
B.
(
2005
).
A MAPK pathway is involved in the control of mitosis after fertilization of the sea urchin egg
.
Dev. Biol.
282
,
192
-
206
.
Zhang
,
G.
,
Li
,
C.
,
Li
,
Q.
,
Li
,
B.
,
Larkin
,
D. M.
,
Lee
,
C.
,
Storz
,
J. F.
,
Antunes
,
A.
,
Greenwold
,
M. J.
,
Meredith
,
R. W.
, et al. 
(
2014a
).
Comparative genomics reveals insights into avian genome evolution and adaptation
.
Science
346
,
1311
-
1320
.
Zhang
,
Y.
,
Duan
,
X.
,
Cao
,
R.
,
Liu
,
H.-L.
,
Cui
,
X.-S.
,
Kim
,
N.-H.
,
Rui
,
R.
and
Sun
,
S.-C.
(
2014b
).
Small GTPase RhoA regulates cytoskeleton dynamics during porcine oocyte maturation and early embryo development
.
Cell Cycle
13
,
3390
-
3403
.
Zhang
,
G.
,
Rahbek
,
C.
,
Graves
,
G. R.
,
Lei
,
F.
,
Jarvis
,
E. D.
and
Gilbert
,
M. T.
(
2015
).
Genomics: bird sequencing project takes off
.
Nature
522
,
34
.

Competing interests

The authors declare no competing or financial interests.

Supplementary information