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.
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.
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.
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.
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.
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).
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.
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.
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.
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
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
We thank the members of the Jae Yong Han's lab for suggestions and inspirational discussions.
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.
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).
The authors declare no competing or financial interests.