Genetic background of semen parameters in Italian Simmental bulls

Abstract

The aim of this study was to estimate genetic parameters and investigate the genomic background of scrotal circumference and semen parameters in the Italian Simmental bulls. Scrotal circumference, number of normal spermatozoa, ejaculate volume, spermatozoa motility, and total spermatozoa were measured on 622 bulls, of which 603 had genotypes for 42,141 SNP. Variance components of scrotal circumference were estimated with an animal model that included the fixed effects of birth year, animal age, and measurement method, and the random effects of day of measurement and animal. In the model for the other traits, the scrotal circumference was added as a covariate to account for its influence on the semen parameters. A genome-wide association study was carried out using the ssGBLUP-approach. Heritabilities ranged from 0.07 ± 0.05 (spermatozoa motility) to 0.50 ± 0.14 (scrotal circumference). A total of 13 SNP passed the Bonferroni correction threshold and the number of significantly associated markers ranged from 1 (ejaculate volume and spermatozoa motility) to 5 (total spermatozoa). Genes already associated with reproduction parameters were retrieved close to the significant SNP. Results of the present study gave preliminary insights about the genetic determinism of semen quality in Italian Simmental bulls.

Highlights

• Low to moderate heritabilities were estimated for semen traits.

• Few markers were associated with the phenotypes, suggesting their polygenic determinism.

• Candidate genes already associated in literature with sperm traits were found.

Keywords:

• Genome-wide association study
• fertility
• single-step genomic BLUP

Introduction

The Italian Simmental (IS) is dual-purpose cattle breed mainly farmed in small herds located in the mountainous area of North-eastern Italy. The breeding program is handled by the Italian Simmental Breeders Association (ANAPRI). The beef aptitude is evaluated in approximately 250 young bulls every year through the performance test (Cesarani et al. 2020), where average daily gain, body size, muscularity, and feet and legs are evaluated. At the end of the performance test, the bulls are used for progeny testing for milk production traits through artificial insemination (the top 15% of the animals, about 20–30 bulls per year) or natural service (the top 16% to 25%); the remaining 60% are slaughtered. The composition of the dual-purpose selection index is: 44% milk production, 24% meat production, 19.5% morphology, and 12.5% fitness. Even if the IS breed is generally characterised by good fertility (e.g. quite short calving to conception interval and calving interval, low non-return rate), the genetic selection for milk production could reduce reproductive performances, because of the negative genetic correlation between these two traits (e.g. Walsh et al. 2011; Ma et al. 2019). In order to prevent or limit the reduction of fertility, the latter has been included in breeding programs worldwide (Philipsson and Lindhé 2003; Miglior et al. 2005). The genetic improvement of fertility is usually based on female phenotypes (such as conception or pregnancy rates, calving interval, and oestrus index) some of which are difficult to measure routinely at population level. Fertility problems could, however, also be due to male subfertility (Han and Peñagaricano 2016). Male fertility parameters are now routinely evaluated in the IS at the end of the PT. Thus, they could be quite easily incorporated into the breeding program. However, an important preliminary step is the assessment of the genetic variability and the genetic determinism of these traits. The aim of this work was to estimate the heritability of some male fertility traits in the IS breed and to investigate their genomic background using a genome-wide association study.

Materials and methods

Data

The study was carried out on a sample of 622 IS bulls intended for both artificial insemination or natural service. Available phenotypes were: scrotal circumference (SC), number of normal spermatozoa (NS), ejaculate volume (EV), spermatozoa motility (SM), total spermatozoa (TS). Since data were recorded in routine measurements at the end of the Performance Test, the Animal Care and Use Committee approval was unnecessary.

These male fertility traits were phenotypically analysed in Corte Pause et al. (2022). Briefly, at the end of the performance test, the SC of each bull is measured transversely using a tape (Reliabull Scrotal Tape, Lane Manufacturing Inc., Denver, CO, USA) according to Hopkins and Spitzer (1997).

At the ned of the PT, the top ranked bulls were examined for semen quality as outlined by the manual of the SFT (Society of Theriogenology; Koziol and Armstrong 2018). The volume (mL) of the ejaculation was measured using a 15 mL collection tube (code 62.554.502, Sarsted s.r.l., Italy) connected to the artificial vagina . Total motility (defined as the percentage of spermatozoa that exhibit motility of any form) and progressive motility (defined as the percentage of spermatozoa that are moving in a rapid linear fashion) were estimated visually, by an experienced veterinarian in charge of the Bull breeding soundness evaluation at the genetic centre, from 10 microscopic fields at 37 °C using a phase-contrast microscope at 200X magnification (CH 2, Olympus, Japan) equipped with a heated stage (HT 50, Minitube, Germany). For the motility evaluation, raw semen was first diluted to 1:60 in sterile 0.9% NaCl solution. Spermatozoa motility was evaluated after placing a 15 µL drop of semen on a pre-warmed glass slide, covered by a 20 × 20 mm coverslip. The spermatozoa concentration was measured using a photometer (Accucell, IMV-Tecnologies, France) by diluting 40 µl of semen in 3960 µl of 0.9% NaCl, and the total number of spermatozoa was calculated by multiplying the ejaculate volume for sperm concentration. For the morphological evaluation of spermatozoa, a smear was prepared by mixing 7 µL of raw semen with 90 µL of eosin-nigrosin solution. Morphological characteristics of spermatozoa (percentage of normal and abnormal spermatozoa) were assessed using bright-field microscopy, at 1000X magnification (Orthoplan, Leitz, Germany), by counting at least 200 sperm cells as defined by Barth and Oko (1989). . Medium density SNP genotypes were available for 603 bulls. After quality control, 42,141 SNP mapped on the 29 autosomal chromosomes of the ARS-UCD 1.3 assembly were retained for the analysis. Phenotypes, pedigree, and genotypes were provided by ANAPRI.

Heritability estimation

Variance components were estimated using a pedigree-based animal model. The model used for the scrotal circumference was: [1] $$SC=YoB+Age+Method+Date+animal+e$$[1] where: YoB was the fixed effect of the year of birth of the bull (9 levels); Age was the covariate of the age of the bull at the measurement day (average of 427 ± 19 days, from 384 to 499); Method was the fixed effect of the measurement method (two different types of tape); Date was the random effect of the measurement day (95 levels); animal was the random effect of the animal (3,044 total animals in the relationship matrix); e was the random residual effect.

The animal model adopted for the semen parameters was the following: [2] $$y=YoB+Age+SC+Date+animal+e$$[2] where YoB, Age, and Date are described above, and SC was the covariate of the scrotal circumference, that was included because of its strong influence on the semen parameters.

In both models, the animal effect was distributed as $\mathrm{N}(0,\mathbf{A}{\mathrm{\sigma }}_{\mathrm{a}}^2)$ where A is the pedigree relationship matrix and ${\sigma }_a^2$ is the variance associated with the additive genetic effect. The random effect of the date was distributed as $\mathrm{N}(0,\mathrm{I}{\mathrm{\sigma }}_{\mathrm{date}}^2)$ where I is an identity matrix and ${\mathrm{\sigma }}_{\mathrm{date}}^2$ is the variance associated with the random effect of the measurement day. Finally, the random error term was distributed as $\mathrm{N}(0,\mathrm{I}{\mathrm{\sigma }}_{\mathrm{e}}^2)$ where I is an identity matrix and ${\mathrm{\sigma }}_{\mathrm{e}}^2$ is the variance associated with the residual error. As suggested by several studies (e.g. Lourenco et al. 2014; Pocrnic et al. 2017; Cesarani et al. 2021), the pedigree was traced back 3 generations from the animals with phenotypes and genotypes (total of 3,044 animals).

Variance components were estimated via a Gibbs sampling algorithm implemented in GIBBSF90+ program (Misztal et al. 2014). From a single chain with 100,000 samples, the first 10,000 were discarded as burn-in, and one sample every 10 was saved to compute means and standard deviations of the posterior distributions using POSTGIBBSF90 (Misztal et al. 2014). Heritability (h²) across date, i.e. not considering in the denominator the variance associated with the date of analysis, was computed as: $$h^2={\sigma }_a^2/({\sigma }_a^2+{\sigma }_e^2)$$ with the variances described above.

Genome-wide association study and gene mapping

The estimated variance components were used in the genome-wide association study (GWAS) carried out using single-step genomic approach (ssGBLUP; Legarra et al. 2009), where the A inverse is replaced by the H inverse (Aguilar et al. 2010), which is a relationship matrix that combines both pedigree (A) and the genomic relationships (G). Once the (G)EBV are estimated with the ssGBLUP, they can be backsolved to marker’s effects because of the equivalence between GBLUP and SNP-BLUP (Wang et al. 2012). Moreover, the individual SNP variances were estimated. Models [1] and [2] were used for the GWAS, with the animal effect distributed as $\mathrm{N}(0,\mathbf{H}{\mathrm{\sigma }}_{\mathrm{a}}^2)$ where H is described above.

The p-values for all SNP were computed based on their prediction error variance according to Aguilar et al. (2019) using the POSTGSF90 software (Misztal et al. 2014). A marker was considered significantly associated with p < 0.05. Moreover, the Bonferroni correction was applied as the ratio between alpha and the number of independent blocks, which was computed using the –blocks flag of PLINK (v. 1.9; Purcell et al. 2007). Using the NCBI online database (https://www.ncbi.nlm.nih.gov/genome/gdv/) genes mapped near each significant SNP (± 250 kb; Cesarani et al. 2019; Manca et al. 2020) were retrieved.

Results and discussion

Variance components estimation

Heritability for SC was of the same magnitude as in a report (0.46 ± 0.07) of a composite cattle population from Nebraska (Russell et al. 2021). Siqueira et al. (2012) reported h² of 0.45 ± 0.02 and 0.63 ± 0.05 for SC in Nellore bulls, offspring of super-precocious or precocious females, respectively. Lower h² were estimated for SC in Mexican Charolais-Charbray (0.18 ± 0.02; Martínez-Velázquez et al. 2020) and in Red Angus (0.32 ± 0.02; McAllister et al. 2011) cattle breeds. Carvalho Filho et al. (2020) analysed the genetic parameter for fertility traits in Nellore bulls and they reported h² of 0.47 ± 0.12 for SC and 0.07 ± 0.08 for progressive sperm motility, which is very close to the h² estimated in the present work for spermatozoa motility (0.07 ± 0.05). The heritability estimated for SM (Table 1) was larger than in a report in Austrian Simmental (0.04 ± 0.01; Gredler et al. 2007), a population closely related to the IS. In the same study the authors reported h² of 0.18 ± 0.02 and 0.22 ± 0.02 for ejaculate volume and total number of spermatozoa, respectively. Similar h² were reported for ejaculate volume in a multi-breed cattle population (0.20 ± 0.04) by Berry et al. (2019) and in Holsten bulls (0.22 ± 0.05) by Druet et al. (2009), whereas lower h² (0.09 ± 0.08) was estimated in Hereford bulls (Kealey et al. 2006). Heritability of 0.09 ± 0.04 was estimated for the number of spermatozoa in Holstein bulls (Druet et al. 2009), which is close to the estimate of the present study (0.11 ± 0.08).

Table 1. Heritability of the investigated traits.

	Variance components
Trait	Sampling Date	Additive	Residual	Heritability
Scrotal circumference	0.90 ± 0.24	1.90 ± 0.60	1.89 ± 0.51	0.50 ± 0.14
Normal spermatozoa	37.54 ± 10.32	29.85 ± 15.97	127.57 ± 16.41	0.16 ± 0.10
Ejaculate volume	0.07 ± 0.05	0.55 ± 0.20	1.19 ± 0.18	0.32 ± 0.11
Spermatozoa motility	38.94 ± 14.94	21.22 ± 19.31	283.58 ± 24.45	0.07 ± 0.05
Total spermatozoa	0.48 ± 0.17	0.37 ± 0.29	3.01 ± 0.32	0.11 ± 0.08

The h² estimates reported in literature for the analysed traits are quite heterogeneous, mainly because the male fertility traits are characterised by a larger variability compared to the female ones, as pointed out by Berry et al. (2019). Typically, heritability for semen traits as well as other reproduction traits is lower (e.g. Ma et al. 2019), compared to morphological ones such as SC, which usually has moderate to high heritability. Moreover, these lower heritabilities estimated for the semen traits (i.e. normal spermatozoa, total spermatozoa and their motility, and ejaculate volume) compared to the SC was expected because the former are largely influenced by the environmental factors (Carvalho Filho et al. 2020). The quite large impact of the environment on the semen traits was considered in the present study by fitting the date of sampling as random effect in the model. The sampling date explained from 4% (EV) to 19% (SC and NS) of the total variance.

Genome-wide association study

Figure 1 shows the Manhattan plot of the GWAS analyses. A total of 11,630 SNP was significant without the Bonferroni correction (Table 2). The Bonferroni correction was implemented in the GWAS procedure to account for multiple comparisons and to decrease the false positive rate (Johnson et al. 2010). When the sample size is small (as the case of the present study) and the correction for the multitude of comparisons is too conservative (as the case of the Bonferroni correction) the rate of false negative could increase (Johnson et al. 2010). The Bonferroni correction is considered overly conservative also because close SNP are in linkage disequilibrium and, thus, they are correlated and tend to be inherited together in blocks that are not independent. In order to prevent possible over-corrections when applying Bonferroni (i.e. increase in the false negative rate), we used the number of independent blocks (N = 866) instead of the number of SNP in the denominator of the Bonferroni formula. After applying this Bonferroni correction,13 SNP were found to be significant in the GWAS analysis (dots above the red line in the Figure 1): the number of markers significantly associated with each trait ranged from 1 for EBV and SM, to 5 fort TS, respectively.

Figure 1. Manhattan plot of the genome-wide association study results. The blue line indicates the significance threshold (0.05), whereas the red line indicates the Bonferroni threshold (0.05 divided by the number of blocks).

Table 2. Result of the genome-wide association study and markers found to be significant with or without the Bonferroni correction¹.

	SNP		Bonferroni correction^a
Trait^b	p < 0.05	Variance explained	BTA	SNP	Range	Variance explained
SC	1,970	26.12%	4	Hapmap41721-BTA-101476	79,878,883–80,378,883	0.05%
			4	BTB-00208407	109,870,572–110,370,572
NS	1,881	25.56%	2	BTA-110520-no-rs	12,405,403–12,905,403	0.00%
			25	ARS-BFGL-NGS-25988	34,064,708–34,564,708
EV	2,026	27.58%	10	BTB-01434227	50,304,831–50,804,831	0.08%
SM	1,864	25.36%	18	BTA-42525-no-rs	16,067,150–16,567,150	0.02%
TS	2,004	27.53%	5	Hapmap51306-BTA-74962	107,777,461–108,277,461	0.07%
			10	ARS-BFGL-NGS-104551	47,420,717–47,920,717
			11	BTA-85510-no-rs	14,540,793–15,040,793
			19	Hapmap60163-rs29015084	21,306,972–21,806,972
			29	ARS-BFGL-NGS-29984	49,198,377–49,698,377

^aBonferroni = 0.05/866 (number of blocks)

^bBTA: Bos Taurus autosome; SC: scrotal circumference; NS: normal spermatozoa; EV: ejaculate volume; SM: spermatozoa motility; SNP: single nucleotide polymorphism; TS: total spermatozoa

Table 2 also shows the percentage of the trait variance explained by the significant SNP before and after the Bonferroni correction, respectively. The variances ranged from 25.36% (NS) to 27.58% (EV) when considering the SNP before Bonferroni, and from 0.00% (NS) to EV (0.08%) after the correction. The latter is interesting because just one SNP (located on BTA10) can explain slightly less than 1% of the ejaculate volume; this could have interesting implications for the weighted version of the ssGBLUP (Zhang et al. 2016; Alvarenga et al. 2020).

Gene mapping

A total of 57 genes were mapped in the significant regions (Supplementary Table 1), some of which have been previously reported to be related to the sperm traits. GWAS on semen traits have been carried out in several cattle breeds (Peñagaricano et al. 2012; Puglisi et al. 2016; Sweett et al. 2020), as well as on other species such as sheep (Serrano et al. 2021), goats (Talouarn et al. 2020), and pigs (Diniz et al. 2014; Marques et al. 2018), in order to decipher their genomic backgrounds.

A GWAS carried out on Italian Holstein bull sperm quality (Puglisi et al. 2016) found a significant region on BTA18 at 16.77 Mb, which is very close to the region highlighted in the present study for SM (BTA18, 16.07–16.57 Mb). Cai et al. (2017) analysed the male infertility in cattle yak and reported two genes found also in the present study: the CNTNAP2, located in the region on BTA4 significant for SC, and the CCL24, located in the region on BTA25 significant for NS. The upregulation of CCL24 in the scrotal skin could indicate a subclinical infection (González-Barrio et al. 2021). The ZP3 gene, located on chromosome 25 and associated with NS, has sperm-binding activity (Kanai et al. 2007).

Among the genes found in the significant regions, the PHKB located on BTA18 is of particular interest. The PHKB has been associated to spermatozoa motility, in pathways related to calcium homeostasis (Mahé et al. 2022); the relationship between calcium and several sperm functions is well known (Bernecic et al. 2019).

The MDH2 gene, which was found in the present study in a region associated with NS, was listed among the overexpressed genes in high-fertile bull spermatozoa (Aslam et al. 2018). This gene is involved in the final steps of the Krebs cycle, and its reduction may negatively impact the motility, capacitation, hyper-activation, and fertilising ability of spermatozoa. Verma et al. (2014) reported the GCNT3 gene as potentially associated with sperm fertility in buffalo.

Two genes mapped in the regions on BTA5 associated with TS were previously associated with bull fertility traits: the WNT5B, reported to be associated with semen volume, number of sperm and sperm motility in a Thai multibreed dairy population (Sarakul et al. 2018), and the ADIPOR2, involved in the organisation of sperm structural and functional traits (Kasimanickam et al. 2013). In a GWAS for male fertility carried out in Holstein-Friesian bulls, the SPAST gene was reported to be associated with poor sperm motility (Hering et al. 2014). This gene was found in the present study among the candidate genes for TS. The abundance of GOSR1 gene, located on chromosome 19, in spermatozoal transcripts was associated with higher fertility dairy sires by Card et al. (2017). The CD81 gene has implications in the intactness of the acrosome of spermatozoa (Jankovicova et al. 2016).

The POR gene, which is expressed exclusively in the late stages of spermatogenesis, has been associated with spermatids and spermatozoa in humans (Keber et al. 2013), whereas IGF-2, H19, and UPK3B genes were associated with sperm traits and maturation of gametes in mice (Yu et al. 2003; Kuriyama et al. 2017).

Interestingly, some of the genes found in the present study were previously found to be associated with female fertility in cattle. Killeen et al. (2016) compared the endometrial gene expression in high and low fertility heifers, and they reported three genes found also in this study: TMIGD1, associated with TS, and GTF2A2 and GCNT3 associated with EV, respectively (Supplementary Table 1). The YWHAG gene, found in the present study associated with NS, is involved in the spontaneous abortions in cattle (Oliver et al. 2019) and it has also implications with intramammary infections (Chen et al. 2015). Finally, Guarini et al. (2019) analysed the reproductive disorders in Canadian Holstein, and they reported three genes potentially related to cystic ovaries in cattle that were associated with NS in the present study: SRRM3, DTX2 and UPK3B genes.

Conclusions

Results of the present work highlighted the possibility of including male fertility parameters in the breeding program of the Italian Simmental cattle breed. The traits seem to be polygenic, with a small proportion of the genetic variance explained by many markers. A few markers were associated with the male fertility traits, and genes in these significant regions are of interest for future studies to further understand the biological pathways that affect male fertility traits.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

Restrictions apply for the datasets used for this study since they belong to the Italian Simmental Breeders Association (ANAPRI), and they were used under agreement for this research.

Correction Statement

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