Comparative Transcriptomic Analysis of Spermatozoa From High- and Low-Fertile Crossbred Bulls: Implications for Fertility Prediction

Crossbred bulls produced by crossing Bos taurus and Bos indicus suffer with high incidence of infertility/subfertility problems; however, the etiology remains poorly understood. The uncertain predictability and the inability of semen evaluation techniques to maintain constant correlation with fertility demand for alternate methods for bull fertility prediction. Therefore, in this study, the global differential gene expression between high- and low-fertile crossbred bull sperm was assessed using a high-throughput RNA sequencing technique with the aim to identify transcripts associated with crossbred bull fertility. Crossbred bull sperm contained transcripts for 13,563 genes, in which 2,093 were unique to high-fertile and 5,454 were unique to low-fertile bulls. After normalization of data, a total of 776 transcripts were detected, in which 84 and 168 transcripts were unique to high-fertile and low-fertile bulls, respectively. A total of 176 transcripts were upregulated (fold change > 1) and 209 were downregulated (<1) in low-fertile bulls. Gene ontology analysis identified that the sperm transcripts involved in the oxidative phosphorylation pathway and biological process such as multicellular organism development, spermatogenesis, and in utero embryonic development were downregulated in low-fertile crossbred bull sperm. Sperm transcripts upregulated and unique to low-fertile bulls were majorly involved in translation (biological process) and ribosomal pathway. With the use of RT-qPCR, selected sperm transcripts (n = 12) were validated in crossbred bulls (n = 12) with different fertility ratings and found that the transcriptional abundance of ZNF706, CRISP2, TNP2, and TNP1 genes was significantly (p < 0.05) lower in low-fertile bulls than high-fertile bulls and was positively (p < 0.05) correlated with conception rate. It is inferred that impaired oxidative phosphorylation could be the predominant reason for low fertility in crossbred bulls and that transcriptional abundance of ZNF706, CRISP2, TNP2, and TNP1 genes could serve as potential biomarkers for fertility in crossbred bulls.


INTRODUCTION
Male fertility is of great importance in dairy cattle breeding industry because semen from a single bull is utilized to breed several thousands of females. The most accurate method for testing the bull fertility is insemination of many fertile females, but this method is time-consuming and expensive, and only a lesser number of males can be tested at any given time (Gillan et al., 2008;Kumaresan et al., 2017). Therefore, it is obligatory to depend on semen fertility prediction techniques (Rodriguez-Martinez, 2003). Available semen evaluation techniques and their foretelling ability of bull fertility are highly variable (Rodríguez-Martínez, 2006;Rodriguez-Martinez and Barth, 2007). Existing semen evaluation assays estimate only few structural attributes of spermatozoa rather than their functional attributes that are having considerable correlation with fertility (Graham, 2001;Kumaresan et al., 2017), which limits the accuracy of bull fertility prediction using these assays. Sperm should comprehend and express many vital attributes in a temporal manner to successfully fertilize the oocyte (Gillan et al., 2008). Bull semen consists of heterogeneous cohorts of subpopulations, as a result of different spermatogenic waves (Rodríguez-Martínez, 2006), and therefore, all the spermatozoa in a given ejaculate are not uniform in terms of their fertilizing capacity. Therefore, it is essential to understand the sperm molecular differences between high-and low-fertile males so that fertility signature molecules could be identified for development of bull fertility prediction tools.
Increasing evidences indicate that males produced by crossing Bos taurus males with Bos indicus females suffer with high rates of infertility/subfertility; however, the etiology remains poorly understood. Reportedly, a greater proportion of crossbred bulls were culled due to infertility/subfertility and poor semen quality; the average ejaculate rejection rate was around 55% in crossbred bulls (Tyagi et al., 2000;Mukhopadhyay et al., 2010;Vijetha et al., 2014a,b). Therefore, our team has been working toward the aim of decoding the reason behind infertility/subfertility in crossbred bulls, and we found the differences in the proportion of Sertoli cells (Tripathi et al., 2015), metabolomic profile of spermatozoa (Saraf et al., 2020), proteomic profile of seminal plasma , spermatozoa , spermatogenic and Sertoli cells (Tripathi et al., 2014), and transcriptomic details of testicular tissue (Elango et al., 2020) between crossbred and purebred cattle. However, transcriptomics of crossbred bull sperm is underexplored, although recent studies revealed the relationship of sperm transcripts with sperm function and fertility of purebred bulls (Card et al., 2017;Ren et al., 2017;Burl et al., 2018;Dhawan et al., 2018). Very recently, we reported the transcriptomic profile of crossbred bull spermatozoa that suggested possible implications of transcriptomic variations on semen quality and fertility . Spermatozoa contain both mature and immature RNAs that are represented as a series of exonic, intronic, and intergenic sequences when mapped back to the genome (Selvaraju et al., 2018). A wide variety of coding and non-coding RNA molecules are having important roles in zygote development (Lalancette et al., 2008;Das et al., 2013), post-fertilization events (Selvaraju et al., 2018), early embryonic development , epigenetic trans generation inherence (Jodar et al., 2013;Bohacek and Mansuy, 2015;Rando, 2016), and placental development (Selvaraju et al., 2018). Although earlier studies reported that the sire conception rates (CRs) were associated with various levels of mRNAs in Holstein bulls (Card et al., 2017;Dhawan et al., 2018), significant variations in the sperm mRNA transcripts between different breeds (Selvaraju et al., 2017) limit the possibility of having a universal sperm transcript panel for fertility prediction across the breeds, which highlight the need for identification of breed-specific fertility-associated transcripts. Nevertheless, information on fertility-associated sperm transcripts is lacking in crossbred bulls. With this background, in this study, a comprehensive analysis of global differential gene expression between spermatozoa from high-and low-fertile crossbred bulls was carried out using high-throughput RNA sequencing technique with the aim to identify the most relevant molecules for fertility prediction and to understand the reason/pathway behind crossbred bull infertility.

MATERIALS AND METHODS
The present study was conducted at Theriogenology Laboratory, ICAR-National Dairy Research Institute, Karnal and Southern Regional Station of ICAR-NDRI, Adugodi, Bengaluru. The study protocol was duly approved by the Institute Animal Ethics Committee (CPCSEA/IAEC/LA/SRS-ICAR-NDRI-2019/No. 09) and performed in accordance with relevant guidelines and regulations.

Sperm Sample
Cryopreserved semen straws of Holstein Friesian crossbred bulls (n = 12) of known fertility status were procured from Kerala Livestock Development Board, Mattupetty, Kerala, India. The CRs of the experimental bulls are shown in the Supplementary Figure 1. CR was calculated based on the number of animals conceived out of total number of animals inseminated (up to three inseminations). The effect of non-genetic factors on CR was studied using least squares analysis of variance for unequal and non-orthogonal data as described previously by Singh et al. (2016). The adjusted CR was used for the calculation of bull fertility. Bulls with CR Mean + 1 standard deviation were considered as high fertile and Mean -1 standard deviation was considered as low fertile. Spermatozoa from high-(HF) and lowfertile (LF) bulls (n = 4 samples) were individually subjected to high-throughput RNA-seq analysis as detailed below.

Spermatozoa RNA Isolation and cDNA Synthesis
For RNA isolation, pure sperm fraction obtained using 90-45% discontinuous Percoll gradient was used. This procedure eliminated contaminating substances like epithelial cells and semen extender. Total RNA was isolated from frozen sperm using TRIzol (Ambion, Thermo Fisher Scientific, United States) as described by Parthipan et al. (2015) with minor modifications. RNA quantification was done using NanoDrop (ND-1000, Thermo Fisher Scientific, United States). RNA samples with 260/280 ratio of 1.7-2.0 were selected for cDNA synthesis (reverse transcription), which was done using a combination of oligo (dT) and random hexamers with an initial concentration of 50-100 ng of total RNA from each crossbred bull spermatozoa sample using the RevertAid First Strand cDNA Synthesis Kit (Thermo Fisher Scientific, United States, Catalog number K1622) based on the manufacturer's instructions of 20 µl final volume. The cDNA samples were stored at −20 • C until further processing.

Transcriptome Library Preparation
Total RNA (1 µg) was used to enrich mRNA using NEB Magnetic mRNA Isolation Kit (Illumina, United States). The transcriptome library was prepared using NEB ultra II RNA library prep kit (Illumina, United States) and sequenced using Illumina NextSeq 500 (Illumina, United States) paired-end technology. The enriched mRNA was fragmented (approximately 200 bp) using fragmentation buffer. Random hexamer primers were then added and hybridized to complementary RNA sequences. These short fragments were used as templates to synthesize the first strand of cDNA using reverse transcriptase and dNTP. The DNA-RNA hybrids produced during first strand cDNA synthesis are converted into full-length double-stranded cDNAs using RNase H and Escherichia coli DNA polymerase I. The second strand of cDNA was synthesized using second strand enzyme mix and buffer. The double-stranded cDNA fragments were purified using 1.8× AMPure beads. The purified doublestranded cDNA was end repaired to ensure that each molecule was free of overhangs and has 5 phosphates and 3 hydroxyls before the adaptor ligation. The adaptor ligated DNA was then purified using AMPure beads and enriched with specific primers, compatible for sequencing on to the Illumina platforms. The final enriched library was purified and quantified by Qubit R Fluorometer, and the size was analyzed by a bio-analyzer.

Sperm RNA Sequencing and Data Analysis
The cDNA of high-fertile (n = 2) and low-fertile (n = 2) bulls were sequenced using Illumina NextSeq 500 sequencing system (Sandor R Lifesciences Pvt. Ltd. Banjara Hills, Hyderabad, India) to generate paired-end 76 bp reads. The sequence analysis was done using the online server tool Galaxy 1 (Shannon et al., 2003). Raw data were generated from the four samples, read quality was checked using FastQC (Galaxy version 0.72) program, and the reads were then processed with Cutadapt tool (Galaxy Version 1.18) (Martin, 2011). Processing includes removal of adapter (AGATCGGAAGA) sequence, length trimming (>15 bp) and quality trimming (30 phred score). With the use of HISAT2 (Galaxy Version 2.1.0+galaxy4) (Kim et al., 2015), all the four sample processed reads were aligned to the bovine genome (Bos taurus UMD 3.1.94/Btau8), and the samples were sorted with aligned sequences using Samtools (Galaxy Version 2.0.2) (Li, 2011). The mapped and properly paired sequence to the reference genome was calculated based on tabular descriptive statistics dataset tool Flagstat (Galaxy Version 2.0.1) (Li, 2011). 1 https://usegalaxy.org/ With the use of the tool Cufflink (Galaxy Version 2.2.1.2) (Trapnell et al., 2010), the presence of individual transcripts and their expression levels were expressed as fragments per kilobase of exon per million fragments mapped (FPKM). The data of cufflink assemblies were then merged using Cuffmerge (Galaxy Version 2.2.1.2). Between high-and low-fertile groups, significant changes in transcript expression, splicing, and promoter were studied using Cuffdiff (Galaxy Version 2.2.1.5).
The transcripts expressed between the high-and low-fertile groups were categorized as differentially expressed transcripts based on log2 fold change and FPKM value. Transcripts classified based on log2 fold change > +1 (upregulated in low fertile), >−1 (downregulated in low fertile), and between > ±1 and < ± 1 (neutral). The transcripts with FPKM present only in highfertile group were considered as unique to high-fertile bulls, while transcripts with FPKM present only in low-fertile group were classified as unique to low-fertile bulls. All the raw data were uploaded and available in the National Center for Biotechnology Information (NCBI) Sequence Read Archive (SRA) database under PRJNA516089 2 . The total number of sperm transcripts expressed between high-and low-fertile populations was plotted using online tool Venny (Version 2.1.0). Selected differentially expressed spermatozoa transcripts were generated as heat map using Clustvis (r programming) (Metsalu and Vilo, 2015).

Gene Ontology and Functional Pathway Analysis
Gene ontology (GO) classification of sperm transcripts was done using Uniport 3 and the Database for Annotation, Visualization, and Integrated Discovery (DAVID) Bioinformatics Resources (v6.8) 4 and categorized as molecular function (MF), biological process (BP), cellular component (CC), and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway. The top 10 BPs, CCs, and MFs were plotted as a donut pie chart using Highcharts 5 . Pathway enrichment was carried out using Custer Profiler and enriched KEGG functions. Interaction of genes and detailed network analysis of combined GO categories and pathway analysis was performed using ClueGo (Version 2.5.4) and Cluepedia (Version 1.5.4) plugins in the open source Cytoscape (Version 3.7.1) 6 platform (Locatelli et al., 2019). A network of interactions between related genes were obtained in the form of different layouts. Every analysis was performed with B. taurus as background.

−2.06
Capacitation and sire fertility Kasimanickam et al., 2013;Kadivar et al., 2016 10.   (NGS). Primer designing was done using web-based software PRIMER-3 across an exon-exon junctions in order to eliminate contaminating genomic DNA amplification. The annealing temperatures of primers for the selected genes were optimized using PCR (Prima-96plus, Himedia). The cDNAs prepared from different bull semen samples were subjected to RT-qPCR experiment using Insta Q96 Plus Real Time Machine PCR system (HiMedia, India) in a 20 µl reaction comprising 2 µl of cDNA, 0.5 µl of (10 pmol/µl) forward and reverse primers, and 10 µl of Maxima SYBR Green/ROX qPCR master mix 2×. The thermal cycling conditions consisted of initial denaturation at 95 • C for 10 min, followed by 40 cycles of 95 • C for 15 s, 60 • C for 30 s, and 72 • C for 30 s. The primer sequence of selected genes for validation using RT-qPCR is given in Supplementary  Table 1. All the reactions were performed in duplicates, and qPCR amplification of selected genes with their desired product sizes was confirmed by 2% agarose gel electrophoresis. Relative gene expression levels were determined using 2 − Ct method (Schmittgen and Livak, 2008), where Ct = Ct target − Ct internal reference. The fold change was calculated as the mean expression in low-fertile bulls in comparison with high-fertile bulls. GAPDH served as the internal reference gene. Initially, we evaluated the expression stability of 10 commonly used housekeeping genes (Supplementary Table 2) and found that GAPDH was the most stable internal control gene (ICG) based on the analysis using Genorm, Normfinder, Delta Ct, and the comprehensive ranking by RefFinder (Supplementary  Figure 2). The differences in the relative expression of genes between high-and low-fertile groups were assessed for statistical significance using Mann-Whitney U-test (SPSS, 22.0, IBM, United States). The difference was considered as significant when p < 0.05.

Sperm Transcriptome Profile of Highand Low-Fertile Bulls
With the use of Illumina Next Seq-500 RNA sequencing, the total raw data generated for the study were 115 million reads. Processing of raw data resulted in 97 million reads, which was mapped against the Bos taurus genome. The crossbred bull spermatozoa contained transcripts for 13,563 genes; highand low-fertile bull spermatozoa contained transcripts for 8,109 and 11,470 genes, respectively. Of the total transcripts detected, 6,016 transcripts were common between the high-and lowfertile populations, 2,093 transcripts were unique to the highfertile population, and 5,454 transcripts were unique to the low-fertile population (Figure 1). Since there are differences in the total number of reads in each sample and a normalized count is required to compare the samples, we normalized the data using Cufflink tool, which generates normalized read counts in FPKM. By this process, the read count matrix was transformed to allow meaningful comparisons of counts across samples. The existence of specific transcripts and the abundance of mRNA transcripts were represented as FPKM using the Cufflink tool. The processed reads were mapped to the B. taurus  reference genome. Read mapping to the exonic region was used to measure the relative abundances of mRNA transcripts. Therefore, the normalized expression in FPKM was determined based on the number of reads mapped. After normalization of data, a total of 776 transcripts were detected, in which 524 transcripts were common to both high-and low-fertile bulls, 84 sperm transcripts were unique to high-fertile bulls, and 168 transcripts were unique to low-fertile bulls (Figure 2), while 176 transcripts were upregulated (fold change > 1) and 209 were downregulated (fold change < 1) in low-fertile bulls. The top 20 differentially expressed sperm transcripts were plotted using heat map (Figure 3), in which ribosomal proteins (RPS8 and RP14) and TPT1 were highly upregulated, whereas ZNF706 was highly downregulated in low-fertile bulls. The top 10 transcripts unique to the high-and low-fertile populations after excluding the non-coding RNAs are listed in Tables 2, 3. A complex pool of coding and non-coding RNAs was also observed in the sperm transcriptome. Among the total transcripts (13,563 transcripts), coding RNAs, non-coding RNAs, pseudogenes, and processed pseudogene were 13,002, 88, 375, and 98, respectively. Among the non-coding RNAs, miscellaneous RNA (misc_RNA), small nuclear RNA (snRNA), small nucleolar RNA (snoRNA), and ribosomal RNA (rRNA) were 48, 21, 12, and 7, respectively. After normalization, a total of 776 transcripts were detected, in which coding RNA (known/categorized), coding RNA (unknown/uncategorized), and non-coding RNA were 585, 148, and 43, respectively. Among the non-coding RNAs, misc_RNA, snRNA, snoRNA, and rRNA were 19, 12, 6, and 6, respectively. A total of 61 protein coding ribosomal mRNAs were observed among the total transcripts after normalization.  Table 3). GO analysis of 209 (178 functionally annotated) sperm transcripts downregulated in the low-fertile population revealed their involvement in 8 MFs, 17 BPs, and 12 CCs (the top 10 in each GO category are shown in Figure 5) and 8 KEGG pathways (Supplementary Table 4). The 524 (436 functionally annotated) sperm transcripts common to the high-and low-fertile populations revealed their involvement in 29 MF, 71 BP, and 39 CC categories and 15 KEGG pathways. The 84 sperm transcripts (50 functionally annotated) unique to the high-fertile population revealed their

Real-Time Expression Analysis of Selected Genes
The fold change in the expression of the 12 genes (TPT1, PFN1, ZNF706, CRISP2, MDB4, TNP2, ADIPOR1, TNP1, IQCF1, RACK1, TMSB10, and TSSK6) between high-and low-fertile bulls is given in Figure 8. Results of RT-qPCR expression analysis revealed that all the genes followed the same trend of expression as observed in NGS, except for RACK1 gene. To understand the variability in the expression of the 12 genes, the geometric mean with 95% confidence interval was calculated for high-and low-fertile bulls, and the results are shown in Supplementary Figure 7. Real-time expression analysis indicated that four genes (ZNF706, CRISP2, TNP2, and TNP1) were significantly (p < 0.05) downregulated in low-fertile bull spermatozoa. Correlation analysis revealed that expression levels of ZNF706, CRISP2, TNP2, and TNP1 genes were positively and significantly (p < 0.05) related to bull fertility (Table 4).

DISCUSSION
The uncertain predictability and the inability of the currently available semen evaluation techniques to maintain constant correlation with bull fertility made us to look for a sperm transcript-based alternative technique for crossbred bull fertility prediction. In that way, earlier studies on humans and livestock ascertained the role of sperm RNA in spermatogenesis, sperm function, and early embryonic and extra embryonic development (Card et al., 2013;Selvaraju et al., 2018). Therefore, we carried out comparative transcriptomic profiling of high-and lowfertile crossbred bull spermatozoa using high-throughput RNA  sequencing technique to decipher the molecular soothsayers for bull fertility and to understand the reason/pathway behind crossbred bull infertility.
The present study revealed the transcripts for 13,814 genes in crossbred bull sperm, which is similar to the earlier studies by and Selvaraju et al. (2017) and Raval et al. (2019), but Card et al. (2013) reported only 6,166 transcripts in bovine sperm. These variations might be due to season of sample collection (Godia et al., 2019), state of spermatozoa (fresh or frozen) Card et al., 2017;Singh et al., 2019), method of RNA isolation (Parthipan et al., 2015;Fraser et al., 2017), integrity of sperm RNA, RNA sequencing instrument (Selvaraju et al., 2017), and library preparation methods (Mao et al., 2014). A high proportion of protein coding ribosomal RNAs was detected as unique and upregulated in low-fertile sperm than high-fertile sperm. Ribosomal RNAs were reportedly packed during the process of spermatogenesis (Garrido et al., 2009;Montjean et al., 2012) and essential for the protein synthesis and sperm motility (Bansal et al., 2015), but their implication in crossbred infertility needs to be studied further.
The highly downregulated sperm transcripts in low-fertile crossbred bull were ZNF706, PICK1, LUZP1, ANKRD9, RUNDC3A, LYRM4, FAM71F1, and EPOP. The zinc finger protein (ZNF706) may be related to zinc that is essential for male fertility and to zinc-containing metalloenzymes that are associated with sperm functions (Kerns et al., 2018). Protein interacting with PRKCA 1 (PICK1) is necessary for cytoskeletal organization (Selvaraju et al., 2017) and acrosome formation (Chen et al., 2016), and its deletion leads to globozoospermia (Liu et al., 2010). Its downregulation in our study is in accordance with Singh et al. (2019). Ankyrin repeat domain 9 (ANKRD9) is involved in lipid metabolism (Wang et al., 2009) and downregulated in asthenozoospermic humans (Jodar et al., 2011). RUN domain-containing 3A (RUNDC34) and LYR motifcontaining 4 (LYRM4) are involved in protein biosynthesis. A family with sequence similarity 71 member F1 (FAM71F1) was detected in Leydig cells and downregulated in azoospermic males (Malcher et al., 2013). Elongin BC and polycomb repressive complex 2-associated protein (EPOP) is a novel gene that has not been reported in bull spermatozoa.
In NGS analysis, we observed that the top sperm transcripts unique to high-fertile bulls were TSSK6, C12H13orf46, FABP3, and IQCF1. However, in RT-qPCR analysis, we found that TSSK6 and IQCF1 were expressed in both high-and low-fertile bull spermatozoa, but the level of expression was lower in low-fertile bulls. With the available knowledge, although it is difficult to explain this paradoxical situation, in NGS data processing, it may be possible that reads of a particular transcript might be discarded during the process of removal of low-quality reads. Testis-specific serine kinase 6 (TSSK6) is involved in protein phosphorylation, sperm chromatin condensation (Selvaraju et al., 2018), sperm motility (Bissonnette et al., 2009;Mondal et al., FIGURE 7 | Oxidative phosphorylation pathway with sperm transcripts downregulated in low-fertile bull sperm. FIGURE 8 | Fold change in expression of selected genes in high-and low-fertile bulls (HF, high-fertile bulls; LF, low-fertile bulls). *p < 0.05; **p < 0.01. 2013), and gamete fusion (Sosnik et al., 2009). The function of chromosome 12 C13orf46 homolog (C12H13orf46) is unknown. Fatty acid-binding protein 3 (FABP3) plays a role in remodeling of member polar lipids in spermatogenesis (Furuhashi and Hotamisligil, 2008), but in contrast to our results, it was highly upregulated in poor motile crossbred semen (Yathish et al., 2017). IQ motif-containing F1 (IQCF1) is localized in the acrosome and involved in tyrosine phosphorylation, motility, capacitation, and acrosome reaction (Fang et al., 2015). The top sperm transcripts unique to low-fertile bulls include ribosomal proteins (RPL37, RPS11, RPS12, RPL13A, RPS3, RPS27, RPL31, RPL30, and RPL32) and thymosin beta 10 (TMSB10). TMSB10 is abundant in the embryos of highfertile sires (Kropp et al., 2017) and is possibly involved in sperm capacitation, fertilization (Selvaraju et al., 2017), and cellular remodeling during trophoblast adhesion (Cammas et al., 2005). Sperm transcripts upregulated and unique to low-fertile bulls were involved in a translation as a BP similar to the reports by Card et al. (2017), Selvaraju et al. (2017), and Raval et al. (2019), and ribosomal pathway. This along with the abovementioned upregulation and abundance of ribosomal proteins in low-fertile crossbred bulls is collectively indicating the possible changes in the translation machinery in low-fertile crossbred bulls. Previous studies described the importance of ribosomal RNAs for the sperm function (Platts et al., 2007;Bansal et al., 2015); however, the exact role of ribosomal RNAs in male fertility is not yet understood.  The important finding of this study is that the sperm transcripts involved in the multicellular organism development (QKI, ODF1, TNP1, PRM2, CFDP1, TNP2, ODF2, SPEM1, and MEA1), spermatogenesis (ODF1, BCL2L11, PRM2, TNP2, ODF2, SPEM1, and MEA1), cell differentiation (QKI, ODF1, TNP2, ODF2, SPEM1, and MEA1), in utero embryonic development (YBX1, UBE2B, BCL2L11, MYH10, and RBBP6), and oxidative phosphorylation pathway (MT-ATP6, ND1, MT-ND2, MT-ND4, ND5, MT-CYB, COX1, MT-CO2, and COX3), which are extremely vital, are downregulated in low-fertile crossbred bull sperm. Oxidative phosphorylation is essential for the sperm to synthesize ATP and produce energy in all mammals (Garrett et al., 2008;Storey, 2008). Alterations in the oxidative phosphorylation process lead to altered sperm function (Terrell et al., 2011) and asthenozoospermia (Rawe et al., 2006;Pelliccione et al., 2011). Therefore, the impaired oxidative phosphorylation could be the predominant contributing factor for crossbred bull infertility.

CONCLUSION
In conclusion, our study established the transcriptomic profile on high-and low-fertile crossbred bull spermatozoa using highthroughput RNA sequencing. These RNAs might have a titanic role on spermatogenesis, post-spermatogenic events, sperm functions, fertilization, and early embryonic development. We identified that spermatozoa from low-fertile bulls had an altered expression of genes involved in oxidative phosphorylation, sperm functions, and embryonic development. The impaired function of oxidative phosphorylation could be the predominant reason for crossbred bull infertility; and significant downregulation of ZNF706, CRISP2, TNP2, and TNP1 genes indicates that they could serve as potential biomarkers for fertility in crossbred bulls.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi.nlm. nih.gov/, PRJNA516089.

ETHICS STATEMENT
The animal study was reviewed and approved by Institute Animal Ethics Committee, SRS of ICAR-NDRI.

AUTHOR CONTRIBUTIONS
MAP: methodology, experiment, writing -original draft, and data curation. AK: conceptualization, project administration, supervision, writing -review, and editing. JE, PN, and AS: methodology and data curation. MS: data curation and bioinformatic analysis. EK: methodology and writing -original draft. TD: formal analysis, writing -review, and editing. All authors contributed to the article and approved the submitted version.

FUNDING
The work was supported by the Bill and Melinda Gates Foundation (grant no. OPP1154401) under the project "Molecular markers for improving reproduction in cattle and buffaloes."