Comprehensive Transcriptome Analysis of mRNA Expression Patterns of Early Embryo Development in Goat under Hypoxic and Normoxic Conditions

Simple Summary Oxygen plays a vital role in the development of early embryos, no matter whether it is too high or low, it will adversely affect the early embryo development, but the mechanisms involved in these effects are still unclear. RNA-seq was performed to compare 8-cell-stage and blastocyst-stage goat embryos under hypoxic and normoxic conditions, the mRNA expression mechanisms of 8-cell- and blastocyst-stage embryos were systematically analyzed under hypoxic and normoxic conditions. Functional enrichment analysis indicated that these differentially expressed genes (DEGs) were mainly related to biological processes and function regulation. In conclusion, we can infer that oxidative stress regulates early embryo development by affecting the expression of zygotic genes and transcription factors, and those stress genes play a potential role in adaptation to normoxic environments in goat embryos. Abstract It has been reported that hypoxic environments were more suitable for the in vitro development of mammalian embryos, but the underlying mechanisms were still unclear. In the present study, RNA-seq was performed to compare 8-cell-stage and blastocyst-stage goat embryos under hypoxic and normoxic conditions; zygotes were checked at 72 and 168 h to 8-cell stage (L8C) and blastocyst stage (LM) in hypoxic conditions and 8-cell stage (H8C) and blastocyst stage (HM) in normoxic conditions. In the H8C and L8C groups, 399 DEGs were identified, including 348 up- and 51 down-regulated DEGs. In the HM and LM groups, 1710 DEGs were identified, including 1516 up- and 194 down-regulated DEGs. The expression levels of zygotic genes, transcription factors, and maternal genes, such as WEE2, GDF9, HSP70.1, BTG4, and UBE2S showed significant changes. Functional enrichment analysis indicated that these DEGs were mainly related to biological processes and function regulation. In addition, combined with the pathway–gene interaction network and protein–protein interaction network, twenty-two of the hub genes were identified and they are mainly involved in energy metabolism, immune stress response, cell cycle, receptor binding, and signal transduction pathways. The present study provides comprehensive insights into the effects of oxidative stress on early embryo development in goats.


Introduction
Goats are widely distributed and diverse domestic animals because of their high social value and short gestation period; they have become an important species in current biological research and applications [1]. Assisted reproductive technologies such as superovulation, in vitro fertilization, in vitro embryo culture, and embryo transfer have been widely promoted and applied in goats and sheep, and they play a significant role in improving livestock productivity [2]. In vitro embryo culture has always been the top priority of our research; it directly determines the development of the individual. In addition, the study of goat embryos in vitro has laid the foundation for human-assisted reproductive technology, which is of great significance for the promotion of biological research.
Oxidative stress refers to a state in which the oxidative and antioxidative effects in the body are unbalanced, leading to excessive reactive oxygen species (ROS) accumulation; oxidative stress participates in a variety of key biological processes [3][4][5]. Recent reports suggested that the in vitro culture of embryos significantly increases the accumulation of ROS in embryonic cells, and the harmful effects of oxidative stress on embryonic development are caused by the excessive increase of ROS [6]. Compared with embryos developed in vivo, embryonic cells developed in vitro were more prone to apoptosis [7]. It has been reported that apoptosis of embryos cultured in vitro occurs at a specific stage, which may be caused by the influence of certain genes on a certain pathway, for example, apoptosis occurs in human 2-8-cell embryos [8], mice 1-cell embryos [9], and cattle 8-16-cell embryos [10]. In early embryos of mice, it was found that high oxygen concentration in the culture environment would accumulate a large amount of ROS and easily increase the incidence of apoptosis, while the hypoxic environment was more suitable for the in vitro development of mammalian embryos [11].
The culture environment of cells is the crucial factor that induces cell apoptosis and leads to embryonic development stagnation. Oxygen plays an important role in in vitro embryo culture, and the oxygen concentration directly affects the embryo development rate and quality [12]. It has been shown that early embryo culture can develop normally under atmospheric oxygen concentrations [13]. However, studies have found that the oxygen concentration in body tissues was generally 2-8% of the volume fraction in the hypoxic environment, and the fallopian tubes of most mammals were a hypoxic environment [14,15]. Studies have shown that reducing the oxygen concentration in the incubator was more conducive to embryo development in humans [16], cattle [17], mice [18], and swine [19]. In subsequent studies, 5% oxygen concentrations (hypoxic) were also used in vitro, resulting in a higher blastocyst rate and greater benefits before and after implantation compared with those of 20% oxygen concentrations (normoxia) [20].
In this study, we aimed to identify critical functional genes and pathways involved in goat embryos the oxidative stress adaptation. The mRNA expression mechanisms of 8-celland blastocyst-stage embryos were systematically analyzed under hypoxic and normoxic conditions. The analysis of these data will lay a preliminary foundation for further research on the biological mechanism of oxidative stress on embryonic development.

Materials and Methods
The goats were supplied by the Jinsheng Goat Breeding Technology Development Co, Ltd. Haimen City, Jiangsu Province, China. The trial was conducted at the Goat R&D Center in Haimen City, Jiangsu Province in December 2019. All biological reagents were sourced from Sigma-Aldrich (St. Louis, MO, USA), unless noted otherwise. This experiment was approved by the Animal Care (SYXK2011-0036) and Use Committee of Nanjing Agricultural University, China.

Collection and Culture of Embryos
Collection and culture of embryos was completed with reference to previous research methods [21]. Briefly, eight-month-old female goats were estrus synchronized, superovulated, and then mated with 2-year-old bucks. Progesterone sponges were intravaginally implanted in goats for 11 days, followed by the administration of 100 IU PG at the time of sponge removal, then donors were mated at 36 h and 48 h after sponge removal. The zygotes of goats within 24 h of mating were removed from the fallopian tubes by surgical method, washed three times with the embryo culture medium in vitro, divided into two groups, and transferred into the pre-balanced BO-IVC medium (IVC Bioscience, Falmouth, United Kingdom) supplemented with 10% fetal bovine se-rum (FBS). They were cultured at the 38.5 • C at a saturated humidity and at different gas composition conditions (5% CO 2 , 5% O 2 , 90% N 2 and 5% CO 2 , 20% O 2 , 75% N 2 ,). The one-cell embryos were cultured and divided into the hypoxia group and the normoxia group, which were checked at 72 and 168 h to 8-cells stage (L8C) and blastocyst stage (LM) under hypoxia and to 8-cell stage (H8C) and blastocyst stage (HM) under normoxia. Five 8-cell-stage or blastocyst-stage embryos were randomly absorbed by the mouth pipette as a sample, each group had 3 replicates, which were directly lysed for RNA sequencing.

RNA Library Construction and Sequencing
The Smart-seq2 method was used for cDNA synthesis as in a previous study [22]. The initiation and terminal connectors of the cDNA sequence were employed as primers for PCR amplification, and the obtained double-stranded cDNA was sufficient for library construction and sequencing. Incubation was performed with the fragmented cDNA using dsDNA fragmentase (NEB, M0348S) at 37 • C for 30 min. The fragmented cDNA at the size of 150-300 bp was screened by the provided sample purification beads for library construction. Illumina Novaseq™ 6000 was used for sequencing after the library was qualified; the sequencing read length was 2 × 150 bp.

Read Alignment and Differential Expression Analysis
The sequencing data were aligned to goat genome ARS1 using HISAT2 (version 2.0) with default parameters. SAMtools (version 0.19) was used to sort and index the reads, and only reads with a unique mapping location in the genome were retained for further analysis. Using Htseq with product count data for each transcript, the expression of genes was normalized to fragments per kilobase of exon model per million mapped reads (FPKM). EdgeR was used for differential analysis of StringTie (version 1.3.0) assembled and quantified genes, the threshold of significant difference was|log2 foldchange|≥ 1, p < 0.05). Graphic displays of the results of the differential expression were created using R software (version 3.2.5).

Functional Annotation of Differentially Expressed Genes
GO function and KEGG pathway enrichment analyses were performed with reference to previous studies [23]. GO and KEGG terms with FDR adjusted p < 0.05 were considered statistically significant.

Whole Transcriptome Amplification of Single-Cell Embryo and Quantitative Real-Time PCR
Single cell whole transcriptome amplification was performed according to the instructions of the CellAmp Whole Transcriptome Amplification Kit (Takara, Dalian, China). In brief, after putting 8-cell-or blastocyst-stage embryos directly into the lysate buffer to completely lyse, we add edMgCl 2 and RT Enzyme Mix to synthesize the first strand of cDNA. Then, dATP and TdT enzymes were added to the above reaction solution for an A tail reaction. Second-strand cDNA was synthesized by combining dNTP mix and Ex TaqHot Start, and the cDNA was amplified using Ex TaqHot, which was diluted for subsequent experiments.
Quantitative Real-Time PCR (qRT-PCR) was implemented using a previous researched method [24]. The relative mRNA expression levels of target genes were calculated using the comparative 2 −∆∆Ct method after normalization to GAPDH control. The PCR primers were designed using NCBI Primer Blast online software, synthesized by Tsingke (Nanjing, China), and are shown in Table S1.

Cell Staining
The Hoechst staining method was performed according to previous research [25]. Embryos were washed thrice with phosphate-buffered saline (PBS) containing 0.3% polyvinyl pyrrolidone and then fixed in 4% paraformaldehyde for 1 h, washed thrice again and transferred to 0.5% TritonX-100 for 20 min. After blockade with 1% bovine serum albumin for 1 h at room temperature they were treated with Hoechst 33,342 for 10 min. Finally, embryos were mounted on glass slides and observed under a confocal laser scanning microscope (LSM 700 META; Zeiss, Jena, Germany).

Statistical and Data Analysis
The data for cell counts were normalized by arcsine transformation and tested by Wilcoxon test. The data obtained for cleavage rate, blastocyst rate, and qRT-PCR were statistically analyzed by using the Student's t-test. All the statistical analysis was performed using SPSS software (version 24.0, SPSS) and * indicates a significant difference (p < 0.05), while ** means the difference was extremely significant (p < 0.01). The data were expressed as mean ± SE.

Effects of Oxygen Concentration on Embryonic Growth and Development
As shown in Figure 1, the zygotes were checked at 72 and 168 h at the 8-cell stage and blastocyst stage, respectively. The cell counts were calculated by staining the blastocysts under normoxic and hypoxic conditions with Hoechst 33,342. The results showed that the cleavage rates (87.30 ± 4.39 vs. 86.67 ± 5.77) were not significantly different under hypoxic conditions compared to normoxic conditions. However, hypoxia significantly (p < 0.05) improved the blastocyst rate (51.44 ± 3.67 vs. 26.85 ± 1.60) and increased the cell count (106 ± 10.44 vs. 74.67 ± 2.52) as compared to normoxia.

Summary of Transcriptomic Profiles of Goat Embryos
Our RNA sequencing data generated a total of 75.37 GB of clean data, in which 84.26-89.59% were mapped to the goat reference genome, and about 60% of the reads were uniquely mapped to the genomes in the mapped reads (Table S2). In the reference genome, alignment region distributions of 88.87-93.58% were located in annotated exons and 5.13-7.86% in introns ( Figure S1). The known annotation information in the goat database was counted ( Figure S2  The images (f,g) were captured with a microscope (LSM 700 META; Zeiss, Jena, Germany). The cleavage rate (d) was not significantly different under hypoxic conditions compared to that of normoxic conditions, but the blastocyst rate (h) and cell count (i) were significantly increased compared to those in the normoxic conditions. Embryonic development data were collected under hypoxic (n = 71) and normoxic (n = 60) conditions. * p < 0.05, ** p < 0.01.

Summary of Transcriptomic Profiles of Goat Embryos
Our RNA sequencing data generated a total of 75.37 GB of clean data, in which 84.26-89.59% were mapped to the goat reference genome, and about 60% of the reads were uniquely mapped to the genomes in the mapped reads (Table S2). In the reference genome, alignment region distributions of 88.87-93.58% were located in annotated exons and 5.13-7.86% in introns ( Figure S1). The known annotation information in the goat database was counted ( Figure S2  The images (f,g) were captured with a microscope (LSM 700 META; Zeiss, Jena, Germany). The cleavage rate (j) was not significantly different under hypoxic conditions compared to that of normoxic conditions, but the blastocyst rate (h) and cell count (i) were significantly increased compared to those in the normoxic conditions. Embryonic development data were collected under hypoxic (n = 71) and normoxic (n = 60) conditions. * p < 0.05, ** p < 0.01.

Identification and Functional Analysis of DEGs in 8-Cell-Stage Goat Embryo
Comparing H8C with L8C, a total of 399 DEGs were identified, including 348 upregulated and 51 down-regulated DEGs (Figure 3a, Table S3). The average FPKM of the DEGs in each sample in H8C was significantly increased compared to that of the samples of L8C (Figure 3b). To further study the role of these DEGs in 8-cell-and blastocyst-stage embryos under oxidative stress, GO and KEGG analyses were implemented to evaluate the function of the DEGs. The DEGs were enriched in GO terms including immune response, neutrophil chemotaxis, nuclear nucleosome, and chemokine activity (Figure 3c, Table S4). KEGG pathway enrichment analysis mainly contained the phagosome, toll-like receptor signaling pathway, cytokine-cytokine receptor interaction, and chemokine signaling pathway, (Figure 3d, Table S5). Similar results were found in the GSEA (Figure 3e).

Identification and Functional Analysis of DEGs in 8-Cell-Stage Goat Embryo
Comparing H8C with L8C, a total of 399 DEGs were identified, including 348 up-regulated and 51 down-regulated DEGs (Figure 3a, Table S3). The average FPKM of the DEGs in each sample in H8C was significantly increased compared to that of the samples of L8C (Figure 3b). To further study the role of these DEGs in 8-cell-and blastocyst-stage embryos under oxidative stress, GO and KEGG analyses were implemented to evaluate the function of the DEGs. The DEGs were enriched in GO terms including immune response, neutrophil chemotaxis, nuclear nucleosome, and chemokine activity (Figure 3c, Table S4). KEGG pathway enrichment analysis mainly contained the phagosome, toll-like receptor signaling pathway, cytokine-cytokine receptor interaction, and chemokine signaling pathway, (Figure 3d, Table S5). Similar results were found in the GSEA (Figure 3e).

Identification and Functional Analysis of DEGs in Blastocyst-Stage Goat Embryo
Comparing HM with LM, a total of 1710 DEGs were identified, including 1515 upregulated and 194 down-regulated DEGs (Figure 4a, Table S6). The average FPKM of the DEGs in each sample in HM was significantly increased compared to that of the samples of LM (Figure 4b). To further study the role of these DEGs in blastocyst-stage embryos under oxidative stress, GO and KEGG analyses were implemented to evaluate the function of the DEGs. GO enrichment analysis indicated that the DEGs were enriched in the regulation of cell proliferation; extracellular exosome, cytoplasm, and phospholipase inhibitor activity; and other correlated biological processes (Figure 4c, Table S7). The KEGG pathway mainly included the Hippo signaling pathway, thyroid hormone synthesis, PI3K-Akt signaling pathway, and apoptosis ( Figure 4d, Table S8). Similar results were found in the GSEA (Figure 4e).

Identification and Functional Analysis of DEGs in Blastocyst-Stage Goat Embryo
Comparing HM with LM, a total of 1710 DEGs were identified, including 1515 upregulated and 194 down-regulated DEGs (Figure 4a, Table S6). The average FPKM of the DEGs in each sample in HM was significantly increased compared to that of the samples of LM (Figure 4b). To further study the role of these DEGs in blastocyst-stage embryos under oxidative stress, GO and KEGG analyses were implemented to evaluate the function of the DEGs. GO enrichment analysis indicated that the DEGs were enriched in the regulation of cell proliferation; extracellular exosome, cytoplasm, and phospholipase inhibitor activity; and other correlated biological processes (Figure 4c, Table S7). The KEGG pathway mainly included the Hippo signaling pathway, thyroid hormone synthesis, PI3K-Akt signaling pathway, and apoptosis ( Figure 4d, Table S8). Similar results were found in the GSEA (Figure 4e).

Identification and Functional Analysis of Common DEGs between H8C and L8C Groups and HM and LM Groups
A total of 66 DEGs were identified in both groups of embryos (Figure 5a, Table S9). The average FPKM in each group of embryos under normoxia was significantly increased when compared to those under hypoxia (Figure 5b). The cluster analysis of 66 DEGs is shown in the Figure 5c. KEGG enrichment analysis was implemented and revealed that these DEGs were principally enriched in pathways related to cellular processes; environmental information processing; genetic information processes; human diseases; metabolism and organismal systems, such as cellular community, signal molecules, and interaction; folding, sorting, and degradation; immune disease; nucleotide metabolism; and immune system (Figure 5d, Table S10).

Identification and Functional Analysis of Common DEGs between H8C and L8C Groups and HM and LM Groups
A total of 66 DEGs were identified in both groups of embryos (Figure 5a, Table S9). The average FPKM in each group of embryos under normoxia was significantly increased when compared to those under hypoxia (Figure 5b). The cluster analysis of 66 DEGs is shown in the Figure 5c. KEGG enrichment analysis was implemented and revealed that these DEGs were principally enriched in pathways related to cellular processes; environ-mental information processing; genetic information processes; human diseases; metabolism and organismal systems, such as cellular community, signal molecules, and interaction; folding, sorting, and degradation; immune disease; nucleotide metabolism; and immune system (Figure 5d, Table S10).

Analysis of the Protein-Protein Interaction (PPI) and Pathway-Gene Interaction Networks of DEGs
In HM and LM groups, the string online database and Cytoscape (version 3.7.2) tool were utilized to establish the PPI network, which was made up of 184 nodes (FPKM ≥ 2 and fold change ≥ 2) and 265 edges ( Figure S3a, Table S11). The MOCODE tool was used to screen the most important modules in the PPI network, including 6 nodes and 14 edges ( Figure S3b), such as PLK2, WEE1, CKS1B, CCNA1, CDC25B, and CDKN1A. KEGG enrichment analysis of the DEGs in these modules associated them with the FOXO signaling pathway, cell cycle, P53 signaling pathway, and PI3K-Akt signaling pathway (Table S12).
The pathway-gene interaction network diagram was established to further analyze the effects of different oxygen concentrations on embryonic development. In HM and LM groups, ten KEGG pathways associated with oxidative stress were screened, and a total of 113 DEGs (FPKM > 1) participated in these pathways, the main genes involved in these pathways were FAS, GAPD1, GADD45A, ZBTB33, CREB3, MDM2, CDKN1A, CDC25B,

Analysis of the Protein-Protein Interaction (PPI) and Pathway-Gene Interaction Networks of DEGs
In HM and LM groups, the string online database and Cytoscape (version 3.7.2) tool were utilized to establish the PPI network, which was made up of 184 nodes (FPKM ≥ 2 and fold change ≥ 2) and 265 edges ( Figure S3a, Table S11). The MOCODE tool was used to screen the most important modules in the PPI network, including 6 nodes and 14 edges ( Figure S3b), such as PLK2, WEE1, CKS1B, CCNA1, CDC25B, and CDKN1A. KEGG enrichment analysis of the DEGs in these modules associated them with the FOXO signaling pathway, cell cycle, P53 signaling pathway, and PI3K-Akt signaling pathway (Table S12).
The pathway-gene interaction network diagram was established to further analyze the effects of different oxygen concentrations on embryonic development. In HM and LM groups, ten KEGG pathways associated with oxidative stress were screened, and a total of 113 DEGs (FPKM > 1) participated in these pathways, the main genes involved in these pathways were FAS, GAPD1, GADD45A, ZBTB33, CREB3, MDM2, CDKN1A, CDC25B, and IL6. (Figure 6). In addition, seven DEGs (HSP70.1, HSPB1, FOXO1, FABP3, PPP1R1A, GPX6, MAPK13) were related to energy metabolism, stress responses, and metabolism. The enriched KEGG mainly contained the AMPK signaling pathway, FOXO signaling pathway, and Rap1 signaling pathway. In summary, the analysis results showed that sixteen DEGs were identified to play a potential role in adaptation to normoxic environments in goat embryos (Table S12). and IL6. (Figure 6). In addition, seven DEGs (HSP70.1, HSPB1, FOXO1, FABP3, PPP1R1A, GPX6, MAPK13) were related to energy metabolism, stress responses, and metabolism. The enriched KEGG mainly contained the AMPK signaling pathway, FOXO signaling pathway, and Rap1 signaling pathway. In summary, the analysis results showed that sixteen DEGs were identified to play a potential role in adaptation to normoxic environments in goat embryos (Table S12). Figure 6. The pathway-gene interaction network diagram was established using ten KEGG pathways and 113 DEGs; the red is the pathway and the green is the gene name.

Key Genes Associated with Oxidative Stress Affecting Embryonic Development and Validation of Sequencing Data
The normal development of an embryo depends on the activation of zygotic genome activation (ZGA), and transcription factors are involved in regulating the activation of ZGA. Eighteen DEGs related to embryonic development were screened and used for heat map analysis in HM and LM ( Figure 7a). As shown in Figure 7b, the expression levels of zygotic genes, BTG4 and HSP70.1, transcription factors TEAD3, ZNF596, ZNF879, ZNF385B, ZAR1L, NPM2, ZFP36L1, SOX6, LITAF, NANOG, and ZBTB33, and maternal genes WEE1, WEE2, and GDF9 were significantly increased, while the expression levels of zygotic genes UBE2S and EIF4EBP1 were significantly decreased. Oxidative stress may affect early embryo development by regulating the expression of these genes.
To further assess the reliability of RNA sequencing, CKS1B, ALDH1A1, UBB, SUMO2, and OXT were randomly selected by qRT-PCR to validate the relative expression level of sequencing in the embryo. The relative expression level of qRT-PCR had the same trend as that of RNA-seq, which showed that sequencing data were efficient and reliable ( Figure S4). Figure 6. The pathway-gene interaction network diagram was established using ten KEGG pathways and 113 DEGs; the red is the pathway and the green is the gene name.

Key Genes Associated with Oxidative Stress Affecting Embryonic Development and Validation of Sequencing Data
The normal development of an embryo depends on the activation of zygotic genome activation (ZGA), and transcription factors are involved in regulating the activation of ZGA. Eighteen DEGs related to embryonic development were screened and used for heat map analysis in HM and LM ( Figure 7a). As shown in Figure 7b, the expression levels of zygotic genes, BTG4 and HSP70.1, transcription factors TEAD3, ZNF596, ZNF879, ZNF385B, ZAR1L, NPM2, ZFP36L1, SOX6, LITAF, NANOG, and ZBTB33, and maternal genes WEE1, WEE2, and GDF9 were significantly increased, while the expression levels of zygotic genes UBE2S and EIF4EBP1 were significantly decreased. Oxidative stress may affect early embryo development by regulating the expression of these genes.
To further assess the reliability of RNA sequencing, CKS1B, ALDH1A1, UBB, SUMO2, and OXT were randomly selected by qRT-PCR to validate the relative expression level of sequencing in the embryo. The relative expression level of qRT-PCR had the same trend as that of RNA-seq, which showed that sequencing data were efficient and reliable ( Figure S4).

Discussion
This experiment showed that different oxygen concentrations had significant effects on embryo development and identified the key genes and pathways involved in oxidative stress, which further revealed the regulatory mechanism of oxidative stress. It lays a preliminary foundation for us to further study the influence of culture environment on early embryo development.
ZGA is the genome-wide gene activation in the zygote, which plays an extremely critical role in early embryonic development [26]. In this experiment, 8-cell-stage and blastocyst-stage embryos were selected for sequencing analysis to explore the effects of different oxygen concentrations on the development of early goat embryos. Some studies have shown that the ZGA period of goats and sheep both occurs in the 8-cell stage [21,27]. During the ZGA period, maternal mRNA and protein were gradually eliminated, and zygotic genes were activated to participate in the regulation of early embryonic development [28].
We found that the number of DEGs in the HM and LM groups was far higher than that in the H8C and L8C groups. A large number of transcription factors were activated after the early embryonic ZGA occur, this result may be related to the ZGA period of goat embryos starting at the 8-cell stage. In addition, we analyzed the difference in the rate of blastocysts under different oxygen concentrations, and the results showed that the expression levels of maternal genes, zygotic genes, and transcription factors were significantly different in the two groups. These genes play an essential role in early embryonic development. Studies have shown that the maternal gene WEE1 was essential for maintaining genomic integrity in mammals; lack of WEE1 can lead to DNA damage and chromosome aneuploidy, which plays an important role in the regulation of early embryonic development [29]. Previous studies found that the transcription level of GDF9 was also different in various stages of embryonic development [30,31]. In mouse oocytes, BTG4 plays the role of the meiotic cell cycle coupling maternal-to-zygotic transition (MZT) permitting factors, while females lacking BTG4 will produce other morphologically normal oocytes, but they are unable to reproduce due to early developmental arrest [32]. As a new effector of SOX2 protein degradation, UBE2S regulates the stability of SOX2 and participates in the maintenance of embryonic stem cell pluripotency [33]. The synthesis of NPM2 maternal

Discussion
This experiment showed that different oxygen concentrations had significant effects on embryo development and identified the key genes and pathways involved in oxidative stress, which further revealed the regulatory mechanism of oxidative stress. It lays a preliminary foundation for us to further study the influence of culture environment on early embryo development.
ZGA is the genome-wide gene activation in the zygote, which plays an extremely critical role in early embryonic development [26]. In this experiment, 8-cell-stage and blastocyststage embryos were selected for sequencing analysis to explore the effects of different oxygen concentrations on the development of early goat embryos. Some studies have shown that the ZGA period of goats and sheep both occurs in the 8-cell stage [21,27]. During the ZGA period, maternal mRNA and protein were gradually eliminated, and zygotic genes were activated to participate in the regulation of early embryonic development [28].
We found that the number of DEGs in the HM and LM groups was far higher than that in the H8C and L8C groups. A large number of transcription factors were activated after the early embryonic ZGA occur, this result may be related to the ZGA period of goat embryos starting at the 8-cell stage. In addition, we analyzed the difference in the rate of blastocysts under different oxygen concentrations, and the results showed that the expression levels of maternal genes, zygotic genes, and transcription factors were significantly different in the two groups. These genes play an essential role in early embryonic development. Studies have shown that the maternal gene WEE1 was essential for maintaining genomic integrity in mammals; lack of WEE1 can lead to DNA damage and chromosome aneuploidy, which plays an important role in the regulation of early embryonic development [29]. Previous studies found that the transcription level of GDF9 was also different in various stages of embryonic development [30,31]. In mouse oocytes, BTG4 plays the role of the meiotic cell cycle coupling maternal-to-zygotic transition (MZT) permitting factors, while females lacking BTG4 will produce other morphologically normal oocytes, but they are unable to reproduce due to early developmental arrest [32]. As a new effector of SOX2 protein degradation, UBE2S regulates the stability of SOX2 and participates in the maintenance of embryonic stem cell pluripotency [33]. The synthesis of NPM2 maternal mRNA was a necessary condition for the activation of the zygotic genome and is essential for normal development after the blastocyst stage. NPM2 knockout mouse embryos rarely develop to the 2-cell stage due to the lack of complexes necessary for transcription [34]. NANOG as a transcription factor that is necessary for the self-renewal of embryonic stem cells and maintenance of pluripotency. The number of vegetative ectoderms and the total number of blastocysts in embryos is reduced with NANOG expression [35]. According to the above, the difference in blastocyst rates between the two groups may be closely related to the differential expression of these genes. The differential expression of zygote genes, maternal genes, and transcription factors was involved in regulating the development of early embryos.
It can be seen from the experiment results that the blastocyst rate of early embryos cultured under hypoxia were greater than those cultured under normoxia, suggesting that the number of arrested embryos under normoxic conditions was more than that of hypoxic conditions after cleavage, so we can infer that differences in embryo development were caused by different oxygen concentrations, and hypoxic culture was more conducive to embryo development. Studies have shown that embryos cultured under normoxia can accelerate the formation of ROS in the environment and produce more ROS than embryos cultured under hypoxia [36,37]. Our results were consistent with previous studies that showed that excessive ROS can affect maternal-to-zygotic transition and embryonic development retardation [38]. It also has been proven that hypoxia was more conducive to improving early embryonic development [39,40]. In addition, the effect of OS on the in vitro culture of sheep embryos has also been reported; L-carnitine supplementation during in vitro maturation reduces oxidative stress-induced embryo toxicity by decreasing intracellular ROS that in turn alters transcript levels of antioxidant enzymes and improves developmental potential of oocytes and embryos [41]. By studying the molecular mechanism of OS in embryo development, we are committed to exploring the suitable culture environment for embryo development, which has a positive significance for promoting the development of in vitro-assisted reproductive technology.
Previous studies have found that oxidative stress significantly alters the gene expression of a variety of cells to adapt to the culture environment [42]. In HM and LM groups, a total of 1710 differential genes were screened and identified: FOXO1, GPX6, GADD45A, HSP70.1 were screened as hub genes and it has been confirmed that these genes were related to the regulation of oxidative stress [43][44][45][46]. Previous studies have shown that the silencing of FABP3 can cause an excessive production of intracellular ROS and lead to mitochondrial dysfunction [47]. Functional enrichment analysis demonstrated that these DEGs were mainly related to biological processes and function regulation. For example, CDC25B, CCNA1, WEE1, PLK2, and CKS1B are related to cell cycle [48][49][50][51][52], while HSPB1, ZBTB33, BMP4, and PDGFA are involved in regulating embryonic growth and development [53][54][55][56].
Through the analysis of the enrichment of 8-cell stage differential genes using GO and KEGG under different oxygen concentrations, we found that the effect of ROS on 8-cell-stage embryos was mainly enriched in the biological process of immune response, inflammatory response, and disease occurrence. Previous studies have shown that oxidative stress can cause inflammation and immune responses in cells [57]. When the embryo developed to the blastocyst stage through zygotic genome activation, the main effect of ROS on the blastocyst-stage embryo was to promote angiogenesis, energy metabolism, hormone synthesis, cell apoptosis, and other biological processes, similar results were obtained in porcine granulosa cells undergoing oxidative studies [58]. A total of 66 DEGs were identified in both groups of embryos; the expression of these DEGs was not related to the stage of the embryo, but was caused by the different culture environment, which shows that oxidative stress is involved in the regulation of early embryonic development. We can infer that both 8-cell-stage embryos and blastocyst-stage embryos adapt to the culture environment by regulating organic defense mechanisms in response to oxidative stress.
Studies have shown that multiple signaling pathways were induced when oxidative stress regulated cellular biological processes, such as Wnt [59], TGF-β [3], FOXO [60], Hippo [61], PI3K-Akt [62], and MAPK signaling pathways [63]. In this experiment, the AMPK signaling pathway, p53 signaling pathway, and Rap1 signaling pathway were also highly enriched in the process of oxidative stress in the early development of goat embryos. Based on the above results, we hypothesize that goat embryos may adapt to oxidative stress through energy metabolism, immune stress response, cell cycle changes, receptor binding, signal transduction pathways, and expression of stress genes.

Conclusions
In summary, differentially expressed RNA profiles were constructed by RNA-seq technology in this experiment and demonstrated that 8-cell-stage-and blastocyst-stage embryos in goats have significant changes in gene expression under different oxygen concentrations. Oxidative stress regulates early embryo development by affecting the expression of zygotic genes and transcription factors, and those stress genes play a potential role in adaptation to normoxic environments in goat embryos. These results provide new insights for further understanding of the regulation mechanism of oxidative response in early embryo development of goats.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/biology10050381/s1, Figure S1. The reference genome alignment region distribution; Figure S2. The known annotation information in the species database was counted; Figure S3. The PPI network established using 132 up-regulated DEGs; Figure S4. Validation of sequencing data. Table S1. The PCR primers were designed using NCBI Primer Blast online software.