Orchestrating the expression levels of sperm mRNAs reveals CCDC174 as an important determinant of semen quality and bull fertility

ABSTRACT

Bulls with acceptable semen quality vary in actual field fertility and this can be elucidated by studying the expression levels of mRNAs in the sperm. The present study aimed at assessing the variations in the sperm gene expression levels of PRM1, CCDC174, RPL36A, TMCO2, SWI5 and OIT3 in bulls differing in fertility status. Frozen semen samples from Holstein-Friesian bulls were classified into high-fertile (n = 8, average field conception rate = 46.1 ± 0.51, p < 0.001) and sub-fertile (n = 7, average field conception rate = 39.4 ± 0.69) groups. In the post-thaw semen samples, sperm kinematics, structural and functional membrane integrities, mitochondrial membrane potential and chromatin distribution were analyzed. The sperm total RNA was subjected to gene expression studies by Real-Time PCR. Multivariate regression analysis was performed using gene expression levels and conception rates. The sperm functional attributes did not differ significantly between the groups. The relative mRNA levels (fold change) of CCDC174 (6.20), RPL36A (4.66), SWI5 (1.86) and OIT3 (1.53) were higher in high-fertile bulls. Further, the expression level of the CCDC174 gene was significantly (p = 0.02) up-regulated in high-fertile bulls. The fertility prediction multivariate model with genes, CCDC174, RPL36A, TMCO2 and OIT3 had the maximum coefficient of determination (R² = 0.68) with the field conception rate. This model had 93.3% bull fertility prediction accuracy with 100% sensitivity and 87.5% specificity. The study suggests that the expression level of CCDC174 can be used as a potential marker for assessing bull fertility.

KEYWORDS:

• Sperm
• semen quality
• conception rate
• sperm mRNA expression
• bull fertility

Introduction

In the dairy industry, artificial insemination (AI) is widely adopted across the globe, but the average conception rate is only 35–40% in India (20th Livestock census, Government of India 2019), USA (Schefers et al. 2010; Norman et al. 2017) and UK (Dobson et al. 2008). A low reproductive efficiency leads to a decrease in milk yield resulting in considerable economic loss to the farmers, which necessitates intense research efforts for augmenting fertility in cows. The field conception rate with AI is influenced by many factors including genetic, micro-environmental and managemental, nutritional and health status of both male and female animals (Walsh et al. 2011; Crowe et al. 2018). Apart from these, one of the major reasons for the reduced conception rate is the insufficiency of quality dairy bulls for breeding (Binsila et al. 2017; Amann et al. 2018). Presently, the breeding bulls are selected based on the breeding soundness evaluation and sperm attributes such as motility, structural membrane integrity, functional membrane integrity and acrosome integrity. These sperm attributes have been reported to have poor correlation with the conception rate and mostly the results were inconsistent (Rodriguez-Martinez 2003; Selvaraju et al. 2008; Somashekar et al. 2015). The subtle changes in the sperm attributes between the high-fertile and the sub-fertile bulls cannot be distinguished by employing the routine semen analysis (Vasan 2011) calling for probing at the molecular level (Binsila et al. 2017).

The sperm carries biomolecules such as RNAs, proteins and metabolites apart from DNA that are used to assess the success of the spermatogenesis and pregnancy establishment (Jodar et al. 2015; Parthipan et al. 2017). Studying the gene expression pattern in sperm is an essential element for understanding the regulated processes at the cellular and molecular level (Ostermeier et al. 2002; Selvaraju et al. 2017) and can be used as a biomarker for assessing fertility status (Jodar et al. 2015). Hence, by analyzing the differences in the expression levels of sperm genes between bulls of varying  fertility rates, the functional competence of the sperm can be assessed.

A study from our group observed that the expression levels of PRM1, CCDC174, RPL36A, TMCO2, SWI5 and OIT3 genes were abundant in the sperm of Holstein Friesian bulls (Selvaraju et al. 2017 and unpublished observations). Their mRNAs have different functions, for instance, the PRM1 is involved in DNA compaction (Aoki et al. 2006) and SWI5 regulates DNA repair (Akamatsu et al. 2010) during spermatogenesis. The testis-specific sperm binding gene, TMCO2 with spliceosome activity might regulate acrosome biogenesis during spermiogenesis (Kaneko et al. 2020). The coiled-coil domain-containing family of protein, CCDC87 has been suggested to have a role in the initiation of motility and acrosomal reaction (Wang et al. 2018). The ribosomal proteins have been reported to regulate the mitochondrial function of the sperm (Bansal et al. 2015) and translation (Naz 1998; de Mateo et al. 2011). The OIT3 has been reported to regulate uric acid metabolism (Yan et al. 2012) might protect sperm from damaging reactive oxygen species (ROS) (Banihani 2018). So far, the expression patterns of these genes that are involved in different functions have not been studied in relation to bull fertility assessment.

Since fertilization is a multifaceted process and its success is governed by many factors, a single test or an attribute will not suffice to classify the bull to be either high-fertile or sub-fertile. Consequently, the fertility diagnostic test with a single attribute or potential marker may be less efficient because of the sub-optimal sensitivity and specificity. A rising trend is to identify the potential markers of prime importance (Gerszten and Wang 2008; Yang et al. 2018) and combine them to deploy as a diagnostic test to improve the accuracy (Mazzara et al. 2017). Therefore, the present study aims (i) to explore the expression levels of these genes to consider as a potential marker(s) of sperm function and (ii) to develop a multivariate regression model for the identification and selection of fertile bulls for the assisted reproductive technologies.

Results

Sperm functional parameters

Based on the average conception rate, the bulls were classified into high-fertile (n = 8) and sub-fertile (n = 7). There was a significant (p < 0.001) difference in conception rate between these two groups. Even then, the post-thaw semen analysis revealed no significant (p > 0.05) differences in the sperm kinematics (Table 1), structural membrane integrity, acrosome integrity, sperm subpopulation positive for both functional membrane and acrosomal integrities, mitochondrial membrane potential and chromatin distribution between the groups (Supplementary Figure 1).

Table 1. Test of significance computed between the high-fertile (n = 8) and sub-fertile bulls (n = 7) for each of the sperm functional parameters analyzed.

Sperm functional parameters	High-fertile (n = 8)	Sub-fertile (n = 7)	p value
Total motile spermatozoa (%)	74.90 ± 3.61	69.07 ± 4.60	0.33
Progressive motile spermatozoa (%)	51.80 ± 4.45	45.13 ± 4.40	0.31
Fast progressive motile spermatozoa (%)	32.23 ± 4.90	32.02 ± 3.67	0.99
Curvilinear velocity (μm/s)	82.80 ± 5.20	93.59 ± 4.94	0.16
Straight-line velocity (μm/s)	47.09 ± 4.37	57.30 ± 5.13	0.15
Average path velocity (μm/s)	41.76 ± 8.58	57.81 ± 8.85	0.22
Structural membrane integrity (%)	66.71 ± 3.06	63.27 ± 2.30	0.45
Acrosome integrity (%)	79.64 ± 4.48	74.17 ± 4.15	0.34
Functional membrane and acrosome integrities (%)	34.44 ± 3.52	30.60 ± 3.41	0.46
Mitochondrial membrane potential (%)	64.48 ± 3.24	61.32 ± 3.40	0.53
Chromatin distribution (%)	91.21 ± 1.06	91.40 ± 2.27	0.99

The sperm functional parameters did not differ significantly between the groups. The values were expressed as mean ± SEM.

Selection of housekeeping gene

The housekeeping genes were expected to have expression stability across all the bulls studied, which was determined by the pair with minimum ΔCt and less average deviation in ΔCt values. The ΔCt has been computed among the selected genes and the GAPDH and RPL23 had less variation were considered as stably expressed genes. The mean standard deviation of ΔCt for GAPDH, RPL23 and PRM1 were 1.329, 1.326 and 1.66, respectively. Since the PRM1 had maximum variation in the standard deviation of the ΔCt values (Figure 1A), the gene was not considered as a housekeeping gene. Even though GAPDH and RPL23 did not show much difference in the average standard deviation, variability in Ct was very less with RPL23, which was evident from the inter-quartile region (Figure 1B). As revealed by the NormFinder analysis, RPL23 had the lowest stability value of 0.6 and hence was chosen as housekeeping gene for the present study (Supplementary Table 1).

Figure 1. Selection of a housekeeping gene for sperm gene expression study. (A): Comparison of the ΔCt among each pair of the housekeeping genes studied: GAPDH vs RPL23, RPL23 vs PRM1 and GAPDH vs PRM1. The variation in the ΔCt of the pair GAPDH vs RPL23 was found to be low. (B): Comparison of the differences in the distribution of Ct values between GAPDH and RPL23. The standard deviation was less with RPL23 evident from the inter-quartile region. Therefore, the RPL23 was selected as the housekeeping gene for the present study

Relative mRNA expression

The genes chosen for the present analysis belong to different functional classes and thus were expected to minimize the chances of them being co-regulated. The RPL23 normalized expression values for the genes, SWI5, RPL36A, CCDC174 and OIT3 were higher in high-fertile as compared to sub-fertile bulls (Figure 2). Among these genes, the expression levels of only CCDC174 were found to have a significant difference (p = 0.02) between the two groups (Figure 3A). The expression levels (fold change) of all the genes, except for PRM1(−1.23) were up-regulated in the high-fertile bulls (SWI5:1.86, RPL36A:4.66, CCDC174:6.20, OIT3:1.53 and TMCO2:1.21) as compared to sub-fertile bulls (fold change = 1).

Figure 2. Changes in the expression levels (ΔCt) of the genes between the two fertility groups. The changes in the Ct values as compared to the housekeeping gene (RPL23) were used for the present analysis. The expression levels of CCDC174 was significantly (p = 0.02) up-regulated in the high-fertile (HF; n = 8) as compared to the sub-fertile (SF; n = 7) bulls. The expression levels of CCDC174 as compared to other genes SWI5, TMCO2, PRM1, RPL36A and OIT3 may influence bull fertility

Figure 3. Univariate regression model for the prediction of bull fertility. (A): The relative mRNA expression levels of CCDC174 in the high-fertile (HF) was 6.2 fold (p = 0.02) higher than the sub-fertile bulls. (B) Univariate analysis of CCDC174 as a prognostic marker employing SPSS software revealed that CCDC174 expression levels had a sensitivity of 87.5% and an accuracy of 86.7% for predicting fertility rate. The area under the curve is 0.83 at the chosen cut-off of 1.12 where the likelihood ratio is 6.12. The analysis suggests that the expression levels of CCDC174 could be a potential marker for predicting bull fertility

Correlation of mRNA expression levels (-ΔCt) with the conception rate

The mRNA expression levels of CCDC174 (r = 0.5) had a significant correlation (p = 0.05) with the conception rate, but the expression levels of other genes viz SWI5 (r = 0.3) and RPL36A (r = 0.3) had only medium correlation. It was also observed that the expression levels of OIT3 (r = 0.06), TMCO2 (r = −0.1) and PRM1 (r = −0.2) had very trivial correlation with the conception rate.

Univariate analysis

The fertility prognostic capability of the expression levels of each of these genes was evaluated using the ROC curve analysis (Table 2). The CCDC174 expression levels had a significant (p = 0.002) influence on the fertility prediction efficiency. The AUC for CCDC174 was 0.83 with a prediction accuracy of 86.7% at a likelihood ratio of 6.12 (Figure 3B). The sensitivity and specificity of the CCDC174 in fertility prediction were 87.5% and 85.7% respectively, at the chosen cut-off value of 1.12. Similarly, for SWI5, RPL36A and OIT3 genes, the prognostic efficiencies were evaluated by ROC curve analysis at the chosen cut-off with a maximum likelihood ratio, but the accuracies of prediction were approximately 60% with the AUC ranging 0.50–0.66 and the likelihood ratio of 2.28–3.42.

Table 2. Univariate analysis of the expression levels of genes, SWI5, TMCO2, PRM1, RPL36A, CCDC174 and OIT3 as a prognostic marker for predicting bull fertility.

Genes	AUC	p value	Sensitivity (%)	Specificity (%)	Cutoff	Likelihood Ratio	Accuracy (%)
SWI5 (S)	0.50	1.00	28.6	87.5	11.3	2.28	60.0
TMCO2 (T)	0.66	0.29	42.9	87.5	2.83	3.42	66.7
PRM1 (P)	0.51	0.90	28.6	87.5	3.47	2.28	66.7
RPL36A (R)	0.64	0.35	28.6	87.5	4.66	2.28	60.0
CCDC174 (C)	0.83	0.002	87.5	85.7	1.12	6.12	86.7
OIT3 (O)	0.58	0.56	28.6	87.5	1.96	2.28	60.0

The expression levels of the CCDC174 gene had significant (p = 0.002) influence on bull fertility prediction sensitivity (87.5%) and accuracy (86.67%).

Multivariate analysis

Multivariate linear regression models were developed by analyzing the possible combinations of these genes expression levels, which would be able to predict the fertility rate with better accuracy. Among the 51 models evaluated, the predicted fertility values generated from 16 models had a significant correlation with the field conception rate. The models with two genes, CCDC174 + TMCO2 had 80% prediction accuracy with the conception rate (Figure 4). The models with three genes, CCDC174 + RPL36A + TMCO2 had a maximum prediction accuracy of 86.6% (Table 3). The models with four genes, CCDC174 + RPL36A + TMCO2 + OIT3 further improved the prediction accuracy to 93.3%. This model had a maximum coefficient of determination (R² = 0.68) with the highest sensitivity (100%) and specificity (87.5%) at the chosen cutoff. The predicted fertility value of this model had a significant correlation (p = 0.001) with the field conception rate. Furthermore, two significant models with five genes had the same prediction accuracy, sensitivity and specificity as that of the CCDC174 + RPL36A + TMCO2 + OIT3 with a likelihood ratio of 8.

Table 3. Development of the multivariate regression models for the prediction of bull fertility

Gene combination models	AUC	p value	Sensitivity (%)	Specificity (%)	Cutoff	Likelihood Ratio	Accuracy (%)
C + T	0.87	0.01	71.4	87.5	42.2	5.74	80.0
C + R + T	0.91	0.01	85.7	87.5	42.8	6.80	86.6
C + R + T + O	0.98	0.001	100.0	87.5	43.6	8.00	93.3
C + R + T + O + P	0.98	0.001	100.0	87.5	43.8	8.00	93.3
C + R + T + O + S	1.00	0.001	100.0	87.5	43.9	8.00	93.3

Overall, 17 models were found to have a significant influence on bull fertility prediction. The multivariate ROC curve analysis of the significant models revealed, C + R + T + O, C + R + T + O + P and C + R + T + O + S were found to have maximum sensitivity and accuracy of prediction. Among these three models, four genes combination model (C + R + T + O) was sufficient to attain maximum bull fertility prediction efficiency. C, R, T, O, P and S denotes the expression levels of genes CCDC174, RPL36A, TMCO2, OIT3, PRM1 and SWI5, respectively.

Figure 4. Multivariate regression models having a significant correlation with the conception rate. Multivariate analysis was carried out employing SPSS software. Progressive inclusion of expression levels of genes in multivariate models improved fertility prediction accuracy as revealed based on ROC curve analysis. The present model suggested that the expression levels of four genes, CCDC174, TMCO2, RPL36A and OIT3 are sufficient to attain maximum sensitivity. The predicted fertility rate of C + R + T + O, C + R + T + O + P and C + R + T + O + S models had the maximum correlation (r = 0.8) with the field conception rate C, R, T, O, P and S denotes the genes CCDC174, RPL36A, TMCO2, OIT3, PRM1 and SWI5, respectively

Protein-protein interaction analysis

The protein-protein interaction analysis carried out for the respective genes revealed that these proteins interact and co-regulate many other proteins influencing sperm function and fertility. For instance, CCDC174 interacts with EIF4A3 and PPP1R35 (Figure 5A). SWI5 along with SFR1, XRCC2 and RAD51 might promote DNA repair activities (Figure 5B). RPL23 accompanied by a group of other RPLs (Figure 5C) involved in the translational machinery. OIT3 may act through its co-localized protein, ABHD2 for successful capacitation (Figure 5D). TMCO2 with DGCR14 and SRRM2 (Figure 5E) and PRM1 in association with TNP1 and TNP2 might regulate spermatid development (Figure 5F).

Figure 5. Interaction analysis among the fertility genes. The corresponding proteins of the studied genes CCDC174, SWI5, RPL36A, OIT3, TMCO2 and PRM1 did not interact among each other. CCDC174 in interaction with EIF4A3 and PPP1R35 may regulate translation-dependent mechanism (A). The SWI5 protein along with SFR1, XRCC2 and RAD51 promote DNA repair activities (B), whereas PRM1 in interaction with TNP1 and TNP2 modulate DNA packaging during spermatid development (F). A group of RPL proteins including RPL36A regulate the translation process (C). OIT3 along with co-expressed protein, ABHD2 may counteract the damaging effects of the ROS, nitrogen species and toxins (D). TMCO2 in interaction with DGCR14 and SRRM2 may regulate acrosome biogenesis and spliceosome activities (E)

Discussion

The present study was conducted to evaluate some of the abundantly expressed spermatozoal mRNAs in predicting bull fertility. Though the fertility rates differed significantly, the sperm functional parameters did not vary between the groups suggesting that the sperm functional attributes are not sufficient to identify the sub-fertile bulls (Selvaraju et al. 2008; Somashekar et al. 2015). This warrants studies at the molecular level to assess the subtle differences that are influencing the sperm function and fertility. The gene expression studies in sperm are invariably associated with many challenges including the selection of housekeeping genes. The choice of housekeeping gene is an important factor for calculating the relative gene expression levels in any cell/tissue. The limited gene expression studies on sperm suggested that the commonly reported housekeeping genes for other cells or tissues were not appropriate for the sperm (Lalancette et al. 2008; Feugang et al. 2010). In some studies, ACTB (Chen et al. 2015), RPLP0 (Arangasamy et al. 2011) and RPL23 (Parthipan et al. 2017) were used as housekeeping genes for expression studies in bull sperm. In the present study, though GAPDH and RPL23 had a similar average standard deviation in ΔCt, considering the compactness of the inter-quartile region and stability values, RPL23 was chosen as a housekeeping gene . This is similar to the findings of the earlier study in neat semen from our laboratory (Parthipan et al. 2017).

The relative mRNA expression levels of SWI5, RPL36A, CCDC174, OIT3 and TMCO2 in sperm were up-regulated and importantly, CCDC174 was found to be 6.2 fold higher in high-fertile as compared to sub-fertile bulls. Further, the expression levels of the CCDC174 had a significant correlation (p = 0.02) with the conception rate. These differentially expressed genes may have a role during spermatogenesis and sperm function. To the best of our knowledge, the involvement of sperm-borne CCDC174 in sperm fertility regulation and especially on improving conception rate has not been reported in any species. The CCDC174 along with EIF4A3 is a major component of the exon junction complex, which is involved in RNA processing, translation and non-sense mediated decay in the neurons (Volodarsky et al. 2015). The CCDC174 together with PPP1R35 might regulate epididymal sperm maturation and motility through testis/sperm enriched protein phosphatase interacting proteins (Goswami et al. 2019). Earlier studies revealed that CCDC genes are regulating sperm fertility. For example, CCDC87 knockout mice were found to have reduced initial motility and progesterone-induced acrosome reaction (Wang et al. 2018). The mutations in other CCDC genes, CCDC39 and CCDC40 perturbed 9 + 2 microtubule cilia and sperm motility leading to the radial spoke defect (Antony et al. 2013). Since the transcriptional and translational activities are limited in sperm, the sperm containing CCDC174 may also probably be translated in the oocyte upon fertilization to support early embryonic development.

SWI5 is a DNA repair gene essential during cell division and also for early embryonic development (Akamatsu et al. 2010; Menezo et al. 2010). The SWI5 protein in interaction with SFR1, XRCC2 and RAD51 might promote DNA repair activities during homologous recombination, which is suggested to be a critical process for the prevention of infertility (Yuan and Chen 2011). Though we did not perform an experiment to explore the DNA repair activities in the present study, a positive association was found between the expression levels of the SWI5 gene with field conception rate. The bulls with an efficient DNA repair mechanism might produce sperm that maintain the DNA quality during various stresses and favor the fertilization process.

RPL36A, a component of the 60S ribosome was 4.66-fold up-regulated in the high-fertile bulls. The biological process of the RPL36A is chordate embryonic development (Provost et al. 2013) and this gene along with other RPLs might regulate the translational mechanism during spermatogenesis and fertilization (Naz 1998; de Mateo et al. 2011).

In our study, OIT3 was up-regulated in the high-fertile sperm, but its role in sperm fertility has not been documented. It can be speculated that the OIT3 may improve bull fertility by regulating sperm function either directly or indirectly. OIT3 is known to have a zona pellucida binding domain, which may directly facilitate the binding of sperm to zona pellucida. The indirect role of this gene in regulating sperm fertility might be through enhancing capacitation function through its co-localized gene, ABHD2. In addition, the OIT3 is involved in the urate homeostasis (Yan et al. 2012) and such a process is essential in the sperm to counteract the damaging effects of the ROS, nitrogen species and toxins (Banihani 2018). Another up-regulated gene in the current study was TMCO2, which is reported to affect acrosomal biogenesis (Kaneko et al. 2020). TMCO2 in interaction with other proteins such as DGCR14 and SRRM2 might regulate spliceosome activities. Any defects in spliceosome activity in testis lead to non-obstructive azoospermic male infertility (Wu et al. 2016). The string network analysis revealed that the co-expressed gene for TMCO2 is OXCT2A, which is a key enzyme involved in ketone metabolism and plays an important role in the energy metabolism of epididymal spermatozoa (Suryawanshi et al. 2013) suggesting that TMCO2 might have a role in the sperm maturational process. Though PRM1 is involved in DNA packaging during spermatid development, the expression levels of PRM1 did not differ between the two groups indicating that the sperm carry leftover mRNAs that had roles The mean fertility rate (%, p <during spermatogenesis (Aoki et al. 2006).

We subjected the expression levels of all the genes either independently or in combination for developing the fertility prognostic model. The selected genes for the present study did not have any direct interaction among each other. However, the predictive ability of the expression levels of the gene, CCDC174 to classify the bulls either into high-fertile or sub-fertile was efficient as the AUC, sensitivity, specificity and accuracy of prediction were found to be high in the univariate prediction method. In the multivariate analysis, all the significant combination models invariably had the CCDC174 gene as an important predictor. The sensitivity and specificity increased upon the addition of other genes in the multivariate model. However, in this study, the model generated by including all the six genes does not have any significant improvement in prediction capability over the model with five genes. From this, it is evident that the right combination of the genes has to be selected for improving the overall prediction accuracy (Gerszten and Wang 2008). Accordingly, the present study suggests that the model with CCDC174, RPL36A, TMCO2 and OIT3 genes could identify the high-fertile bulls with 100% accuracy at the chosen cut-off wherein the likelihood ratio was 8. Though a likelihood ratio of 10 and above was suggested as a strong model (Gerszten and Wang 2008; Nadeau et al. 2013), the present panel is a moderate evidence model for the prediction of fertility rate. Overall, the model comprised of various genes regulating different facets of spermatogenesis to fertilization (Figure 6). The fact that the TMCO2 regulates germ cell development as evidenced from the literatures and protein-protein interactions. The CCDC174 may influence sperm maturational events, sperm functional attributes, and RNA processing and translation activities. OIT3 can balance the metabolic activity of the spermatozoa during stress and promote zona binding activity of the spermatozoa. The RPL36A may synchronize the embryonic development through its translational activity. The present gene expression study in relation to field conception rate provided insight for understanding the mechanism behind selective retention of certain transcripts during spermatogenesis and their importance in improving field fertility.

Figure 6. The fertility prediction model developed by incorporating the sperm genes involved in fertility regulation. The resultant model included four genes for predicting bull fertility. They are, CCDC174 influences sperm maturational events, sperm functional attributes, and RNA processing and translation activities; TMCO2 is associated with germ cell development; OIT3 balances the metabolic activity of the spermatozoa during stress and to promote zona binding activity of the spermatozoa and RPL36A through translational activity may synchronize the embryonic development. Thus, the expression levels of these genes are predictors of bull fertility

The present study suggests that the CCDC174 could be a potential marker for assessing the fertility status of the bulls. The combined expression levels of four genes, TMCO2, RPL36A, CCDC174 and OIT3 may be used for the identification and selection of the potential high-fertile bulls for the AI program.

Materials and methods

Semen sample collection

Frozen semen samples from 15 Holstein Friesian bulls were procured from the Nandini Sperm Station, Karnataka Milk Federation, Hessarghatta, Bengaluru, India. The frozen-semen samples were stored at −196°C until further analyses. The average fertility rate (%) of the contemporary bulls (n = 74) maintained in this Sperm Station was 42.87 ± 0.4. The selected bulls for this study were classified into high-fertile (n = 8) and sub-fertile (n = 7) based on the fertility rate (%) obtained from approximately 2000 inseminations per bull. The mean fertility rate (%, p < 0.001) for the high-fertile bulls was +3.25 (46.12 ± 0.51%), whereas for the sub-fertile bulls it was −3.48 (39.39 ± 0.69%). The semen from all these bulls was used for carrying out inseminations in a particular agro-climatic region by a group of inseminators. The cows with apparently normal estrous symptoms were inseminated without considering the age and parity of cows. Since large numbers of inseminations were carried out to obtain field conception rates, the influence of environmental, managemental and cow effects on the bull fertility were considered to be minimal. The sperm functional parameters were analyzed from the same batches of frozen semen samples, which went to the field.

Sperm functional parameters

The frozen semen samples were thawed in a water bath at 37°C for 30 s and then subjected to the sperm functional parameters analysis.

Kinematics

The kinematic parameters of the frozen-thawed semen samples were analyzed using computer-assisted semen analyzer (CASA, Sperm class analyzer, Microptic, Spain). The thawed semen samples were diluted five times (~5 million sperm/ml) in 1X phosphate-buffered saline (PBS, pH 7.3) and 10 μl of the diluted semen sample was analyzed in CASA by capturing 10 homogenous fields/sample (25 frames/second) as described previously (Selvaraju et al. 2012). The set parameters in CASA were: cell area (range 10 to 70 μm²); sperm with VCL < 10, 10 to 25, 25 to 50 and >50 μm/s were classified as static, slow, medium and rapid motile, respectively; progressive motile sperm, > 70% straightness (STR); circular motile sperm, < 50% linearity (LIN); curvilinear velocity (VCL, μm/s); the average path velocity of the sperm head along its actual trajectory); straight-line velocity (VSL, μm/s: the average path velocity of the sperm head along a straight line from its first to the last position); average path velocity (VAP, μm/s: the velocity of the sperm head along its average path); linearity (LIN %: the ratio between VSL and VCL); straightness (STR %: the ratio between VSL and VAP).

Structural membrane integrity

The structural membrane integrity was assessed by mixing the frozen-thawed semen sample (10 μl) with 10 μl of eosin (5%) and 10 μl nigrosin (10%) (Selvaraju et al. 2008). The semen-dye mixture was smeared onto a clean glass slide and air-dried. A minimum of 400 sperm was counted under the 100x objective in the phase-contrast microscope (Nikon Eclipse 80i, Nikon, Japan). Sperm that did not take the eosin dye were considered as structural membrane intact and those that stained with eosin dye were considered as structural membrane integrity lost. The percentage of sperm with structural membrane integrity was calculated considering the total number of cells counted.

Functional membrane integrity

Functional membrane and acrosome integrities were assessed using the hypoosmotic swelling-Giemsa (HOS-G) test. The semen samples (50 μl) were incubated in 450 μl of hypo-osmotic (150 mOsm) and 450 μl of iso-osmotic (300 mOsm; control) solutions at 37°C for 30 min. After incubation, the semen samples were smeared onto a clean glass slide, air-dried and then stained with Giemsa (Selvaraju et al. 2008). A minimum of 400 sperm was counted under the oil immersion objective (100x) in a phase-contrast microscope (Nikon Eclipse 80i, Nikon, Japan). The sperm subpopulation positive for functional membrane intactness (FMI) as observed from the hairpin bend in the principal piece of the tail and acrosome intactness in 150 mOsm was calculated by subtracting the corresponding population in 300 mOsm [FMI (%) = FMI at 150 mOsm (%) – FMI at 300 mOsm (%)]. The sperm acrosome intactness in 300 mOsm was considered for calculating the percentage of acrosomal integrity.

Mitochondrial membrane potential

The sperm mitochondrial membrane potential was assessed using cationic dye, 5, 5ʹ,6,6ʹ-tetrachloro −1,1ʹ,3,3ʹ-tetraethyl benzimidazolyl carbocyanine iodide (JC-1). The semen samples (10 μl) were added with 1.53 mM of JC-1 stain (50 μl) dissolved in phosphate-buffered saline (pH 7.3) and incubated for 30 min in dark. After incubation, the semen samples were washed in PBS and smeared onto the glass slide, air-dried and counted under oil immersion (100x) in the fluorescent microscope with an excitation filter of 510–560 nm and a barrier filter of 505 nm (Nikon Eclipse 80i, Nikon, Japan). A minimum of 400 sperm per sample was counted and the sperm having an orange colored mid-piece (J-aggregates in the mitochondria) were considered to have high mitochondrial membrane potential (Selvaraju et al. 2008).

Chromatin distribution

The sperm chromatin distribution was assessed based on Feulgen’s staining (Selvaraju et al. 2008). Briefly, the sperm samples were smeared onto glass slides and treated with 5 N HCl for 30 min followed by washing and staining with Schiff’s reagent. At least 200 sperm per sample were observed under the 100x objective in a phase-contrast microscope (Nikon Eclipse 80i, Nikon, Japan). The sperm with evenly distributed chromatin was considered normal and the sperm with pyriform shaped, unevenly distributed chromatin and vacuole in the sperm head was considered abnormal. The percentage of sperm with normal chromatin distribution was calculated for the total number of sperm observed.

Gene expression studies

Total RNA isolation and cDNA synthesis

The total RNA from bovine sperm was isolated by employing double lysis followed by a membrane-based kit method (Goodrich et al. 2007) with minor modifications as standardized in the Reproductive Physiology Laboratory of ICAR-NIANP, Bangalore (Parthipan et al. 2015). Briefly, frozen-thawed semen samples were washed in Bovipure (Nidacon, Sweden) solution to remove the egg yolk and other contaminating cells. The sperm (30–40 × 10⁶ million cells) were subjected to double lysis followed by eluting the RNA using the membrane-based kit (PureLink RNA mini kit, Invitrogen, USA). The isolated RNA samples had a 260/280 ratio of 1.88 ± 0.01 and 260/230 ratio of 1.86 ± 0.10 (Supplementary Table 2) as assessed using a spectrophotometer (NanoDrop, ND-1000, Thermo Scientific, USA). The isolated RNA was subjected to DNase digestion (TURBO DNase, Ambion, USA) and quantified using a fluorometer (Qubit 4.0, Invitrogen, USA). An equal concentration of RNA from each semen sample was reverse transcribed into cDNA using the SuperScript IV (Thermo Scientific, USA) in 20 μl reactions using OligodT primer. The cDNA was stored at −20°C until qPCR analysis.

Real-time PCR

The expression levels of mRNAs were determined using the SYBR Green Mastermix (TB Green Premix Ex Taq II, Takara Bio, Japan) and 1.2ng of cDNA per reaction in real-time PCR (Step One Plus, Applied Biosystems, USA). Each sample was checked with an intron spanning primer (PRM1) for genomic DNA contamination. The contamination of RNAs with other cells was also checked using primers for PTPRC (leukocytes), CDH1 (epithelial cells) and KIT (germ cells) genes. The samples free from genomic DNA and other cells’ RNA contaminants were selected for the gene expression studies. The primers used for the genes SWI5, TMCO2, RPL36A, CCDC174 and OIT3 (Table 4) were checked for their specificity using melt curve and agarose gel electrophoresis. The efficiency of each primer was calculated to be 90 to 110%. The cycle conditions used were: initial denaturation of  95°C for 30 s; 40 cycles of  95°C for 3 s and 60°C for 1 min and with the default melt curve settings of the instrument. The Cq data were acquired and analyzed using StepOne software (v2.2.2).

Table 4. The list of primers used for the assessment of contamination and gene expression studies

S.No	Gene ID	Primer Sequence (5ʹ to 3ʹ)	Primer Length (bp)	Product Size (bp)	NCBI Accession Number
1	PRM1*	F – AAGATGTCGCAGACGAAGGAG	21	222	NM_174156.2
		R – GTGGCATTGTTCGTTAGCAGG	21
2	KIT**	F – GAATAGCTGGCATCAGGGTG	20	224	AF263827.1
		R – CCAGATCCACATTCTCTCCATC	22
3	PTPRC**	F – TGGACGAAATTGCATCCCTCAGGA	24	237	NM_174156.2
		R – TGGTCAGGACGTTTACAGCTCACA	24
4	CDH1**	F – CTGCATTCCTGGCTTTGGTG	20	171	NM_001002763.1
		R – GTAAGCACGCCATCTGTGTG	20
5	GAPDH	F – CTGAGGACCAGGTTGTCTCCTG	22	141	NM_001034034.1
		R – CCCTGTTGCTGTAGCCAAATTC	22
6	RPL23	F – CAGCGGTGGTAATTCGACAAC	21	116	NM_001035014.2
		R – GGCGGAACCTTTCATCTCG	19
7	SWI5	F – GAGCTCACAGCCCTTTGTCC	20	138	XM_010810563.1
		R – GACTCGGCTTGGGCTCTCTA	20
8	TMCO2	F – ACTACTGCAGCTGCTCTGACTG	22	122	NM_001077077.2
		R – CCACCGCACACACAAGGAAA	20
9	RPL36A	F – CGCTCTCGCAAACATGGTGA	20	166	NM_001034316.2
		R – CTGCCCACCATAGCCACTCT	20
10	CCDC174	F – GGGTACTGGTCGAAGAGGCA	20	112	NM_001075578.1
		R – CCTGGCTGCTCAGAGACTCA	20
11	OIT3	F – TATGAGCCGCAGTCACGGAA	20	103	NM_001046064.1
		R – TCGCAGCTTGAGAGTAGGCA	20

*The assessment of genomic DNA contaminations was studied using intron spanning PRM1. **The contaminating RNA from germ cells, leukocytes and epithelial cells were assessed using KIT, PTPRC and CDH1, respectively.

Selection of housekeeping gene and relative quantification of mRNA expression levels

For the gene expression studies in sperm, the housekeeping gene has not been well-established and hence, in the present study, three probable housekeeping genes such as GAPDH, RPL23 and PRM1 expression levels were compared across all the samples. The gene pair with the lowest ΔCt and the gene with the less average deviation in ΔCt values across the samples was chosen as the housekeeping gene (Silver et al. 2006). The validity of this analysis was also verified using NormFinder software (Andersen et al. 2004). The fertility rate was correlated with -ΔCt values (Ct of target gene – Ct of RPL23) of target genes in each bull. The relative quantification of mRNA expression levels was estimated using 2^−ΔΔCt method (Livak and Schmittgen 2001).

Development of fertility prediction model

The fertility prediction models were developed using SPSS software. As a first step, the predictive power of the expression levels of each gene was examined in the univariate analysis. Subsequently, multivariate regression analysis was performed by combining the expression levels of each of these genes to predict the ability of each of the combined models to classify a bull based on fertility status. Linear regression models were developed with the conception rate as an independent variable and ΔCt of the genes as a dependent variable. The resultant predicted values calculated with the model-fit statistics for each model were then correlated with the conception rate. The receiver operating characteristic (ROC) curve analysis was performed to evaluate the sensitivity and specificity for every point in the analysis to evaluate the diagnostic efficiency of a predictor for classifying the bulls into either as high-fertile or sub-fertile group. The area under the curve (AUC), sensitivity (%), specificity (%), the accuracy of prediction (%) was determined at the chosen cut-off with maximum likelihood ratio value. The percentage of sensitivity is the number of high-fertile bulls that are correctly classified. The percentage of specificity is the number of sub-fertile bulls that are correctly classified. The percent accuracy is the closeness of the predicted value to that of the true conception rate.

Statistical and bioinformatic analyses

The percentage values were log-transformed before subjecting to statistical analyses using IBM SPSS statistics 20 and GraphPad Prism 6. The difference in gene expression levels between the high-fertile and the sub-fertile groups were analyzed using the Student’s t-test. The Pearson correlation analysis was employed to assess the relationship among the sperm functional parameters, genes expression pattern and conception rate. All the values were presented as mean ± SEM. The significance was set at p < 0.05. The protein-protein interaction analysis was studied using STRING (version 11.0, https://string-db.org/) with Homo sapiens as background to assess the influence of associated or co-expressed proteins on the sperm fertility regulation.

Ethical approval for animal study

The experiment was conducted with the approval of the Institutional animal ethical committee (IAEC).

Author contributions

Conceptualization and funding acquisition: SS; Designed the experiment: SS, DS, and LR; Conducted the experiment: DS, LR, LM and SSA; Analyzed the data and drafted the manuscript: SS, DS and MS; All authors reviewed and finalized the manuscript.

Abbreviations

Table

AI:	artificial insemination;
PRM1:	Protamine 1;
CCDC174:	coiled coil domain containing 174;
RPL36A:	ribosomal protein L36a;
OIT3:	onco-protein induced transcript 3;
TMCO2:	Transmembrane and coiled coil domains 2;
SWI5:	SWI5 recombination repair homolog;
ROS:	reactive oxygen species;
GAPDH:	Glyceraldehyde 3- phosphate dehydrogenase;
RPL23:	ribosomal protein L23;
ACTB:	actin beta;
RPLP0:	ribosomal protein P0;
PPV:	predicted pregnancy value;
CASA:	Computer assisted sperm analysis;
LIN:	Linearity index;
STR:	Straightness index;
VCL:	Curvilinear velocity;
VSL:	Straight-line velocity;
VAP:	Average path velocity;
HOS-G:	Hyposmotic swelling-Giemsa;
HPAP:	Hos positive and acrosome positive;
HPAN:	Hos positive and acrosome negative;
HNAP:	Hos negative and acrosome positive;
HNAN:	Hos negative and acrosome negative;
JC-1:	5, 5ʹ,6,6ʹ-tetrachloro −1,1ʹ,3,3ʹ-tetraethyl benzimidazolyl carbocyanine iodide;
PTPRC:	Protein tyrosine phosphatase receptor type C;
CDH1:	Cadherin 1;
KIT:	KIT proto-oncogene, receptor tyrosine kinase;
ROC:	Receiver operating characteristic;
AUC:	Area under the ROC curve;
EIF4A3:	Eukaryotic translation initiation factor 4A3;
PPP1R35:	Protein phosphatase 1 Regulatory Subunit 35;
ABHD2:	Abhydrolase domain containing 2;
OXCT2A:	3-Oxoacid Co-A Transferase 2A.

Acknowledgments

The authors are grateful to Dr. Raghavendra Bhatta, Director, ICAR-NIANP, Bengaluru, India for his critical technical inputs and necessary facilities to carry out this work. The authors acknowledge Dr. J.P. Ravindra Former Head, Animal Physiology Division, ICAR-NIANP, Bangalore for proofreading of the manuscript. The authors also acknowledge Dr. S. Parthipan, Senior Research Fellow, ICAR-NIANP, Bangalore for providing technical input.

Disclosure statement

All the authors declare no conflict of interest.

Supplementary material

Supplemental data for this article can be accessed here.

Correction Statement

This article has been republished with minor changes. These changes do not impact the academic content of the article.

Additional information

Funding

This research was carried out under the ICAR-National Fellow Project funded by the Indian Council of Agricultural Research, Ministry of Agriculture, Government of India. Dr. S. Selvaraju is supported by the ICAR-National Fellow Project, ICAR, Ministry of Agriculture, Government of India