Since, we haven’t interacted with the alive animals in this study and all the data were provided by the other researchers or were obtained from the public databases, therefore the approval by the animal ethics approval was not required.

In the present study, omics data in the levels of genome (SNP array genotypes and whole genome sequences data) and transcriptome (RNA-seq data) were used to identify potential lethal genes in Holstein dairy cattle. There are two approaches to study the early embryo lethal variants and genome regions, including variant- and haplotype-based detection. In this work, we applied population genomic methods to identify the lethal variants and genes following the procedure proposed by Todd et al. (2020). The flow chart of the consecutive steps of processing the data and applied approaches in this study is shown in Figure 1 .

The SNP genotypes consisted of ≥ 617885 markers originated from three Holstein cattle populations (n=3117) genotyped by Illumina BovineHD 770K SNP BeadChip platform. Animals with more than 5% missing genotypes were excluded (n=180). SNPs with unknown position or located on sex chromosomes were filtered out. Also, only SNPs with ≥ 95% genotype call rate and minor allele frequency (MAF) ≥ 0.05 were retained. The quality control step was performed using PLINK v 1.9. Following the quality control step, the genotype information (> 559K SNPs) for 2937 individuals were available for further analysis. The number of samples and SNPs in raw data and after quality control step are shown in Tables 1 , 2 , respectively. The distribution of the SNPs over the chromosomes is summarized in Supplementary Table S1 .

The whole genome sequence (WGS) data were obtained from Alemu et al. (2021), which is publicly free available. The data included 13,037,955 markers for 743 Holstein cattle. The whole process of generating the final Variant Call Format (VCF) file including DNA extraction, library preparation, reads alignment to the reference genome (Bos Taurus UMD 3.1), base quality calibration, variant calling and variant quality score recalibration have been described by Kadri et al. (2016). We selected 12,735,685 bi-allelic SNP located on autosomal chromosomes. After filtering out the SNPs with MAF (≥ 0.05), 4,279,775 SNPs for 743 sample remained for further analysis.

Genome-wide SNP data for each population was scanned for lethal SNPs which were defined as those: (1) are significantly deviated from the Hardy-Weinberg equilibrium, tested using test (p <7*10^-8) using PLINK v 1.9, (2) lack of one of the homozygous genotypes, and (3) with ≥5 expected homozygote cows for the minor allele. Since the technical and systematic issues may affect the allele frequency, we considered only SNPs which were common in at least two populations, to avoid any unbiasedness. The same procedure was applied on WGS data for detecting the lethal variants in this dataset. In order to identify the lethal genes underlying the embryonic mortality in cattle, the genic variants were mapped on the last genome assembly (ARS-UCD1.2 Bos taurus). To achieve high-confidence results, we applied a narrow filtration on the genes harboring the resulted SNPs. Therefore, only those genes meeting all the following criteria were considered as potential candidate lethal genes: (1) harboring SNPs common in at least two populations, or genes which harbor a SNP from the chip array genotypes and ≥5 SNPs from the WGS data, (2) are annotated to the functions related to embryo mortality and development, and (3) have been expressed constantly over the early embryo stages.

The Bovine HD arrays and WGS positions, which were initially based on UMD3.1 genome assembly, were converted to ARS-UCD1.2 reference genome using UCSC liftOver online tool (https://genome.ucsc.edu/cgi-bin/hgLiftOver) with the default options. The genes harboring potential candidate lethal SNPs were explored using Ensembl Biomart online tool and the Ensembl genes 111 database. Also, the the Ensembl Variant Effect Predictor (VEP) online tool was used to check the mutation type within the gene (synonymous or non-synonymous). The comparative genomics tool in Ensembl was used to determine the orthologous genes in placental mammals’ species (n=103) and all species available on Ensemble database (n=214). Enrichment analyses of gene ontology terms, including biological process (BP), cellular component (CC), and molecular function (MF) categories were conducted using various bioinformatics platforms, ensemble (https://www.ensembl.org/Bos_taurus/Info/Index), Genecards (https://www.genecards.org/), NCBI (https://www.ncbi.nlm.nih.gov/), UniProt (https://www.uniprot.org/) and DAVID (https://david.ncifcrf.gov/) tools.

In order to explore the expression patterns of the genes of interest across the various stages of embryonic development, the publicly free available transcriptome data on NCBI were used. Two transcriptome data including RNA-seq of embryonic tissues (trophoblast, endoderm, and mesoderm, collected on 14, 15, and 18 days of embryo development; under project PRJNA230971) and RNA-seq of fetal stem cells (under project PRJNA718333) were investigated. These datasets had been generated employing standard protocols on the Hiseq 2000 platform, facilitated by the FastTrack service at Illumina. The gene expression analysis of all annotated genes was conducted using CLC Genomics Workbench (version: 8.5.1) QIAGEN Bioinformatics (https://www.qiagenbioinformatics.com/), including quality control of reads, filtering out unwanted reads, removing adapter contamination, and mapping reads to the Bos taurus ARS-UCD1.2 reference genome. Expression levels were quantified in the RPKM scale. Additionally, we utilized NGS-based transcriptome data of 10 studies available on cattle Genotype-Tissue Expression (cGTEx, http://cgtex.roslin.ed.ac.uk/) atlas to assess the expression level of the genes, measured by RNA-seq, in different tissues over various developmental stages of embryo. The accession number of these transcriptome data and sample size are given in Supplementary Table S2 .

We also assessed the overlaps of the final gene list with identified QTLs underlying mortality related traits. To do this, the position of genes was inquired in Bovine QTL database available online at AnimalQTLdb (accessed in September 2024: https://www.animalgenome.org/cgi-bin/QTLdb/BT/index).

Par-3 Family Cell Polarity Regulator (PARD3), as member of the PARD family, plays crucial roles in asymmetric cell division and directing polarized cell growth (Cui et al., 2022). Studies have highlighted its involvement in key signaling pathways such as MAPK and Hippo, contributing to follicle growth and development in sheep (Wang et al., 2023a), morula formation during mouse embryonic development (Pfeffer, 2018), and the polarization process of mammalian embryonic blastomeres (Negrón-Pérez and Hansen, 2018). Additionally, research focused on carcass traits and gene networks has underscored PARD3's significant role in lipid and carbohydrate metabolism (Zhang et al., 2020). Notably, Cole et al. (2011) through a genome-wide analysis of predicted transmissibility (PTA) identified 5 QTLs (QTL IDs in AnimalQTLdb: 46882, 46852, 46866, 46835, and 46838), overlapping with the PARD3 gene, to be associated with stillbirth in Holstein cattle.

Bromodomain Adjacent Homology Domain-Containing Protein 1 (BAHD1) encodes a pivotal protein for chromatin organization and gene silencing, thereby contributing to cell differentiation and maintenance of homeostasis in mammals (Zhao et al., 2016). The Bromo Adjacent Homology (BAH) domain is essential for its function, facilitating gene silencing by mediating interactions with specific combinations of histones, such as H3K27me3 (Currie et al., 2024). Disruption of this interaction through point mutation leads to chromatin remodeling and increased histone acetylation at Polycomb gene targets. Mice with dysfunctional mutations in this pathway exhibits severe embryonic instability, thereby underscoring its critical role in normal embryo development (Fan et al., 2021).

In humans, Bierne et al. (2009) proposed that BAHD1 acts as a silencer by linking heterochromatin factors (e.g., HP1, MBD1, and HDAC5) to DNA-bound transcription factors (e.g., SP1), thereby functioning as a tumor suppressor through the silencing of cancer-related genes (e.g., IGF2). Studies in mice have shown that BAHD1 plays a crucial role in fetal growth and placental development. Lakisic et al. (2016) demonstrated that BAHD1-/- embryos exhibited structural changes in the placenta, including reduced weight and smaller placentas, with notable reductions in the junctional and cardio-fetal areas essential for nutrient exchange between the mother and fetus. Additionally, there was a decrease in the number of trophoblast glycogen (GC) cells, which is crucial for providing glucose for photo-placenta development.

Furthermore, Salilew-Wondim et al. (2018) suggested that BAHD1 is involved in DNA methylation in bovine embryos due to exposure conditions. They studied three groups of bovine embryos produced in vivo but cultured in vitro at different stages (2-, 8-, and 16-cell) and measured gene methylation and expression at the blastocyst stage. The authors reported that BAHD1 was one of the 28 (2%) genes differentially methylated and expressed in all groups. Additionally, BAHD1 exhibited a negative correlation between DNA methylation patterns and expression level, showing hyper-methylation but downregulation of mRNA expression in all blastocyst groups, likely due to suboptimal culture conditions.

Endoplasmic Reticulum-Golgi Intermediate Compartment Protein 2 (ERGIC2) is fundamental in early secretion pathways, particularly in the trafficking pathway responsible for transporting proteins, lipids, and other newly translated molecules. Located in the intermediate compartment of the endoplasmic reticulum-Golgi apparatus, ERGIC2 plays a crucial role in modulating ER-to-Golgi transport (Moreau et al., 2011). Additionally, Yu et al. (2014) reported that ERGIC2 is essential for efficient Wnt protein secretion and signaling. In their study, knocking down ERGIC2 in Xenopus through morpholino injection into two-cell embryos resulted in lethality before gastrulation, highlighting its indispensable role in embryonic growth (Yu et al., 2014).

Frizzled class receptor 3 (FZD3) plays a pivotal role in orchestrating various cellular processes, embryonic and postnatal development, particularly in the nervous and integumentary systems (Millar et al., 1999; Deardorff et al., 2001; Stuebner et al., 2010). Seigfried et al. (2017) stated that deficiency in FZD3 leads to severe developmental abnormalities, including impaired regulation of angiogenesis and disrupted neural tube closure. Deletion or absence of FZD3 has resulted in embryonic lethality and impedes neural tube closure (Kemp et al., 2007). Moreover, FZD3 has been reported to be involved in midbrain morphogenesis and cell migration (Stuebner et al., 2010), highlighting its crucial role in embryonic development, particularly in neural tube closure and organogenesis. The FZD3 regulates critical processes during embryonic development and the pre-implantation uterine environment in various species, including cattle, sheep, mice, and humans. High expression of FZD3 is observed in early gastrulation-stage bovine embryos, indicating its crucial role in embryonic tissue patterning and differentiation (Pfeffer et al., 2017). In general, FZD3 acts as a critical component of the Wnt signaling pathway, wherein Wnt proteins bind to Frizzled receptors (e.g., FZD3). It is essential for embryonic development and tissue homeostasis (Pan et al., 2022) and regulates various aspects of organogenesis, including nephrogenesis processes such as differentiation, proliferation, and morphogenesis (Cizelsky et al., 2014). In cattle, FZD3 participates in these pathways to regulate embryonic development and angiogenesis under hypoxic conditions (Verma et al., 2018).

IQ Motif Containing Nucleoporin (IQCN) is required to facilitate proper acrosome attachment during spermatogenesis (Dai et al., 2022). Dysfunction leads to fertilization failure, highlighting its significance in reproductive processes (Dai et al., 2022; Wang et al., 2023b). Wang et al. (2023b) reported that homozygous frameshift variants of this gene cause dysfunction in sperm, leading to male infertility in humans and mice. This dysfunctionality is mediated by abnormalities in the sperm oocyte activating factor PLCζ and head deformity. In a study on mice, Dai et al. (2022) showed that Iqcn-knockout (Iqcn-/-) mice can survive and grow, but are not fertile due to abnormal acrosome structures. From this study, it can be deduced that this gene does not directly cause embryo lethality. Since this gene is protein-coding and located on an autosome (chromosome 7), it likely has functions in ova as well. Moreover, the effect of the homozygous status of this gene on ova function, similar to its effect on sperm, remains unclear. To our knowledge, there are no published reports in this regard; therefore, further investigations into the function of this gene in bovine semen, ova, zygotes, and embryos are highly recommended.

Prokineticin receptor 1 (PROK1) plays a pivotal role in regulating endometrial angiogenesis and secretory function during the estrous cycle and early gestation. Its heightened expression is linked with monocyte recruitment, crucial for successful pregnancy establishment (Baryla et al., 2023). PROK1's involvement spans various species, as evidenced by its increased expression during the mid-luteal phase in pigs, humans, and cattle. Notably, its expression dynamics differ in the late luteal phase, suggesting a role in monocyte recruitment to the regressing corpus luteum (Baryla et al., 2023). In Spanish beef cattle, PROK1 is emerged as a fertility-related gene (González-Rodríguez et al., 2016).

During implantation and early placentation, PROK1 mRNA expression undergoes dynamic changes across species. In humans and mice, mRNA levels peak during early pregnancy, whereas in pigs, upregulation occurs in trophoblasts during implantation and early placentation (Hoffmann et al., 2007).

Phosphatidylcholine Transfer Protein (PCTP), also known as StARD2, is an important member of the steroidogenic acute regulatory protein-related transfer (START) domain superfamily (Kanno et al., 2007a, b). This family encompasses proteins involved in lipid transfer and metabolism, with PCTP facilitating the transfer of phospholipids, especially phosphatidylcholine (PtdCho), between cellular membranes. The evolutionarily conserved START domain, to which PCTP belongs, underscores its functional importance across species (Schrick et al., 2004). STARD2, with its minimal START domain structure, functions as a cytosolic phosphatidylcholine transfer protein and crucial for phospholipid trafficking between membranes (Soccio et al., 2002). Studies have highlighted PCTP's significance in various biological processes. The PCTP gene overlaps with quantitative trait loci (QTLs) associated with stillbirth reported by Cole et al. (2011) in Holstein dairy cattle, indicating its role in reproductive health. Additionally, PCTP's expression throughout mice embryo development stages suggests its involvement in cell growth and differentiation (Geijtenbeek et al., 1996).

The SH3 Domain-Containing GRB2-Like Protein 1 (SH3GLB1) gene encodes a multifunctional protein, also known as endophilin B1, which drives membrane curvature. This protein plays a role in apoptosis, a process involving dying cells characterized by cytoplasmic shrinkage, membrane blebbing, chromatin condensation, and fragmentation into membrane-bound vesicles or apoptotic bodies (Peterson et al., 2015). Additionally, SH3GLB1 contributes to regulating molecular pathways associated with spermatogenesis and male fertility (Henderson et al., 2011; Peterson et al., 2015). SH3GLB1 activity is crucial in the early development of the auditory, cardiovascular, and muscular systems in zebrafish. Therefore, abnormal SH3GLB1 function might lead to imbalanced mitophagy (Gao et al., 2023). Furthermore, SH3GLB1 acts as a pivotal tumor suppressor, essential for preventing chromosomal instability and suppressing resistance to apoptosis (Takahashi et al., 2013). Deficiency in SH3GLB1 could result in embryonic lethality and/or misregulated lymphoma cell homing through up-regulation of the expression of the Mcl-1 gene (Moshnikova et al., 2006).

Ras Association Domain Family Member 5 (RASSF5), also known as novel Ras effector 1 (NORE1), belongs to the Ras-association domain family and serves as a crucial regulator of cellular homeostasis. It exerts its functions through growth suppression and regulation of apoptosis (Moshnikova et al., 2006). RASSF5 plays multifaceted roles in cellular processes and is thus act as an essential player in maintaining cellular homeostasis. Its predominant cytoplasmic localization allows it to interact with cytoskeletal proteins, thereby participating in microtubule stabilization and potentially regulating mitosis. Furthermore, RASSF5's interaction with KRAS highlights its unique role in Ras-induced pro-apoptotic pathways, distinguishing it from other RASSF family members. The existence of multiple isoforms of RASSF5, generated through alternative splicing or promoter usage, underscores its diverse functions in cellular regulation, including cellular motility, cell cycle control, apoptosis, stability of microtubules, and adaptive immune defense (Ehrkamp et al., 2013).

Park et al. (2010) observed significant protection from TNF-α-mediated apoptosis of the liver in mice due to the inactivation of RASSF5. They concluded that RASSF5 acts as a tumor suppressor and plays a direct role in the activation of the proapoptotic kinase Mst1 after TNF-α stimulation. Moreover, RASSF5 and RASSF1 are homologous and play similar roles in regulating early embryonic differentiation in mice.