Single nucleotide polymorphisms cumulating to genetic variation for fertility in crossbred (Bos taurus × Bos indicus) bull spermatozoa

Abstract Spermatozoa from high-fertile (HF) and low-fertile (LF) breeding bulls were subjected to high-throughput next-generation sequencing to identify important Single nucleotide polymorphisms (SNPs) and novel variants associated with fertility. A total of 77,038 genome-wide SNPs were identified, among which, 10,788 were novel variants. A total of 42,290 and 34,748 variants were recorded with 6115 and 4673 novel variants in in HF and LF bulls, respectively. Higher number of SNPs were identified in HF compared to LF bulls. GO analysis of filtered genes with significant variations in HF bulls indicated their involvement in oxidative phosphorylation and metabolic pathways. GO analysis of filtered genes with significant variation in LF bulls revealed their involvement in Ca2++ ion binding, structural constituent of ribosome, and biological processes like translation and ribosomal small subunit assembly. The study identified SNPs in candidate genes including TPT1, BOLA-DRA, CD74, RPS17, RPS28, RPS29, RPL14, RPL13, and RPS27A, which are linked to sperm functionality, survival, oxidative stress, and bull fertility. The identified SNPs could be used in selection of bulls for high fertility and the variation in these genes could be established as an explanation for the fertility differences in bulls upon validation in large number of bulls.


Introduction
The decline in the reproductive efficiency of dairy cattle has become a challenging problem worldwide, and it may cause colossal economic loss to the dairy industry. 1,2Although both dam and sire contribute to reproductive success, much of the research on dairy cattle has been directed toward female fertility.In contrast, bull fertility has received much less consideration, and the potential contributions of the bull have been largely ignored.Recent research envisages that a significant percentage of reproductive failure in dairy cattle is attributed to infertility or subfertility in bulls. 3his emphasizes the fact that bull fertility cannot be disregarded in breeding schemes aimed at improving the reproductive performance of dairy cattle. 4Thus, the economic impacts of bull fertility greatly influence the AI industry. 5A recent review emphasized the need for developing markers for bull fertility, especially in case of crossbred bulls. 6Crossbreeding of indigenous cattle (Bos indicus) with improved (Bos taurus) breeds gained momentum and economic relevance in several countries to increase milk production.While production performance of the crossbred offspring is high due to hybrid vigor, they suffer from a high incidence of reproductive problems, especially the infertility problems in crossbred males as compared to purebred bulls. 6urrently, breeding bulls are evaluated based on a breeding soundness exam (BSE) that is composed of preliminary inspection of basic semen quality attributes such as sperm morphology, sperm concentration and sperm motility 7 combined with a phenotypic evaluation of bulls' features.Despite great efforts to evaluate bulls using BSE, bull fertility is deemed suboptimal under field conditions, with a conception rate varying from 20 to 40%. 6][9][10] The source or etiology for the potential differences in bull fertility may be genetic or epigenetic. 11An increasing research evidences suggests that genetic factors could explain part to the differences in fertility among the sires and evaluation of service sire/bull needs much attention. 12Genetic variants with moderate to significant effects on fertility can be identified by sequencing the genomes of high-fertile (HF) and LF and subfertile or infertile sires. 11n the recent past, there has been much acceleration in the application of high-throughput sequencing technologies like next-generation sequencing (NGS) platforms for the rapid screening of thousands of molecules that are related to bull fertility.The main advantages of using these sequencing platforms is that they are simple, scalable, more accurate and give highest confirmed output; their sensitivity of screening is very high and they take little time for screening large number of samples.The only disadvantage could be the cost involved and the requirement of specialized software for analyzing the data.4][15][16][17][18] Identifying genomic regions, and preferably individual genes responsible for genetic variation in bull fertility will enhance the understanding of biological pathways involving this trait and may point to opportunities for improving sire fertility via selective breeding.Genetic estimates of fertility can be improved by studying genome-wide associations and by employing single nucleotide polymorphism (SNP) arrays.4][25][26][27] Many bulls are now genotyped for thousands of SNP markers used in genomic evaluation of breeding value. 28ovine whole genome sequence information, 29 combined with estimated effects of many SNP can readily facilitate localization of specific genes or gene regions related to traits of interest.1][32] Detection of candidate genes may provide valuable information for better understanding and identifying causal polymorphisms that can lead to more efficient genetic improvement using SNP-assisted selection.The previous use of SNPs for fertility studies has been limited to a few markers, and their implication in male infertility has not yet been fully elucidated.Hence, in the current study, we aimed at identifying the various genetic variants significantly related to fertility by thorough screening of the SNPs in the sperm transcripts.

Ethical approval statement
The present study was carried out at the Theriogenology Laboratory, Southern Regional Station of ICAR-National Dairy Research Institute, Bengaluru, India.All the experiments and procedures performed according to the guidelines laid down and as approved by Institute Animal Ethical Committee (CPCSEA/ IAEC/LA/SRS-ICAR-NDRI-2019/No.04) in the current study.

Experimental animals and samples
Twelve experimental bulls (Holstein Friesian crossbred with 50-75% exotic inheritance; n ¼ 12) were selected out of 50 bulls based on their known fertility status, and these bulls are routinely used for artificial breeding program.Based on their field conception rates (CR), the bulls were classified either into HF or LF (n ¼ 6 each).CR for each bull was calculated based on the proportion of pregnant animals to the total number of inseminations per bull (up to third insemination).Furthermore, CR was adjusted to nongenetic parameters according to the procedure described previously by Singh et al. 33 Bulls with CR Mean þ 1 standard deviation were considered as HF and Mean À 1 standard deviation was considered as LF.Cryopreserved spermatozoa from the experimental bulls (six HF and six LF) were used for RNA-Seq analysis.

Extraction of RNA and cDNA synthesis
For RNA extraction, semen was purified using Percoll gradient (90-45% discontinuous) centrifugation to eliminate the semen extender and epithelial cell contamination.Total RNA was extracted from frozen sperm using TRIzol (Ambion, Thermo Fisher Scientific, United States) 17 with minor modifications.Quantification of RNA was done by NanoDrop (ND-1000, Thermo Fisher Scientific, United States).RNA samples with good quality (260/280 ratio of 1.85-2.0)were selected for cDNA synthesis (reverse transcription).cDNA synthesis 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 ll final volume.The cDNA samples were stored at À20 C until further processing.Two representative cDNA samples were made from the six HF bulls by pooling equal quantities of sperm cDNA from three HF bulls each.Similarly, two representative cDNA samples were prepared from six LF bulls.All the samples were subjected to Illumina Nextseq-500 sequencing system for analysis.

Library construction and sequencing
The transcriptome library was produced with the NEB Ultra II RNA library prepkit (Illumina, United States) and sequenced with Illumina NextSeq 500 paired-end technology (Illumina, United States).The NEB Magnetic mRNA Isolation Kit was used to enrich mRNA from total RNA (1 g) (Illumina, United States).Fragmentation buffer was used to fragment the enriched mRNA (about 200 bp).Random hexamer primers were then added and hybridized to complementary RNA sequences.The first strand of cDNA was synthesized by using fragments as templates and done by 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. Second strand enzyme mix and buffer were used for obtaining the second strand of cDNA.1.8AMPure beads were used for purifying the double-stranded cDNA fragments.The adaptor-ligated DNA was subsequently purified using AMPure beads and enriched with particular primers that were compatible with Illumina platforms for sequencing.A QubitV R Fluorometer was used to quantify the final enriched library, and a bio-analyzer was used to determine its size.

Data analysis
The sequence data quality was assessed using the FastQC version 0.10.0. 34The raw data generated in fastq format was formatted by trimming (removing low quality based below À q 30), removing the adapters and length filtering (minimum of 15 bp reads are retained) using the tool 'Cutadapt v2.8'. 35urthermore, the data were aligned against the Bos taurus reference genome (ARS-UCD1.2-GCA_002263795.2) using 'HISAT2 v2.2.1' program that maps the processed data to the reference genome resulting .bamfiles.For variant identification analysis, the number of mutational events were quantified using 'bcftool v1.10'.Identified variants were annotated using the VEP (Ensembl Variant Effect Predictor) and SnpEff program.

Variant calling, annotation, and analysis
The aligned samples and the reference genome sequence were used for variant (SNPs) calling.We further processed the .bamfiles for the identification and acquisition of variants by employing multiple bioinformatic tools 'Samtools v1.9' and 'bcftool v1.10' 36 for variant calling and 'Samtools v1.9' programme 'mpileup' function was used for generating a pileup file that provides a single-base resolution of coverage, ie, coverage of mapped reads to the reference.The identified variants were filtered using the 'VCFtools v0.1.8'and significant variants were obtained with a minimum quality of Phred score Q30 and minimum coverage at each base depth of 20 using the 'varFilter' function and filtered vcf file was generated. 37The SNPs discovered were annotated against the reference bovine genome with general feature file (gff) database and functional effect were predicted with the help of VEP and SnpEff 38 tool.VEP reports Sorting intolerant from tolerant (SIFT) and gives scores for the nonsynonymous variants in the dataset.A common threshold is to classify variants with SIFT score <0.05 as potentially deleterious.The SNPs obtained after the variant calling were compared with the cattle SNPs in dbSNP (build 150) of NCBI to identify already reported known SNPs and unreported novel SNPs.The variants identified in the study were mapped to genomic regions harboring the major candidate genes identified through association studies for production traits in dairy cattle. 39

Gene ontology analysis and functional annotation
Gene ontology analysis for the critical genes with differential observation of SNPs related to fertility and for the novel variants and KEGG pathway analysis was carried using DAVID Bio-informatics Resources 6.8 (Laboratory of Human Retrovirology and Immuno-informatics, USA) for important genes involved in biological and molecular processes.KEGG pathway enrichment analysis and visualization for the critical genes were carried out using R package 'clusterProfiler'.Interaction network analysis of combined GO categories and pathway analysis among the SNPs was performed using ClueGo (Version 2.5.4,Integrative Cancer Immunology, Jerome Galon) and Cluepedia (Version 1.5.4)plugins in the open source Cytoscape (Version 3.7.1,National Institute of General Medical Sciences (NIGMS), USA) platform 40 and obtained as different network layouts.

Total variants in crossbred bull spermatozoa
With Illumina, Next Seq-500 RNA sequencing, the total raw data of about 115 million reads were generated.Furthermore, raw data processing resulted in 97 million reads, which were mapped against the Bos taurus reference genome (ARS-UCD1.2-GCA_002263795.2).A total of 77,036 SNPs were identified to be present in the crossbred bull spermatozoa samples with a minimum read depth of 20, and among them, 10,269 were observed to be novel variants.All these variants were found to have a minimum quality of q30 and coverage of 20.A total of 42,290 and 34,748 variants were recorded in HF and LF bulls, respectively, among which 6115 and 4673 were found to be novel in HF and LF groups of bulls, respectively.After filtering, 92 and 64 significant genes (ie, genes having the SNPs with mutational event whose impact is high, sorting intolerant from tolerant (SIFT) score is deleterious and consequence is a missense variant) were observed in HF and LF bulls sample groups.The Venn diagram showing the total and novel variants were represented in Fig. 1A,B, respectively.The total variants identified in each sample were annotated using the VEP database (Ensemble Variant Effect Predictor database).The details of total variants, existing and novel variants observed in HF and LF group bulls were presented in Table 1.A Circos plot was also generated for the global distribution of SNP variants, novel variants along the cytogenetic location (Fig. 2A,B).The genes with variants and genes with novel variants were also identified and presented in a Venn diagram in Fig. 3A,B.
The coding region consequences were enumerated in the HF and LF bulls and are classified as synonymous variants, missense variants, frameshift variants, stop retained, stop loss, and stop gained.The data for the same is presented in Fig. 4A,B.In the identified SNPs, silent mutations were most frequent among the functional class, followed by missense mutations and nonsense mutations (Table 2).Classifying SNPs based on the genomic region, most of the SNPs were found in intron specific regions followed by intergenic regions (Supplementary Figure 1).The missense to silent ratio was also presented in Table 3.The Ts/Tv ratio was found to be 1.89, comprising of 52,340 transitions and 25,613 transversions (Table 3).The number of G/C and T/G substitutions was relatively lower than other types of substitutions (Table 4).

Gene ontology (GO) analysis
GO analysis was carried for the specific genes with novel mutations (1127) in HF groups and specific genes with novel mutations (783) genes in LF group.For the HF group, genes with novel sperm-associated variants revealed their involvement in 27 molecular functions (MFs), 79 biological processes (BPs), and 32 cellular components (CCs (the top 10 in each GO category are shown in Fig. 5 (Supplementary Table 1) and 59 KEGG pathways (Supplementary Table 2).Similarly, LF group genes with novel sperm-associated variants revealed their involvement in 21 MFs, 37  BPs, and 7 CCs (the top 10 in each GO category are shown in Fig. 5 (Supplementary Table 3) and 12 KEGG pathways (Supplementary Table 4).
GO analysis of 92 filtered genes with significant variation in the spermatozoa of HF bulls revealed their involvement in 7 MFs, 3 BPs, and 12 CCs, and 3   KEGG pathways.The genes listed were involved in various important molecular functions like Zn and Mg binding, phospholipid-translocating ATPase activity and biological processes like ATP synthesis coupled proton transport and complement activation, classical pathway.The filtered genes were also found to be involved in various important pathways like oxidative phosphorylation and metabolic pathways.GO analysis of 64 filtered genes with significant variation in the spermatozoa of LF bulls revealed their involvement in 4 MFs, 3 BPs, 8 CCs categories, and 2 KEGG pathways.The genes listed involved various important molecular functions like Ca 2þþ ion binding, structural constituent of the ribosome and biological process like translation, ribosomal small subunit assembly.

Network and pathway analysis
Pathway enrichment of genes in the HF group with novel sperm variants revealed that they were involved in important pathways like cAMP signaling pathway (26 gene counts), Calcium signaling pathway (24 gene counts), and top ten enriched pathways are shown in Fig. 6A.In the LF group, the pathway analysis of genes with novel variants were involved in the ribosome (47 gene counts), regulation of actin cytoskeleton (18 gene counts), Endocytosis (16 gene counts), and Ubiquitin mediated proteolysis (12 gene counts), and the top 10 pathways are shown in Fig. 6B.Network analysis using Cluego was carried for important filtered genes of HF and LF groups.Network analysis of filtered genes in the HF group revealed their involvement involving mitochondrial respiratory chain complex, apical plasma membrane, phospholipid translocating ATPase activity and adenylate cyclase binding (Fig. 7A).And that of LF group network analysis of important filtered genes revealed their role in important processes and pathways like small ribosomal unit, MHC protein complex, response to hydrogen peroxide and cysteine type endopeptidase regulator activity.(Fig. 7B).

Discussion
The ability to identify bulls with greater reproductive performance would significantly improve the efficiency of the cattle industry.Identifying genomic regions, particularly individual genes associated with fertility traits and reproduction might improve the reliability of genomic estimates and are important for enhancing sire fertility via selective breeding.Bulls' fertility-related traits are moderately heritable (0.05-0.22) but significantly influenced by genetics. 41herefore, improving our understanding about the genes and epigenetic modifications that contribute to bull fertility could improve reproduction success in crossbred bulls. 41,42Genome-wide association studies (GWAS) have been effective in applying genetic markers, such as SNP markers, to determine genomic regions associated with economically important phenotypes such as fertility.Hence in the current study, we for the first time made an attempt to elucidate the various existing SNPs and novel variants of SNPs in HF and LF bull spermatozoa.
In the current study, we identified a total of 77,038 variants and upon processing we could find 10,788 variants as novel.We could identify 42,290 and 34,748 SNPs in HF and LF bull spermatozoa respectively.We detected a higher number of SNPs in the HF group bull spermatozoa than the LF group bull spermatozoa.The same trend is also observed in Genes with missense, genes with high impact and genes with SIFT deleterious.We also observed more number of silent mutations and less number of nonsense and missense mutations in the spermatozoa of LF bulls in comparison to HF bull spermatozoa.The  genome sequencing has revealed that there are more variants present in the HF bull spermatozoa genome than in the LF bull spermatozoa genome reflecting the difference in the fertility of these two groups.More number of variants observed in HF bulls compared to LF bulls might be explained based on the fact that the  genes associated with fertility traits are polymorphic and are under continuous selection pressure; therefore may not be genetically stable but could be of more evolutionary significance.VEP analysis of spermatozoa of LF bulls revealed some potential transcripts to be stop gained (RPS27A, NRP2, TNP2, CD74, and ZNF280B), stop lost (LYZL6, MANSC4) and deleterious (RPL13, RPS28, TNP2, JUNB, HMGB4, TTLL5, ZFN280B, LIMD1, OIT3, and ORI3E12).
The genes that are stop gained (abnormal cessation of translation due to premature termination codon) plays a vital role in spermatogenesis, chromatin integrity of the spermatozoa and its blastocyst formation blastocyst development and blastocyst growth, placentation and embryogenesis fertilizing ability.There was a deleterious effect on some of the ribosomal genes.Differential expression of the ribosomal proteins could affect the orchestration of ribosomes in mitochondria, which may further lead to dysfunction of spermatozoa mitochondria, which further could compromise the fertilization process in LF group bulls.RPS27A, TNP2 play a crucial role in the process of spermatogenesis and integrity of the spermatozoa.In RPS27A gene, a mutation was observed at chromosome number 11 at position 37970665 where C > T converted.The impact of this mutation was found to be high, and it has a stop gained effect, which will lead to premature termination of the codon and might cause the protein to be abnormally untranslated.Earlier studies indicated that dysregulated expression of RPS27A in spermatozoa had been linked to lower conception rates. 43ibosomal RNAs are said to be compacted at spermatogenesis 44,45 and are necessary for protein synthesis and sperm motility; 46 however, the role of mutations in these ribosomal coding regions in bull infertility has not yet been investigated.ZNF280B gene is required for the successful binding of spermatozoa to the oocyte.These genes were stop gained, indicating their protein formation is hampered, which may be the reason for compromised fertility in LF group bulls.
The important genes filtered with significant mutation corresponding to the respective groups ie, 92 genes for HF group and 64 for LF group and were further considered for the functional annotation and analysis.The KEGG pathway analysis of these genes revealed to be associated with Ribosome, and antigen processing and presentation are the two pathways.In the ribosomal pathway, RPS17, RPS28, RPS29, RPL14, RPL13, and RPS27A genes had significant mutations.RPL14 gene has two missense mutations in BTA22 (chr22) at location 13295907 bp where C > G and other mutations at position 13295942 bp where C > T. Ribosomal protein L14 (RPL14) was found to be prevalent in the testes and is reported to be associated with low fertility in crossbred bulls. 17,47,48Similarly, RPS28 gene has variation at Chr 7 at position 16968899 bp and T > C, a missense variation with a deleterious effect (SIFT value of 0).RPS17 also has a missense variation at Chr 21 at position 22841433 bp.RPL13 is seen to have missense variation, and the mutation also has a deleterious effect (SIFT value of 0).The RPL13 gene has variation at BTA18 at 14490777 bp where A > G. SNP determined with deleterious effect were reported to change the function of a protein and hence contribute to a genetic disorder. 49O analysis of the RPL13 gene revealed its role in mRNAs' stabilization and translational regulation during spermatogenesis. 50The deleterious effect of SNP in this gene may cause dysregulation to the spermatogenesis process.In the RPS29 gene, in chromosome number 10, we found a missense mutation from G > C at 26779309 bp.Genetic variation in RPS29, RPS27A, and RPL13 were reported to be resulted in a low pregnancy rate in humans. 43GO analysis of RPS29 revealed its molecular function in binding to the Zn ion.Zinc is required for the survival of germ cells and the development of spermatogenesis.It also protects the sperm flagellum from early oxidation, as it will attach to sulfhydryl groups of outer dense fiber protein cysteine groups and hence plays a pivotal role in the reduction of oxidation of sperm membranes. 51All these proteins have a significant role in producing healthy and active spermatozoa.Mutations to these important genes might have contributed to the defects in the process of spermatogenesis in LF group bulls.
Another pathway in which the majority of the mutated genes were involved was antigen processing and presentation pathway, which is neither well studied nor described earlier.CD74, PSME1, BOLA-DRA are the three mutated genes in this pathway.CD74 was reported to play a critical role as cytokines contribute to testicular function and maintain male reproductive health. 52This gene was found to have a mutation which is stop gained at position 61755486 bp in chromosome 7, where C > A transition and two missense mutations at Chr 7 (at positions 61750597 bp where C > T; at 61752414 bp where T > A).GO analysis revealed this gene has beta-amyloid binding and cytokine receptor activity as its molecular function and, positive regulation of cytokine-mediated signaling pathway as a biological process.CD74 was previously reported to have a vital role in controlling the antigen presentation for immune responses in chickens 33 and cyprinid fish, 53 and it is highly upregulated in LF. 17 Proteasome activator subunit 1 (PSME1) was found to have a mutation at chromosome number 10 (bta10) at position 21012919 bp where there is a transition of A > C.This gene was involved in various important biological processes like positive regulation of endopeptidase activity, antigen processing and presentation of exogenous antigen, and regulation of catabolic proteasome protein.PSME1 is also involved in epididymis secretory sperm binding protein and are essential for promoting sperm motility. 54,55BOLA-DRA, a Major histocompatibility complex, class II, DR alpha (BOLA-DRA) gene was observed to have two missense mutations in Chr 23 at 25840067 bp, 25840771 bp and both transitions were found to be G > A. This gene was also reported to have a role in membrane integrity and sperm morphology. 56Another important gene observed to have missense mutation is the tumor protein, translationally controlled 1 (TPT1) gene.TPT gene has a major role in apoptosis, cellular differentiation, and sperm functions. 57This gene has missense variation at chr 12 (bta12) at position 15475714 bp where A is replaced with T. It is reportedly abundant in the sperm of humans 58 and chickens 33 but downregulated in oligozoospermic men. 45etwork analysis of different GO terms and KEGG pathway interaction analysis revealed that there is mutation in the intronic regions of the genes related to the NADH dehydrogenase activity (ND3, ND5), Oxidoreductase activity (ND3, ND5), adenylate cyclase binding (ADRB2, AKAP12) and phospholipid translocating activity (ATP8B1, ATP9A) in the spermatozoa of HF bulls, which is not having a deleterious effect and this mutation is moderate and silence accepted.In the LF bull spermatozoa mutation in important genes were involved in the ribosomal pathway (RPL8, RPS14, RPL18, RPS27A, RPS11, RPS10, RPL19, RPS12, RPS9, RPS7, RPL21, RPS8, RPL23, RPS5, RPS26, RPS25, RPS28, RPS27, RPL27A, and RPL37A), regulation of actin cytoskeleton (APC2, AKAP4, PXN, ARPC4, ACTB, ACTG1, FGF16, and PFN1), spliceosome (HSPA8, PRPF38B, DDX5, TRA2A) and ubiquitin mediated proteolysis (UBE2Q2, UBE2E3, TRIM37, TRIM32, and BIRC3).
Profilin 1 (PFN1) modulates actin and is involved in oocyte maturation, fertilization, embryo development 59 and spermatogenesis. 47Ribosomes are actively involved in protein translation in spermatozoa. 60ibosomal protein S8 (RPS8) is found abundant in spermatozoa of humans, 58 whereas Ribosomal protein L14 (RPL14) is abundant in the testis of Bos taurus and Bos indicus bulls.47,48 ACTB is expressed in spermatozoa and is distributed in the acrosomal and postacrosomal regions of ejaculated spermatozoa where it is involved in membrane changes during acrosome reaction with an important implication on sperm function.61 AKAP4 is one of the major components of the sperm fibrous sheath, a structure known to be involved in sperm motility.The molecular chaperone HSPA8 plays a key role in remodeling the sperm surface during capacitation. Itis established that mice lacking HSPA gene are infertile and spermatozoa that fail to interact with the zona pellucida of the oocyte consistently lack HSPA2 protein expression.HSPA8 was also reported to have a role in protecting spermatozoa in the oviduct. 62 AKPs can bind other protein kinases, protein phosphatases, ion channels, and small GTP binding proteins and thus serve as platforms for the integration of cAMP and other signaling pathways.63,64 AKAP proteins are reported to be key molecules in the biochemical machinery regulating the functions of flagella and cilia.Spermatozoa from mice lacking AKAP4 failed to show progressive motility and the male mice were infertile.65 These genes play a vital role in spermatogenesis, microtubule formation in the spermatozoa, and the process of ubiquitination.These genes were found to have high impact mutation, which hampers the important processes which might be the probable cause for the low fertility in LF bulls.

Conclusion
In conclusion, the findings of the present study showed that RNA-Seq is a potential approach for genomic profiling variants and mutational events in relation to fertility.Using RNA-seq technique, candidate genes for high and low fertility characteristics of interest in bulls were determined and discovered the functional assessment of genes with substantial mutations related to bull fertility.Specifically, SNPs in candidate genes including TPT1, BOLA-DRA, CD74, RPS17, RPS28, RPS29, RPL14, RPL13, and RPS27A, which are linked to sperm functionality, survival, and protection from oxidative stress and fertilizing potential, were related to bull fertility.The variation in these genes with the highly significant genic region might be associated with low fertility in breeding bulls.

Figure 1 .
Figure 1.(A) Total and (B) novel variants observed in HF and LF crossbreed bull spermatozoa.

Figure 2 .
Figure 2. Circos plot was also generated for the global distribution of SNP variants, novel variants along the cytogenetic location.The chromosomes were represented in the outermost, next GC% for the chromosome.(A) Two circle highlights all the variants observed outer being HF (second circle from inside) and next is LF (first circle from inside).(B) Two circle highlights novel variants observed outer being HF (second circle from inside) and next is LF (first circle from inside).

Figure 3 .
Figure 3. Venn diagram of HF and LF.(A) Genes with the novel variants.(B) Important filtered gene with consequence as missense variant, impact is high, and SIFT is deleterious.

Figure 4 .
Figure 4. Summary statistics for coding consequences.(A) HF sample and (B) LF sample.

Figure 5 .
Figure 5. Functional annotation of various gene ontology terms (biological processes, molecular functions and cellular components) of novel variants observed in high-and low-fertile groups of bulls.

Figure 7 .
Figure 7. Interaction network analysis of various biological processes, molecular functions, and pathways expressed by the novel variants.(A) High-fertile bull group.(B) Low-fertile groups.

Table 1 .
Summary statistics of total SNPs identified.

Table 2 .
Summary statistics for important filtered genes with their functional class.

Table 3 .
List of total type of mutations observed in crossbred bull spermatozoa.

Table 4 .
Nucleotide changes in identified SNPs of crossbred bull spermatozoa.