Autosomal recessive loci contribute significantly to quantitative variation of male fertility in a dairy cattle population

Background Cattle are ideally suited to investigate the genetics of male fertility. Semen from individual bulls is used for thousands of artificial inseminations for which the fertilization success is monitored. Results from the breeding soundness examination and repeated observations of semen quality complement the fertility evaluation for each bull. Results In a cohort of 3881 Brown Swiss bulls that had genotypes at 683,609 SNPs, we reveal four novel recessive QTL for male fertility on BTA1, 18, 25, and 26 using haplotype-based association testing. A QTL for bull fertility on BTA1 is also associated with sperm head shape anomalies. All other QTL are not associated with any of the semen quality traits investigated. We perform complementary fine-mapping approaches using publicly available transcriptomes as well as whole-genome sequencing data of 125 Brown Swiss bulls to reveal candidate causal variants. We show that missense or nonsense variants in SPATA16, VWA3A, ENSBTAG00000006717 and ENSBTAG00000019919 are in linkage disequilibrium with the QTL. Using whole-genome sequence data, we detect strong association (P = 4.83 × 10− 12) of a missense variant (p.Ile193Met) in SPATA16 with male fertility. However, non-coding variants exhibit stronger association at all QTL suggesting that variants in regulatory regions contribute to variation in bull fertility. Conclusion Our findings in a dairy cattle population provide evidence that recessive variants may contribute substantially to quantitative variation in male fertility in mammals. Detecting causal variants that underpin variation in male fertility remains difficult because the most strongly associated variants reside in poorly annotated non-coding regions. Supplementary Information The online version contains supplementary material available at 10.1186/s12864-021-07523-3.


Background
Male fertility is a complex trait that is determined by genetic and non-genetic sources of variation. Because of its low heritability, a large sample size is required to investigate the genetic architecture of male fertility. Large cohorts of males with repeated measurements for reproductive traits are not available in many species including humans, but are accessible in cattle where semen samples from individual bulls are used for thousands of artificial inseminations [1]. Because these data facilitate the disentanglement of male and female factors contributing to establishing pregnancy [2], the fertility of bulls can be quantified objectively.
Semen quality has a large effect on insemination success [3][4][5][6][7]. On entering the semen collection center, the semen quality of each bull is examined as part of an indepth breeding soundness evaluation [8]. The semen quality of artificial insemination bulls varies over time, consequently the fertilization rates may differ substantially between ejaculates. In order to ensure high and uniform insemination success, parameters such as ejaculate volume, sperm concentration, sperm head and tail morphology, and sperm motility are examined in all ejaculates immediately after semen collection. Ejaculates that fulfill predefined quality requirements are diluted and filled into doses that contain between 15 and 25 million sperm, depending on semen quality. After subsequent cryopreservation, at least the frozen-thawed sperm motility is checked. Computer-assisted sperm analysis or flow cytometric assays can add additional information on the semen quality [9]. The combination of multiple semen quality parameters may enable predicting bull fertility to a certain extent [10][11][12][13][14].
Phenotypes describing male reproductive traits and dense microarray-derived genotypes for many artificial insemination sires enable genome-wide association testing. Quantitative trait loci (QTL) for bull fertility have been detected in Holstein [15][16][17][18][19] and Jersey [20] cattle. Nevertheless, these studies did not attempt to reveal the causal variants underlying the QTL with whole-genome sequence information. Moreover, except in [19], the sample sizes were relatively small providing low statistical power to detect QTL for a trait with low to moderate heritability like bull fertility. Recent investigations provide growing evidence that non-additive effects contribute to inherited variation in bull fertility [16,[21][22][23][24]. However, the underlying genetic mechanisms remain largely unknown.
Here, we apply haplotype-based association testing to detect QTL for bull fertility and semen quality in a large mapping cohort of bulls from the Brown Swiss cattle breed that had dense microarray-derived genotypes. Our association studies revealed five recessive QTL for bull fertility. Using two complementary sequence-based fine-mapping approaches and transcriptome data, we identify candidate causal variants in coding sequences of genes that are expressed in the male reproductive tract.

Results
To detect QTL for bull fertility, we considered 3736 Brown Swiss (BSW) bulls that were kept at semen collection centers in Switzerland, Germany and Austria. Bull fertility was quantified based on artificial insemination success adjusted for confounding factors (see Methods). The bulls had partially imputed genotypes at 589,791 autosomal SNPs with minor allele frequency greater than 0.5%. Genomic restricted maximum likelihood estimation indicated that the autosomal markers explained 0.10 ± 0.02% of the phenotypic variance in bull fertility.
Recessive QTL for bull fertility are located on chromosomes 1, 6, 18, 25 and 26 Haplotype-based association tests that were based on an additive mode of inheritance revealed three QTL located on chromosomes 1, 6, and 26 for bull fertility that exceeded the Bonferroni-corrected genome-wide significance threshold of 2.68 × 10 − 7 (Table 1, Fig. 1a).
We previously showed that the QTL for bull fertility at BTA6 is in linkage disequilibrium with the recessive BTA6:58,373,887 T-allele (rs474302732) that activates cryptic splicing in WDR19, resulting in reduced semen quality [22]. In order to investigate the association of the BTA6:58,373,887 T-allele in our mapping cohort, we imputed the genotypes of the 3736 bulls to the whole-genome sequence level using a breed-specific reference panel of 368 sequenced cattle. A sequence-based association analysis that was based on a recessive mode of inheritance confirmed that rs474302732 is strongly associated (P = 2.72 × 10 − 45 ) with bull fertility. The association with bull fertility was stronger for rs474302732 than the most significantly associated haplotype (Table  1). Four non-coding variants (57,900,948 bp, 57,950,075 bp, 57,408,389 bp, 57,810,181 bp) had slightly lower P values (P = 5.36 × 10 − 46 -3.81 × 10 − 47 ) than rs474302732. Due to its large effect on male fertility, we subsequently fixed the top haplotype of the BTA6 QTL as a covariate in the association model.
When the top haplotype of the BTA6 QTL was fixed as a covariate in the model, the association signal at BTA6 disappeared (Fig. 1b). In the conditional analysis, the P value of the top-associated haplotype at the BTA1 QTL did not meet the Bonferroni-corrected significance threshold by a small margin (Fig. 1b). The QTL on BTA26 remained significant (P = 9.96 × 10 − 11 ).
Next, we repeated the haplotype-based association test assuming a recessive mode of inheritance. The association with bull fertility was more pronounced (i.e., the P values of the top haplotypes were lower) at the QTL on chromosomes 1, 6 and 26 (Table 1, Fig. 1c & d). Moreover, the recessive model revealed two QTL at chromosomes 18 and 25 that were not detected using the additive model. Visual inspection of the haplotype effects corroborated that the four novel QTL detected on chromosomes 1, 18, 25, and 26 act in a recessive manner (Fig. 1e, f, g, h). Compared to non-carrier and heterozygous haplotype carrier bulls, the fertility of bulls that carry the BTA1, BTA18, BTA25 and BTA26 top haplotype in the homozygous state is reduced by 0.80 ± 0.10, 0.28 ± 0.05, 1.09 ± 0.20 and 0.52 ± 0.06 phenotypic standard deviations, respectively. A linear regression analysis conditional on the top 10 principal components and the top haplotype at the BTA6 QTL confirmed that fertility does not differ between heterozygous and non-carrier bulls at the four novel QTL (P BTA1 = 0.29, P BTA18 = 0.61, Between 24 and 434 out of the 3736 bulls were homozygous carriers of fertility-decreasing haplotypes at the four novel and one previously detected QTL. Across the five QTL, 2827, 604, 259, 44, and 2 bulls were homozygous carriers of 0, 1, 2, 3, and 4 top haplotypes, respectively. In order to quantify the impact of the five recessive QTL on insemination success in the BSW population, we investigated the 56-day non-return rate (NRR56; mean = 65.06 ± 4.53%) in cows inseminated with semen from 1322 BSW bulls, estimated from at least 300 inseminations per bull. The NRR56 is the proportion of cows that is not re-inseminated within a 56day interval after the first insemination and is considered as a proxy for the insemination success. The mean NRR56 is 0.5 phenotypic standard deviations higher (P = 2.75 × 10 − 8 , two-tailed t-test) in 1000 bulls that do not carry any of the top haplotypes in the homozygous state than in 322 bulls that are homozygous for at least one of the five recessive QTL ( In order to unravel the mechanism through which the QTL impacts male fertility, we analysed routinely collected semen quality data that were available for 32 homozygous, 372 heterozygous and 498 non-carrier bulls. The BTA1 QTL was not associated with ejaculate volume (P = 0.07), sperm concentration (P = 0.19), sperm motility (P = 0.10), sperm head morphology (P = 0.53) or sperm tail morphology (P = 0.19, Supplementary Figure  1) assessed in fresh semen. However, the insemination straws contained more (+ 1.57 million, P = 4.78 × 10 − 6 , Supplementary Figure 1f) sperm in homozygous than heterozygous and non-carrier bulls, which indicates an attempt to compensate for low sperm quality by providing additional sperm per straw, although the association with the number of sperm per straw did not meet the Bonferroni-corrected significance threshold (P Bonf = 2.68 × 10 − 7 ) by a small margin. Moreover, more ejaculates were discarded due to insufficient semen quality from homozygous than heterozygous and non-carrier bulls (14.8 vs. 6.5%).
The association pattern of the BTA1 top haplotype was puzzling. Homozygosity for the top haplotype compromises bull fertility but not routinely recorded fresh semen quality. However, the insemination doses of homozygous bulls contain more sperm per straw. In order to further examine this apparent contradiction, we studied 1577 sperm morphology evaluations that were collected as part of the andrological examination of 575 prospective artificial insemination bulls. During the andrological examination, sperm morphology of at least one ejaculate per bull is systematically evaluated, before the semen is used for artificial inseminations. This evaluation is repeated either until the semen quality meets predefined quality thresholds or the bull is rejected for artificial insemination due to insufficient semen quality. Sperm morphology evaluations were provided for between 1 and 13 (median = 2) ejaculates of 575 Swiss BSW bulls. From these reports, we derived 21 sperm quality phenotypes (see Methods).
Genome-wide haplotype-based association testing revealed that the QTL at BTA1 affects four (out of 21) sperm morphology features. Bulls that are homozygous for the top haplotype produce ejaculates that contain less normal spermatozoa (− 10.6%, P = 1.72 × 10 − 10 ), twice as many sperm with non-compensatory defects (+ 5.4%, P = 9.07 × 10 − 16 ), 1.5 times more sperm with major defects (+ 8.9%, P = 2.55 × 10 − 11 ) and twice as many sperm with head shape anomalies (+ 5.6%, P = 3.98 × 10 − 17 ) (Fig. 2, Supplementary Figure 1 h, i, l, & m). These data suggest that an increased proportion of sperm with abnormal head morphology compromises bull fertility, as no other sperm morphology feature was affected by the BTA1 QTL. The number of sperm morphology evaluations was also higher (P = 1.63 × 10 − 10 ) for homozygous than heterozygous and noncarrier bulls (Supplementary Figure 1 g). Eventually, only 21 out of 34 homozygous bulls (61.7%) compared to 520 out of 561 (92.7%) heterozygous and non-carrier bulls passed the sperm morphology examination and produced ejaculates that are suitable for breeding (≥ 65% normal sperm and ≤ 20% non-compensatory defects). Sperm morphology did not differ between heterozygous haplotype carriers and non-carrier bulls, corroborating recessive inheritance.
We examined and imaged fresh sperm from two bulls homozygous for the top haplotype using oil-immersion phase-contrast light microscopy. At semen collection, the bulls were 644 and 504 days old. Some spermatozoa had morphological anomalies (i.e., rounder and shorter heads, heads with rounder frontal part (pear-shape /pyriform), narrower head base (tapered), and heads with abnormal contour (uneven shaped)) that are classified as head shape anomalies in the sperm morphology evaluation ( Supplementary Figures 2 & 3). However, the vast majority of the spermatozoa (85 and 84.6%) were normal. Acrosome defects were not apparent.
A missense variant in and non-coding variants upstream SPATA16 are associated with the BTA1 QTL In order to detect candidate causal variants for the increased number of abnormal sperm and reduced fertility associated with the BTA1 QTL, we applied a two-step fine-mapping approach. First, we used whole-genome sequence data of 125 BSW bulls to identify variants that are compatible with the inheritance of the top haplotype. The average fold coverage of the sequenced bulls was 10.1 ± 4.3-fold. The BTA1 QTL top haplotype was determined for the sequenced bulls based on their microarray-derived genotypes. Second, we imputed whole-genome sequence variant genotypes for the 3736 BSW bulls of the mapping cohort using 368 reference animals to perform an association study between imputed sequence variant genotypes and bull fertility.
One and 40 sequenced bulls carried the top haplotype in the homozygous and heterozygous state, respectively. Of 44,948 sequence variants that were polymorphic within a window encompassing 3 Mb on either side of Fig. 2 Effect of the BTA1 top haplotype on sperm morphology. Boxplots representing the sperm morphology of non-carrier, heterozygous, and homozygous (haplotype status 0, 1, and 2) bulls. The numbers above the boxplots indicate the number of bulls in the respective group. Compared to non-carrier and heterozygous bulls, a the number of sperm morphology examinations is twice as high (5.32 ± 3.25 vs. 2.49 ± 2.50), b the amount of normal spermatozoa is~10 percentage points lower (66.19 ± 10.43% vs. 76.45 ± 8.51%), c the amount of sperm with noncompensatory defects is~5 percentage points higher (10.61 ± 6.28% vs. 5.79 ± 3.19%), d the amount of sperm with major defects is~8 percentage points higher (24.08 ± 8.99% vs. 15.70 ± 6.82%), e and the amount of spermatozoa with head shape anomalies is~5 percentage points higher (10.49 ± 6.30% vs. 5.42 ± 3.09%) in homozygous bulls the top haplotype, we detected 764 variants between 93, 614,265 and 96,742,540 bp that were compatible with recessive inheritance, i.e. these variants were heterozygous in haplotype carriers and homozygous in the bull that carried the haplotype in the homozygous state. The sequence coverage did not differ between heterozygous, homozygous and non-carrier bulls within that interval.
Of the 764 compatible variants, 505 were also significantly associated with bull fertility at the genome-wide Bonferroni-corrected significance threshold of 3.88 × 10 − 9 (Fig. 3a). The significantly associated variants clustered between SPATA16 encoding spermatogenesis associated protein 16 and NLGN1 encoding neuroligin 1 (Fig. 3a). The most significantly associated variant (BTA1:93,972,058G > A, rs379712951) is in an intron of NLGN1, 373 kb upstream the translation start site of SPATA16, within the top window from the haplotypebased association study. The P value is slightly lower for the BTA1:93,972,058G > A variant than for the most significantly associated haplotype (P = 3.00 × 10 − 17 vs. 4.16 × 10 − 17 ). We did not consider NLGN1 as a candidate gene for an impaired bull fertility, as it is not notably expressed in adult bull testis (TPM < 1). However, we considered SPATA16 as a positional and functional candidate gene for the BTA1 QTL, because it is testes-specific expressed in cattle (http://cattlegeneatlas. roslin.ed.ac.uk) and human (https://gtexportal.org/ home/gene/SPATA16). Moreover, SPATA16 mRNA is highly abundant in testis tissue of adult bulls (TPM = 295) (Fig. 3b).
Two of the significantly associated variants compatible with recessive inheritance were in coding regions: a missense variant in SPATA16 (rs440830663 at BTA1:94,396,804A > G, ENSBTAP00000053460.3: p.Ile193Met, P = 4.83 × 10 − 12 ) and a synonymous variant in GHSR encoding growth hormone secretagogue receptor (rs714884352, BTA1:95,029,804C > T, ENSBTAP00000014446.5:p.228Leu, P = 3.7 × 10 − 9 ). Pathogenic alleles of human SPATA16 are associated with globozoospermia, i.e., round-headed, often acrosome lacking spermatozoa [25,26]. The p.Ile193Met variant compatible with recessive inheritance resides in an evolutionarily conserved tetratricopeptide repeat domain of SPATA16 (Fig. 3c). The isoleucine at position 193 is on the surface of the protein (Fig. 3d). According to protein structure modelling with Mis-sense3D, a methionine at position 193 does not alter the tertiary structure of SPATA16. However, it is predicted to be deleterious to SPATA16 function (SIFT score: 0.03, PolyPhen-2 score: 0.998). The P value is lower and the effect on male fertility (− 0.80 vs. -0.57) is larger for the most significantly associated intergenic than the p.Ile193Met variant in SPATA16.
We performed haplotype and sequence-based association studies conditioned on either the top associated haplotype, the SPATA16:p.Ile193Met variant (rs440830663) or intronic BTA1:93,972,058G > A, (rs379712951) to disentangle the QTL (Supplementary Figure 4a,b,c). Conditioning on the top associated haplotype removed the signal in the haplotype study, but not in the sequence-based study (Supplementary Figure 4a). When the SPATA16:p.Ile193Met variant was fixed as covariate in the haplotype and sequence-based association model (Supplementary Figure  4b), the original top haplotype and the variants upstream SPATA16 were still associated with bull fertility, albeit not at the Bonferroni-corrected significance threshold. The P values of the most significantly associated variant (BTA1: 93,972,058 bp) and of the top haplotype from the conditional analysis were 1.52 × 10 − 6 and 1.74 × 10 − 6 . When the association analysis was conditioned on the BTA1:93,972, 058 bp variant, i.e., the most significantly associated variant from the sequence-based association study, the QTL signal was absent in both association studies ( Supplementary Figure 4c). In our set of partially imputed sequence variants, BTA1:94,396,804A > G and BTA1:93,972,058G > A were in high linkage disequilibrium (r 2 = 0.858). ENSBTAG00000006717 encoding the nondescript ATP-binding cassette sub-family A member 3-like protein is highly expressed in testis tissue of adult bulls (127 TPM, Fig. 4b) and shows a testis-specific expression in cattle [27] (http://cattlegeneatlas.roslin.ed.ac.uk). The frameshift partly truncates an evolutionarily conserved domain "ATPase associated with a variety of cellular activities (AAA)". VWA3A encoding von Willebrand factor A domain-containing protein 3A is expressed in testis tissue of adult bulls at 17 TPM. The stop mutation truncates VWA3A by 58%. The von Willebrand factor A domain-containing protein 3A is a motile ciliaassociated protein and might contribute to the beating movement of the flagella [28].

Coding variants in
An association study between imputed sequence variants and bull fertility revealed no variants exceeding the Bonferroni-corrected significance threshold (3.88 × 10 − 9 ) at the BTA25 QTL (Fig. 4a). However, 38 variants that did not meet the significance threshold by a small margin (P < 7.66 × 10 -7 ) were between 18,963,390 and 20, 033,799 bp. The most significantly associated variant (P = 4.56 × 10 − 8 ) was at 19,675,705 bp in an intron of PDZD9 encoding PDZ domain containing 9. The P values were considerably larger for the frameshift and nonsense variant in ENSBTAG00000006717 (P = 7.98 × 10 − 3 ) and VWA3A (P = 5.11 × 10 − 5 ), respectively.  46. The two coding variants in ENSBTAG00000019919 were in complete linkage disequilibrium and highly significantly associated with bull fertility (P = 1.3 × 10 − 15 ) in the sequence-based association study (Fig. 5a). The p.Asn576Lys-variant is predicted to be deleterious (SIFT: 0.02, PolyPhen-2: 0.921) to protein function and the similarity to the amino acid sequence of other mammalian species is weak in Clustal Omega protein alignment (Fig. 5c). ENSBTAG00000019919 is expressed at 21 TPM in adult bull testis (Fig. 5b). In Chlamydomonas reinhardtii, a CFAP46 homolog is essential for normal flagellar motility [29].
The association study with imputed sequence variant genotypes revealed 18 variants that were stronger associated with bull fertility than the p.Asn576Lys-variant. The two variants with the lowest P value (P = 1.5 × 10 − 24 ) were BTA26:50145932C > T (rs720936782) and BTA26: Fig. 4 Detailed view of a QTL for bull fertility on BTA25. a Manhattan plot representing the association of imputed sequence variants (diamonds) and haplotypes (bars) with bull fertility. Red framed bars represent haplotypes that met the Bonferroni-corrected significance threshold. The orange and blue diamonds represent the nonsense variant (p.Arg505Ter) in VWA3A and frameshift variant (p.Ile1167LeufsTer) in ENSBTAG00000006717. Symbols over the Manhattan plot indicate genes and transcripts that are annotated at the QTL region. Genes and transcripts with testis-biased expression, male and female reproductive tract specific expression, and testis and brain specific expression, are indicated in blue, violet and green colour, respectively. b Transcript abundance quantified in TPM in testis tissue of adult (blue) and newborn cattle. Genes and transcripts with testis-biased expression, male and female reproductive tract specific expression, and testis and brain specific expression, respectively, are indicated in blue, violet and green colour, respectively 50178052C > T (rs210673393), that are 4375 bp downstream of the stop codon of SYCE1 encoding synaptonemal complex central element protein 1 and 2171 bp upstream of the start codon of CYP2E1 encoding cytochrome P450, family 2, subfamily E, polypeptide (Fig. 5). SYCE1 is expressed at 167 TPM in adult bull testis, whereas CYPE1 is not expressed (TPM < 1) (Fig. 5b). However, both variants were not compatible with the inheritance of the top haplotype.
In order to disentangle the QTL, we performed haplotype and sequence-based association studies conditioned on either the top associated haplotype, the p.Asn576Lysvariant in ENSBTAG00000019919 (BTA26:50,850, 915C > G, rs378141069), or the BTA26:50145932C > T (rs720936782) variant downstream SYCE1 (Supplementary Figure 5a,b,c). The conditional analyses indicated that the top haplotype and the p.Asn576Lys-variant capture the QTL variance only partially, because the rs720936782 variant was still associated with bull fertility (P = 1.1 × 10 − 9 and 8.5 × 10 − 9 ), albeit when conditioned on rs378141069 not at the stringent Bonferronicorrected significance threshold (3.88 × 10 − 9 ). However, the QTL was absent when the analyses were conditioned on rs720936782 (Supplementary Figure 5c). In our set of the imputed sequence variant genotypes, rs378141069 and rs720936782 were in moderate linkage disequilibrium (r 2 = 0.598).

Discussion
In a mapping cohort of 3736 BSW bulls from which semen was collected at semen collection centers under standardized conditions, only 10% of the phenotypic variation of bull fertility was explained by autosomal variants. This estimate suggests that the additive genetic heritability of bull fertility is low in BSW cattle, which agrees with findings in other cattle populations [30][31][32][33][34]. It is possible that we underestimated the genomic heritability of bull fertility, because the genomic relationship matrix was constructed from autosomal SNPs only. Variants on the sex chromosomes may also contribute to inherited variation in bull fertility [35,36]. Our mapping cohort provided sufficient statistical power to detect large-effect QTL for a complex trait with an estimated narrow-sense heritability of 10% [37]. However, we did not detect additive QTL. Therefore, we conclude that the additive genetic architecture of bull fertility is polygenic and determined by QTL with small to moderate effects.
Our haplotype-based association analysis identified five recessive QTL with moderate frequencies that, taken together, had a substantial effect on bull fertility in the BSW cattle population. Homozygous haplotype carriers are fertile, but their non-return rate is reduced in artificial inseminations. Our data did not allow investigating bull fertility in natural matings. Phenotypic manifestations of the five recessive QTL detected in our study are less detrimental than effects arising from pathogenic alleles in CCDC189, ARMC3 and TMEM95 that lead to male sterility in cattle [38][39][40]. Yet, the hypothetical exclusion of bulls carrying at least one of the five detected top haplotypes in the homozygous state from breeding would significantly increase the 56-day non-return rate in cows and heifers. Thus, the recessive QTL detected in our study contribute to quantitative variation in bull fertility. Similar findings in Fleckvieh, Jersey, and Holstein bulls suggest that the genetic architecture of male reproductive traits is shaped by non-additive effects [16,20,40]. Repeatability estimates for traits describing male reproductive performance in cattle are substantially higher than narrow-sense heritability estimates in different cattle breeds including Brown Swiss [22,41,42]. In addition to additive effects, the repeatability also includes permanent environmental and non-additive sources of variation. Our study shows that non-additive QTL for bull fertility can be detected in large mapping cohorts using genome-wide association testing. The fact that non-additive effects contribute substantially to the genetic variation of male reproductive traits has implications in efforts to predict bull fertility from dense molecular markers. While improving the fertility of bulls is not a breeding objective in most cattle improvement programs, the genome-based prediction of bull fertility might facilitate implementing compensation strategies for bulls for which the anticipated insemination success is low. Statistical models that consider only additive effects have low predictive power for bull fertility [8,35], whereas methods that also include non-additive effects may improve the accuracy of genomic predictions [19]. Our findings suggest likewise that non-additive effects should be considered in order to improve the accuracy of genomic prediction for male reproductive traits.
The sequence-based fine-mapping of the QTL revealed candidate causal variants in coding sequences of genes that are either testis-specific expressed in cattle or contribute to variation of male fertility in species other than cattle. The top haplotype on BTA6 is in strong linkage disequilibrium with a likely causal variant that activates cryptic splicing in WDR19 [22]. A variant in human WDR19 is associated with morphological abnormalities of sperm and impaired male fertility [43]. The top haplotypes at the BTA1, BTA25 and BTA26 QTL encompass putatively deleterious variants in genes (SPATA16, VWA3A, ENSBTAG00000006717 and ENSBTAG00000019919) that are expressed in the bovine male reproductive tract. Furthermore, SPATA16 and ENSBTAG00000006717 are specifically expressed in the male reproductive tract (http://cattlegeneatlas.roslin.ed.ac. uk) and the produced proteins SPATA16 and ABCA14like are prevalent in spermatozoa [44,45]. Pathogenic variants in SPATA16 and ENSBTAG00000019919 orthologs have been associated with either impaired male fertility in humans or impaired movement of flagella in a flagellated green alga [25,29]. VWA3A is associated with motile cilia in humans [28].
Candidate causal variants for recessive traits are frequently prioritized if they are compatible with the segregation of the top haplotype [22,[46][47][48]. Strictly applying this criterion, the missense and nonsense variants in the coding sequences of SPATA16, VW3A3, ENSBTAG00000006717 and ENSBTAG00000019919 qualify as compelling candidate causal variants for impaired bull fertility; these variants are in linkage disequilibrium with the respective top haplotype and predicted to be deleterious to protein function. Our screen for variants that are compatible with the Fig. 6 Detailed view of a QTL for bull fertility at BTA18. a Manhattan plot representing the association of imputed sequence variants (diamonds) and haplotypes (bars) located on chromosome 18 between 32 and 38.5 Mb. Red framed bars and diamonds represent haplotypes and variants that meet the Bonferroni-corrected significance threshold. Symbols over the Manhattan plot indicate genes and transcripts that are annotated at the QTL region. Blue, violet, and green colour indicates genes and transcripts with testis-biased expression, male and female reproductive tract specific expression, and testis and brain specific expression. b Expression of genes within the segment 34 to 36.5 Mb in adult (blue) and newborn (black) bull testis in TPM. Blue, violet, and green labels indicate genes with testis-biased expression, male and female reproductive tract specific expression, and testis and brain specific expression inheritance of the top haplotype takes genotyping errors into account. However, due to the low number of sequenced bulls that were homozygous haplotype carriers, it is possible that some variants were misclassified. Thus, we performed a sequence based association study as a complementary fine-mapping approach. An association study with imputed sequence variant genotypes showed that the candidate causal variants are also associated with bull fertility to a certain degree. Yet, non-coding variants were stronger associated with bull fertility than the putatively functional coding variants at all QTL. It is possible that the true causal variants reside in regulatory regions, but are in linkage disequilibrium with the coding variants prioritized in our study. Such a pattern could also indicate multiple trait-associated alleles within a QTL [49].
The p.Asn576Lys variant in ENSBTAG00000019919 is a particularly compelling candidate causal variant for the BTA26 QTL because it was (i) in linkage disequilibrium with the top haplotype, (ii) highly significantly associated (P = 1.3 × 10 − 15 ) in the sequence based association study, (iii) a lysine at position 576 is predicted to be deleterious to protein function, and (iv) the expression of ENSBTAG00000019919 is restricted to the male and female reproductive tract. However, the sequence-based association study also revealed highly significantly associated variants nearby SYCE1 that were not in complete linkage disequilibrium with the top haplotype. In fact, an association analysis conditional on the p.Asn576Lys variant did not fully account for the BTA26 QTL variation. The non-coding variant downstream SYCE1 captured the QTL variance better. SYCE1 is also a candidate gene for impaired male fertility, because it is testis-specific expressed at very high levels (TPM > 150). Loss of function of SYCE1 leads to spermatogenic arrest in human and mice [50][51][52]. The Synaptonemal Complex Central Element Protein 1 is part of the central element of the synaptonemal complex that is required for chromosome synapsis in prophase I of meiosis [53]. Errors in the meiotic recombination can lead to sperm aneuploidy [54] and thus to embryonic loss [55,56]. It seems plausible that sequence variation affecting SYCE1 and the p.Asn576Lys variant in ENSBTAG00000019919 contribute to the BTA26 QTL variation.
It is worth mentioning, that the true causal variants may not necessarily be the most significantly associated variants in association studies with imputed sequence variant genotypes due to imputation errors or sampling effects [57][58][59]. Interestingly, also the BTA6:58373887C > T variant (rs474302732) that had been reported as a putatively causal variant for the BTA6 QTL [22], was not the top variant in our association study between imputed sequence variant genotypes and bull fertility. Thus, we consider the coding variants in SPATA16, VWA3A, ENSBTAG00000006717 and ENSBTAG00000019919 as candidate causal variants for impaired male reproductive performance, but further functional investigations and replication studies in independent populations are required to corroborate their association with bull fertility.
Because our association studies considered both bull fertility and semen quality, we were able to reveal the mechanisms through which the BTA1 QTL contributes to variation in fertility. Bulls that are homozygous for the top haplotype encompassing SPATA16 are less fertile likely because their ejaculates contain an increased number of sperm with morphological abnormalities of the sperm head. Abnormal head morphology is negatively correlated with NRR56 [7,60]. However, semen quality of homozygous bulls is not as severely compromised as, e.g., in males with globozoospermia [25,26] and spermatogenic arrest [61], respectively, due to loss-of-function alleles in human and murine SPATA16. It is plausible that the SPATA16:p.Ile193Met variant compromises semen quality in the homozygous state, as the methionine at position 193 is predicted to be deleterious to SPATA16 function. However, we can not preclude that non-coding variants in linkage disequilibrium with p.Ile193Met contribute to the BTA1 QTL variation. In any case, the isoleucine residue at position 193 of SPATA16 is not essential for fertilization, because homozygous bulls are fertile. Sperm head shape abnormalities were only detected in fixed and therefore immotile spermatozoa of homozygous haplotype carriers. Because the morphological abnormalities of the sperm heads are subtle, they remain unnoticed in fresh and motile spermatozoa. In fact, the majority of sperm produced from homozygous haplotype carriers is normal, but the overall proportion of normal sperm is reduced. Many bulls homozygous for the BTA1 QTL are judged to be not suitable for breeding due to insufficient semen quality before their semen is used for artificial inseminations. Thus, the actual effect of the BTA1 QTL on bull fertility is likely more pronounced than observed in our study that was based on insemination success from ejaculates meeting minimum quality requirements.
We found no clues from the analysis of semen quality data why fertility is reduced in bulls that are homozygous carriers of the top haplotypes at the BTA18, BTA25, and BTA26 QTL. Thus, our study confirms indirectly that routine semen quality parameters, once bulls producing ejaculates with clearly insufficient quality are removed, are not sufficient to anticipate bull fertility [62]. Computer-assisted and flow-cytometric sperm analyses may provide further insights into mechanisms that underpin phenotypic variation in bull fertility [11]. Moreover, the phenotype bull fertility does not enable differentiating between fertilization failure and early embryonic losses, because it is based on the 56-day non-return rate. It is thus possible that fertilization is normal, but the apparently reduced bull fertility is primarily a manifestation of early embryonic losses [63,64].

Conclusion
Recessive variants contribute substantially to quantitative variation in male fertility in Brown Swiss cattle. Our study suggests that non-additive effects are an important source of variation for traits describing male reproductive performance. Of the five fertility-associated QTL, only two QTL on BTA1 and BTA6 are also associated with reduced semen quality. Mechanisms through which the other QTL impact bull fertility remain unknown. Coding variants are in linkage disequilibrium with the QTL, but it remains difficult to determine their causality because many strongly associated variants in linkage equilibrium reside in poorly annotated non-coding regions.

Genotypes of BSW bulls
Microarray-derived genotypes of 33,045 BSW cattle (12, 887 males and 20,158 females) were provided by breeding associations. The genotypes were determined using various versions of Illumina BovineSNP Bead chips that interrogated between 20 K and 777 K SNPs. Quality control including the removal of variants and samples with more than 10% missing genotypes and the exclusion of variants deviating significantly (P < 10 − 8 ) from Hardy-Weinberg proportions was carried out using plink (version 1.9) [65].
Using a reference panel of 1166 cattle (a subset of the 33,045 BSW) that had Illumina BovineHD BeadChip genotypes for 777,962 SNPs, the genotypes for animals typed at lower density were imputed to higher density (HD) and phased using Beagle 5.1 [66] with default settings. The final dataset comprised 33,045 animals and 683,609 SNPs. For association testing, we considered only 3881 bulls for which either bull fertility or semen quality data were available.

Bull fertility in Swiss BSW sires
Bull fertility (as at February 2020) was provided by Swissgenetics for 1382 BSW bulls that were used for between 232 and 15,923 first inseminations with conventional (i.e., semen was not sorted for sex) frozen-thawed semen. Bull fertility was estimated based on the 56-day non-return rate after the first insemination in cows and heifers using a linear mixed model similar to the one proposed by Schaeffer et al. [2]: y ijklmnopq = μ + MO i + PA j + PR k + (BS x BK) lm + TE n + h o + s lp + e ijklmnopq , where y ijklmnopq is either 0 or 1, indicating whether or not a subsequent insemination was recorded within a 56-day interval after the first insemination, μ is the intercept, MO i is the insemination month, PA j is the parity (heifer or cow), PR k is the cost of the semen dose, (BS x BK) lm is a combination of bull's breed x cow's breed, TE n is the insemination technician, h o is the random effect of the herd, s lp is the random effect of the bull expressed as deviation from the average non-return rate, and e ijklmnopq is a random residual term. Bull fertility was subsequently standardized to a mean of 100 ± 12. The resulting value is an objective measurement of bull fertility. Three phenotypic outliers (i.e., bulls for which the fertility was more than 5 standard deviations below average) were not considered for subsequent analyses because such bulls might carry rare genetic conditions [40] that can lead to spurious associations [67]. For the GWAS, we considered the 1130 out of the 1382 Swiss BSW bulls which had HD level genotypes that were partially imputed from low-density micro-array derived genotypes. Their fertility was quantified based on 850,708 first artificial inseminations (minimum 232, median 366.5, mean 752.8, and maximum 15,923 per bull).

Bull fertility in German and Austrian BSW sires
Phenotypes for bull fertility were provided by ZuchtData EDV-Dienstleistungen GmbH, Austria, for 4617 BSW bulls from Germany and Austria (as at December 2017) that were used for 4,267,990 and 1,646,254 inseminations in cows and heifers, respectively. Bull fertility was obtained from routine breeding value estimation for reproductive traits, which are jointly estimated for males and females [68]. The statistical model used to estimate breeding values for female reproductive traits included a fixed effect for the service sire that represents bull fertility as deviation from the population mean. We excluded four bulls with exceptionally poor fertility (see above). For association testing, we considered 2606 bulls from the German and Austrian BSW populations which had HD level genotypes that were partially imputed from low-density micro-array derived genotypes. Those bulls were used in a total of 3,746,202 first inseminations (minimum 100, median 360.5, mean 1437.5, maximum 73,858 per bull).

Compilation of a mapping cohort for association testing of bull fertility
To estimate the genetic correlation of bull fertility between the Swiss and German / Austrian bulls, a bivariate restricted maximum likelihood (REML) analysis was performed using the GCTA software (version 1.93.2) [69]. Using the "reml-bivar" option in GCTA [70], we fitted a genomic relationship matrix that was built using 589,791 autosomal SNPs with minor allele frequency greater than 0.5%. To control for effects attributable to population stratification, we additionally fitted the top 10 principal components of the genomic relationship matrix using the "-qcovar" option [71]. The REML estimation revealed that the same genetic determinants control the trait in both mapping cohorts (r g = 1 ± 0.28). We standardized the phenotypes within the Swiss and German / Austrian cohorts to compile a joint mapping population of 3736 bulls. The proportion of phenotypic variation explained by 589,791 autosomal SNPs with minor allele frequency greater than 0.5% was estimated using GCTA.
Semen quality of Swiss BSW bulls from routine fresh semen assessment Semen quality data were provided by Swissgenetics for 74,945 ejaculates that were collected from 1432 BSW bulls between January 2000 and November 2019. Parameters assessed for fresh ejaculates were semen volume (in ml), sperm concentration (million sperm per ml) quantified using photometric analysis, and sperm motility (percentage of sperm with forward motility) assessed visually using a heated-stage microscope at 200-fold magnification. Moreover, the amount of sperm with head and tail anomalies was classified for each ejaculate with scores ranging from 0 to 3 (0: no or very few anomalies, 1: less than 10% sperm with anomalies, 2: between 10 and 30% sperm with anomalies, 3: more than 30% sperm with anomalies). Ejaculates that fulfilled the requirements for artificial insemination (semen volume above 1 ml, more than 300 million sperm per ml, at least 70% motile sperm, no apparent impurities and no excessive amount of sperm with head and tail abnormalities) were diluted using a Tris-egg yolk based extender, filled in straws containing between 15 and 25 million sperm, and cryoconserved in liquid nitrogen. The raw ejaculate data were filtered according to the parameters listed in Table 2. Response variables for the association study were the mean semen quality values of 902 bulls assessed in at least eight ejaculates per bull.
Sperm morphology data of Swiss BSW bulls from breeding soundness examinations Sperm morphology was provided by Swissgenetics for 1766 ejaculates of 669 BSW bulls that entered the semen collection center between August 2009 and February 2020. Technicians systematically evaluate sperm morphology of prospective breeding bulls as part of the andrological examination in formalin-fixed fresh semen samples under oil immersion light microscope. Between 150 and 200 sperm per ejaculate are classified according to three non-compensatory and 13 compensatory sperm morphology categories that are considered as either minor or major defects (Table 3). For ejaculates that evidently fail the quality assessment, only 50 spermatozoa are classified. For sperm with multiple morphological abnormalities, only the most severe defect is recorded. Bulls producing ejaculates that contain at least 75% normal sperm and not more than 20% sperm with noncompensatory defects are considered suitable for breeding. If these requirements are not fulfilled, the sperm morphology examination is repeated approximately 4 weeks later. If a bull produces ejaculates that contain between 65 and 74% normal sperm and not more than 20% sperm with non-compensatory defects, its ejaculates are processed but the reduced semen quality is compensated with more sperm per insemination dose. Bulls that repeatedly produce ejaculates with less than 65% normal sperm and more than 20% sperm with noncompensatory defects are not used for breeding.
For the association study, we considered three trait categories of the sperm morphology examination: The total number of examinations per bull (n = 595 bulls with 1577 ejaculates). The number of examined ejaculates ranged from 1 to 26 (median 2). The average proportion of sperm for 13 defect categories (Table 3) was quantified for 575 bulls based on 1165 examined ejaculates for which at least 150 spermatozoa were assessed. Between 1 and 13 (median = 1) ejaculates per bull were examined. Aggregate phenotypes were the average proportion of sperm with normal morphology (assessed spermsum of all morphological defects), major defects, minor defects, compensatory defects and noncompensatory defects (Table 3).

Haplotype-based association testing
Haplotype-based association testing was implemented in R using a sliding window-based approach that was previously applied to map complex and binary traits [22,72]. We shifted a window of 50 contiguous SNPs in steps of 15 SNPs along the autosomal haplotypes. Within each sliding window, all haplotypes with a frequency above 1% were tested for association with bull fertility or semen quality using the linear model y ¼ μ þ P 10 j¼1 a j PC j þ bHT þ ck þ e, where y is a vector of phenotypes (see above), μ is the intercept, PC j are the top 10 principal components of the genomic relationship matrix, a, b and c are effects of the principal components, the haplotype (HT) tested and a covariable k (e.g., the recessive BTA6 haplotype associated with male fertility), respectively, and e is a vector of residuals that are assumed to be normally distributed. Haplotypes were tested for association assuming either additive or recessive mode of inheritance. Recessive tests were carried out for haplotypes that were observed in the homozygous state in at least 0.5% of individuals. The genomic relationship matrix was built with plink (version 1.9) [65] using 589,791 autosomal SNPs with minor allele frequency greater than 0.5%. Its principle components were calculated with GCTA (version version 1.92.1) [69]. The genomic inflation factor lambda was calculated as: lambda = median (qchisq ( qchisq (0.5, 1) where p is a vector/ qchisq (0.5, 1) where p is a vector of P values.

Redefining a previously detected fertility-associated haplotype on BTA6
We previously detected a recessive QTL for semen quality and bull fertility on BTA6 [22] that segregates at high frequency in BSW cattle. In order to localize the QTL in our mapping cohort, we performed an association study between bull fertility and 6962 haplotypes (31,362 partially imputed SNPs) of chromosome 6 in 1130 Swiss BSW bulls using a haplotype-based association mapping approach (see above) assuming recessive inheritance. The association study revealed the strongest association (P = 2.4 × 10 − 20 ) with bull fertility for a haplotype located between 57,628,627 and 57,919,676 bp. In association analyses that were conditioned on the BTA6 QTL, we Whole-genome sequencing, sequence variant genotyping and imputation We considered 368 BSW and Original Braunvieh (OBV) cattle that were sequenced at between 3 and 60-fold (median of 11.6) genome coverage using either Illumina HiSeq 2500, Illumina HiSeq 4000, or Illumina Nova-Seq6000 instruments. The raw sequencing reads were filtered using the fastp software [73] with default parameter settings and subsequently aligned to the ARS-UCD1.2 assembly of the bovine genome [74] using the mem-algorithm of the Burrows Wheeler Aligner (BWA) software package [75]. Single nucleotide and short insertion and deletion polymorphisms were detected and genotyped using the HaplotypeCaller, GenomicsDBImport and GenotypeGVCF modules of GATK [76]. Low quality variants were filtered out using GATK's best practice recommendations for hard filtering. Following filtration, 24,614,962 SNPs and 3,981,732 indels were retained for subsequent analyses. Beagle 4.1 [77] phasing and imputation was applied to improve the primary genotype calls from GATK and infer missing genotypes. We used the Beagle 5.1 [66] software with default settings to impute the genotypes for 28,596,694 sequence variants for 3736 bulls from the mapping cohort that had partially imputed genotypes at 683,609 SNPs from a reference panel of 368 sequenced BSW and OBV animals.

Sequence-based association study
Imputed sequence variant genotypes that had allele frequency greater than 0.05 and imputation accuracy (Beagle r 2 ) greater than 0.4 were considered for association testing. The statistical model applied to test for an association between the imputed sequence variants and bull fertility was analogous to the haplotype-based association study. We used the plink (version 1.9) [65] software to fit a linear regression model assuming a recessive mode of inheritance conditional on the top ten principal components of a genomic relationship matrix (see above) and the fertility-associated haplotype on BTA6 (see above).

Prioritization of candidate causal variants
Variants that exceeded the Bonferroni-corrected significance threshold of 3.88 × 10 − 9 (defined based on 12,861, 528 imputed sequence variants with MAF > 0.05 and Beagle r 2 > 0.40) and were compatible with the inheritance of the top haplotype were considered as candidate causal variants. To identify variants compatible with the inheritance of the top haplotype, we considered sequence variant genotypes from 125 BSW cattle with bull fertility phenotypes that had both array-derived genotypes and sequence variant genotypes. Their carrier status for the top haplotype at the QTL was inferred from the array-derived genotypes. We considered sequence variants within ±3 Mb of the top haplotype as positional candidate causal variants. We filtered for variants that had allele frequency greater than 0.8 in animals that carried the top haplotype in the homozygous state, between 0.4 and 0.6 in heterozygous haplotype carriers, and less than 0.05 in non-carrier animals.
Functional consequences of candidate causal variants were predicted based on the Ensembl (release 101, August 2020) annotation of the ARS-UCD1.2 Bos taurus assembly using the Variant Effect Predictor (VEP) software tool [78].

Bioinformatic analyses
The Clustal Omega [79] software was used for multi-species alignment of protein sequences and Missense3D [80] to predict putative effects of missense variants on protein tertiary structure. The three dimensional protein structure was modelled using the SWISS-MODEL platform [81].

Transcriptome analysis
To assess the abundance of transcripts in bovine testis, we downloaded between 47 and 58 million 2 × 150 bp pairedend sequencing reads from the European Nucleotide Archive (ENA, https://www.ebi.ac.uk/ena) that were generated using RNA extracted from testis tissue samples of three mature bulls and three newborn male calves of the Angus beef cattle breed (ENA accession numbers: SAMN09205186-SAMN09205191 [82];). The RNA sequencing reads were pseudo-aligned to an index of the bovine transcriptome (ftp://ftp.ensembl.org/pub/release-98/fasta/bos_taurus/cdna/Bos_taurus.ARS-UCD1.2.cdna. all.fa.gz) and transcript abundance was quantified using the kallisto software [83]. We used the R package tximport [84] to aggregate transcript abundances to the gene level.

Imaging of spermatozoa from bulls homozygous for the haplotype on BTA1
Fresh semen samples were provided by Swissgenetics for two bulls that were homozygous carriers of the fertilityassociated haplotype on BTA1. The haplotype status of the bulls was derived using microarray-derived genotypes. Sperm morphology was assessed in semen samples fixed in buffered formaldehyde saline solution using oil-immersion phase contrast microscopy (Olympus BX50, UplanF1 100×/1.30, Olympus, Wallisellen, Switzerland).