The present study used frozen semen from 12 NR bulls, divided into two groups based on fertility. The bulls were selected from a database of 470 NR bulls, born between 2010 and 2014, all with at least 500 registered first inseminations (AI). Fertility was recorded as the non-return rate (NRR) at 56 days after first AI of the inseminated cows, and varied from 49.3 to 80.5 with an average of 72.5 (s.e. = 3.5). The 12 selected bulls had a record of 661–901 first artificial inseminations and represented the highest and lowest fertility among all registrations. Bulls in the high fertility (HF) group had a NRR of 78.7–80.5 (n = 6), while the low fertility group (LF) had an NRR of 49.3–62.1 (n = 6). The reason for the difference in fertility was not known, and the semen was no longer in commercial use. The semen had passed standard testing requirements performed by Geno SA ¹ , the breeding organization for NR cattle, before commercial use; Macroscopic evaluation, a concentration threshold of 390 million cells per ml, at least 70 and 50% motile spermatozoa pre-freezing and post-thawing, respectively, and a threshold of at least 83–90% morphologically normal spermatozoa depending on the specific deviation (personal communication, Geno SA).

For embryo production in vivo, we used 14 NR heifers in the age span of 12–15 months. In order to reduce individual differences and the maternal effect to a minimum, all animals came from the same genetic line with 28 years of targeted breeding, with high fertility and a low occurrence of clinical mastitis as target traits (Heringstad and Larsgard, 2010). They were free from disease or medical treatments according to their health records for the last 6 months before the sampling, they had no earlier inseminations and at least two visually registered estrus cycles. They were held indoors in the same free-range barn and received the same feeding throughout the study. Their body condition scorings were considered normal with an individual variation of 3.5–4.0, using a scale from 1-5 with increments of 0.25 (Edmonson et al., 1989), modified and adjusted for NR according to Gillund et al. (1999).

The ethical approval for the present study was provided by the Norwegian Food Safety authority with approval ID 11732. The combination of bull and heifer was randomized with block randomization; heifers were listed according to age, and every second animal was appointed a randomly chosen HF bull or LF bull, respectively. Randomization was performed using Sergeant, ESG, 2018, Epitools Epidemiological Calculators, Ausvet ² . Semen from the two bulls with the lowest NRR were appointed to be used for two different heifers.

Embryo production was performed at The Animal Production Experimental Centre, NMBU in Ås, Norway in the spring of 2017. A protocol for synchronization and superovulation was developed for young NR heifers. The animals were synchronized with an intramuscular (i.m) injection of 2 ml cloprostenol 0.25 mg/ml (Estrumat vet., MSD Animal Health, Intervet International B.V., Nederland) twice at a 12-day interval and the following heat was visually detected. On day 9 after the first signs of standing heat, the 4-day administration with decreasing amounts of follicle stimulating hormone (FSH; Follitropin 500 IE and lutinizing hormone 500 IE, Pluset^® vet, Laboraotiros Calier, Barcelona, Spain) started.; Two i.m. injections were given daily (day one 2.0 ml, day two and three 1.5 ml and day four 1.0 ml). On the fourth day, a 2.0 ml i.m. injection of cloprostenol 0.25 mg/ml was administered in the morning and in the evening to induce luteolysis. The heat occurred after 2 days and AI was performed two or three times with 12 h in between depending on length of heat behavior.

Embryo flushing was performed 7 days after first AI. The ovaries were controlled for superovulation response using both manual palpation and rectal ultrasound, followed by an administration of epidural anesthesia. A Foley catheter was used to flush each uterine horn at least five times with ViGRO™ Complete Flush Solution (Vetoquinol, Lure Cedex, France. Previously Bioniche Animal Health, United States). The fluid was retrieved in an embryo collection filter (Emcon, Panningen, Netherland).

The embryos were transferred from the filter to a petri dish with SYNGRO ™ holding solution (Bioniche Animal Health, United States) and evaluated under a magnifying loupe at 40x magnification by three veterinarians according to the IETS guidelines (IETS-manual 3rd edition, IETS bovine in vivo embryo slide set tutorial, 2010). All embryos were then cleaned with phosphate buffered saline, moved to separate sterile mini tubes containing 1U/µl RNAsin in nuclease free water (RNasin Ribonuclease Inhibitor, Promega corporation, WI, United States), instantly frozen in liquid nitrogen, and stored in a −80°C freezer.

Embryos derived from HF sires were referred to as HF embryos, while those from the LF sires were referred to as LF embryos.

Embryos with first grade quality and IETS stage 5 (early blastocyst), 6 (blastocyst) or 7 (expanded blastocyst) from both HF and LF groups were selected for further processing. RNA isolation was performed using the RNAqueous-Micro Total RNA isolation Kit (Thermo-Fisher Scientific, MA, United States) according to the producer’s instructions with some modifications; In order to break the zona pellucida and cell walls, embryos placed in 100 µl of the kit lysis solution were submitted to five cycles of 2 min freeze in liquid nitrogen and 2 min thaw in a 50°C water bath, followed by an incubation at 42°C overnight. The next day 50 µl of ethanol was added and isolation was continued according to the protocol. RNA was eluted twice in 6.0 and 6.5 µl elution solution, and the pooled eluate was treated with DNase1 as described by the producer. Production of cDNA was performed using the SMART-Seq v4 Ultra Low Input RNA Kit for sequencing (TaKaRa Bio Europe, Göteborg, Sweden) according to the producer’s protocol, using 10.5 µl of input RNA and 18 cycles of amplification. The Agencourt AMPure XP kit (Beckman Coulter, IN, United States) was used to purify the amplified cDNA as described in the Smart-Seq kit protocol. The cDNA concentration and quality were measured by Qubit fluorometer using the Qubit dsDNA HS Assay Kit (Thermo Fisher Scientific, MA, United States) and Agilent TapeStation D1000 using High Sensitivity reagents (Agilent, Santa Clara, CA, United States), respectively.

cDNA samples with concentrations between 0.262 ng/µl and 17.1 ng/µl (Supplementary Table S1) and sufficient quality according to the TapeStation profiles (Supplementary Figure S1) were selected for sequencing. In total, cDNA from 24 embryos were sent for sequencing at the Norwegian Sequencing Centre ³ , i.e., 13 embryos from four HF bulls and 11 embryos from three LF bulls. The total distribution between embryo stages 5, 6 and 7 were 1, 12 and 11 embryos, respectively. In the LF group, the embryo stages were 6 (n = 5) and 7 (n = 6), and in the HF group, the stages were 5 (n = 1), 6 (n = 7) and 7 (n = 5) (Supplementary Table S1).

Sequencing libraries from cDNA were prepared using SMARTer ThruPLEX DNA-prep kit (TaKaRa Bio United States Inc., San Jose, CA, United States) using unique indexes. Libraries were pooled and 150 bp paired end sequencing was performed on one lane of HiSeq 4,000 (Illumina, United States). Raw reads were processed using BBDuk (part of BBMap v34.56) (Bushnell, 2014) (parameters: ktrim = r k = 23 mink = 11 hdist = 1 tbo tpe qtrim = r trimq = 15 maq = 15 minlen = 36 forcetrimright = 149) to remove/trim low-quality reads and adapter sequences. Cleaned reads were aligned to the Bos taurus genome (ARS-UCD1.2; ENSEMBL release 95) using hisat2 v2.1.0 (Pertea et al., 2016) and the resulting sam files were converted to bam format using samtools v1.2. Reads mapping to the genes (ARS-UCD1.2; ENSEMBL release 95) were counted using featureCounts v1.4.6-p1 (Liao et al., 2014). Differential expression analysis was performed using DESeq2 v1.22.2 package (Love et al., 2014) in R v3. In brief, the counts were normalized followed by outlier detection (Cook’s distance), dispersion estimation (fittype: parametric) and statistical testing (hypothesis testing: Wald test). Independent filtering was performed which discarded 14,355 genes due to very low count values. Finally, multiple testing was performed using the Benjamini-Hochberg method. Genes with the adj.p-value less than 0.05 were considered to be significantly differentially expressed.

A gene set enrichment analysis of the DEGs was performed using the g:GOSt function in g:profiler version e104_eg51_p15_3922dba (Raudvere et al., 2019), which also integrated results from KEGG Pathway database, WikiPathways and Reactome. Separate analysis was performed for the genes highly expressed in the HF embryos and for the genes highly expressed in the LF embryos. A custom gene list consisting of 12,826 genes that were detected (adj.p-value not equal to “NA”) in the DE-analysis was used as background list for statistical domain scope (Supplementary Table S2). To correct for multiple testing, we used the Benjamini-Hochberg FDR algorithm and a threshold level of 0.05 for significance. For further functional analysis and visualization, we used the online version of Pathview v1.3.2 ⁴ (Luo et al., 2009; Luo and Brouwer, 2013; Luo et al., 2017).

The bovine genome is not as well studied and annotated as the human genome. To further enrich results interpretation, DEGs and the background list for statistical domain were converted to human ENSEMBL orthologs from bovine ENSEMBL IDs, by g:profiler and run though gene ontology analysis in g:profiler and pathway analysis in Pathviews, as described above.

All the heifers responded well to the superovulation protocol with normal size ovaries and no un-ovulated follicles. Sampling was normal in all heifers except for one, where collection from one of the uterine horns was not performed due to practical challenges. In total, 73 embryos were collected from 8 heifers inseminated with 6 LF bulls (=30 embryos) and 4 heifers inseminated with 4 HF bulls (=43 embryos), with an individual distribution of 1–19 embryos per animal. In the LF group, 17 of the 43 embryos collected had developed to the blastocyst stage, compared to 20 out of 30 in the HF group. From two heifers, no embryos could be found at collection. An overview of embryos sequenced from the different sires is given in Table 1. After the collection of embryos, Geno SA conducted an independent investigation of bulls that had previously been included in their breeding program, by analyzing genomic profiles, i.e., SNP data. That process revealed that one bull (K) had a deletion on chromosome 12.

The sequencing resulted in an average of 12.5 million paired end reads per sample and more than 98% passed quality check. Read information and alignment statistics are provided in Supplementary Table S1. Raw sequence reads were uploaded to the NCBI SRA database with bioproject ID number: PRJNA762262.

The distribution of gene expression in the individual blastocysts was visualized by principal component analysis plots (Supplementary Figure S2). Two embryos (number 48 and 66) were shown as outliers and were removed from further analyses. These two samples also had an inferior cDNA quality according to the TapeStation profiles, which justified our decision. Due to this, the lowest cDNA concentration for samples analyzed in DESeq2 was 1.64 ng/µl. Embryos of IETS stages 5, 6 and 7 did not show strong signs of clustering within these groups in any of the three dimensions. The embryo of IETS stage 5 (from a HF bull) did not stand out as outlier or represent the highest or lowest value in any aspect of the PCA. Embryos from LF bulls compared to HF bulls clustered in the second dimension, but not in the first and third dimension.

Among the 14,744 annotated genes that were detected during the analysis, there was a significant difference in the expression of 62 genes; 28 genes had a higher expression in the HF group, and 34 genes a higher expression in the LF group (Log2-foldchange > 0.64 adj.p-value < 0.05), see Table 2.

An overview over the significant pathways and gene ontology (GO)-terms found in g:profiler is provided in Table 3. The full DE-list revealed three significantly affected pathways in Pathview, represented by 9 different genes; terpenoid backbone synthesis (4 genes represented), insulin signaling pathway (3 genes represented) and leukocyte transendothelial migration (2 genes represented).

Genes that were highly expressed in the LF group (n = 34) were associated with 29 GO-terms for biological processes and pathways, most prominently biosynthetic process- and metabolism of sterol, steroids, isoprenoids and cholesterol. Pathway analysis in KEGG, Reactome and WikiPathways through g:profiler showed enrichment of terpenoid backbone synthesis, and metabolism and biosynthesis of cholesterol, steroids and lipids. The individual transcript showing the biggest difference between the groups, with higher expression in the LF group, was mRNA coding for bactericidal permeability increasing protein (gene symbol BPI, log2fold change = 10.897).Another transcript, with higher expression in the LF group was mRNA coding for neutrophil cytosol factor 1 (gene symbol NCF1, log2 fold change = 5.812), a protein engaged in leukocyte transendothelial migration and a NADPH oxidase regulator. The leukocyte transendothelial migration pathway was further represented by the gene claudin 9 (CDLN9, log2 fold change = 3.688), the third most up-regulated transcript in the LF group and coding for cell adhesion molecules (CAMs).

Genes that were highly expressed in the HF group (n = 28) were associated with only one GO-term through Corum in g:profiler, the G protein complex. The biggest foldchange was represented by a non-annotated gene (ENSBTAG00000031825, log2fold change = 8.016) followed by GIMAP4 (log2fold change = 5.850) and HLX (log2fold change = 5.724).

Out of all significant DEGs in the DE-list, nine genes were described as non-annotated genes, of which seven did not have any pathway and GO-terms associated. Three of the non-annotated genes were highly expressed in LF embryos, and six in HF embryos (Table 2). Some of these represented a very high level of difference between the LF and HF groups, with the highest log2fold change of 8.016 and 2.823 in the HF group, and 2.994 and 2.425 in the LF group.

To enrich results interpretation, we used human orthologs, as the human genome is better studied and annotated. Out of the 62 bovine genes from the DE-list, g:profiler found 58 human orthologues (Supplementary Table S3). Only four genes did not have a human orthologue. The Pathview program pathways represented by these orthologues matched the pathways for the bovine genome. Out of the 29 GO-terms or pathways linked to genes highly expressed in the LF group, 27 were enriched in the human orthologs. Out of these 27, 16 terms or pathways were each associated with one or more genes in the human orthologs compared to the original material. In 6 other enriched processes or pathways, the same number of genes represented both the cattle and human ortholog outcome. Fourteen terms or pathways that did not appear in the original material, were shown in the analysis of orthologs. These pathways were strongly related to, or represented a higher hierarchy of, the same pathways that had already been identified. Only two pathways in the human orthologs, Omega-9 FA synthesis and Cuteneous photosensitivity, were not related but these were only represented by two and three genes, respectivley. For the genes that were highly expressed in the HF group, the analysis of human orthologs did not give any GO-terms or pathways.

Embryos derived from LF bulls showed a higher genetic expression corresponding to a more active metabolism. These results are in correlation with earlier literature, proposing an association between a high survival rate in the embryo, and a lower level of metabolism (Leese, 2002; Baumann et al., 2007; Leese et al., 2007). One of the pathways that was highly expressed in the LF group was the terpenoid backbone biosynthesis (Kanehisa and Goto, 2000; Kanehisa, 2019; Kanehisa et al., 2021) (Supplementary Figure S3), which initiates the production of sterol isoprenoids, such as cholesterol, and non-sterol isoprenoids (Buhaescu and Izzedine, 2007; Miziorko, 2011). The products derived from this pathway play an essential role in various cellular processes such as cell growth and differentiation, and cell signaling (Goldstein and Brown, 1990). One of the continuations of this pathway; steroid biosynthesis (Kanehisa and Goto, 2000; Kanehisa, 2019; Kanehisa et al., 2021) (Supplementary Figure S4) was also highly expressed in the LF group. The increased activity in these connected pathways was represented by nine transcripts in the present study. The first three, HMGCS1, HMGCR and PMVK all code for enzymes in the production of mevalonate. HMGCR has been described as a rate limiting enzyme and is of major importance for the entire downstream process (Goldstein and Brown, 1990). Both HMGCR and HMGCS1 underlie well-studied, vast mechanisms of regulation. One of those mechanisms is an end-product feedback system that allows any absence of sterol isoprenoids to activate the transcription of the HMGCR gene through a family of transcription factors called sterol regulatory binding proteins (SREBP) (Brown and Goldstein, 1997; Buhaescu and Izzedine, 2007).

The next two genes represented in these pathways, DID1 and FDPS, encode enzymes that catalyze the further descendance of metabolites towards steroid biosynthesis while four other genes, FDFT1, MSMO1, SC5DL and CYP17A1 encode enzymes that lead the metabolism down to the biosynthesis of cholesterol (and several other steroids) and steroid hormones. One study demonstrated that an increased SREBP activity not only acted on HMGCR and HMGCS1, but also increased the mRNA expression of several enzymes along the entire pathway of cholesterol production (Sakakura et al., 2001). Based on this information, we speculate that the high expression of the mentioned enzymes in the LF embryos, could be a result of any dysfunction of the pathways related to sterols, or any exaggerated degradation or demands of its products. Cholesterol is essential for the developing embryo, as it forms part of the cell membrane, and acts in cell signaling crucial for developmental patterning, in collaboration with the hedgehog gene family (Porter et al., 1996; Roux et al., 2000). It can be toxic in too large quantities and its production entails high metabolic costs for cells to produce. Hence, its production is under strict regulation (Sharpe and Brown, 2013). The complete knockout of the HMGCR gene in mouse embryos resulted in the recovery of morphologically normal blastocyst but no later developmental stages. This suggest that at least some of HMGCR products are essential for development from the blastocyst stage, either prior to implantation, or for the implantation process itself (Ohashi et al., 2003). Equally, CYP17A1 disruption leads to early embryonic lethargy in murine embryos (Bair and Mellon, 2004). Hence, a suboptimal level of cholesterol and/or its precursors in the LF embryos would lower their developmental potential. This marks a difference between LF and HF embryos that could explain at least part of the reason for low field fertility in the LF bulls.

Interestingly, the sterol biosynthetic process and cholesterol pathway have been highlighted in the comparison between morphologically similar in vitro and in vivo produced embryos (Driver et al., 2012). In vitro embryos have a reduced developmental potential from the zygote to blastocyst stage, and a lower success in embryo transfer (Rizos et al., 2008). Driver et al. (2012) performed a transcriptome study in stage 7 blastocysts, where the in vitro group had an increased expression in 11 genes related to the cholesterol pathway. These genes included HMGCS1, HMGCR, PMVK, IDI1 and FDFT1 which are all identical to our findings. The present study only analyzed in vivo produced embryos, but equal to the study by Driver et al. (2012), it compared the transcripts of embryos with a hypothetical difference in developmental potential. The fact that the results of the two studies are in agreement, confirms the central role of the mentioned pathways for successful embryo development.

The leukocyte transendothelial migration pathway (Kanehisa and Goto, 2000; Kanehisa, 2019; Kanehisa et al., 2021) (Supplementary Figure S5) was represented by two of the most highly expressed genes in the LF group. Interestingly, the pathway was also represented by one gene with a high expression in the HF group. Earlier literature pointed to several similarities between leukocyte transendothelial migration and human implantation, stating that both processes use the same mechanisms of adhesion, molecular interaction and migration (Genbacev et al., 2003; Dominguez et al., 2005). Hence, one explanation for the low bull fertility, might be through an effect on the control of implantation. The leukocyte transendothelial migration pathway was represented by NCF1 and CLDN9 in the LF group, and CLDN10 in the HF group. Claudins (CLDN9 and CLDN10) also have a role in embryo development, independent of their role in this pathway. The claudin gene family forms part of tight junctions, which are transmembrane compounds with functions in the maintenance of apical-basal polarity and cell adhesion (Gupta and Ryan, 2010). Tight junctions are crucial for morphogenesis (Furuse and Moriwaki, 2009), and a loss of function-study revealed that some claudins are essential for the formation of the murine blastocyst (Moriwaki et al., 2007). In a review of claudin function in embryogenesis, the authors hypothesized that the combined expression of claudin, or the “claudin signature,” is critical to embryonic tissues (Gupta and Ryan, 2010). The specific importance of a high expression of NCF1 is uncertain. NCF1, also known as p47phox, takes part in the production of reactive oxygen species (ROS) through its role in NADPH oxidase (Babior, 2004). The change in expression of NCF1 might point to a change in the redox state (reviewed by Harvey et al., 2002) in any of the two embryo groups. Changes in the embryo redox state through the limited accumulation of ROS is naturally occurring, enabling developmental progress in the embryo (Dennery, 2007). However, it also controls programmed cell death (Pierce et al., 1991), and in excess, oxidative stress is embryotoxic (Dennery, 2007). Consequently, one possible causative factor of the poorer outcome for the LF embryos could be through a lower competence in the regulation of redox activity. This hypothesis is supported by the findings in one earlier study of the sire’s contribution to embryo development (Kropp et al., 2017).

Another pathway differing between the HF and LF groups was the insulin signaling pathway (Kanehisa and Goto, 2000; Kanehisa, 2019; Kanehisa et al., 2021) (Supplementary Figure S6). Two genes were abundant in the HF group: SOCS1 and BAD. SOCS1 encodes an enzyme which inhibits the action of the insulin receptor (Mooney et al., 2001), hence suppresses the full signaling pathway, which contains PI3K-Akt signaling. PI3K-Akt signaling has a central role in embryo survival, regulating differentiation and cell growth, proliferation, anti-apoptosis and calcium metabolism (Leese and Brison, 2015, p.184). The enzyme encoded by SOCS1 also has several roles in the negative feedback mechanism of cytokine signaling (Krebs and Hilton, 2001). Cytokines are produced by the embryo itself, as well as the female reproductive tracts as a mediator of maternal-embryo communication. They affect a vast range of processes, again related to cell differentiation and cell survival (Leese and Brison, 2015, p.173–174). SOCS1 expression is a product of both interferon-γ and interleukin-4 and the protein encoded by SOCS1 has a negative feedback on these cytokines (Fujimoto and Naka, 2003). Interferon-γ has several important roles in embryo development, but excess production is detrimental (Leese and Brison, 2015, p.193), implying that a well-functioning regulatory mechanism is beneficial for the embryo. BAD, on the other hand, encodes an enzyme that is inhibited by the activation of the insulin signaling pathway, so a higher expression could be a result of the elevated activity in SOCS1. BAD is an antagonist of apoptosis, which is interesting since a higher apoptotic cell ratio indicates a lower developmental competence in the embryo (Maddox-Hyttel et al., 2003). Certainly, the lower degree of apoptosis in the HF group would make a logical explanation for a higher developmental potential. The apoptotic cell ratio is inversely correlated to early cleavage in zygotes (Byrne et al., 1999), which is again positively correlated to bull fertility (Ward et al., 2001). However, one study intended to demonstrate a direct association between bull fertility and apoptotic cell ratio, but failed to do so (Vandaele et al., 2006).

Some of the of the most highly expressed genes in the two groups were not represented in any pathways. GIMAP4, encoding a small GTPase active in the immune system (Heinonen et al., 2015) was highly expressed in the HF embryos. It regulates cytokine secretion in the early human CD4⁺ Th lymphocytes and initiates the secretion of interferon-γ (Heinonen et al., 2015). GIMAP4 is also an important regulator of calcium signaling (Schnell et al., 2006), a process which in recent years has been shown to have several functions in the pre- and peri implantation period (Leese and Brison, 2015, p.158–164). To the authors knowledge, the exact role of GIMAP4 in embryo development has not yet been defined. HLX was another highly expressed gene in the HF embryos. Similar to SOCS1, HLX is also a regulator of cytokines, allowing trophoblast proliferation and the development of the placenta (Rajaraman et al., 2010). This, again, proposes an association between the paternal contribution and the mechanism for implantation in our material.

Another interesting finding in the HF group was the higher expression of POUF51. This gene encodes the transcription factor Oct4, which is essential for pluripotency and the formation of an intracellular matrix (Nichols et al., 1998). Our results could denote that HF embryos are more competent in this matter, and are more likely to develop beyond the blastocyst stage. The highest DE seen in the HF group was of a non-annotated gene; ENSBTAG00000031825. Its homologue C19orf12 (e value 0.0) has recently been shown to be important in neuronal development in zebrafish embryos, as a downregulation of the gene had severe effects on brain morphology and resulted in embryo death before day 7. Its function was suggested to be related to lipid metabolism even though the cellular mechanism is poorly understood (Mignani et al., 2020). The higher expression of this gene in our HF embryos is indeed an interesting finding that could explain differences in bull fertility but requires more investigation.

BPI was the single gene showing the highest DE in the LF group. The gene product is a lipid-transfer protein with the capacity to neutralize endotoxin. In humans, it is produced by neutrophils and the epithelial lining of mucosa as part of an antimicrobial defense mechanism (Schultz and Weiss, 2007). Proteins encoded by BPI and the BPI-like PLUNC genes from the same superfamily, have been found in the seminal plasma of rams (Soleilhavoup et al., 2014; van Tilburg et al., 2020), and the spermatozoa membrane of mice (Zhou et al., 2014) and rodents (Yano et al., 2010), and are hypothesized to have a role in the sperm-oocyte fusion process (Li et al., 2013). Even if BPI were to be identified in the semen of bulls, the finding in the embryos of the present study is likely not a direct result of paternal transcripts deposited to the oocyte at fertilization, since these transcripts start to degrade at EGA and should not be abundant in the analyses of embryonic gene expression at the blastocysts stage (Graf et al., 2014; Jiang et al., 2014). Neither BPI in cattle, nor its human orthologue has been assigned to any pathway, and to the authors knowledge, the role of BPI in embryo development has not been reported in earlier studies. However, one RNA-seq study that compared different stages of in vivo produced cow embryos found that BPIAF1 (a BPI-like PLUNC gene) is a hub gene in blastocysts. This information was validated with literature of human and mice blastocysts (Jiang et al., 2014). Further research is necessary to study the role of BPI in embryo development.

Defining the significance of our interpretations to our findings is challenging, given that cell signaling in embryo development is controlled by a vast number of processes with overlapping actions and shared receptors (Leese and Brison, 2015, p.180). Superovulation could have altered the gene expression of some genes in the current study compared to a normal in vivo produced embryo (Mundim et al., 2009), but this alteration applies to both LF and HF embryos and should not affect the differences between the groups. Moreover, it is not certain that the death of the conceptuses from LF bulls occurs at the blastocyst stage or before implantation, even if most embryo death probably occurs before day 8 after conception (Diskin et al., 2016). Equally, although fertilization failure is not the main problem of non-successful coupling (Sreenan and Diskin, 1986; Diskin et al., 2016), we cannot rule out that the LF bulls in our study might have had a weak fertilization capacity. However, two recent studies on early embryo development in high and low fertility bulls, showed no difference in fertilization rate (Kropp et al., 2017; O'Callaghan et al., 2021), while one showed a difference in the development until day 7 between the two groups (O’Callaghan et al., 2021). In the present study, it is uncertain whether embryos were lost in the in vivo collection process, or whether the embryos that did not develop to the blastocyst stage, had the potential to do so. Equally, although the maternal effect was reduced by using very similar heifers with equal living conditions, each embryo was inevitably affected to some extent by the individual differences in the genetics of the heifers. Another limitation of the study was that we had to choose a mix of embryos of IETS stages 5, 6 and 7, which could have affected the relative expression of some genes. However, despite of this and individual differences, the distribution in the PCA plots supports the argument that the selected embryo stages are sufficiently uniform to study the differences between the HF and LF bulls. Equally, it would have been interesting to include embryos of all qualities in the two groups, and not only the highest quality. However, it is well known that gene expression varies between individuals and if this variation is too large it may obscure potential differences between groups of individuals. Embryos of different quality are expected to differ in expression profiles. Hence, in order to keep the individual variation within the groups to a minimum, only embryos of the highest quality were used.

To our knowledge, this is the first study to investigate the transcriptome of in vivo produced embryos for the influence of paternal field fertility. Comparing our results to a similar study that investigated bull field fertility and embryo transcriptomic profiles in in vitro produced blastocysts (Kropp et al., 2017), we found few evident similarities in the genes or pathways that were differentially expressed. This underlines the importance of studying in vivo produced embryos even though it is a challenging approach. To understand the sire’s effect on the embryo, one needs to study a complex relationship between several factors such as aspects of the semen and spermatozoa, molecular genetics and epigenetics. The bull’s effect on the embryo as reported in the present study, might be caused by either bull DNA, or regulations by proteins, transcriptome or epigenetic factors deposited in the oocyte at fertilization. Regardless of the type of contribution, it is certain that it originates from the spermatozoa. Therefore, further epigenetic investigations of both spermatozoa and the resulting embryos from the same bulls would be highly interesting.

The present study adds important information to the current understanding of the paternal influence on the genetic components in embryo development. Although the field of bull fertility has received clear attention and progress during the last decade, further research is needed to clarify this complex matter, with the goal to find biomarkers that aid the prediction of bull fertility.