Identification of Signaling Pathways for Early Embryonic Lethality and Developmental Retardation in Sephs1−/− Mice

Selenophosphate synthetase 1 (SEPHS1) plays an essential role in cell growth and survival. However, the underlying molecular mechanisms remain unclear. In the present study, the pathways regulated by SEPHS1 during gastrulation were determined by bioinformatical analyses and experimental verification using systemic knockout mice targeting Sephs1. We found that the coagulation system and retinoic acid signaling were most highly affected by SEPHS1 deficiency throughout gastrulation. Gene expression patterns of altered embryo morphogenesis and inhibition of Wnt signaling were predicted with high probability at E6.5. These predictions were verified by structural abnormalities in the dermal layer of Sephs1−/− embryos. At E7.5, organogenesis and activation of prolactin signaling were predicted to be affected by Sephs1 knockout. Delay of head fold formation was observed in the Sephs1−/− embryos. At E8.5, gene expression associated with organ development and insulin-like growth hormone signaling that regulates organ growth during development was altered. Consistent with these observations, various morphological abnormalities of organs and axial rotation failure were observed. We also found that the gene sets related to redox homeostasis and apoptosis were gradually enriched in a time-dependent manner until E8.5. However, DNA damage and apoptosis markers were detected only when the Sephs1−/− embryos aged to E9.5. Our results suggest that SEPHS1 deficiency causes a gradual increase of oxidative stress which changes signaling pathways during gastrulation, and afterwards leads to apoptosis.


Introduction
Selenium (Se) is an essential trace element required in the diet of humans and other forms of life. An adequate amount of selenium is essential to good health. For example, Se is implicated in cancer prevention, antiviral response, boosting the immune system, male reproduction, and embryo development [1][2][3]. In addition to those beneficial effects of Se, it is notable that Se also affects the progression of pregnancy. Specifically, Se levels and the activity of blood glutathione peroxidase were lower than average in women who had experienced miscarriage, premature birth, or preeclampsia [4,5]. Moreover, pregnant women who had low serum Se levels showed a high incidence of spontaneous miscarriage [6].
Se is the only trace element to be specified in the genetic code, and selenocysteine (Sec) is the 21st amino acid, which is incorporated into protein in response to the UGA codon during translation [7]. Selenoproteins contain single or multiple Sec residues in their active sites and are known to carry out important roles, often with the beneficial effects of Se described above [1,2]. One of most common functions of selenoproteins is to protect cells and tissues from oxidative stress by serving as reactive oxygen species (ROS) scavengers [8,9]. During Sec synthesis, selenophosphate serves as the selenium donor. Selenophosphate synthetase (SEPHS) catalyzes the reaction of selenophosphate synthesis, in which inorganic selenium and ATP are used as substrates [10]. There are two types of SEPHSs in higher eukaryotes, SEPHS1 and 2 [11]. Both isotypes have ATP binding and catalytic domains, and there is a high amino acid sequence homology between them. One of the biggest differences between SEPHS1 and 2 is that SEPHS2 contains Sec in the catalytic domain, while SEPHS1 does not. Instead, SEPHS1 contains threonine at the position corresponding to Sec in SEPHS2 [12]. Another feature is that only SEPHS2 has the ability to synthesize selenophosphate [13]. According to data from the International Mouse Phenotyping Consortium, whole-gene deletion of Sephs2 did not show embryonic lethality but showed abnormalities in heart morphology of the knockout fetuses [14]. SEPHS1, however, is required for cell survival and proliferation [15].
The most prominent role of SEPHS1 is that it participates in the regulation of cellular redox homeostasis [3]. In Drosophila, a P-element insertion mutation in Sephs1 (SelD) led to embryonic lethality following the loss of imaginal disc formation [16]. Subsequently, it was demonstrated that the embryonic lethality in Drosophila is mediated by ROS-induced apoptosis [17]. An in vitro study using SL2 cells, an embryonic cell line of Drosophila, showed that SEPHS1 deficiency induced ROS accumulation which in turn led to the inhibition of cell proliferation and glutamine-dependent megamitochondria formation [18]. In addition to studies in Drosophila, systemic knockout mice targeting Sephs1 showed embryonic lethality and complete resorption by E14.5 [19]. Knockdown of Sephs1 in a mouse embryonic cancer cell line (F9) showed the inhibition of cell proliferation by increased levels of ROS, specifically hydrogen peroxide [19]. Deficiency of SEPHS1 in F9 cells also reversed cancer malignancy characteristics such as cell invasion and anchorage independence. The expression levels of glutaredoxin 1 (Glrx1) and glutathione-s-transferase O1 (GstO1), which are involved in redox homeostasis, were significantly decreased by SEPHS1-deficiency in F9 cells.
Oxidative stress causes numerous types of damage during embryonic development by altering cellular macromolecules such as lipids, proteins, and nucleic acids [20]. Consequently, an affected embryo undergoes growth inhibition, development retardation, metabolic dysfunction, and apoptosis [21]. In the case of embryos cultured in vitro, high ROS levels have detrimental effects on embryo growth, but the addition of free radical scavengers recovers cells from the detrimental effects of ROS [22]. The importance of ROS during development was also demonstrated by regulating the expression of antioxidant genes. For example, disruption of thioredoxin 1 (Txn1) resulted in embryo hatching failure and lethality shortly after implantation [23]. Thioredoxin 2 (Txn2) mutation inhibited neural tube formation and induced massive apoptosis at E10.5 [24].
Although it is known that SEPHS1 plays an essential role for cell survival and proliferation, the underlying molecular and biochemical mechanisms have not been fully clarified. To elucidate the role of SEPHS1 during early embryogenesis, we generated systemic knockout mice targeting Sephs1. Pathways and genes regulated by SEPHS1 were predicted using various bioinformatical tools described herein, and the predicted pathways were verified by analyzing the anatomical structure of the developing embryo.

RNA-seq Data Analysis
We previously reported that the systemic knockout targeting Sephs1 in a mouse model resulted in the embryo beginning to show a difference in size at E7.5 and lethality at E11.5 [19].
In this study, the effect of SEPHS1 deficiency on early embryogenesis was analyzed in more detail during this period. Figure 1A shows whole-embryo images obtained by optical microscopy. There are no differences between wild-type (Sephs1 +/+ and Sephs1 +/− ) and the Sephs1 −/− embryo at E6.5. The Sephs1 −/− embryo began to show size differences from wildtype at E7.5, and the difference extended at E8.5, wherein head folds were observed to be less developed than in wild-type. In addition, the Sephs1 −/− embryo does not turn in the final fetal position at E9.5 and is smaller than wild-type. Unlike in the wild-type embryo, optic and otic vesicles were not observed, although the allantois still remained, and three primary brain vesicles (prosencephalon, mesencephalon, and rhombencephalon) were observed as being immature. At E9.5, no vitelline vessel was found in the yolk sac (arrow in Figure S1A). At E10.5, the size of the Sephs1 −/− embryo was dramatically reduced from that at E9.5 (data not shown) and the embryo appeared to be fully absorbed at E11.5 (arrowhead in panel E11.5 of Figure 1A). We used 32,24,34,20,17, and 10 embryos at E6.5, E7.5, E8.5, E9.5, E10.5, and E11.5, respectively (See Materials and Methods for more detailed information). All the embryos that had the same genotype and were prepared at the same embryonic day showed similar phenotypes both in size and in morphology. We previously reported that the systemic knockout targeting Sephs1 in a mouse model resulted in the embryo beginning to show a difference in size at E7.5 and lethality at E11.5 [19]. In this study, the effect of SEPHS1 deficiency on early embryogenesis was analyzed in more detail during this period. Figure 1A shows whole-embryo images obtained by optical microscopy. There are no differences between wild-type (Sephs1 +/+ and Sephs1 +/− ) and the Sephs1 −/− embryo at E6.5. The Sephs1 −/− embryo began to show size differences from wild-type at E7.5, and the difference extended at E8.5, wherein head folds were observed to be less developed than in wild-type. In addition, the Sephs1 −/− embryo does not turn in the final fetal position at E9.5 and is smaller than wild-type. Unlike in the wildtype embryo, optic and otic vesicles were not observed, although the allantois still remained, and three primary brain vesicles (prosencephalon, mesencephalon, and rhombencephalon) were observed as being immature. At E9.5, no vitelline vessel was found in the yolk sac (arrow in Figure S1A). At E10.5, the size of the Sephs1 −/− embryo was dramatically reduced from that at E9.5 (data not shown) and the embryo appeared to be fully absorbed at E11.5 (arrowhead in panel E11.5 of Figure 1A). We used 32,24,34,20,17, and 10 embryos at E6.5, E7.5, E8.5, E9.5, E10.5, and E11.5, respectively (See Materials and Methods for more detailed information). All the embryos that had the same genotype and were prepared at the same embryonic day showed similar phenotypes both in size and in morphology. overlapping between embryonic day E6.5-E8.5. DEG cut-off: max(FPKM) > 1, |Log2(Fold Change)| > 1 and p < 0.01. al, allantois; amc, amniotic cavity; ba, first branchial arch; ec, ectoderm; en, endoderm; epi, epiblast; exe, extra embryonic region; hf, head fold; ht, heart; lb, limb bud; me, mesoderm; mes, mesencephalon; op, optic vesicle; ot, otic vesicle; rho, rhombencephalon; t, tail; tel, telencephalon.
In order to assess the effect of SEPHS1 deficiency on embryonic development and the related signaling pathways responsible for phenotypic changes, RNA-seq was performed using purified RNA from wild-type and Sephs1 −/− embryos at the E6.5, E7.5, and E8.5 stages.
As shown in Figure 1B, Pou5f1 (Oct4) and Brachyury (T) were expressed most abundantly at E6.5 and at E7.5, respectively, and Six3 was expressed only at E8.5 in wild-type embryo. These results indicate that the pooled RNAs isolated from embryos at the same embryonic day were highly homogeneous. In addition, principal component analysis (PCA) was performed to examine the relationship between the read sets. PCA revealed that the differences in developmental stages are much more distinct than genotypic differences between wild-type and Sephs1 knockout ( Figure 1C). Differentially expressed genes (DEGs) were obtained at E6.5, E7.5 and E8.5 with the log 2 (fold change) cut off of ±1 at p < 0.01 (Table S1). There are 21 genes commonly affected by SEPHS1 deficiency at all three stages ( Figure 1D). The functions of the down-regulated genes include those relating to organogenesis (Cntnap2, Gata4, Mmp15, Asrgl1, and Arl6ip5) and cell survival (Hectd3) [25][26][27]. Interestingly, the up-regulated genes showed region-specificity; Fgg, Afp, Trf, and Serpinf2 in the extraembryonic region and Spp2 in the placenta. Notably, Galectin 2 (Lgals2), the most up-regulated DEG, is reported to be an oxidative stress-responsive gene shown to be up-regulated under H 2 O 2 treatment [28]. These results suggest that SEPHS1-deficiency causes developmental retardation and the induction of oxidative stress.

Pathway Analysis of Differentially Expressed Genes
In order to identify the pathways that are regulated by SEPHS1, pathway enrichment analysis using Metascape was performed with DEGs ( Figure 2A), with the log 2 (fold change) cut off of ± 0.5 at p < 0.01. At E6.5, embryo and tissue morphogenesis were significantly affected by Sephs1 knockout. At E7.5 and E8.5, development of differentiated tissues such as 'vascular morphogenesis', 'epithelial cell differentiation', 'mesenchymal development', 'head' development', and 'heart development' were greatly affected. Some of the predicted results in Figure 2A were consistent with the results of Ingenuity Pathway Analysis (IPA) ( Table S2). For example, the genes included in 'LXR/RXR activation' predicted by IPA were also found in the 'Regulation of body fluid level' and 'Plasma lipoprotein assembly, remodeling, and clearance' predicted by the Metascape analysis. These genes participate in retinoic acid (RA) signaling. In addition, most of the genes in the 'Coagulation system' of IPA were included in the Metascape category of 'Hemostasis'.
In order to identify the transcription factors targeting the DEGs during gastrulation, transcription factors that are known to be activated during gastrulation were selected from Transcriptional Regulatory Relationships Unraveled by Sentence-based Text mining (TRRUST) database. Among the gastrulation-specific transcription factors, the transcription factors that contain target genes (of which more than 50% are DEGs) were further selected to identify the transcription factors governing the expression of target DEGs ( Figure 2B). Among the selected transcription factors, 12 transcription factors were DEG themselves (asterisks in Figure 2B). We defined these transcription factor genes as differentially expressed transcription factor genes (DTFGs). Interestingly, DTFGs were apt to be clustered together. The expression of DEGs by non-DTFGs, the transcription factors whose expression was not changed by Sephs1 knockout, may be regulated indirectly by changes in protein stability, phosphorylation, and interaction with other co-regulators.
To analyze the expression pattern of DEGs, hierarchical clustering was performed ( Figure 2C). Expression patterns of DTFGs (asterisks in Figure 2B) were indicated on the right of each cluster to which they belong. The DTFGs in the same cluster in Figure 2B were in the same or neighboring cluster in Figure 2C. For example, Eomes, Tal1, Hnf4α, Gata4, Stat3, and Pitx2 were included in C3 and C4. This suggests that DTFGs are expressed in the same pattern with their target genes.
To identify biological processes that were most highly affected by SEPHS1 deficiency, gene ontology (GO) analysis was performed with DEGs belonging to each cluster ( Figure 2D). As a result, genes in each cluster were found to regulate distinct pathways. For example, the genes in cluster 1 were predicted to regulate pathways participating in axis formation and the genes in cluster 2 in neuron development.  In order to identify the transcription factors targeting the DEGs during gastrulation, transcription factors that are known to be activated during gastrulation were selected Since dermal layers of the embryos at the early gastrula stage determine the cell lineage leading to tissue differentiation, the gene expression pattern and/or morphological feature of each dermal layer are important to organ development at the following stages. We examined the expression levels and locations of dermal layer markers (Sox2, Otx2, Foxa2, and Eomes) selected from the transcription factors shown in Figure 2B.
Sex determining region Y-box 2 (Sox2) encodes a transcription factor that is essential for maintaining pluripotency of undifferentiated embryonic stem cells, but is also known to be expressed specifically in the ectoderm of both embryonic and extraembryonic regions during the gastrulation stage [29]. As shown in Figure 2E, the levels of Sox2 expression were significantly decreased in the Sephs1 −/− embryo compared to wild-type (log 2 (Fold Change) at E6.5 = −0.42). Orthodenticle homeobox 2 (Otx2), which encodes a member of the bicoid subfamily of homeodomain-containing transcription factors, is an embryonic mesoderm-specific gene that plays a key role in nervous system development [30]. The area in which Otx2 was expressed was reduced to the bottom half of the Sephs1 −/− embryo, while the size of the embryo was similar with that of wild-type. The fact that the structure of mesodermal layer was changed by SEPHS1-deficiency suggests that the development of connective tissues such as blood, blood vessels, muscles, and heart will proceed abnormally in the Sephs1 −/− embryos, because these tissues are differentiated from the mesodermal lineage cells. Notably, the expression level of Otx2 was not changed by Sephs1 knockout (log 2 (Fold Change) at E6.5 = −0.22). Forkhead box protein A2 (Foxa2) encodes a protein belonging to a subfamily of the Forkhead box transcription factors and is an endoderm-specific gene [31]. There was no difference in expression levels of Foxa2 and location of expression between wild-type and the Sephs1 −/− embryo. We observed that only the rate of differentiation of the endoderm-originated organs, such as the digestive and respiratory systems, was retarded and these organs appeared to be at the E9.5 stage in the Sephs1 −/− embryo. It appears that SEPHS1-deficiency does not affect the morphology of endodermal lineage tissues during development. Eomesodermin (Eomes), also referred to as T-box brain protein 2 (Tbr2), is a member of the T-box family of transcription factors initially expressed in the extraembryonic ectoderm and is known to play an important role in anterior visceral endoderm formation and epithelial-mesenchymal transition (EMT) [32,33]. In the wild-type embryo, Eomes was expressed throughout the extraembryonic ectoderm. However, in the Sephs1 −/− embryo, the expression area was limited to the bottom of the extraembryonic ectoderm region, wherein the expression levels were slightly decreased (log 2 (Fold Change) at E6.5 = −0.15), suggesting that extraembryonic lineage tissues, such as the yolk sac and placenta, may differentiate abnormally after E6.5. In conclusion, expression levels and locations of the dermal layer marker genes provide evidence that knockout of Sephs1 causes morphological changes as well as the retardation of development during dermal layer formation.

Morphological Changes in the Embryonic Region by Sephs1 Knockout
We confirmed the developmental abnormalities predicted in Figure 2 by examining morphological changes in the Sephs1 −/− embryos using X-ray microscopy (XRM) technology ( Figure 3). The central nervous system of the Sephs1 −/− embryo manifested retarded development and abnormal shape compared to those of the wild-type embryo. At E7.75, head fold (hf) was found in the wild-type embryo (panel (a) of Figure 3A) but not in the Figure 3A). At E8.5 and E9.5, although three primary brain vesicles (prosencephalon (pro) and mesencephalon (ms) and rhombencephalon (rho)) were formed, the size of brain tissue in the Sephs1 −/− embryo was much smaller than that in wild-type, and the neural groove between prosencephalon (pro) and mesencephalon (ms) was not closed (arrowheads in Figure 3B). In the transverse section of the cervical region in E8.5, the closures of the neural grooves in wild-type embryos were progressing, whereas in the Sephs1 −/− embryo, it was open to the outside ( Figure 3C). The neural plate (blue) was closed around somites (yellow) in the wild-type embryo (Video S1), but remained open in the Sephs1 −/− embryo (Video S2). In addition, the development of somites was delayed. Segmentation of somites was clearly observed in the wild-type embryo at E8.5 and E9.5, while somite segmentation did not occur in the Sephs1 −/− embryo until E9.5 (see yellow-colored regions in Figure 3B). Segmented somites (yellow) were straight and parallel with the neural plate (blue), and clearly distinguished from the neural plate (blue) in wild type at E9.5 (Video S3). On the other hand, somites (yellow) in the Sephs1 −/− embryo were severely twisted and had ambiguous boundaries (Video S4).   Heart development was also inhibited by SEPHS1-deficiency. As shown in Figure 3B, the differentiation of heart (ht) in the Sephs1 −/− embryo at E8.5 proceeded slowly compared to wild-type. Heart structure of Sephs1 −/− embryo remained in the cardiac crescent stage, while that in wild-type embryo developed into the heart tube stage showing ventricles and heart septum at E8.5 ( Figure 3C). Notably, the size of the heart in both wild-type-and Sephs1 −/− embryos was similar at E9.5, but the cell density was significantly reduced and the morphology of the heart was irregular in Sephs1 −/− embryos (compare Sephs1 +/+ with Sephs1 −/− of E9.5 in Figure 3C).
The archenteron (ar), which differentiates into foregut (fg) was not observed in the Sepsh1 −/− embryo at E7.75 ( Figure 3A), but observed at E8.5 and E9.5 ( Figure 3C). Hindgut (hg) development seemed to proceed more slowly than foregut formation in the Sephs1 −/− embryo. Hindgut formation was observed only at E9.5 in the Sephs1 −/− embryo ( Figure 3D). These data suggest that gut formation was not significantly affected by Sephs1 knockout and this phenomenon is consistent with the result of FOXA2 expression patterns in Figure 2E.
The process that connects allantois (al) with the chorion (ch) was completed in wildtype at E8.5, but al and ch were not connected in the Sephs1 −/− embryo at E8.5 ( Figure 3E). At E9.5, umbilical (uv) and vitelline vessels (vv, insets (a) and (b) in Figure 3E) were formed in wild-type, but not in the Sephs1 −/− embryo (insets (c) and (d) in Figure 3E). The fact that the formation of vv and uv is inhibited by SEPHS1-deficiency suggests that SEPHS1 plays an essential role in the transport of substances between mother and embryo.
In addition, the total volume of the embryo was measured from the 3-dimensional structure obtained by reconstructing XRM images using a volume calculation software (Dragonfly). The volume ratio between the Sephs1 −/− embryo and the wild-type decreased by approximately 3.2 times and 3.4 times at E8.5 and E9.5, respectively ( Figure 3F). At E8.5, the number of somites (sm) in Sephs1 −/− and wild-type embryos was 3 and 8.2, respectively. The gap in the somite number between Sephs1 −/− embryo and wild-type embryos was increased at E9.5 (5.3 somites in the Sephs1 −/− embryo and 21 somites in the wild-type embryo). The heart volume of the wild-type embryo at E8.5 was 10.1 times greater than that of the Sephs1 −/− embryo. However, at E9.5, the heart volume in the Sephs1 −/− embryo became abnormally enlarged ( Figure 3C). Unlike the wild-type embryo, the Sephs1 −/− embryo did not exhibit axial rotation at E9.5 and the embryonic axis was oriented in the same direction as at E8.5 (compare the E9.5 wild-type with the E9.5 Sephs1 −/− embryos in Figure 3B). No further axial rotation was observed at E10.5 (data not shown) suggesting that axial rotation stopped at E8.5.

Pattern Analysis of DEGs Expressed in Extraembryonic Region
Since biological processes (for example, 'placenta development' in Figure 2D) occurring in the extraembryonic region were also predicted, we hypothesized that development of the extraembryonic region also would proceed abnormally by SEPHS1 deficiency. To test if there is any abnormality of development in the extraembryonic region in the Sephs1 −/− embryo, pathways involved in the development of extraembryonic region were identified and the pathway-related morphological changes were examined. We first selected DEGs among targets of extraembryonic region-specific DTFGs listed in Figure 2C. Four DTFGs expressed in the extraembryonic region and 70 of their target DEGs were identified. GO analysis was performed for the identified genes ( Figure 4A). As a result, four pathways (Retinoid metabolism and transport (RMT), Vitamin transport (VT), Embryonic morphogenesis (EM), and Epithelial cell differentiation (ECD)) were predicted with high probability (log 10 (p) < −4.9).
identified and the pathway-related morphological changes were examined. We first selected DEGs among targets of extraembryonic region-specific DTFGs listed in Figure 2C. Four DTFGs expressed in the extraembryonic region and 70 of their target DEGs were identified. GO analysis was performed for the identified genes ( Figure 4A). As a result, four pathways (Retinoid metabolism and transport (RMT), Vitamin transport (VT), Embryonic morphogenesis (EM), and Epithelial cell differentiation (ECD)) were predicted with high probability (log10(p) < −4.9).  Hnf4α and target genes are predicted to participate mainly in 'Retinoid metabolism and transport' and 'Vitamin transport' (RMT and VT in Figure 4A). Most of the genes in 'Vitamin transport' are included in 'Retinoid metabolism and transport' except Duox2, Duoxa2, Cubn, and Amn. Duox2 is target of retinoic acid signaling and Duoxa2 is a maturation factor of Duox2 [34]. Cubn, which is the cobalamin (vitamin B12) receptor gene, has an inhibitory function against retinoic acid signaling and Amn encodes the protein that facilitates uptake of vitamin B12. Therefore, genes in 'Vitamin transport' are related with retinoic acid signaling and can be categorized into the same group with 'Retinoid metabolism and transport'. It should be noted that GO analysis using total DEGs by IPA also predicted the retinoic acid signaling with the highest probability (Table S2). Retinoic acid is a morphogen participating in axis formation, embryo growth, and cell fate determination during early embryogenesis [35]. In addition, Hnf4α and its target genes are expressed in extraembryonic endoderm which develops into the yolk sac at a later stage suggesting that HNF4α plays an essential role in yolk sac development and its function via retinoic acid signaling.
Gata4 and Eomes, and their target genes did not show enrichment in any specific pathway. Instead, they commonly participate in the pathways regulated by both Hnf4α and Tfap2c. It should be noted that Gata4 is an activator of Hnf4α and is expressed more widely than Hnf4α [36].
Tfap2c and target genes are predicted to participate mainly in embryo morphogenesis and epithelial cell differentiation (EM and ECD in Figure 4A). During organogenesis, embryonic morphogenesis and epithelial cell differentiation should occur together. Therefore, embryo morphogenesis and epithelial cell differentiation can be categorized into organogenesis. In addition, both embryo morphogenesis and epithelial cell differentiation share BMP4 as a key protein suggesting BMP4 is used as a common morphogen for embryo morphogenesis and epithelial cell differentiation. Recently, it was reported that Tfap2c, which is expressed in extraembryonic ectoderm, plays an essential role in the development of trophoblast and placenta [37]. Chorioallantoic fusion is a process in trophoblast differentiation. As described in the previous section, we found chorioallantoic fusion did not occur in the Sephs1 −/− embryo ( Figure 3E) suggesting that Tfap2c and its target genes inhibit chorioallantoic fusion in the SEPHS1-deficient embryo.
Since Hnf4α is known to be expressed in extraembryonic endoderm, we examined the region where this gene was expressed during gastrulation. At E6.5, HNF4α was distributed more widely in the extraembryonic region of the Sephs1 −/− embryo than that of wildtype, and at E9.5, the cell density of HNF4α-expressing yolk sac endoderm was lower in the Sephs1 −/− embryo than that of the wild-type ( Figure 4B). The yolk sac endoderm morphology was concavo-convex in the wild-type but became flat in the Sephs1 −/− embryo. Another significant difference in the structure of the yolk sac between wild-type and Sephs1 −/− embryo was that the yolk sac endoderm and mesoderm were separated in the Sephs1 −/− embryo at E9.5. XRM image showed that vitelline vessels (vv), which are normally generated in the yolk sac at E9.5 stage, were not found in the Sephs1 −/− embryo (compare panels (a) and (b) of Figure 4C). Separation of the yolk sac endoderm and mesoderm was also observed in the Sephs1 −/− embryo (panel (d) of Figure 4C). At the same time point, expression of CD31, which is known as a hemato-endothelial progenitor marker [38], was not detected in the Sephs1 −/− embryo, while it was detected in the wildtype embryo suggesting that both vitelline vessel and blood progenitors were not formed due to SEPHS1 deficiency ( Figure 4D).
It should be noted that SEPHS1, which is known to be expressed in all tissues, was expressed only in yolk sac endoderm, not yolk sac mesoderm in the wild-type embryo ( Figure 4B) suggesting that SEPHS1 is expressed cell-type specifically in the same tissue and its deficiency affects mainly the function of yolk sac endoderm at late gastrulation.

Pathway Prediction through Protein-Protein Interaction and Gene Set Enrichment Analysis
In order to identify the signaling pathway and molecular mechanism that participate in gastrulation, additional analyses were performed. By analyzing protein-protein interaction (PPI) of DEGs at E6.5, we could obtain three modules with high probability (log 10 (p) < −9.5); Wnt signaling pathway, Glutathione metabolism, and the post-translational protein phosphorylation pathway ( Figure 5A).

Pathway Prediction through Protein-Protein Interaction and Gene Set Enrichment Analysis
In order to identify the signaling pathway and molecular mechanism that participate in gastrulation, additional analyses were performed. By analyzing protein-protein interaction (PPI) of DEGs at E6.5, we could obtain three modules with high probability (log10(p) < −9.5); Wnt signaling pathway, Glutathione metabolism, and the post-translational protein phosphorylation pathway ( Figure 5A).  Wnt signaling regulates mainly morphogenesis, such as growth in the early embryo, axis formation, and pattern formation [39]. Wnt signaling was inhibited in Sephs1 −/− embryos at E6.5, suggesting that pathways determining embryonic morphology were affected by Sephs1 knockout before the organogenesis stage. Glutathione is one of the molecules that plays an essential role in maintaining redox homeostasis. Imbalance in glutathione metabolism by Sephs1 knockout leads to the accumulation of ROS that causes oxidative stress in the cell [19]. Prediction of post-translational protein phosphorylation with high probability suggests that signaling pathways are actively regulated, since almost all signaling pathways are regulated through phosphorylation of proteins participating in each cognate pathway. These results of PPI analysis suggest that SEPHS1 deficiency causes oxidative stress by disrupting redox homeostasis, and that the oxidative stress will primarily affect Wnt signaling.
Besides GO analysis, another useful method for pathway prediction is Gene Set Enrichment Analysis (GSEA). Because GO analysis uses only a gene list to predict pathways or processes without giving weight according to fold change value of each gene, we cannot determine to what extent a specific gene contributes to those pathways or processes. The advantage of GO analysis is that one can identify all possible pathways or processes enriched by DEGs. On the other hand, with GSEA methods, we can determine how a specific gene contributes to a specific pathway or process, since DEGs are ranked according to their fold change value and the ranked DEGs are used to calculate enrichment score for a pathway or process of interest [40].
GSEA with total genes using the Hallmark Gene Sets from the Molecular Signatures Database revealed that 'Reactive Oxygen Species' was predicted with high significance where normalized enrichment score (NES) and false discovery rate (FDR) were 1.50 and 0.02, respectively ( Figure 5B). There are 28 genes within the leading-edge subset (LES) of 'Reactive Oxygen Species' (Table S3). The genes include Glrx1, Prdx2, Prdx6, Txnrd1, Txnrd2, and Nqo1, and these genes participate in redox homeostasis or ROS generation. We then examined the gene number within LES and NES of 'Reactive Oxygen Species' in more detail at each embryonic day. The gene numbers within LES were 25, 27, and 31 at E6.5, E7.5, and E8.5, respectively (Table S3). In addition, the NES of 'Reactive Oxygen Species' at E6.5, E7.5, and E8.5 were 0.56, 1.0, and 1.42, respectively ( Figure S2B). These results suggest that ROS levels in the embryo are gradually increased. DNA damage, such as formation of γH2AX and 8-oxoguanine, is commonly used as an oxidative stress marker.
In order to determine whether DNA damage due to oxidative stress occurred in the Sephs1 knockout, immunohistochemistry (IHC) using an antibody against 8-oxoguanine and γH2AX (phosphorylated H2AX) was performed. Neither formation of γH2AX nor of 8-oxoguanine were detected from E6.5 to E8.5 ( Figure S2C) in Sephs1 +/− or Sephs1 −/− embryos. However, both 8-oxoguanine and γH2AX were generated only in the Sephs1 −/− embryo at E9.5 ( Figure 5C). These results indicate that oxidative stress was too mild to cause DNA damage in Sephs1 −/− embryos by E8.5, but was strong enough to cause DNA damage at E9.5.
In addition to 'Reactive Oxygen Species', 'Apoptosis' was also predicted with high significance where NES and FDR were 1.6 and 0.00, respectively ( Figure 5D). The leadingedge subset of 'Apoptosis' consists of 63 genes (Table S4) including Sqstm1, Pdcd4, Smad7, Hspb1, Faslg, Bax, Madd, Bmf, and Hgf, which are known to participate in apoptotic signaling. The gene numbers within LES and the NESs were gradually increased as embryogenesis proceeded ( Figure S2B and Table S4). Unexpectedly, we could not detect the activated form of caspase-3 in SEPHS1-deficient embryos until E8 ( Figure S2D). However, the activated form of caspase-3 was detected in the Sephs1 −/− embryo at E9.5, suggesting that SEPHS1 deficiency induced cell death through apoptosis in the embryo ( Figure 5D). Notably, genes such as Gadd45b, cdkn1a, and Btg2, which are related to cell proliferation, were also included in the apoptosis pathway. Therefore, we assumed that cell proliferation was also inhibited by Sephs1 knockout. PCNA is used as the most reliable marker for evaluating cell proliferation. As shown in Figure 5E, the level of PCNA was significantly decreased in the Sephs1 −/− embryo. These results strongly suggest that SEPHS1-deficiency causes inhibition of cell proliferation followed by apoptosis.

Discussion
SEPHS1 regulates various cellular functions such as redox homeostasis, and thus is known to be essential for embryo survival and growth. However, the detailed mechanisms of which genes or pathways are controlled by Sephs1 during development and how they result in embryonic lethality have not been fully elucidated. Through transcriptome analysis, the pathways that are responsible for morphological defects in Sephs1 −/− embryo were predicted. To predict pathways affected by Sephs1 knockout, a combination of various bioinformatic tools (AmiGO and Metascape analysis, IPA analysis and GSEA) with various databases for protein-protein interaction (String and BioGrid), transcription factor (TRRUST and ChIP-Atlas) and cell type (Mouse Gastrulation Atlas and MGI) were applied.
Bioinformatic analyses suggested several interesting features of the effect of SEPHS1deficiency on the development at early embryonal stages. Throughout the gastrulation stage, coagulation was predicted as the most highly affected pathway (log 10 (p) = −11.9). Genes included in this GO term are implicated in cell migration and adhesion which are the most prominent features of cells involved in the dermal layer and in axis formation [41]. Among signaling pathways participating in gastrulation, retinoic acid signaling was predicted to be activated throughout gastrulation with high significance (log 10 (p) = −8.9). Retinoic acid is a morphogen, and its signaling pathway is known to regulate axis formation and cell fate determination [35]. The upregulation of Rbp4, which mediates the transport of retinoic acid into the cell suggests that the levels of intracellular retinoic acids were increased by Sephs1 knockout. In addition, we found that the targets of retinoic acid pathway were also increased in their expression levels, suggesting that retinoic acid signaling was hyperactivated. Hyperactivation of retinoic acid signaling, due to imbalance of retinoic acid concentration in the cells, may adversely affect embryonal axis formation and cell fate determination. We observed various abnormalities during gastrulation and early organogenesis of SEPHS1-deficient embryo, including the absence of axial rotation.
PPI analysis produced other interesting findings. Wnt, prolactin and insulin-like growth factor (IGF) signaling pathways were predicted to be affected by Sephs1 knockout at E6.5, E7.5, and E8.5, respectively ( Figure 5A and Figure S2A). Wnt signaling is known to be involved in regulation of growth in early embryo development, axis formation, and pattern formation [39,42]. Since the expression of Wnt3a, Wnt2b, Wnt8a, and Fzd3 were significantly decreased in Sephs1 −/− embryo at E6.5, Wnt signaling seemed to be inhibited. Among the direct targets of Wnt signaling listed in the database 'the Wnt' [43], the expression of 85% of the targets was down-regulated in the Sephs1 −/− embryos. At E6.5, Sephs1 −/− embryos showed changes in the expression level and location of dermal layer markers, such as Sox, Otx2, and Eomes, which may lead to the retardation of embryonic growth and axis formation at E7.5 (see Figures 1A and 2E). Prolactin is known to play an important role in the trophoblast growth and the development and differentiation of neural crest cells in neurulation stage, but its role in the gastrulation stage has not been clearly identified [44]. Since the expression of prl2c2, prl2c3, prl4a1, and prl5a1 was increased in the Sephs1 −/− embryos, we hypothesized that prolactin signaling was hyperactivated. Among the target genes of the prolactin pathway in KEGG DB (map04917; Prolactin signaling pathway), the expression of Socs3, Elf5, Prlr, and Slc2a2 were increased by more than 1.5-fold in the Sephs1 −/− embryos at E7.5. Interestingly, although not at E7.5, abnormal development of the central nervous system was observed at E8.5 ( Figure 3B,C). One of the important roles of prolactin is to control the development of the nervous system [44]. IGF signaling is known to be involved in promoting growth of embryos and organ development [45]. The expression of IGFBP3 and other proteins including proteases that binds to IGFBP3 were up-regulated in the Sephs1 −/− embryos ( Figure S2A). Binding of IGF to IGFBP3 and binding of IGFBP3 to proteases inhibit IGF signaling [46]. Therefore, SEPHS1-deficiency seems to inhibit IGF signaling during late gastrulation and early neurulation. The Sephs1 −/− embryos showed growth retardation, structural brain abnormalities and disrupted cardiac development at E8.5 and these phenotypes are consistent with those described in another study [45]. These results suggest that SEPHS1-deficiency affects embryo morphogenesis by regulating Wnt and retinoic acid signaling and then organogenesis by regulating retinoic acid, prolactin, and IGF signal pathways in combination.
It should be noted that most of the genes responsible for reception of retinoic acid signaling (apolipoprotein and retinoic acid binding protein genes) are expressed in the extraembryonic region specifically, but its target genes are expressed both in the embryonic and extraembryonic regions. This finding is consistent with the experimental results. For example, both vitelline vessel formation and chorioallantoic fusion did not occur in the Sephs1 −/− embryo (Figures 3 and 4). Abnormalities in the yolk sac and placental development will inhibit nutrient uptake and waste disposal from the embryo.
Taking the results from both the bioinformatic analyses and experimental evidence into consideration, we propose the role of SEPHS1 during early mouse embryogenesis as follows ( Figure 6). SEPHS1 regulates intracellular ROS levels through regulating the redox homeostasis system. The deficiency of SEPHS1 causes disruption of the redox homeostasis system and leads to a gradual increase in ROS levels. Low levels of ROS will cause mild oxidative stress, leading to abnormalities in signaling pathways. During early gastrulation, abnormalities in embryonic morphogenesis, such as dermal layer and axis formation, occur presumably through Wnt and retinoic acid signaling. Abnormalities in organogenesis occur presumably through prolactin and retinoic acid signaling, followed by insulin-like growth hormone and retinoic acid signaling during mid-and late-gastrulation, respectively. ROS accumulate sufficiently to cause cell death through DNA damage at E9.5, followed by embryonic death. Dead embryos then undergo resorption by E14.5 [19]. Since the signaling pathways were predicted using bioinformatic tools, future molecular studies will be needed to further validate our results.
IGFBP3 and binding of IGFBP3 to proteases inhibit IGF signaling [46]. Therefore, SEPHS1deficiency seems to inhibit IGF signaling during late gastrulation and early neurulation. The Sephs1 −/− embryos showed growth retardation, structural brain abnormalities and disrupted cardiac development at E8.5 and these phenotypes are consistent with those described in another study [45]. These results suggest that SEPHS1-deficiency affects embryo morphogenesis by regulating Wnt and retinoic acid signaling and then organogenesis by regulating retinoic acid, prolactin, and IGF signal pathways in combination.
It should be noted that most of the genes responsible for reception of retinoic acid signaling (apolipoprotein and retinoic acid binding protein genes) are expressed in the extraembryonic region specifically, but its target genes are expressed both in the embryonic and extraembryonic regions. This finding is consistent with the experimental results. For example, both vitelline vessel formation and chorioallantoic fusion did not occur in the Sephs1 −/− embryo (Figures 3 and 4). Abnormalities in the yolk sac and placental development will inhibit nutrient uptake and waste disposal from the embryo.
Taking the results from both the bioinformatic analyses and experimental evidence into consideration, we propose the role of SEPHS1 during early mouse embryogenesis as follows ( Figure 6). SEPHS1 regulates intracellular ROS levels through regulating the redox homeostasis system. The deficiency of SEPHS1 causes disruption of the redox homeostasis system and leads to a gradual increase in ROS levels. Low levels of ROS will cause mild oxidative stress, leading to abnormalities in signaling pathways. During early gastrulation, abnormalities in embryonic morphogenesis, such as dermal layer and axis formation, occur presumably through Wnt and retinoic acid signaling. Abnormalities in organogenesis occur presumably through prolactin and retinoic acid signaling, followed by insulinlike growth hormone and retinoic acid signaling during mid-and late-gastrulation, respectively. ROS accumulate sufficiently to cause cell death through DNA damage at E9.5, followed by embryonic death. Dead embryos then undergo resorption by E14.5 [19]. Since the signaling pathways were predicted using bioinformatic tools, future molecular studies will be needed to further validate our results.  In this study, we elucidated how SEPHS1-deficiency induces developmental abnormality and embryonic lethality. Our results provide evidence that SEPHS1-deficiency is one of the contributors of natural miscarriages, and that the levels of SEPHS1 can be used as a marker for diagnosis or prognosis of natural miscarriages.

Materials
Antibodies against SEPHS1, γH2AX, 8-oxo-guanine, POU5F1, CFL488-conjugated mouse IgG, and CFL488-conjugated rabbit IgG were purchased from Santa Cruz Biotech by StringTie, gene counts by HTseq-count and differential expression analysis by Ballgown and DESeq2 packages in R. Differential gene expression was evaluated using max (FPKM) > 1, p < 0.01 and |log 2 (Fold Change)| > 1. To reduce the sample-to-sample systematic bias that may affect the interpretation, the data were calibrated by Trimmed Mean of M-values (TMM) normalization and estimating the size factor using count data 'edgeR' in R.

Pathway Analyses and Transcription Factor Prediction
Ingenuity Pathway Analysis (IPA) was carried out using the canonical pathway module in IPA and Gene Set Enrichment Analysis (GSEA) using pre-ranked modules. Functional enrichment analysis was carried out using the specified gene lists by PANTHER or Metascape. Transcription factor analysis was performed using TRRUST and PaGenBase modules in Metascape. Protein-protein interaction analysis was carried out using BioGrid and String modules in Metascape.
Manual transcription factor analysis was performed using TRRUST database. All transcription factors and their target lists were recruited from the TRRUST database, and then filtered into the 'gastrulation' GO category. The ratio of DEGs among these target genes was calculated. Transcription factors were filtered with max (FPKM) > 10, the number of target gene (s) in 'gastrulation' GO category > 1 and DEGs in target gene (%) > 50%.
Hierarchical clustering was performed using GenePattern. Distance was calculated by Euclidean distance or Pearson correlation, and clustered by pairwise complete-linkage method. Heatmap was illustrated by Treeview.
4.8. X-ray Microscopy Imaging 4.8.1. Sample Preparation Embryos were prepared, stained, and imaged as described previously with slight modification [62]. Embryos were dissected within decidua to fully preserve embryo to extra-embryonic association without disturbing the orientation or morphology. Embryos were fixed in 4% paraformaldehyde overnight at 4 • C and were washed with 1X PBS three times for 1 h each. Fixed samples were immersed in 0.1 N (v/v) Lugol's solution at room temperature for varying time according to the size of the embryo, 3-5 days. After staining, samples were mounted individually within a microcentrifuge tube filled with 0.5% w/v agarose and imaged immediately.

Image Acquisition
The raw data for 3D imaging of the samples were acquired via Xradia 620 versa (Carl Zeiss, Oberkochen, Baden-Württemberg, Germany). Each data set was acquired with the Xray source at 60 kV and 142 µA with a 0.5 mm aluminum attenuation filter. The acquisition time/isotropic voxel is 2 h/1.3 voxel to 4 h/0.8 voxel. Each sample was rotated 360 • along the anterior-posterior axis, and a projection image at 2016 × 1344 pixels was generated every 0.3 • at an average of 4 images. Acquired projection images were reconstructed and analyzed via Dragonfly (Version: 2021.1; ORS, Montreal, Quebec, Canada) software.

Histological Analysis
Embryos were dissected and fixed with 4% paraformaldehyde overnight at 4 • C. Fixed embryos were washed in running tap water for at least 6 hr and then paraffin embedded. The paraffin embedding process was performed automatically by Tissue processor (Leica, Wetzlar, Germany) following the routine overnight protocol. Paraffin blocks were sliced to 5 µm thickness, and then the slices attached to slide glasses coated with 3-aminopropyl triethoxysilane. After overnight drying at 37 • C on a slide warmer, the sample slides were used for further analysis. Hematoxylin and eosin staining was processed automatically by Autostainer (Leica).

Statistical Analysis
Statistical analyses were performed using an unpaired Student's t-test or one-way ANOVA with GraphPad PRISM 7.0. Statistical analysis for DEGs were performed by exact test using edgeR. A value of p < 0.01 was considered significant.