Accurate prediction of bull fertility is essential for genetic improvement, herd sustainability and economic gains. Conventional breeding soundness examinations provide basic insights into semen quality by assessing semen volume, sperm concentration, motility, morphology and general physical health. However, despite passing these evaluations, some bulls still exhibit suboptimal fertility. This discrepancy suggests that conventional semen evaluation fails to capture the underlying molecular and functional determinants of sperm competence, including metabolic activity, mitochondrial function, DNA integrity and regulatory RNA molecules (Fortes et al., 2014; Takeda et al., 2015). Consequently, there is an increasing need to identify reliable molecular and cellular biomarkers that more accurately reflect the fertilizing potential of spermatozoa.
       
In recent years, high-throughput omics technologies, particularly transcriptomics, have become powerful tools in unravelling the molecular basis of fertility (Lowe et al., 2017). Transcriptomic profiling of spermatozoa has revealed that, although mature sperm are transcriptionally inert, they carry a variety of coding and non-coding RNAs that reflect testicular function, spermatogenesis and epididymal maturation which have been shown to reflect fertility status in mammals, including cattle (Raval et al., 2019; Somashekar et al., 2017). These sperm RNAs not only play a role in fertilization and early embryonic development but also serve as molecular fingerprints that can discriminate between high-and low-fertility bulls. Differential gene expression related to energy metabolism, capacitation, acrosome reaction and stress response has been associated with variations in male fertility, indicating that transcriptomic analysis offers a promising approach to identifying fertility-associated biomarkers beyond traditional semen evaluation (Wang et al., 2009). Further transcriptomic data alone do not directly reflect the functional competence of spermatozoa and therefore require complementary functional assays to interpret their biological significance in relation to fertilizing ability.
       
Complementing transcriptomic profiling, flow cytometry has emerged as a critical tool for assessing sperm function at the single-cell level. Unlike traditional microscopy, flow cytometry allows rapid, quantitative and multiparametric evaluation of sperm populations with high reproducibility. Various flow cytometric assays have been standardized to assess key determinants of sperm function, including plasma membrane integrity, acrosomal integrity (AcI), mitochondrial membrane potential (MMP) and DNA fragmentation. These parameters are closely associated with sperm viability, capacitation potential and fertilization capacity in cattle. Despite its advantages, flow cytometric analysis primarily evaluates phenotypic sperm characteristics and does not reveal the underlying molecular mechanisms responsible for functional differences among bulls. Therefore, combining functional sperm assays with molecular profiling approaches such as transcriptomics may provide a more comprehensive understanding of fertility-related variation. Recent evidence also highlights the importance of long non-coding RNAs (lncRNAs), with RNA-seq analyses identifying several testis-expressed lncRNAs associated with sexual maturation and spermatogenesis in bovines, suggesting their potential role as fertility-related molecular biomarkers (Han et al., 2026).
       
While much of the research has focused on transcriptomics or sperm functional assays in isolation, integrated studies combining both approaches remain limited in Bos indicus cattle. Such an integrated approach may help bridge the gap between molecular signatures and sperm functional competence, thereby improving the predictive accuracy of fertility assessment. Therefore, the present study aimed to integrate sperm transcriptomic profiling with flow cytometric evaluation of sperm viability, acrosomal integrity, mitochondrial membrane potential and DNA fragmentation to identify molecular and functional indicators associated with fertility differences between high- and low-fertility Kankrej bulls.

Experimental animals and sampling
 
The study was performed from June 2023 to April 2024 on mature, healthy pedigree-breeding Kankrej bulls maintained at the Livestock Research Station, Kamdhenu University, Sardarkrushinagar, Dantiwada, Dist. Bansakantha, Gujarat, India. All bulls were in good health, under uniform veterinary care and maintained in good sanitary conditions. High (n = 4) and low-fertility (n = 4) Kankrej bulls each with five ejaculates were selected based on retrospective conception rate. Bulls were grouped based on a conception rate: a) mean of + 1 standard deviation as a high fertility (n = 4; 52.81±1.3%) and b) a mean of -1 standard deviation as low fertility (n = 4; 39.46±5.3%). The conception rate was calculated considering minimum 300 AI in the field. The samples for frozen semen doses (FSDs) were collected for the study considering the above data.
 
Evaluation of frozen semen
 
After thawing, frozen semen was evaluated by conventional methods for sperm viability (SV), Hypo-osmotic swelling test (HOST), AcI, sperm mucus penetration test (SMPT) and intracytoplasmic morphological sperm evaluation (IMSE). The IMSE test was evaluated as per method described by Patel et al., (2022). During the evaluation, the semen sample was thawed by placing it in a water bath with the temperature at 37°C for 30 seconds. 
 
Transcriptomic profiling of frozen semen
 
Spermatozoa RNA extraction and cDNA synthesis
 
Total RNA was purified using a hybrid protocol, which included lysis with a 2 ml syringe and the trizol + kit (RNeasy Plus Mini Kit, Qiagen, UK). Based on the results of agarose gel electrophoresis in terms of gDNA (genomic DNA) elimination, the final protocol was followed with the modified protocol of Card et al., (2013). Briefly, 2.5 ml of frozen semen were washed with phosphate buffer saline (PBS), followed by sperm enrichment with 4 ml of DMEM media. The final pellet was resuspended in 1 ml Trizol having 45 μl of β-mercaptoethanol, lysed with a 2 ml syringe needle for 10-15 times on ice, followed by adding 300 μl chloroform and centrifuged at 12,000 g for 10 min at 4°C. The phase separation step was followed twice. The uppermost aqueous phase containing total RNA+ gDNA was passed to gDNA spin column and centrifuged at 12,000 g for 30 sec. In flow through, 5 μl (10 pmol) of DNase I (RNase-Free DNase Set, Qiagen, UK) was added and incubated at 37°C for 15 min followed by incubation at 70°C for 5 min. In that suspension, 650 μl of 100% ethanol was added and allowed 2 min incubation at room temperature. Total RNA was purified as per the manufacturer’s instructions (RNeasy plus micro kit, Qiagen, UK). Total RNA was quantified on Qubit 4.0 Fluorometer (Invitrogen, USA) using RNA high sensitivity kit (Invitrogen, USA). To check the isolated RNA is free from gDNA cell specific primer (PRM1) were used to carry out PCR and ran on 1.8% agarose gel. The complementary DNA (cDNA) preparation from the total RNA was done by high-capacity cDNA reverse transcription kit (Applied biosystems, thermo scientific).
 
Library preparation and sequencing
 
The RNA-Seq libraries were prepared using the NEBNext® single cell/low input RNA library prep kit for Illumina® (Illumina, USA). Libraries were stored at -20°C until further use. Library were assessed for quality and quantity on a qiaxcel advanced bioanalyzer. This library was sequenced using miseq platform. (Illumina, USA).
 
Bioinformatics
 
A total of 15.02 Gbps of data was generated from Illumina MiSeq. The quality (Q30) of the data was assessed using fastp (v0.23.4) with default parameters. QC filtered data was mapped to the refseq genome (GCA_002263795.4) using STAR2.7.11a and read counts were obtained with feature counts (v2.0.6). The gene expression was calculated by linna voom method. Gene ontology analysis was performed with IDEP 2.0 (Integrated differential expression and pathway analysis) software.
 
Sperm apoptotic rate by TUNEL assay
 
Sperm apoptotic rate was evaluated using in situ cell death detection kit, terminal deoxynucleotidyl transferase dUTP nick end labeling assay (TUNEL assay) and fluorescein (Roche diagnostics GmbH, Germany, cat. no. 11-684-795-910) as per manufacturer instructions. Apoptotic rate was calculated by total no of sperm having fragmented DNA visualised by bright green fluorescence out of 100 sperms counted.
 
Flowcytometric analysis of frozen semen
 
The SV, AcI and MMP were recorded using a BD FACS Aria Fusion flow cytometer (Becton dickinson, franklin lakes, United States) equipped with fluorescent probes excited by a 20-mW argon ion laser (488 nm). Forward-scatter (FSC) vs. side-scatter (SSC) plots were used to separate sperm cells from debris. Non-sperm events were excluded from further analysis. Fluorescence detection was set with two photomultiplier tubes: detector FL1 (green: 530/30 nm) and detector FL3 (red: 616/23 nm). A total of 30,000 events per sample for each bull were analyzed at a flow rate of 200 cells/s. Data were acquired and analyzed using BD FACSDiva software.
 
Sperm viability
 
The SV was assessed using BD cell viability kit (catalog no: 349483, Becton dickinson, franklin lakes, United States). Briefly, 2 μl of thiazole orange and 2 μl of propidium iodide (PI) was added to 200 μl of sperm suspension and incubated for 5 min, followed by analysis via flow cytometer. Sample excitation was done by a blue laser (488 nm). The emitted fluorescence of thiazole orange was observed in a FL1 band-pass filter (530/30 nm, green fluorescence) and PI fluorescence was measured using a FL3 band-pass filter (616/23 nm). The resultant populations were divided into live, dead and moribund sub-populations. Three sperm populations were detected on the FL1/FL3 Pseudocolor plot: viable (green), dead (red) and moribund (double-stained) spermatozoa.
 
Acrosomal integrity
 
The AcI was assessed using FITC-PNA (Fluorescein isothiocyanate-labeled Pisum sativum agglutinin) and PI combination as per the protocol described by Singh et al., (2016) with some modification. Briefly, 5 μl of 0.025 mg/ml FITC-PNA was added to approximately 2 million spermatozoa and incubated at 37°C for 10 min, then 2 μl of 1 mg/ml PI was added and incubated for 2 min in dark and assessed by flow cytometer. The samples were excited by a blue laser (488 nm). The emitted fluorescence by FITC was observed in a FL1 bandpass filter (530/30 nm) and PI fluorescence was measured using a FL3 band-pass filter (616/23 nm). The resultant subpopulations were divided into four quarters i.e. live acrosome intact sperm, live acrosome reacted sperm, dead acrosome intact sperm and dead acrosome reacted sperm.
 
Mitochondrial membrane potential
 
The MMP was assessed using JC-1 (5,5,6,6′-tetrachloro-1,1′,3,3′ tetraethylbenzimi-dazoylcarbocyanine iodide) stain as per the procedure described by Nag et al., (2021) with modifications. Briefly, 2 μL of 3 mM JC-1 was added to approximately 2 million spermatozoa and incubated at 37°C for 30 minutes and analysed in a flow cytometer. The fluorescence was measured in FL1 (530/30nm) and FL2 (585/42 nm) channels following excitation with a blue laser. The subpopulation of high MMP and low MMP w-as analysed.
 
Sample size calculation and Statistical analysis
 
Considering bulls as the experimental unit (n = 4 per group), the study design has 80% power to detect large effect sizes (Cohen’s d ≥2.0) at α = 0.05. If ejaculates are treated as replicates (n = 20 per group), the design provides 80% power to detect medium-to-large effect sizes (Cohen’s d ≥0.7). Given the exploratory nature and high-throughput molecular analyses, the current sample size is adequate for identifying major fertility-associated molecular signatures.
       
Descriptive statistics were performed for the data obtained for various semen parameters i.e. microscopic (SV, AcI, HOST test, SMPT and IMSE), flowcytometric (SV, AcI and MMP) and sperm apoptotic rate using semen of high and low fertility (n = 4 each, 5 ejaculates) Kankrej bulls. Semen parameters were estimated considering avoiding pseudoreplication, bull was treated as the experimental unit, with repeated ejaculates analysed using linear mixed models. Fertility group (high vs. low) was specified as a fixed effect and bull identity as a random effect. Data were checked for normality (Shapiro-Wilk test) and homogeneity of variance (Levene’s test). Variables expressed as percentages were arcsine square root transformed before analysis. Model residuals were examined to confirm assumptions. Estimated marginal means ± standard errors (SE) were reported and pairwise comparisons were adjusted. Differences were considered statistically significant at p<0.05. All statistical analyses of collected data were performed using SPSS 26 (IBM, Bangaluru, India).

Recent molecular tools such as flow cytometric phenotyping, apoptotic profiling and transcriptomic analyses enable linking molecular signatures with functional semen traits. While such approaches are well explored in Bos taurus, corresponding information in Bos indicus bulls remains limited (Kumaresan et al., 2017; Talluri et al., 2022). To our knowledge, this is the first study integrating conventional semen evaluation, flow cytometry, TUNEL assay and sperm transcriptomics in indigenous Kankrej cattle and correlating these molecular traits with field fertility outcomes.
 
Sperm functional attributes in bulls with different fertility ratings
 
The average seminal attributes and their statistical significance between high-and low-fertility bulls are presented in Table 1. The semen attributes SV, AcI, SMPD and HOST response were significantly higher (P<0.01) in high-fertility bulls (62.45±0.46%, 69.55±0.52%, 29.15±0.47 mm and 59.35±0.58%, respectively) compared to low-fertility bulls (51.45±0.40%, 55.55± 0.48%, 15.10±0.41 mm and 46.25±0.57%). According to the grading system of Vanderzwalmen et al., (2008), semen from high-fertility bulls fell under the good grade category (20-30 mm), whereas low-fertility bull semen was graded as medium (12-20 mm).


       
IMSE analysis revealed a significantly higher proportion of Grade I spermatozoa in high-fertility bulls (82.15±0.12%) compared to low-fertility bulls (77.45±0.13%). Conversely, Grade IV abnormalities were significantly higher in low-fertility bulls (13.05±0.12%) than in high-fertility bulls (8.55±0.10%), while Grades II and III did not differ significantly. These observations corroborate earlier findings in Kankrej bulls (Patel et al., 2022).
       
In contrast to the present study, Dogan et al., (2013) did not observe significant differences in sperm viability between fertility groups. However, Yániz et al. (2021) reported significantly higher acrosome integrity in high-fertility bulls, consistent with our findings. Variability in post-thaw sperm kinematics and fertility relationships may arise from differences in CASA settings, thawing protocols, sperm concentration and extender composition (Karunakaran and Devanathan, 2017). Similar semen characteristics were reported earlier in frozen semen of Kankrej bulls irrespective of fertility status, aligning closely with the present low-fertility group (Patel et al., 2022).
 
Transcriptomic profiling of frozen semen
 
RNA quantification and library quality control
 
PCR analysis confirmed the absence of genomic DNA contamination and library concentrations ranged from 3.76 to 12.90 ng/μl (mean: 8.33±0.66 ng/μl), with an average fragment size of 263 bp. The relatively low RNA yield observed is consistent with earlier reports showing sperm RNA content ranging from 2-100 fg per cell across species (Sahoo et al., 2021). RNA yield is influenced by species, semen type, extraction protocol, sperm concentration and enrichment media (Tiwari et al., 2023).

Sequencing quality and alignment
 
After quality filtering, 97.66% of reads were retained, yielding an average of 0.77 million reads per sample. Reads were aligned to the Bos taurus reference genome (GCA_ 002263795.4), with an average mapping rate of 78.05%. A total of 9,113 transcripts were detected, comparable to earlier reports in bovine sperm (Selvaraju et al., 2017; Raval et al., 2019), though higher than those reported by Card et al., (2013). Variations among studies may be attributed to seasonal effects, semen status (fresh vs frozen), RNA integrity, sequencing platform and library preparation methods (Mao et al., 2013).
 
Differential gene expression
 
Gene ontology analysis identified 36 differentially expressed transcripts between high-and low-fertility bulls. Among these, 26 genes were upregulated (>1 log2 fold change) and 10 were downregulated (<-1 log2 fold change) in high-fertility bulls (Fig 1; Table 2). Compared to the present findings, Prakash et al., (2021) reported a higher number of differentially expressed transcripts, while Karuthadurai et al., (2022) observed extensive transcriptomic alterations in poor-quality semen.



​
 
Functional enrichment and pathway analysis
 
Pathway enrichment analysis of upregulated genes in high-fertility bulls revealed significant involvement in metabolic pathways, lysosome function, carbon metabolism, biosynthesis of nucleotide sugars and amino acid biosynthesis. Key upregulated transcripts included CRYZ, PLD3, TKT, PGAP3, TBCC and PSMC1. Phospholipase D (PLD) is known to regulate sperm capacitation and hyperactivation (Brener et al., 2003) and PLD1 localization in the acrosomal region suggests a role in acrosome reaction (Garbi et al., 2000).
       
Post-translational modifications such as ubiquitination and phosphorylation play critical roles in sperm function (Hunter 2007). The PSMC1 gene, involved in protein ubiquitination, has been linked to DNA integrity and sperm quality (Zhang et al., 2021). Its upregulation in high-fertility bulls may explain the superior seminal attributes and reduced DNA damage observed in this group. Downregulated transcripts in high-fertility bulls included AGBL4, TXNIP and FNIP2. Notably, AGBL4 has been associated with teratozoospermia in humans (Han et al., 2021), suggesting its potential as a negative fertility marker, although further validation in bovines is required.
       
Spermatozoal mRNA profiles reflect transcriptional events during spermatogenesis and may indicate inefficiencies affecting sperm quality (Platts et al., 2007). The enrichment of metabolic pathways observed here aligns with evidence that mature sperm metabolize exogenous substrates to regulate processes such as motility, capacitation and acrosome reaction (Odet et al., 2013). Supporting this, Ozbek et al., (2021) demonstrated differential abundance of metabolites like GABA and benzoic acid between high-and low-fertility bulls.
       
Genes upregulated in high-fertility bulls were associated with mitochondrial and membrane components, including translocase of the outer mitochondrial membrane complexes, which facilitate protein transport and support mitochondrial function. This may explain the higher MMP observed in high-fertility bulls. Ribosomal RNAs, known to persist in sperm, may also contribute to mitochondrial protein synthesis during capacitation (Zhao et al., 2009), although their precise roles remain to be elucidated.
 
Flow cytometric evaluation of frozen semen
 
High-fertility bulls exhibited significantly higher SV, AcI and MMP compared to low-fertility bulls (P<0.001; Fig 2) which is in accordance with previous reports (Kumaresan et al., 2017; Saraf et al., 2021; Talluri et al., 2022). Lower sperm viability in low-fertility bulls likely reflects compromised membrane integrity, permitting propidium iodide uptake. Similar associations between sperm viability and bull fertility have been widely reported (Singh et al., 2016; Turri et al., 2021).


       
The proportion of live sperm with intact acrosomes was significantly higher in high-fertility bulls, supporting their enhanced fertilizing capacity. Comparable findings were reported in holstein-friesian bulls by Talluri et al., (2022) and others (Garbi et al., 2000; Yániz et al., 2021).
 
Sperm apoptotic rate by TUNEL assay
 
The TUNEL assay revealed a significantly lower apoptotic index in high-fertility bulls (4.90±0.33%) compared to low-fertility bulls (7.60±0.30%; P<0.001; Fig 3). These results agree with earlier studies demonstrating a negative correlation between sperm DNA fragmentation and fertility (Dogan et al., 2013; Kumaresan et al., 2017). A recent meta-analysis (Abah et al., 2025) further confirmed higher sperm DNA fragmentation in low-fertility bulls, emphasizing the value of incorporating DNA integrity assays into breeding soundness evaluations. However, standardized thresholds for sperm DNA fragmentation remain to be established.

Overall, high-fertility Kankrej bulls demonstrated superior sperm functional attributes, including higher SV, AcI, HOST response, SMPD and lower nuclear abnormalities. Transcriptomic profiling identified 36 differentially expressed genes, with enrichment in metabolic and mitochondrial pathways and key upregulated genes (PLD3, PSMC1 and TKT) associated with capacitation, motility and DNA integrity. Flow cytometry confirmed higher MMP and reduced apoptosis in high-fertility bulls. Collectively, these findings establish that integrating flow cytometric evaluation, TUNEL-based DNA integrity assessment and sperm transcriptomic profiling provides a robust framework for predicting bull fertility and advancing genetic improvement programs in indigenous cattle breeds such as Kankrej.

We express our sincere gratitude to Livestock Research Station, Kamdhenu University, Datiwada, Banaskantha, Gujarat, India for providing semen for evaluation. We are also grateful to Gujarat Biotechnology Research Centre, Department of Science and Technology, Gujarat, India for providing environment of experiment specially Flowcytometry and Miseq sequencer.
 
Compliance with ethical standards
 
Informed consent
 
All animal procedures for experiments were approved by the Committee of Experimental Animal Care and Handling Techniques were approved by the University of Animal Care Committee.
 
Availability of data and materials
 
The RNA-seq dataset generated in this study, including raw FASTQ files and processed count matrices, has been deposited in the NCBI Sequence Read Archive (SRA) under BioProject accession number PRJNA 1331512.
       
Flow cytometry (FCS) files: Representative raw cytometry files and gating strategies have been deposited in Zenodo (https://doi.org/10.5281/zenodo.17197655).

The authors declare that there are no conflicts of interest.