The Transcriptional Landscape of the Early Embryonic Developmental Biology in Hulunbuir Short-Tailed Sheep

Ethics statement
 
All sheep experimental procedures and protocols were approved and authorized by the animal care and use committee of Inner Mongolia Agricultural University in this study (License NO. SYXK, Inner Mongolia, 2016-0017).
 
Samples
 
Hulunbuir short tailed sheep (Fig 1a) were selected from Ewenki Autonomous Banner, Hulunbuir City, Inner Mongolia Autonomous Region.In the spring of 2021, it will be transported to Siziwang Banner Inner Mongolia Mengyuan Meat Seed Co., Ltd. for test and the test period is one and a half years. Synchronous estrus treatment was given for Hulunbuir short tailed sheep: Synchronous estrus treatment was performed by progesterone-releasing  intravaginal  device (PRID).
 
Methods
 
Implantation of PRID (PRIDs contain progestin drugs such as 30 mg flugestone or 50mg medroxyprogesterone acetate) sponge suppository. At the time of the implantation of suppository, a layer of penicillin powder and glycerol was applied on the outside of the sterilized catheter. Then, sponge suppository was sent into the deep vagina of the ewe near the cervix with a catheter and the thin thread was led to outside the vulva to remove the suppository at the end of treatment. The sponge suppository was placed for 12 days to achieve the purpose of Synchronous estrus. PMSG (330 IU / animal) was injected intramuscularly before suppository removal. Next, the sponge suppository was removed by pulling a thin thread and meanwhile, intramuscular injection of PG was also given to achieve a better effect of synchronous estrus. In this experiment, the estrus rate of ewes reached more than 90% after synchronous estrus treatment. The ewes were mated 36 hours after synchronous estrus treatment.
       
The complete uteri were surgically removed from the ewes on days 16 and 25 after mating, stored in tissue storage solution and brought back to the laboratory. The excess fetal membrane tissue was removed under a stereomicroscope. E16 (Fig 1b) and E25 (Fig 1c) samples were taken out and photographed. Finally, the complete embryo samples were stored in cryovials containing 2ml tissue preservation solution and three embryo samples were taken out respectively. Lianchuan Biotechnology Co., Ltd. is responsible for sequencing.
 

       
RNA extraction and library construction

Total RNA was isolated and purified using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) following the manufacturer’s procedure. The RNA amount and purity of each sample was quantified using NanoDrop ND-1000 (NanoDrop, Wilmington, DE, USA). The RNA integrity was assessed by Bioanalyzer 2100 (Agilent, CA, USA) with RIN number >7.0 and confirmed by electrophoresis with denaturing agarose gel. Poly (A) RNA is purified from 1 μg total RNA using Dynabeads Oligo (dT) 25-61005 (Thermo Fisher, CA, USA) using two rounds of purification. Then the poly(A) RNA was fragmented into small pieces using Magnesium RNA Fragmentation Module (NEB, cat.e6150, USA) under 94°C 5-7min. Then the cleaved RNA fragments were reverse-transcribed to create the cDNA by SuperScript™ II Reverse Transcriptase (Invitrogen, cat. 1896649, USA), which were next used to synthesise U-labeled second-stranded DNAs with E. coli DNA polymerase I (NEB, cat.m0209, USA), RNase H (NEB, cat.m0297, USA) and dUTP Solution (Thermo Fisher, cat.R0133, USA. An A-base is then added to the blunt ends of each strand, preparing them for ligation to the indexed adapters. Each adapter contains a T-base overhang for ligating the adapter to the A-tailed fragmented DNA. Single-or dual-index adapters are ligated to the fragments and size selection was performed with AMPureXP beads. After the heat-labile UDG enzyme (NEB, cat. m0280, USA) treatment of the U-labeled second-stranded DNAs, the ligated products are amplified with PCR by the following conditions: initial denaturation at 95°C for 3 min; 8 cycles of denaturation at 98°C for 15 sec, annealing at 60°C for 15 sec and extension at 72°C for 30 sec; and then final extension at 72°C for 5 min. The average insert size for the final cDNA library was 300±50 bp. At last, we performed the 2x150 bp paired-end sequencing (PE150) on an Illumina Novaseq™ 6000 (LC-Bio Technology CO., Ltd., Hangzhou, China) following the vendor’s recommended protocol.
 
Bioinformatics analysis of RNA-seq
 
Fastp software (https://github.com/OpenGene/fastp) were used to remove the reads that contained adaptor contamination, low quality bases and undetermined bases with default parameter. Then sequence quality was also verified using fastp. We used HISAT2 (https://ccb.jhu.edu/software/hisat2) to map reads to the reference genome of Ovis aries ARS-UI_Ramb_v2.0. The mapped reads of each sample were assembled using StringTie (https://ccb.jhu.edu/software/stringtie) with default parameters. Then, all transcriptomes from all samples were merged to reconstruct a comprehensive transcriptome using gffcompare (https://github.com/gpertea/gffcompare/). After the final transcriptome was generated, StringTie and was used to estimate the expression levels of all transcripts. StringTie was used to perform expression level for mRNAs by calculating FPKM.
   
The differentially expressed mRNAs were selected with fold change > 2 or fold change < 0.5 and with parametric F-test comparing nested linear models (p value < 0.05) by R package edgeR (https://bioconductor.org/packages/release/bioc/html/edgeR.html).

RNA-sequencing data
 
From the whole transcriptome of conceptuses samples (n = 3 embryos of Hulunbuir Short-tailed Sheep and n = 3 embryos of Ujumqin Sheep), an average of 52.41 million reads/sample was generated. After quality control, an average of 50.7 million reads/sample was kept. Moreover, an average of 87.36% of the reads were mapped against the sheep reference genome (Ovis aries, ARS-UI_Ramb_v2.0). The sequence data were deposited in the Sequence Read Archive of the NCBI databases under Bioproject number GSE186602 (https://www.ncbi.nlm.nih.gov/geo/info/linking.html,published on October 22, 2023).
 
Differentially expressed genes between E16 and E25
 
In the comparison of DEGs between E16 and E25, there were 4348 genes with significant differences (P<0.05), 1921 up-regulated genes in E16 and 2427 up-regulated genes in E25 (Fig 2).
 

 
Functional analysis
 
Through GO enrichment analysis of 1921 up-regulated genes (P<0.05) in E16, it was found that these genes could be enriched to 1165 GO terms, including 809 biological processes; 224 molecular functions; 132 cellular components (Fig 3a). In the enriched GO terms, several crucial for developmental biology were highlighted, such as endoplasmic reticulum, apical plasma membrane, membrane, etc. Through GO annotation analysis of 2427 up-regulated genes in E25, it was found that these genes could be annotated to 500 GO terms, including 327 biological processes; 103 molecular functions; 70 cellular components (Fig 3b). In the enriched GO terms, several crucial for developmental biology were highlighted, such as embryonic skeletal system morphogenesis, sequence-specific DNA binding, etc.
 

 
KEGG enrichment analysis of 1921 up-regulated genes (P<0.05) in E16 showed that these genes were enriched into 61 signaling pathways (Fig 4a), among which, the pathways were found, such as E25 receptor interaction, PI3K-Akt signaling pathway, Hippo signaling pathway and MAPK signaling pathway. KEGG enrichment analysis of 2427 up-regulated genes inE25 showed that these genes were enriched into 33 signaling pathways (Fig 4b). Among these signaling pathways, these signaling pathways were found, such as complement and coagulation cascades, PPAR, myocardial contraction, amino acid synthesis and metabolism.
 

       
Due to ethical reasons, many questions about human embryos cannot be answered directly. No animal model can truly reproduce the process of human embryonic development and however, pregnant sheep have been widely used to study maternal-fetal interaction. Therefore, it is necessary to have a thorough understanding of the development events of sheep embryos (Barry and Anthony, 2008; Yu and Guan, 2013). In addition, the research on early embryonic development in sheep will help us better regulate the gestation and reproductive process of sheep, which is also crucial for economic development.
 
E16
 
KEGG pathway annotation showed that DEGs are significantly enriched in ECM- receptor interaction pathway, focal adhesion, protein digestion and absorption pathway and some cancer pathways, which seemed to indicate that many DEGs and signaling pathways are involved in the occurrence of some cancer (Bao et al., 2019). This was due to the detection of a large number of epithelial cells with mesenchymal characteristics during organogenesis, which is similar to the characteristics of epithelial/mesenchymal cells during tumorigenesis (Dong et al., 2018). Collagen (COL) gene-involved ECM-receptor interaction, focal adhesion, as well as protein digestion and absorption pathways play a key role in cell proliferation and metastasis in organogenesis (Zhong et al., 2022) which is consistent with our analysis results. In the experimental data of this study, FN1 was significantly up-regulated and involved in the regulation of ECM receptor interaction, focal adhesion, cancer and actin cytoskeleton pathways (Wang et al., 2013). Integrins are heterodimeric transmembrane adhesion receptors composed of α- and α-subunits, whose members promote ECM-to-cell interactions and participate in cell adhesion and signal tansconduction, such as tissue development, tissue differentiation and immune responses (Zhu et al., 2020). Lama protein family is chiefly distributed in endothelial cells and some epithelial basement membrane (Shan et al., 2015), which mainly mediates cell adhesion and migration and plays an important role in tissue and organ formation by interacting with other ECM components during embryonic development (Jaluria et al., 2008). In addition, it was found in this study that PDGFRB and VEGFA were also enriched in focal adhesion. Platelet-derived growth factors (PDGFs) belong to a family of the growth factor peptides, which bind to their receptors (PDGFRA and PDGFRB) and stimulate cell growth and migration (Zhang et al., 2012; Gong et al., 2012). At the same time, it was also found in this study that the up-regulated genes in E16 were enriched in Axon guidance pathway, including semaphorins family members SEMA4G, SEMA3G, SEMA4A, SEMA3E and SEMA6D, which can affect such processes as axon bundling, branching and synaptic formation, as well as plexin PLXNA2 and PLXNA4 (Moskowitz and Lo, 2003), in the semaphorin receptor family. Hippo signaling pathway is an evolutionarily conserved signaling one, which can control the size of organs from flies to humans (Zhao et al., 2011). However, in the enrichment results of this study, the canonical transcription regulator YAP1 (Piccolo et al., 2014), the factor DLG11 (Enomoto and Igaki, 2011) that inhibits the entry of  YAP/TAZ complex into the nucleus and the transcription factors TEAD1 and TEAD4 (Lin et al., 2017) were found in Hippo pathway. In addition, it was noticed that in genes enriched in Hippo pathway, there were also TGFB1, TGFB2, TGFBR1, TGFBR2, BMP5, BMPR2, WNT2 and WNT11. Studies have shown that Hippo signaling pathway, as an important pathway regulating organogenesis during embryonic development, is not only related to Wnt signaling pathway, but also interacts with the TGF-β signaling pathway (Zhao et al., 2010), which is consistent with the results of this study. BMP2, BMP5, TGFB1 and TGFB2 as important ligands in the TGF-β signaling pathway, as well as BMPR2, TGFBR1 and TGFBR2 as indispensable receptors in the pathway, affect the osteogenic differentiation and apoptosis of cells via the BMP signaling and the classical TGF-âsignaling pathways, respectively (Lochab and Extavour, 2017).
 
E25
 
During the E25 period, the following pathways were significantly enriched. The complement system is canonically considered to be the main effector in innate immunity. It is a germline-encoded system composed of more than 50 circulating and membrane- bound proteins, most circulating proteins are produced in the liver, the core component of the complement system is C3 (Morgan and Gasque, 1997; Merle et al., 2015). The local production of complement proteins appears to be sufficient to generate humoral immune responses and it is the main source of complement in immunological privileged sites (e.g. the brain and eyes) (Barnum, 1995). In addition, studies also showed that complement plays a role in angiogenesis and synapse formation, especially an integral role in mouse development (Haynes et al., 2013). At the E25 stage, the significant enrichment of complement and coagulation cascades indicates that the embryo at this stage is constructing the innate immune system and meanwhile the occurrence of blood vessels, synapses, brain, eyes and liver is thus confirmed.
       
The ATP-binding cassette (ABC) transporters constitute one of the largest known protein families (Liu, 2019). In the immune system, antigen processing-related transporters TAP1 and TAP2 belong to subclass B of ABC transporters (ABCB2 and ABCB3, respectively), which are crucial molecules for antigen processing and loading onto major histocompatibility complex I (MHC I) (Eggensperger and Tampé, 2015). In addition, outside the immune system, bile acid production in the liver depends on several different ABC transporters, including ABCB4 (Morita and Terada, 2014). The generation of immature T cells occurs in the thymus, during which a subset of T cell receptors (TCRs) have the ability to transmit appropriate signals. The intensity of TCR signals largely depends on the composition and structure of plasma membrane (Wang et al., 2016). Therefore, the changes in plasma membrane lipid composition directly affect TCR-mediated signal transduction, as well as the development of thymocytes. The gene ABCC8 encodes the SUR1 subunit of the pancreatic β-cell ATP-sensitive K+ (K ATP) channel and affects pancreatic β-cell formation (Kapoor et al., 2013). The successful enrichment of ABC transporters shows the occurrence process of pancreatic β-cells, T-cells and liver.
               
In addition, in this study, the PPAR signaling pathway was also received attention. PPARs have three subtypes (PPAR alpha, beta/delta and gamma), which exhibit different expression patterns in vertebrates. PPAR alpha plays a role in circulation or the clearance of cellular lipid by regulating the expression of genes involved in lipid metabolism in liver and skeletal muscle. PPAR beta/delta is involved in lipid oxidation and cell proliferation. PPAR gamma promotes adipocyte differentiation to enhance blood glucose uptake (Ramirez et al., 2018). Retinoid × receptor (RXR) may exert an effect on the development and function of brain and nervous system. In the development of the nervous system, PL of neural membrane is rich in long-chain polyunsaturated fatty acids (PUFA), such as docosahexaenoic acid (DHA) and the brain-FABP (BFABP; FABP7) is mainly detected in the prenatal brain with a relatively low level in the adult brain. Hepatic carnitine palmitoyltransferase 1 (CPT1A) is the primary isoform (Choudhary et al., 2022) which is found in the liver, spleen, kidneys, lungs, intestines, pancreas, brain and ovaries. Brain carnitine palmitoy ltransferase 1 (CPT1C) is mainly expressed in the brain. The successful enrichment of the two classes of genes in PPAR signaling pathway indicates that fat metabolism in liver and brain is relatively vigorous at the E25 stage.

This study provides an exploratory landscape of transcriptional changes related to morphogenesis and organogenesis during early embryonic development in Hulunbuir short tailed sheep. The up-regulated genes in E16, which reflected that the transcriptional regulation of cell interaction and organogenesis. Additionally, the up-regulated genes in E25, which reflected that the organ and body structures are more complex at this stage of development. Our findings provide an exploratory and informative transcriptional landscape of sheep organogenesis, which might contribute to further studies addressing specific developmental events in sheep and in other mammals.

This work was supported by the National Natural Science Foundation of China (No. 31860689); Inner Mongolia Autonomous Region Science and Technology Plan of China (No.2019GG241), (No.2020ZD0007); Inner Mongolia Autonomous Region Science and Technology Achievements Transformation Project (No.2020CG0078); Key Technology Project of Inner Mongolia Autonomous Region (No. 2021GG0062).