Differential transcriptome dynamics during the onset of conceptus elongation and between female and male porcine embryos

Background Porcine embryos undergo rapid differentiation and expansion between Days 8 and 12 before attaching to the maternal uterine epithelial surface after Day 13. It is known that maternal recognition of pregnancy and successful implantation are driven by mutual interactions between the elongated conceptus and the maternal endometrium. While most of the genes involved in regulation of embryo development are located on autosomal chromosomes, gene expression on sex chromosomes is modulating development through sex-specific transcription. To gain more insights into the dynamic transcriptome of preimplantation embryos at the onset of elongation and into X-linked gene expression, RNA-seq analyses were performed for single female and male porcine embryos collected on Days 8, 10, and 12 of pregnancy. Results A high number of genes were differentially expressed across the developmental stages (2174 and 3275 for Days 8 vs 10, and 10 vs 12, respectively). The majority of differentially expressed genes (DEGs) were involved in embryo elongation, development, and embryo-maternal interaction. Interestingly, a number of DEGs was found with respect to embryo sex (137, 37, and 56 on Days 8, 10 and 12, respectively). At Day 8, most of these DEGs were X-linked (96). Strikingly, the number of DEGs encoded on the X chromosome dramatically decreased from Day 10 to Day 12. Conclusions The obtained results deepen the understanding about temporary transcriptomic changes in porcine embryos during the phase of conceptus elongation, meanwhile reveal dynamic compensation of X chromosome in the female and distinct transcriptional differences between female and male embryos. Electronic supplementary material The online version of this article (10.1186/s12864-019-6044-z) contains supplementary material, which is available to authorized users.


Background
Porcine embryos undergo a morphological transformation before attaching to the uterine surface, which is supposed to maximize the intimate embryo-maternal communication. This is essential for maternal recognition of pregnancy and proper embryonic development during the preimplantation period, and thus a prerequisite for a successful pregnancy [1]. The transformation of morphology is likewise essential for embryo survival as the increased contact surface area improves the nutrient exchange between the conceptus and the uterus [2]. Embryos in the morula stage are transported into the anterior tips of the uterine horns around 60-72 h after estrus, thereafter developing into blastocysts until Day 5 of pregnancy. The blastocysts start hatching around Day 8 [3]. After hatching, porcine embryos continue to grow and form a sphere of around 2 to 6 mm in diameter on Day 10. Finally, the spherical blastocysts are transformed to tubular and then a long and thin filament approximately 100-150 mm in length on Day 12 of pregnancy [4]. Unlike the ruminant trophoblast-growth-driven elongation, the elongation from spherical to filamentous form of the pig conceptus is mainly contributed by comprehensive cellular migration and reorganization [5,6].
Along with embryo elongation, estrogens secreted from the conceptuses on Days 11 and 12 initiate the most pronounced phase of maternal recognition of pregnancy [7,8]. A number of genes from the uterine epithelium mediating cell growth, adhesion, transcription, transport, as well as prostaglandins, amino acids, and glucose synthesis are stimulated by estrogens [8,9]. Besides, the cellular signaling pathways and their corresponding receptors located in the conceptuses were found to play a role in modulating cell proliferation, movement, adhesion, and trophectoderm cells survival during preimplantation and implantation periods [3]. The developing conceptus also contributes to the accumulation of prostaglandin F2α (PGF) and prostaglandin E 2 (PGE) in the uterine luminal content [5,10,11]. According to previous studies, porcine conceptus interleukin one beta (IL1B) mRNA and protein abundance is increasing during the transformation from the tubular to the filamentous stage, reaching a maximum at the highly elongating period [12]. It is known that IL1B binds to a large number of receptors and its antagonists, together with an accessory protein, thereby regulating innate immunity and inflammation [13,14]. Furthermore, interleukin 1 receptor (IL1R1) in the endometrium activated by conceptus IL1B2 triggers a cellular signaling pathway cascade through extracellular signal-regulated kinase 1/2 (ERK1/2), mitogen activated kinase-like protein (MAPK) [15] and nuclear factors kappa-B (NFKB) [16]. Then, the NFKB activation mediates numerous biological networks including cytokines, chemokines, and prostaglandin-endoperoxide synthase 2 (PTGS2) expression [17]. However, the dynamic changes of gene expression during the preimplantation period in the conceptus are still poorly understood.
Genes located on the sex chromosomes can modulate the genome machinery resulting in differences in the development of female and male embryos. However, the sole determinant in phenotypic differences is different chromosome dosages during embryo preimplantation development [18]. In females, most genes on one X chromosome are silenced because of X-chromosome inactivation (XCI) during early development [19]. As the model of dosage compensation during preimplantation development is proved by previous studies, a reversible dynamic X inactivation may cause X-linked gene upregulation in female embryos [20]. The molecular events that cause differences between male and female embryos were observed in in vitro culture experiments including embryo developmental speed, blastocyst cell number, and metabolism [21]. In a global gene expression study of bovine Day 7 embryos, 193 X-linked transcripts were upregulated in female compared to male embryos, suggesting that XCI is partially achieved at the blastocyst stage [21]. To date, this is not clear in porcine embryos.
Considering the limited knowledge in pigs in terms of transcriptome changes of preimplantation embryos during the onset of elongation and respective differences between female and male embryos, this study was designed to analyze global gene expression in individual porcine embryos during the preimplantation period (Day 8, 10, and 12 of pregnancy) with the aims of (i) characterizing the dynamics of gene expression profiles Fig. 1 Principal component analysis of all genes expression patterns of embryos. The multidimensional scaling plot displays that embryos from the same stage are grouped together. Samples from the same stage with the same symbol (square, round, triangle shape indicate samples of Day 8, 10, and 12, respectively). Red and blue colors indicate female and male embryos, respectively and the involved functions and thereby (ii) investigating X-linked gene expression between female and male embryos.

RNA sequencing of porcine embryo samples
At least two embryos (one female and one male) were selected from the same sow in each stage and they were mixed within the same treatment. The embryos size collected on Days 10 and 12 ranged from 2.2 to 3 mm, and from 45 to 200 mm, respectively. RNA-seq libraries were prepared from 30 individual embryos (5 female and 5 male embryos each on Days 8, 10, and 12). The number of raw reads (100-bp single-end reads) ranged from 19 to 41 million per sample. From these, 18 to 39 million clean reads were obtained after quality filtering and adapter clip. All data used in this study have been included in the article and its Additional files. The sequence data (GSE113366) is available at National Center for Biotechnology Information (NCBI) Gene Expression Omnibus. The principal component analysis for all detected genes (13103) in these embryos ( Fig. 1) revealed that embryos with the same days of pregnancy were grouped together. It was notable that embryo from Day 12 showed a higher distance with embryos from Days 8 and 10, besides the variation between embryos on Day 10 was highest compared with the rest embryos.
Differentially expressed genes across the stages during the preimplantation phase Based on the analysis of DEG (FDR 0.1%, log2 fold change > 1), the total number of 2174 DEG were found between Days 8 and 10, and even more genes (3275 DEG) were differentially expressed between Days 10 and 12 (Fig. 2a). A number of 1072 DEG were shared by these two comparisons and many of them shared a similar expression profile. The DEG of each comparison were combined (4377 DEG) for analysis of typical expression profiles across the developmental stages. Hierarchical cluster analysis provided an overview of the expression profiles for all the combined DEG (Fig. 2b). To have a closer look into the regulation of these DEGs, the genes with similar expressions at the same stage were identified by Self-Organizing Tree Algorithm (SOTA) analysis (Fig. 3), which and Day 10 vs Day 12 are shown for genes with significant differences in gene expression (FDR < 0.1% and │log2 fold change│ > 1). The color scale is from −2 (blue, lower than mean) to 2 (red, higher than mean). Each row represents 1 gene, each column 1 sample revealed the expression of almost half of the genes was similar on Days 8 and 10 but increased on Day 12 (genes of cluster 1 and 5). Inversely, a number of genes showed increased expression on Days 8 and 10 but decreased expression on Day 12 (cluster 3). Furthermore, clusters of genes were found that showed increased expression only on Day 8 (contained in cluster 3) and only on Day 10 (cluster 2), respectively. In cluster 4, the gene expression dramatically increased on Day 10, and kept a stable upregulated expression level on Day 12. The DEGs across the three stages are shown in Additional file 1: Table S1.

Functional category analysis of gene expression profile clusters
A functional annotation clustering analysis was performed with the Database for Annotation, Visualization, and Integrated Discovery (DAVID) tool for the individual SOTA clusters. Highly significant functional annotation clusters (geometric mean of member categories' P-values < 0.001 which was calculated based on Fisher's exact test) for each of the 5 SOTA expression profile clusters are shown in Table 1. Strong enrichment of functional categories related to cell communication, cell motility and migration, cell adhesion and junction, cytoskeleton organization, embryo development, apoptosis, and immune response was found for genes specifically upregulated on Day 12 (clusters 1 and 5). For genes up-regulated genes on Day 10 (cluster 2), lipid biosynthetic and metabolic process, transport, developmental and homeostasis were enriched. Genes upregulated on Days 8 and 10 (cluster 3) revealed overrepresented functional categories such as vacuole, lipid metabolic and biosynthetic process, cell-cell junction, and ion transport. For cluster 4, containing genes with increased expression on Days 10 and 12 functional categories involved in development and morphogenesis, cell-cell signaling, and apoptotic process were overrepresented. More detailed information is shown in Additional file 2: Table S2.
An overview of a network of enriched functional classification categories for the obtained DEG of the SOTA clusters is shown in Fig. 4. The DEG from these five clusters were analyzed with the online tool ToppCluster for Gene Ontology (GO) and pathway analysis. The most overrepresented categories including response to estrogen, cell migration, actin filament-based process, vasculature development, as well as animal organ morphogenesis were shared in most of the SOTA expression profile clusters (Fig. 4). In addition, there were also functional categories specifically overrepresented in cluster 1, e.g., interferon signaling, extracellular matrix organization, cytokine signaling in immune system, heparin binding, growth factor binding, and cell adhesion molecule binding. Cluster 3 was enriched for transport and cluster 4 in actin cytoskeleton organization, developmental growth, epithelial tube formation, and cytoskeletal protein binding. In addition, neuron development was specifically enriched in cluster 5.
Top 5% functional annotations of upregulated genes from the embryos on Day 12 To explore more details of the upregulated genes on Day 12, our data were subjected to the online tool DAVID (Gene Ontology (GO) categories and KEGG pathways) to show the systematic functional classification for upregulated genes in the embryos from Day 12   compared with Day 10. The functional categories were selected based on the condition of FDR < 1%, then filtered by a score combining FDR and fold enrichment. The overrepresented GO terms and pathways with the top 5% scores were selected and the results were shown in Table 2.
The categories related to 'morphogenesis' and 'development' were disclosed to be the main functions on Day 12.
Signaling pathways related to maternal recognition of pregnancy on Day 12 In the embryos of Day 12, a number of 164 DEGs assigned to a selection of particularly interesting pathways and processes, e.g., estrogen signaling, steroid hormone biosynthesis, prostaglandin (PG) metabolism, signaling and transport are shown in Additional file 3:  Fig. 5a and b, respectively.

Comparison of DEGs in the elongating conceptus between pig and sheep
Data of dynamic transcriptome changes from ovoid to a filamentous conceptus in sheep [22] have been adopted to compare with our findings (from Day 10 to Day 12 in pigs, and from Day 12 to Day 14 in sheep). The results show that the majority of the DEGs (1829) were identified as upregulated specifically in porcine embryos (Fig. 6), and 1025 gene were identified as downregulated. A number of 1093 downregulated genes and 795 upregulated genes were specifically expressed in sheep embryos. Besides, 133 upregulated and 77 downregulated genes were commonly expressed during the elongation in both species, while 209 genes had different regulations in pig and sheep embryos. The detailed information of these DEGs is depicted in Additional file 4: Table S4.
Analysis of expression of genes located on the X chromosome in the three embryo stages At an FDR of 5%, a number of 137, 37, and 56 differentially DEG between female and male embryos were obtained on Days 8, 10, and 12, respectively. A hierarchical cluster analysis of these DEG is shown in Fig. 7a, b, and c. Gene Ontology (GO) functional classification network of clusters. All significant differential genes (human Entrez Gene IDs) from clusters were used as input for the ToppCluster. The following databases were used, i.e. "biological process", "cellular component", "molecular function" and pathway. Finally, the data were uploaded in Cytoscape 3.6.0 to modify the network. Nodes were colored based on specificity: red nodes specific for different cluster; nodes for the three GO functions and pathway were in different colors With respect to higher or lower expression of DEG in female compared to male embryos, 107 upregulated and 30 downregulated, 23 upregulated and 14 downregulated, and 31 upregulated and 25 downregulated genes were identified on Days 8, 10, and 12, respectively (Additional file 5: Table S5). Regarding chromosomal location of these DEG, a number of 96, 19, and 11 X-linked DEG were found on Days 8, 10, and 12, respectively (Fig. 8a). The number of DEG located on the X chromosome decreased on Day 12 compared with Day 8 embryos (20% vs 70% of the DEG). Four X-linked genes, namely B-cell receptor-associated protein 31 (BCAP31), family with sequence similarity 3 member A (FAM3A), ribosomal protein L10 (RPL10), and tafazzin  5 The networks of DEGs involved in estrogen (a) and prostaglandin (b) signaling pathways. The sources of co-expression, physical interactions, co-localization, pathway, shared protein domains, and genetic interactions were selected to weight the network. One edge indicates one source, and one node means one gene for the female and male embryos are shown with significant differences in gene expression (FDR < 5%). Color scale is from − 2 (blue, lower than mean) to 2 (red, higher than mean). Each row represents one gene, each column represents one embryo (TAZ), were upregulated in female embryos on all 3 days (Additional file 5: Table S5). All DEG located on the X chromosome were upregulated except three genes, namely family with sequence similarity 155 member B (FAM155B), collagen type IV alpha 5 chain (COL4A5) and spindlin-2B (LOC100526148), that were downregulated in female embryos on Day 8 (Fig. 8b). More detailed information about the DEG between female and male embryos for each stage is shown in Additional file 5: Table S5.

Functional annotations for sex-based DEGs and X-linked DEGs in the embryos
The functional classification based on the DEGs from the comparison between female and male embryos was identified at each stage (Table 3). A number of 22 functional categories were identified on Day 8. Among these categories, genes involved in "phosphatidylinositol-3,5-bisphosphate 3-phosphatase activity" had the highest enrichment. Interestingly, on Days 10 and 12, only one category, the "translational initiation", was enriched in both two stages. X-linked DEGs were analyzed to reveal the related functional categories ( Table 4). The majority of functions were similar as the categories from the sex-based DEGs on Day 8. However, no functional categories with significant fold enrichment from the X-linked genes were identified on days 10 and 12.

Discussion
Prior to attaching to the maternal uterine epithelial surface, porcine embryos normally undergo programmed differentiation and expansion. This reflects the rapid morphological changes between Days 8 and 12, which results from various hormonal and gene expression changes. The transcriptome of the embryo either in autosomal or sex chromosomes is hypothesized to be dynamic, but only little evidence supports this hypothesis. By using RNA-Seq, the present study analyzed embryos sampled at three critical time points and revealed a range of phase-specific gene expressions regulating embryo elongation, development and communication with the maternal uterus. These events are critical for successful implantation in pigs.

Genes related to conceptuses elongation (cell growth, cell movement, and cellular remolding)
The current results revealed that the elongating conceptuses under investigation displayed various dynamic gene expression changes most highly expressed on Day 12, especially relating to cell growth, cell movement as well as cellular remolding.
A number of growth factors were found to have increased gene expression on Day 12, including insulin growth factor (IGF) family, fibroblast growth factor (FGF) family, transforming growth factor beta 1 (TGFβ1), transforming growth factor beta 2 (TGFβ2), platelet derived growth factor receptor alpha (PDGFRA), insulin like growth factor binding protein 5 (IGFBP5), and fibroblast growth factor receptor 2 (FGFR2). These results confirm a related study in mouse, where a range of growth factors were produced by the preimplantation embryo and the reproductive tract, and many of the respective receptors were detected on the embryo surface [23]. The TGFβ was regarded as an essential factor in Fig. 8 Venn diagram for X-linked differentially expressed genes (a) and expression profile of X-linked differentially expressed genes (b). Green, blue, and pink filed in Venn graph indicate X-linked differentially expressed genes on Day 8, 10, and 12 of the embryos, respectively. The red and blue plots in the scatter mean the up and down-regulated genes on the X chromosome in female embryos, respectively modulating extracellular matrix (ECM) protein expression and the accumulation of ECM proteins results in increasing cell migration [24]. Expression of TGFβ1 and TGFβ3 were observed to be unregulated in porcine conceptuses on Day 10 compared with Day 8, and mRNA of TGFβ1, TGFβ2, and TGFβI were still increasing except TGFβ3 was decreasing from Day 10 to Day 12. From our findings, the increasing TGFβ may contribute to some extent to morphological changes of the embryo via migrating cells and remolding tissues. Growth differentiation factor 6 (GDF6) was observed with the highest expression on Day 12. Recently, GDF6 gene was detected in several distinct embryonic locations in mouse, which facilitates the skeletal and soft tissues formation [25]. We found placenta expressed transcript 1 (PLET1) upregulated on Day 12, which is in line with the previous result that PLET1 was found with high expression in elongated conceptuses [26], indicating a strong activity of trophoblast differentiation [27].
Coordination between assembly and disassembly of actin filaments is a key factor to provide a driving force that initiates the cell movement in animals [28]. It is reported that the cellular motility in vitro can be reconstituted by a core set proteins including actin, actin related protein2/3 (Arp2/3) complex, actin depolymerizing factor (ADF)/cofilin, profilin, and capping protein [29]. Consistently, the current study observed that genes coding actin proteins, actin binding LIM protein 1 (ABLIM1), actin related protein 2/3 complex subunit 1B (ARPC1B), cofilin 2 (CFL2), profilin 1 (PFN1), and capping actin protein gelsolin like (CAPG) were differently expressed from Days 8 to 12. With reference to porcine embryos, the filamentous actin cytoskeleton reassembly is rather outstanding in the early pregnant stage implying that compaction and blastogenesis are activated during the preimplantation period [30]. Herein, we report that the gene coding for cytoskeleton associated protein 4 (CKAP4) revealed highest expression on Day 12 suggesting that stage-specific changes in actin organization occur between Days 8 and 12 of pregnancy. This is consistent with the formation of cytoskeletal elements occurring during mouse embryo preimplantation stage [31]. On the other hand, the low expression of genes coding for two major cytoskeletal proteins (desmin (DES) and vimentin (VIM)) was observed. The expression of VIM was greater in Day 12 conceptuses compared to Day 8 and 10 conceptuses, which is in line with previously observed VIM protein patterns indicating that mesodermal differentiation and migration activities are greater during the filamentous stage [32][33][34]. Though DES and VIM proteins were confirmed to be expressed dynamically, their transcripts were kept a low level based on the present results. Though DES and VIM regulating filaments were found not to be essential for embryo development, the cytokeratin formation is an initial step for the epithelium differentiation [31]. A DES and VIM co-expression may exist to make up the insufficient filament protein production in developing embryonic tissue [35]. The majority of genes regulating the cell differentiation, growth and movement were expressed highly on Day 12. As indicated in Fig. 2b, the total number of activated genes increased gradually along developmental progression suggesting that Day 12 of pregnancy is a critical period with more pathways and functions involved in embryo elongation. Meanwhile, the expressions of some genes reached the peak on Day 10, like KRT, MDK that are responsible for blastocyst formation [36] and embryo survival [37].

Genes related to embryo development
It is known that the transformation of morphology is essential for embryo survival as the increased placental surface area can improve the nutrient exchange between the conceptus and the uterus [2]. Lipids are not only components of cellular membrane and cytoplasm but also of great importance in producing energy that is essential for proper embryo development [38]. Genes related with lipid biosynthetic process, such as lipoprotein (LPL) and acyl-CoA synthetase long chain family member 3 (ACSL3), were found to be expressed in the present study, especially fatty acid binding protein 3 (FABP3) and ELOVL fatty acid elongase 6 (ELOVL6), which reached a peak on Day 10. As reported previously, the content of lipids in pig embryos was markedly decreased in the late blastocyst stage [38], which proves that the energy for embryo development from Day 8 to Day 10 is provided by lipid synthesis and metabolism. As the embryo achieves a more advanced development from Day 10 onwards, genes involved in the category of energy production were rationally downregulated on Day 12.
Significant different gene expression of KRT8 and KRT18 were found between bovine morula and blastocyst embryos indicating their roles in blastocyst formation and embryo implantation [39]. The latter, KRT18 mRNA and protein were confirmed to be differentially expressed in bovine embryo development and it was regarded as a marker for blastocyst formation [36]. The present findings were in line with a previous report that expression of KRT18 mRNA was decreased in the filamentous conceptus compared with the ovoid conceptus, implying that KRTs play a role in trophectoderm development [40]. The neurite growth-promoting factor, midkine (MDK) protein is a developmentally regulated heparin binding cytokine that is induced by retinoic acid [41,42]. The MDK transcripts reached at a peak in the Day 10 conceptus of our study suggesting that it may be associated with extra-embryonic tissue development from Day 8 to Day 10 [40] rather than the tissue remodeling from Day 10 to Day 12.
The genes under similar regulations between porcine and ovine elongating embryos (as shown in Fig. 6), revealing the conserved DEGs in two species. While the opposite regulated genes from these two species may likely indicate the distinct pathways that facilitate embryo elongation in pig different from sheep. The genes involved in nutrient requirements displayed stage-specific expression. For example, Day 10 embryos required more energy reflected by the high gene expression coding for lipid and energy synthesis. In addition, novel and stage-specific genes transcripts were discovered regarding estrogen and embryo development.

Genes related to embryo-maternal interaction
Seven integrin alpha genes (ITGA2, ITGA3, ITGA5, ITGA6, ITGA9, ITGA10, and ITGAL) and five integrin beta genes (ITGB2, ITGB3, ITGB5, ITGB6, and ITGB8) were notably differentially expressed. A higher abundance of ITGA3, ITGB6, and ITGB8 was identified on Day 10 compared with Day 8. However, expression level of nine integrin genes (including ITGA2, ITGA3, ITGA6, ITGA10, ITGAL, ITGB2, ITGB3, ITGB5, and ITGB6) higher abundant at Day 12 compared to Day 10. Integrins have been recognized as cell surface adhesion receptors [43]. According to the network of ligand receptor-mediated multicellular signaling from the FAN-TOM5 database of human cells [44], the ligand of integrin has been located in the endometrium on Day 12 of pregnancy in pigs [45,46]. The expression of integrins with the detection of their ligands indicates the role of cell adhesion taking place in the maternal-embryo crosstalk, especially on Day 12, a time with an intensive embryomaternal interaction.
A number of genes coding gap junction protein (GJA1, GJB1, GJB3, GJB5) were highly expressed, especially on Day 12. Gap junction is a complex of proteins locating distinctly on the plasma membrane of bordering cells to establish the intercellular connection [47]. The gene coding desmosome associated protein, PNN, was likewise among the highly expressed transcripts in the present study. Desmosomes were reduced in the mouse uterine epithelium during the preimplantation revealing that the reduction in desmosome adhesion contributes to embryonic-penetration in luminal epithelium [48]. In contrast to the invasive implantation in mouse, a high expression of PNN was observed in all three stages of development under investigation, which leads us to assume that there was no penetrating process due to the specific non-invasive implantation in pigs. In our study, genes coding tight junction proteins (TJP1, TJP2), and TJP3 were upregulated, which are produced by embryonic cells during preimplantation period and act as a tight junction for polarized transporting, and intercellular integrity and signaling [49].
The expression of IL1B in the embryo during the preimplantation was reported with dynamic changes, however, the pattern currently seems to be controversial. Ross et al. reported that the expression of IL1B was increasing from Day 10 to Day 12 along with the embryo elongation [12]. Another related study showed that the level of IL1B mRNA of porcine conceptus was decreased on Day 12 when compared with Day 11 [50]. Similar to the later, we noticed that transcripts of IL1B mRNA reached a peak on Day 10 compared to Days 8 and 12. It has been shown that IL1 was induced by estrogen and progesterone, while high levels of estrogen and progesterone could inversely inhibit IL1 secretion [51,52]. Therefore, we propose that from Days 8 to 10, the increasing secretion of estrogen by the elongating conceptuses most likely stimulates IL1B expression. However, the large amount of accumulated estrogen on Day 12 thereafter inhibit ILIB transcription. IL1B has a unique form in pig conceptuses, namely IL1B2, which was highly expressed on Day 10. Though interleukin 1 receptor type 1 (IL1R1) and interleukin 1 receptor accessory protein (IL1RAP), receptors for the IL1B-signaling system, were found low expression from Day 8 to 12. The IL1B of the embryo together with IL1R1 and IL1RAP initiate the IL1B-signaling system plays an important role in implantation by regulating expression of endometrial genes for prostaglandin (PG) synthesis [53].
Along with morphological changes, expression patterns indicate the intimate dialogue between endometrium and embryo. Over time, the communication intensity increases indicated by the increased expressions of various genes involved in the cell surface adhesion and gap junction.

The maternal recognition of pregnancy on Day 12
The rapidly elongated embryo secretes large amount of estrogen on Day 12, which initiates the maternal recognition of pregnancy in pigs. Expression of estrogen receptor has been described in the endometrium of sows during early pregnancy [54]. Cytochrome P450 family members are involved in estrogen synthesis and metabolism [55,56]. In the current study, a number of genes coding specific CYP members, i.e. CYP11A1, CYP17A1, CYP19A1, CYP2C9, CYP2C18, CYP2S1, and CYP4F22 were identified as DEGs specifically on Day 12. Hydroxysteroid 17-beta dehydrogenase 1 (HSD17B1) catalyzes cholesterol to estrogen in endometriosis [57], and is massively expressed in the trophectoderm rather than in the embryonic disc of the porcine conceptus tissue on the day 12 [58]. Compared with Day 10, the expression of HSD17B1 was decreased on Day 12, suggesting that HSD17B1 secreted from the trophectoderm was not the most important component in the estrogen synthesis. Phospholipase A2 family members (PLA2s), the key enzymes for the release of PG precursor molecules [59], were upregulated in the embryos on Day 12. In contrast, endometrial PLA2s were described as downregulated on Day 12 in pigs [46]. The differential expression of PLA2s in embryos and endometrial cells point towards the regulation of the corpus luteum maintenance via PG. The transporter of PGs, the solute carrier organic anion transport family member 2A1 (SLCO2A1) mRNA was upregulated in the embryos on Day 12 compared with Day 10. The large amount of estrogen secreted from the embryos could stimulate the endometrium to switch the secretion of PG from endocrine to exocrine. Taken together, the genes involved in estrogen secretion and regulation of PG synthesizing, metabolizing and transporting proteins revealed a complex pattern of maternal recognition of pregnancy.

DEGs between female and male
Dosage compensation for X-linked genes occurs by inactivating one X chromosome in the female during embryonic development, the event of which is highly conserved in the mammal and essential for embryogenesis [60]. The current findings revealed a number of differentially expressed genes between female and male pig embryos, namely 137, 37, and 56 on Day 8, 10, and 12, respectively. Distributing the gene location on the chromosomes, 100, 19, and 11 genes were X-linked DEG on the respective days. In bovine blastocysts, more DEG, namely 382 genes, were detected between female and male embryos, from which 218 genes were X-linked [21]. In previous study, around 600 DEG were discovered in between female versus male mouse blastocysts using transgenesis and microarray analyses [61]. Thus, less X-linked genes were detected in pigs than mouse and cows. Sex chromatin as a part of the X chromosome was first detected at the blastocyst stage with 45 cells in porcine embryos, which was an approximate guide to the presence of an inactive chromosome [62]. This could also be a reason that XCI started partially before Day 8 in porcine embryos. Additionally, the majority of the upregulated genes in female embryos were X-linked on Day 8, similar to the results observed in bovine [21]. Despite that parts of the X-linked genes were not differentially expressed before Day 8, our results still provide the evidence that 81 X-linked gene expression are changed between male and female embryos from Day 8 to Day 10, and eight from Day 10 to Day 12. Therefore, the data allow the assumption that XCI probably proceeds gradually and causes the decreasing of X-linked DEG from Day 8 to Day 12.
Prior to the initiation of XCI, all X-linked genes should be present as 2-fold dose in females [63]. In fact, most X-linked genes determined here were less than two-fold expressed in female embryos on Day 8 compared to males. The housekeeping X-linked gene glucose-6-phosphate dehydrogenase (G6PD) has been considered as a potential candidate involved in sex difference [64]. In our findings, G6PD transcripts were high abundant, but only 1.5-fold in female embryos on Day 8 compared to males. Moreover, G6PD was not differentially expressed on Days 10 and 12. Previous studies reported that inactivation of paternal X-linked genes was often incomplete, which is associated with the specific locus, species, and tissue, respectively [65,66]. In addition, it is known for many years that sex chromatin and late replication of a X chromosome do not occur in all cells of the embryo at the same time [67]. This is likely the reason for most genes on X chromosomes being differentially expressed on Day 8 and no difference observed on Days 10 and 12. Functional annotations for the sex-based DEGs and X-linked DEGs revealed that the majority of functional categories were attributed by the DEGs on the X chromosome rather than the autosome on Day 8. Neither on Day 10 nor Day 12, significant functional categories were detected, indicating that the compensation on the X chromosome has balanced the difference between female and male embryos. Overall, based on the comparison of X-linked genes between female and male preimplantation embryos, the number of DEG on X chromosome were decreasing along with the progression of pregnancy. These finding suggested that the compensation between female and male has occurred before Day 8, and would persist until Day 12 with most X-linked genes expressed with similarity between male and female embryos.

Conclusion
In conclusion, we identified comprehensive transcriptome changes associated with embryo elongation, development, and embryo-maternal interaction during the preimplantation period. A number of biological processes and pathways with temporal changes were revealed governing the embryonic cell movement and remodeling to form the elongated embryo. Genes involved in cell communication and adhesion were highly expressed on Day 12, which indicates the increased interaction between the mother and embryo. Bioinformatics analyses of gene expression between female and male embryos showed that a number of X-linked genes were differentially expressed on Day 8. These DEGs disappeared gradually along with the embryo elongation on Days 10 and 12. The latter findings suggest that dynamic changes of transcriptome on the X chromosome may reveal the changed dosage compensation between sexes before embryo implantation in pigs.

Methods
The animal experiment was conducted in Freising, Germany and approved by the District Government of Upper Bavaria and were in accordance with the accepted standards of humane animal care in Germany.

Animal experiment
The animals were owned and kept at the Research station "Thalhausen" of the Technical University of Munich, Freising, Germany. Twelve German Landrace × Pietrain crossbred gilts were cycle synchronized using Altrenogest ReguMate® for twelve days, then Intergonan® (PMSG) at 750 iU was applied once in the following evening, and Ovogest® (human chorion gonadotropin) at 750 iU was applied once 3.5 days later. The next day (day 0), all animals were inseminated with sperm of the same single Pietrain boar. The animals were slaughtered in a commercial slaughterhouse. Four gilts were slaughtered by stabbing in the neck for bloodletting after anesthesia each on Days 8, 10, and 12 post insemination, respectively. The respective gilts were randomly assigned to the date of slaughtering. The reproductive tracts were collected immediately after slaughter.

Embryo collection
The complete embryos (including embryonic disc and trophectoderm cells) were flushed from the uterus using 10 ml of phosphate-buffered saline (PBS, pH 7.4) per horn, and then were transferred into a petri-dish containing PBS for morphological observation. All embryos were then individually shock-frozen in liquid nitrogen immediately and stored at − 80°C for further DNA and RNA analyses.

DNA sex determination and RNA-seq library preparation
All embryos underwent total DNA and RNA extraction using the AllPrep DNA/RNA Micro Kit (Qiagen, Valencia, CA, USA) according to the manufacturer's recommendation. The integrity and the quantity of the RNA were assessed on the Agilent 2100 Bioanalyzer (Agilent Technologies, Waldbronn, Germany). The DNA samples were used for sex determination by measuring relative gene expression with quantitative real-time PCR (qPCR). Gene-specific primers (SRY and Histone genes) were designed with online tool NCBI Primer-BLAST (see Additional file 6: Table S6). The amplification products were obtained from qPCR with reaction mixture including 2 × FastStart Universal SYBR Green Master 10 μl, DNA 1 μl, forward primer 0.6 μl, reverse primer 0.6 μl, and add water to 20 μl. The following qPCR program: 95°C for 10 min, amplification for 40 cycles at 95°C for 15 s and 60°C for 60 s, and ΔCt value method was used for relative quantification.
We then randomly selected 5 embryos per time point and sex of high RNA quality (RIN > 8) for RNA-seq. After RNA extraction, 100 ng total RNA was used as input for the Illumina TruSeq Stranded mRNA library construction (Illumina Inc., San Diego, CA, USA) regarding the manufacturer's recommendation. Briefly, RNA was fragmented and random primers were hybridized for cDNA synthesis. The resulting cDNA was followed by 3′ adenylation and adapter ligation. Finally, PCR amplification was performed with the subsequent protocol (98°C for 30 s; 15 cycles of 98°C for 10 s, 60°C for 30 s, 72°C for 30 s; 72°C for 5 min). Library quality control was performed with Agilent 2100 Bioanalyzer (Agilent Technologies, Waldbronn, Germany), followed by sample pooling. Finally, pooled libraries was sequenced as 100-bp single-end on the HiSeq 2500 (Illumina Inc.) platform.

RNA sequencing data analysis
Analysis of RNA-seq data was performed using a locally installed version of Galaxy [68]. Raw reads were quality trimmed with Trim Galore and 3 bp were removed from the 5′ end of each read. All sequences were mapped to the Sus scrofa genome (version 11.1) from NCBI (ftp://ftp.ncbi. nih.gov/genomes/Sus_scrofa/GFF) and filtered by CPM cutoff. The resulting read count table with CPM was used for statistical analysis in EdgeR to identify differentially expressed genes (DEG) by using GLM robust method [69]. The DEG between the female and male embryos were filtered with a false discovery rate (FDR) at 5%, and DEG between the different stages were filtered with FDR < 0.1% and │log2 fold change│ > 1. Hierarchical cluster analysis was performed for preliminary assessments of the number of gene clusters in MultiExperiment Viewer (MeV) [70], then the list of DEGs from three stages were used for Self Organizing Tree Algorithm (SOTA) clustering. Gene ID of DEGs in each cluster were uploaded into DAVID for gene ontology (GO) and pathway analyses [71]. DEGs from individual clusters were analyzed with the online tool Toppcluster (http://toppcluster.cchmc.org) to generate the overview networks between the GO categories and pathways and the p-Value cutoff was set at 0.05 then improved on Cytoscape to show the final results [72]. The network of the DEGs involved in the estrogen and prostaglandin signaling pathway was analyzed on online tool GeneMANIA.