Advertisement

Identification of genomic regions associated with total and progressive sperm motility in Italian Holstein bulls

Open AccessPublished:November 15, 2022DOI:https://doi.org/10.3168/jds.2021-21700

      ABSTRACT

      Sperm motility is directly related to the ability of sperm to move through the female reproductive tract to reach the ovum. Sperm motility is a complex trait that is influenced by environmental and genetic factors and is associated with male fertility, oocyte penetration rate, and reproductive success of cattle. In this study we carried out a GWAS in Italian Holstein bulls to identify candidate regions and genes associated with variations in progressive and total motility (PM and TM, respectively). After quality control, the final data set consisted of 5,960 records from 949 bulls having semen collected in 10 artificial insemination stations and genotyped at 412,737 SNPs (call rate >95%; minor allele frequency >5%). (Co)variance components were estimated using single trait mixed models, and associations between SNPs and phenotypes were assessed using a genomic BLUP approach. Ten windows that explained the greatest percentage of genetic variance were located on Bos taurus autosomes 1, 2, 4, 6, 7, 23, and 26 for TM and Bos taurus autosomes 1, 2, 4, 6, 8, 16, 23, and 26 for PM. A total of 150 genes for TM and 72 genes for PM were identified within these genomic regions. Gene Ontology enrichment analyses identified significant Gene Ontology terms involved with energy homeostasis, membrane functions, sperm-egg interactions, protection against oxidative stress, olfactory receptors, and immune system. There was significant enrichment of quantitative trait loci for fertility, calving ease, immune response, feed intake, and carcass weight within the candidate windows. These results contribute to understanding the architecture of the genetic control of sperm motility and may aid in the development of strategies to identify subfertile bulls and improve reproductive success.

      Key words

      INTRODUCTION

      In 2020, the world market for AI was valued ∼$3.96 billion with an expected annual growth of ∼5.94% by 2028. Among other factors, the quantity, quality, and sanitary status of the semen will determine the reproductive success of a bull. Although fertility of semen is highly influenced by age, season, nutrition, management, and the interval between ejaculations (
      • Fuerst-Waltl B.
      • Schwarzenbacher H.
      • Perner C.
      • Sölkner J.
      Effects of age and environmental factors on semen production and semen quality of Austrian Simmental bulls.
      ;
      • Suchocki T.
      • Szyda J.
      Genome-wide association study for semen production traits in Holstein-Friesian bulls.
      ), it is also known that genetics and biological processes involved in sperm production and sperm-ovum interactions play a significant role (
      • Druet T.
      • Fritz S.
      • Sellem E.
      • Basso B.
      • Gérard O.
      • Salas-Cortes L.
      • Humblot P.
      • Druart X.
      • Eggen A.
      Estimation of genetic parameters and genome scan for 15 semen characteristics traits of Holstein bulls.
      ;
      • Yin H.
      • Zhou C.
      • Shi S.
      • Fang L.
      • Liu J.
      • Sun D.
      • Jiang L.
      • Zhang S.
      Weighted single-step genome-wide association study of semen traits in Holstein bulls of China.
      ;
      • Sweett H.
      • Fonseca P.A.S.
      • Suárez-Vega A.
      • Livernois A.
      • Miglior F.
      • Cánovas A.
      Genome-wide association study to identify genomic regions and positional candidate genes associated with male fertility in beef cattle.
      ).
      Sperm motility (SM) affects the ability of the sperm to reach and penetrate the zona pellucida, fertilize the ovum, and initiate embryogenesis (
      • Alves M.B.R.
      • Celeghini E.C.C.
      • Belleannée C.
      From sperm motility to sperm-borne microRNA signatures: New approaches to predict male fertility potential.
      ). Sperm motility is routinely measured by insemination centers using computer-assisted sperm analysis. However, it has been shown that SM is a complex phenotype which is affected by the environmental (
      • Druet T.
      • Fritz S.
      • Sellem E.
      • Basso B.
      • Gérard O.
      • Salas-Cortes L.
      • Humblot P.
      • Druart X.
      • Eggen A.
      Estimation of genetic parameters and genome scan for 15 semen characteristics traits of Holstein bulls.
      ) and both additive and nonadditive genetic effects (
      • VanRaden P.M.
      • Olson K.M.
      • Null D.J.
      • Hutchison J.L.
      Harmful recessive effects on fertility detected by absence of homozygous haplotypes.
      ;
      • Peñagaricano F.
      • Weigel K.A.
      • Khatib H.
      Genome-wide association study identifies candidate markers for bull fertility in Holstein dairy cattle.
      ). Consequently, a better understanding of the genomic architecture of semen-related complex traits is expected to help improve fertility through the assisted selective breeding of bulls. The heritability estimates for SM vary widely among studies on dairy and beef cattle, from very low values of 0.01–0.13 (
      • Gredler B.
      • Fuerst C.
      • Fuerst-Waltl B.
      • Schwarzenbacher H.
      • Sölkner J.
      Genetic parameters for semen production traits in Austrian dual-purpose Simmental bulls.
      ;
      • Corbet N.J.
      • Burns B.M.
      • Johnston D.J.
      • Wolcott M.L.
      • Corbet D.H.
      • Venus B.K.
      • Li Y.
      • McGowan M.R.
      • Holroyd R.G.
      Male traits and herd reproductive capability in tropical beef cattle. 2. Genetic parameters of bull traits.
      ;
      • Suchocki T.
      • Szyda J.
      Genome-wide association study for semen production traits in Holstein-Friesian bulls.
      ) to moderate and high values of 0.43–0.60 (
      • Hering D.M.
      • Olenski K.
      • Kaminski S.
      Genome-wide association study for poor sperm motility in Holstein-Friesian bulls.
      ).
      A large number of genes, some having pleiotropic effects on different semen traits and cell types, are involved in the production, maturation, and sperm viability and can affect semen quality and fertility (
      • Ballow D.
      • Meistrich M.L.
      • Matzuk M.
      • Rajkovic A.
      Sohlh1 is essential for spermatogonial differentiation.
      ;
      • Zhang T.
      • Oatley J.
      • Bardwell V.J.
      • Zarkower D.
      DMRT1 is required for mouse spermatogonial stem cell maintenance and replenishment.
      ). Here, we identify novel genomic regions and candidate genes putatively involved in total and progressive motility (PM) in Italian Holstein bulls using a GWAS approach.

      MATERIALS AND METHODS

      Phenotypic and Genotypic Data

      We used 5,960 records of total and progressive SM from 949 Italian Holstein bulls with ages between 1 and 11 yr from 10 AI centers (Figure 1a), collected between 1995 and 2014 (Figure 1b). The number of measurements per sire varied between 3 and 78 for both total motility and PM traits. The procedures for semen collection are standardized across AI centers. Frozen-thawed sperm are routinely analyzed in the collecting centers using a CASA Integrated Visual Optical System (IVOS; CASA System, Hamilton Thorne Inc.) and a standardized protocol to assess SM. Total motility (TM) is the percentage of sperm moving irrespective of direction, and PM is the percentage of sperm that exhibit linear movement. Sperm phenotypic data were retrieved from the database of the Istituto Lazzari Spallanzani, in charge of official sperm quality controls of national bulls in Italy. The DNA extracted from semen was already available at the Biobank of Università Cattolica (see
      • Mancini G.
      • Gargani M.
      • Chillemi G.
      • Nicolazzi E.L.
      • Marsan P.A.
      • Valentini A.
      • Pariset L.
      Signatures of selection in five Italian cattle breeds detected by a 54K SNP panel.
      ).
      Figure thumbnail gr1
      Figure 1(a) Number of bulls collected at each AI center identified by anonymous codes and (b) number of bulls sampled per year of birth.
      Illumina BovineSNP50 genotypes for 636 animals were imputed to high-density (777,962) with Beagle 4.1 (
      • Browning B.L.
      • Browning S.R.
      Genotype imputation with millions of reference samples.
      ) using a reference set of 1,006 animals, genotyped with the Illumina BovineHD array. Both data sets were previously updated to the ARS-UCD1.2 (
      • Rosen B.D.
      • Bickhart D.M.
      • Schnabel R.D.
      • Koren S.
      • Elsik C.G.
      • Tseng E.
      • Rowan T.N.
      • Low W.Y.
      • Zimin A.
      • Couldrey C.
      • Hall R.
      • Li W.
      • Rhie A.
      • Ghurye J.
      • McKay S.D.
      • Thibaud-Nissen F.
      • Hoffman J.
      • Murdoch B.M.
      • Snelling W.M.
      • McDaneld T.G.
      • Hammond J.A.
      • Schwartz J.C.
      • Nandolo W.
      • Hagen D.E.
      • Dreischer C.
      • Schultheiss S.J.
      • Schroeder S.G.
      • Phillippy A.M.
      • Cole J.B.
      • Van Tassell C.P.
      • Liu G.
      • Smith T.P.L.
      • Medrano J.F.
      De novo assembly of the cattle reference genome with single-molecule sequencing.
      ) bovine assembly version using Plink 1.9 (
      • Chang C.C.
      • Chow C.C.
      • Tellier L.C.
      • Vattikuti S.
      • Purcell S.M.
      • Lee J.J.
      Second-generation PLINK: Rising to the challenge of larger and richer datasets.
      ). Quality control excluded indels, SNPs with minor allele frequency lower than 0.05 (
      • Wiggans G.R.
      • Sonstegard T.S.
      • VanRaden P.M.
      • Matukumalli L.K.
      • Schnabel R.D.
      • Taylor J.F.
      • Schenkel F.S.
      • Van Tassell C.P.
      Selection of single-nucleotide polymorphisms and quality of genotypes used in genomic evaluation of dairy cattle in the United States and Canada.
      ), and duplicated markers from both data sets. For the BovineSNP50 data set, only SNPs in common with HD data set were used. Finally, the filtered reference data set, comprising 569,630 SNPs, was phased and missing genotypes imputed through Beagle 5.2 (
      • Browning B.L.
      • Tian X.
      • Zhou Y.
      • Browning S.R.
      Fast two-stage phasing of large-scale sequence data.
      ). Low quality imputed variants, with DR2 (dosage R2, which is the estimated squared correlation between the estimated allele dose and the true allele dose), lower than 0.3 were excluded (
      • Wu D.
      • Dou J.
      • Chai X.
      • Bellis C.
      • Wilm A.
      • Shih C.C.
      • Soon W.W.J.
      • Bertin N.
      • Lin C.B.
      • Khor C.C.
      • DeGiorgio M.
      • Cheng S.
      • Bao L.
      • Karnani N.
      • Hwang W.Y.K.
      • Davila S.
      • Tan P.
      • Shabbir A.
      • Moh A.
      • Tan E.-K.
      • Foo J.N.
      • Goh L.L.
      • Leong K.P.
      • Foo R.S.Y.
      • Lam C.S.P.
      • Richards A.M.
      • Cheng C.-Y.
      • Aung T.
      • Wong T.Y.
      • Ng H.H.
      • Liu J.
      • Wang C.
      • Ackers-Johnson M.A.
      • Aliwarga E.
      • Ban K.H.K.
      • Bertrand D.
      • Chambers J.C.
      • Chan D.L.H.
      • Chan C.X.L.
      • Chee M.L.
      • Chee M.L.
      • Chen P.
      • Chen Y.
      • Chew E.G.Y.
      • Chew W.J.
      • Chiam L.H.Y.
      • Chong J.P.C.
      • Chua I.
      • Cook S.A.
      • Dai W.
      • Dorajoo R.
      • Foo C.-S.
      • Goh R.S.M.
      • Hillmer A.M.
      • Irwan I.D.
      • Jaufeerally F.
      • Javed A.
      • Jeyakani J.
      • Koh J.T.H.
      • Koh J.Y.
      • Krishnaswamy P.
      • Kuan J.L.
      • Kumari N.
      • Lee A.S.
      • Lee S.E.
      • Lee S.
      • Lee Y.L.
      • Leong S.T.
      • Li Z.
      • Li P.Y.
      • Liew J.X.
      • Liew O.W.
      • Lim S.C.
      • Lim W.K.
      • Lim C.W.
      • Lim T.B.
      • Lim C.K.
      • Loh S.Y.
      • Lok A.W.
      • Chin C.W.L.
      • Majithia S.
      • Maurer-Stroh S.
      • Meah W.Y.
      • Mok S.Q.
      • Nargarajan N.
      • Ng P.
      • Ng S.B.
      • Ng Z.
      • Ng J.Y.X.
      • Ng E.
      • Ng S.L.
      • Nusinovici S.
      • Ong C.T.
      • Pan B.
      • Pedergnana V.
      • Poh S.
      • Prabhakar S.
      • Prakash K.M.
      • Quek I.
      • Sabanayagam C.
      • See W.Q.
      • Sia Y.Y.
      • Sim X.
      • Sim W.C.
      • So J.
      • Soon D.K.N.
      • Tai E.S.
      • Tan N.Y.
      • Tan L.C.S.
      • Tan H.C.
      • Tan W.L.W.
      • Tandiono M.
      • Tay A.
      • Thakur S.
      • Tham Y.C.
      • Tiang Z.
      • Toh G.L.-X.
      • Tsai P.K.
      • Veeravalli L.
      • Verma C.S.
      • Wang L.
      • Wang M.R.
      • Wong W.-C.
      • Xie Z.
      • Yeo K.K.
      • Zhang L.
      • Zhai W.
      • Zhao Y.
      Large-scale whole-genome sequencing of three diverse Asian populations in Singapore.
      ;
      • Chen L.
      • Abel H.J.
      • Das I.
      • Larson D.E.
      • Ganel L.
      • Kanchi K.L.
      • Regier A.A.
      • Young E.P.
      • Kang C.J.
      • Scott A.J.
      • Chiang C.
      • Wang X.
      • Lu S.
      • Christ R.
      • Service S.K.
      • Chiang C.W.K.
      • Havulinna A.S.
      • Kuusisto J.
      • Boehnke M.
      • Laakso M.
      • Palotie A.
      • Ripatti S.
      • Freimer N.B.
      • Locke A.E.
      • Stitziel N.O.
      • Hall I.M.
      Association of structural variation with cardiometabolic traits in Finns.
      ), giving a final imputed data set of 417,256 SNPs. As 91 animals were genotyped both with 50k and HD panels, these were removed from the imputed data set and their HD data were considered for the downstream analysis. Therefore, the GWAS analysis was performed using 949 animals with phenotype and genotype information (404 animals genotyped with BovineHD and 545 animals imputed from BovineSNP50). A further quality control was carried out to exclude SNPs with minor allele frequency <0.05. After quality control, 412,737 SNPs and 949 animals were used for further analyses.

      Genome-Wide Association Study

      The genomic BLUP using individuals with repeated measurements were analyzed with a repeatability single trait animal model:
      y = + Zu + Wpe + e,


      where yis a vector of observed phenotypes (TM or PM) for all sires, βis the vector of fixed effects (AI center, contemporary group: year-season, quadratic, and linear effects of age), u is the random vector of additive genetic effects, peis the random vector of permanent environmental effects, and e is the random vector of residual effects. The matrices X, Z, and W are the incidence matrices of fixed effects, additive genetic effects, and permanent environmental effects, respectively. Random effects were assumed to follow a multivariate normal distribution, as follows:
      (upee|σu2,σpe2,σe2)&sim;N[0,(Gσu2000Iσpe2000Iσe2)],


      where σu2, σpe2, and σe2 are the additive genetic variance, the permanent environment variance, and the residual variance, respectively; I is an identity matrix. G is a genomic relationship matrix derived from marker information (
      • VanRaden P.M.
      Efficient methods to compute genomic predictions.
      ). The heritability was estimated as the ratio of additive genetic variance to the total phenotypic variance (σu2 + σpe2 + σe2) and the repeatability was calculated as sum of individual variance (σu2 + σpe2) divided by phenotypic variance.
      The percentage of additive genetic variance explained by the ith SNP window was then calculated as
      Var(ui)σu2×100=Var(j=1xZjuj)σu2×100,


      where ui is the genetic value of the ith genomic region under consideration, σu2 is the total additive genetic variance, x is the total number of adjacent SNPs within the 1 Mb region, Zj is the vector of gene content of the jth SNP for all windows, and uj is the effect of the jth SNP within ith window. GWAS results are reported as the percentage of genetic additive variance explained by windows of 1 Mb and plotted using CMplot package in R v.1.3.1073 (R core team 2020). Additionally, we calculated the P-value of individual SNP according to
      • Aguilar I.
      • Legarra A.
      • Cardoso F.
      • Masuda Y.
      • Lourenco D.
      • Misztal I.
      Frequentist p-values for large-scale-single step genome-wide association, with an application to birth weight in American Angus cattle.
      . All procedures were performed using BLUPF90 family programs (
      • Misztal I.
      • Tsuruta S.
      • Lourenco D.A.L.
      • Aguilar I.
      • Legarra A.
      • Vitezica Z.
      Manual for BLUPF90 family of programs. University of Georgia.
      ). To control false positives, we applied the false discovery rate (FDR) method for multiple testing according to
      • Benjamini Y.
      • Hochberg Y.
      Controlling the false discovery rate: A practical and powerful approach to multiple testing.
      as follows:
      FDR = m × Pmax/n,


      where m is the number of times to be tested, Pmax is the genome-wide significance level empirical P-value of FDR adjusted, and n is the number of significant SNPs at the assigned FDR level (i.e., 0.05).

      Gene Annotation and Functional Enrichment

      Our results indicated that single SNPs with the smallest P-values did not explain a large percentage of genetic variance. Therefore, we used the percentage of variance explained by SNPs within consecutive 1 Mb windows to identify genomic regions of interest for further investigation. Genes within the 10 windows explaining the largest amount of variance were annotated using the GALLO v.0.99 R package (
      • Fonseca P.A.S.
      • Suárez-Vega A.
      • Marras G.
      • Cánovas Á.
      GALLO: An R package for genomic annotation and integration of multiple data sources in livestock for positional candidate loci.
      ) and the ENSEMBL annotation of the Bos taurus ARS-UCD 1.2 genome assembly ( http://bovinegenome.org ). The genes identified from ENSEMBL were collected to perform a functional annotation by DAVID v6.8 ( https://david.ncifcrf.gov/home.jsp ) and REVIGO (
      • Supek F.
      • Bošnjak M.
      • Škunca N.
      • Šmuc T.
      REVIGO summarizes and visualizes long lists of gene ontology terms.
      ) using Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway ( https://www.genome.jp/kegg ) databases for B. taurus. The annotation and enrichment analyses of QTL in adjacent windows (1 Mb) were assessed for the effects of neighboring SNPs, which can be captured due to linkage disequilibrium. The QTL enrichment was evaluated for all the annotated QTL and the P-values were calculated and corrected for multiple testing, using a value of FDR <0.05 to identify the enriched QTL. The annotation and enrichment analyses were performed using the GALLO v.0.99. package and the Animal QTLdb for the ARS-UCD 1.2 bovine genome ( http://bovinegenome.org ).

      RESULTS

      Descriptive Statistics and Genetic Parameters

      The bulls had a mean of 55.33% (12.70 SD) for TM and 43.21% (11.04 SD) for PM. The fixed effect of the contemporary group and the quadratic and linear effect of age (covariable) were highly significant in the model for both traits (P < 0.001). The effect of the AI center was significant (P < 0.01) and the estimates of heritability were low for both TM (0.13 ± 0.0062) and PM (0.10 ± 0.0045). The genetic correlation between TM and PM was high (0.62 ± 0.010) and the repeatability was similar for the 2 traits: 0.31 and 0.33 for TM and PM, respectively (Table 1).
      Table 1Additive genetic (σ^a2), permanent environmental (σ^ep2), and residual (σ^e2) variance and genetic parameters heritability (h^2), genetic correlation (r^g2), and repeatability (± SE) for total motility (TM) and progressive motility (PM) in Italian Holsteins bulls
      Repeatability: r=σa2+σep2/σep2σa2σa2+σep2+σe2.
      Traitσ^a2σ^ep2σ^e2h^2r^g2r
      TM14.96 ± 2.0220.85 ± 2.2778.80 ± 1.130.13 ± 0.00620.62 ± 0.0100.31 ± 0.0027
      PM13.05 ± 2.2824.36 ± 2.5182.43 ± 1.170.10 ± 0.00450.33 ± 0.0019
      1 Repeatability: r=σa2+σep2/σep2σa2σa2+σep2+σe2.

      QTL Regions and Candidate Genes

      As previously noted, our results indicated that the SNPs with the smallest P-values did not explain a large percentage of genetic variance (Figure 23), contrary to expectation. P-values are influenced by sample size (
      • Greenland S.
      • Senn S.J.
      • Rothman K.J.
      • Carlin J.B.
      • Poole C.
      • Goodman S.N.
      • Altman D.G.
      Statistical tests, p values, confidence intervals, and power: A guide to misinterpretations.
      ) and, in particular, by allele frequencies (
      • Aguilar I.
      • Legarra A.
      • Cardoso F.
      • Masuda Y.
      • Lourenco D.
      • Misztal I.
      Frequentist p-values for large-scale-single step genome-wide association, with an application to birth weight in American Angus cattle.
      ).
      Figure thumbnail gr2
      Figure 2Manhattan plot displaying the association of all SNP for progressive motility (PM; P-values of individual SNP effects) and quantile-quantile plots (right), which show the observed distribution of P-values against the expected P-values under the null hypothesis of no association. The horizontal axis shows −log10 transformed expected P-values and the vertical axis indicates −log10 transformed observed P-values.
      Figure thumbnail gr3
      Figure 3Manhattan plot displaying the association of all SNP for progressive motility (PM; P-values of individual SNP effects) and quantile-quantile plots (right), which show the observed distribution of P-values against the expected P-values under the null hypothesis of no association. The horizontal axis shows −log10 transformed expected P-values and the vertical axis indicates −log10 transformed observed P-values.
      The 10 windows that explained the greatest percentage of genetic variance were located on BTA 1, 2, 4, 6, 7, 23, and 26 for TM and BTA 1, 2, 4, 6, 8, 16, 23, and 26 for PM (Figure 4). These regions explained 6.91% and 6.07% of the additive genetic variance for TM and PM, respectively (Table 2). A total of 150 genes for TM and 107 genes for PM were annotated by GALLO within these significant windows (Table 2). Among them, 130 genes for TM and 90 genes for PM were found in the Bos taurus DAVID database.
      Figure thumbnail gr4
      Figure 4Circular Manhattan plot for the proportion of genetic variance explained by nonoverlapping windows of 1 Mb for (a) total motility (TM) and (b) progressive motility (PM) in Italian Holstein Bulls. The scale bar represents the percentage of variance explained by each window. Chr = chromosome.
      Table 2Genomic regions and genes in the 10 windows explaining the highest percentage of genetic variance for total motility (TM) and progressive motility (PM) in Italian Holstein Bulls
      Chr
      Chr = chromosome.
      Window regionTM variance (%)PM variance (%)Gene
      StartEnd
      165,216,44966,211,1800.55
      Regions not identified for the trait; agenes associated only with TM, bgenes associated only PM. *Regions not identified for the trait; agenes associated only with TM, cgenes associated with TM and PM.
      ENSBTAG00000048371, FBXO40, GOLGB1, GTF2E1, HCLS1, HGD, NDUFB4, POLQ, RABL3, STXBP5L
      68,304,71369,303,1890.48
      Regions not identified for the trait; agenes associated only with TM, bgenes associated only PM. *Regions not identified for the trait; agenes associated only with TM, cgenes associated with TM and PM.
      CCDC14, ENSBTAG00000053239, ITGB5, KALRN, ROPN1, UMPS
      102,624,551103,622,4541.291.62E1BHJ0
      148,523,433149,497,3740.36
      Regions not identified for the trait; agenes associated only with TM, bgenes associated only PM. *Regions not identified for the trait; agenes associated only with TM, cgenes associated with TM and PM.
      CHAF1B, CLDN14, DOP1B, HLCS, MORC3, PIGP, RIPPLY3, SIM2, TTC3, VPS26C
      2134,981,487135,980,9970.500.43ARHGEF10L, ARHGEF19, ATP13A2, CPLANE2, CROCC, EPHA2, FBXO42, MFAP2, NECAP2, PADI1, PADI2, PADI3, PADI4, PADI6, RCC2, SDHB, SPATA21, SZRD1
      469,538,51870,533,4610.78CBX3, HNRNPA2B1, NFE2L3, SNX10
      90,956,17191,949,0790.390.44ARF5, FSCN3, GCC1, GRM8, PAX4, SND1, ZNF800
      637,145,52238,144,1950.38
      Regions not identified for the trait; agenes associated only with TM, bgenes associated only PM. *Regions not identified for the trait; agenes associated only with TM, cgenes associated with TM and PM.
      DCAF16, FAM184B, LAP3, LCORL, MED28, NCAPG
      102,406,084103,403,1930.410.42CRMP1,a DMP1,c DSPP,c ENSBTAG00000006508b ENSBTAG00000006509a ENSBTAG00000017829b ENSBTAG00000017830a ENSBTAG00000021730a ENSBTAG00000054995c JAKMIP1,a PPP2R2C,c SPARCL1,c WFS1c
      744,213,44845,212,0960.44
      Regions not identified for the trait; agenes associated only with TM, bgenes associated only PM. *Regions not identified for the trait; agenes associated only with TM, cgenes associated with TM and PM.
      AFF4, BTBD2, CSNK1G2, FSTL4, GDF9, HSPA4, LEAP2, SCAMP4, SHROOM1, UQCRQ, ZCCHC10
      854,914,79855,908,878
      Regions not identified for the trait; agenes associated only with TM, bgenes associated only PM. *Regions not identified for the trait; agenes associated only with TM, cgenes associated with TM and PM.
      0.52TLE4
      1670,458,02471,454,118
      Regions not identified for the trait; agenes associated only with TM, bgenes associated only PM. *Regions not identified for the trait; agenes associated only with TM, cgenes associated with TM and PM.
      0.49ANGEL2, ATF3, BATF3, DTL, FLVCR1, NENF, NSL1, PACC1, PPP2R5A, RPS6KC1, SPATA45, TATDN3, VASH2
      23336,7481,330,072
      Regions not identified for the trait; agenes associated only with TM, bgenes associated only PM. *Regions not identified for the trait; agenes associated only with TM, cgenes associated with TM and PM.
      0.49KHDRBS2
      24,376,52925,372,4860.83
      Regions not identified for the trait; agenes associated only with TM, bgenes associated only PM. *Regions not identified for the trait; agenes associated only with TM, cgenes associated with TM and PM.
      CILK1, EFHC1, ENSBTAG00000025486 ENSBTAG00000055193, FBXO9, GCM1, GSTA1, GSTA2, GSTA3, GSTA4, GSTA5, IL17A, IL17F, MCM3, PAQR8, PKHD1, TMEM14A, TRAM2
      27,974,48828,973,5100.74
      Regions not identified for the trait; agenes associated only with TM, bgenes associated only PM. *Regions not identified for the trait; agenes associated only with TM, cgenes associated with TM and PM.
      ABCF1, ATAT1, BOLA, BOLA-NC1, C23H6orf136, C23H6orf15, CCHCR1, CDSN, DDR1, DHX16 ENSBTAG00000004491 ENSBTAG00000005147 ENSBTAG00000007076 ENSBTAG00000037422 ENSBTAG00000038620 ENSBTAG00000054589 FLOT1, GNL1, GTF2H4, IER3, JSP.1, MDC1, MOG, MRPS18B, NRM, POLR1H, POU5F1, PPP1R10, PPP1R11, PPP1R18, PRR3, PSORS1C2, RNF39, RPP21, SFTA2, TCF19, TRIM10, TRIM15, TRIM26, TRIM31, TRIM39, TRIM40, TUBB, VARS2, ZFP57
      28,014,88229,014,696
      Regions not identified for the trait; agenes associated only with TM, bgenes associated only PM. *Regions not identified for the trait; agenes associated only with TM, cgenes associated with TM and PM.
      0.56ABCF1, ATAT1, BOLA BOLA-NC1, C23H6orf136, C23H6orf15, CDSN, DDR1, DHX16 ENSBTAG00000005147 ENSBTAG00000007075 ENSBTAG00000017240 ENSBTAG00000037422 ENSBTAG00000038619 ENSBTAG00000054440 ENSBTAG00000054589 FLOT1, GABBR1, GNL1, GTF2H4, IER3, JSP.1, MDC1, MOG, MRPS18B, NRM, POLR1H, PPP1R10, PPP1R11, PPP1R18, PRR3, RNF39, RPP21, SFTA2, TRIM10, TRIM15, TRIM26, TRIM31, TRIM39, TRIM40, TUBB OR5V1, OR2H1D, OR2H1, OR10C1, OR11W11, VARS2, ZFP57
      2646,055,58547,048,3420.520.98DOCK1, INSYN2
      1 Chr = chromosome.
      * Regions not identified for the trait; agenes associated only with TM, bgenes associated only PM.*Regions not identified for the trait; agenes associated only with TM, cgenes associated with TM and PM.
      The DAVID annotations and their relative P-values were grouped using REVIGO into biological processes (BP), cellular components (CC), and molecular functions (MF), each comprising 13, 23, and 45 genes for TM and 24, 18, and 61 for PM (Table 3). A total of 42 genes are in common between TM and PM.
      Table 3Gene Ontology terms significantly associated with biological processes (BP), cellular components (CC), and molecular functions (MF) related to total motility (TM) and progressive motility (PM)
      CategoryGroupGO: TermGene nameBTATM FDR
      FDR = false discovery rate.
      PM FDR
      Biological processAlpha-AA biosynthetic processGO: 0036414 Histone citrullinationPADI1, PADI2 PADI6, PADI4 PADI329.74 × 10−91.5 × 10−6
      GO: 0018101 Protein citrullination9.74 × 10−91.5 × 10−7
      GO: 0018195 Peptidyl-arginine modification2.35 × 10−43.8 ×10−4
      GO: 0006749 Glutathione metabolic processGSTA3, GSTA4 GSTA2, GSTA1 GSTA5232.23 × 10−3
      Genome Ontology not identified for the trait.
      Antigen processing and presentationGO: 0002428 antigen processing and presentationJSP.1, BOLA BOLA-NC1230.00160.0046
      Immune system processGO: 0002376 immune system processEPHA2, PADI2 SNX102
      Genome Ontology not identified for the trait.
      0.0192
      BOLA, BOLA-NC1, FLOT1, JSP.1, MOG RNF39, TRIM10 TRIM15, TRIM31 TRIM4023
      DOCK126
      Sensory perception of chemical stimulusGO: 0003008 system processOR12D2, OR10C1, OR2H1, OR5V1, OR12D323
      Genome Ontology not identified for the trait.
      4.8 × 10−10
      GO: 0007606 sensory perception of chemical stimulus
      G protein-coupled receptor signaling pathwayGO: 0051606 detection of stimulusOR12D2, OR12D323
      Genome Ontology not identified for the trait.
      5.1 × 10−6
      GO: 0007186 G protein-coupled receptor signaling pathwayGABBR1, OR12D2, OR10C1, OR2H1 OR5V1, OR12D323
      Genome Ontology not identified for the trait.
      1.3 × 10−8
      Cellular componentsMicrotubule cytoskeletonGO: 0015630 microtubule cytoskeletonCCHCR1, CILK1 EFHC1, MCM3, PKHD1230.0037
      Genome Ontology not identified for the trait.
      CBX3, SNX104
      Genome Ontology not identified for the trait.
      0.0058
      EVC6
      PPP2R5A16
      CROCC20.00370.0058
      CRMP1, JAKMIP16
      ATAT1, FLOT1, TUBB23
      GO: 0072686 mitotic spindleATAT1, EFHC1 PKHD1230.0002
      Genome Ontology not identified for the trait.
      GO: 0005856 cytoskeletonCCHCR1 EFHC1, MCM3, PKHD1230.0320
      Genome Ontology not identified for the trait.
      CBX34
      Genome Ontology not identified for the trait.
      0.0310
      EVC6
      PPP2R5A16
      CPLANE2 CROCC PADI6, RCC220.03200.0310
      FSCN34
      CRMP1, JAKMIP16
      ATAT1, FLOT1 TUBB23
      Protein serine/threonine phosphatase complexGO: 1903293 phosphatase complex GO: 0008287 protein serine/threonine phosphatase complexPPP1R10, PPP1R11 PPP2R2C23
      Genome Ontology not identified for the trait.
      0.0045
      Focal adhesionGO: 0005925 focal adhesionEPHA2, NECAP2 ATP13A220.0713
      Genome Ontology not identified for the trait.
      MDC1, ATAT1 FLOT123
      Molecular FunctionsHydrolase activity acting on carbon-nitrogen (but not peptide) bonsGO: 0004668 protein-arginine deiminase activityPADI2 PADI1 PADI6 PADI4 PADI323.6 × 10−102.91 × 10−9
      GO: 0016810 hydrolase activity acting on carbon-nitrogen (but not peptide) bons3.2 × 10−51.2 × 10−5
      BindingGO: 0043167 ion bindingFBXO40, GTF2E1, HGD, KALRN, MORC111.8 × 10−4
      Genome Ontology not identified for the trait.
      LAP36
      CSNK1G2, FSTL4, HSPA47
      CILK1, EFHC1 GCM1, MCM3 TCF1923
      GO: 0005488 bindingCROCC, MFAP2 RCC22
      Genome Ontology not identified for the trait.
      0.00098
      CBX3, FSCN3 HNRNPA2B1 NFE2L3 PAX4, SND1 SNX104
      CRMP1, DSPP JAKMIP1, WFS16
      TLE48
      PPP2R5A, RPS6KC116
      CDSN, FLOT1 GTF2H4 KHDRBS2 MDC1, MOG NRM, PPP1R11 PPP1R18, RPP2123
      GO: 0043167 ion binding GO: 0005488 bindingATP13A2 CPLANE2, EPHA2 PADI1, PADI2 PADI3, PADI4, PADI6, SDHB, SPATA2121.8 × 10−40.00098
      ARF54
      SPARCL16
      ABCF1, DDR1, DHX16 GNL1, POLR1H PPP1R10, PRR3, RNF39, TRIM10, TRIM15, TRIM31 TRIM40, TUBB, VARS2, ZFP5723
      Glutathione transferase activityGO: 0004364 glutathione transferase activityGSTA5, GSTA3, GSTA4, GSTA2 GSTA3232.3 × 10−8
      Genome Ontology not identified for the trait.
      Olfactory reception activityGO: 0004984 olfactory reception activityOR12D12 OR10C1 OR2H1 OR5V1 OR12D323
      Genome Ontology not identified for the trait.
      9.1 × 10−11
      Molecular transducer activityGO: 0060089 molecular transducer activityEPHA22
      Genome Ontology not identified for the trait.
      6.8 × 10−7
      OR12D12 OR10C1 OR2H1 OR5V1 OR12D323
      1 FDR = false discovery rate.
      * Genome Ontology not identified for the trait.
      The BoLA and JSP.1 genes on BTA 23 were associated with both PM and TM, whereas the OR (olfactory receptor) family genes were only related to PM. The gutathione S-transferase alpha genes, involved with the metabolism of glutathione, which is an antioxidant, were only associated with TM.
      Some QTL within the genomic windows explaining the greatest percentage of genetic variance were associated with reproductive traits: these included BTA 1 and 6 for TM, and BTA 1 and 23 for PM. The enrichment analysis of the QTL (FDR-corrected P < 0.05) for TM identified 28 QTL related to the ease of calving and 6 QTL associated with the immune response on BTA 6 and BTA 23. In addition, QTL involved with feed intake (7), carcass weight (24), and fat thickness (10) were associated with TM on BTA 6. For PM, 2 QTL related to reproductive efficiency and with milk yield were identified on BTA 23 (Table 4).
      Table 4Quantitative trait loci annotated and significant in the enrichment analysis for total (TM) and progressive (PM) motility
      TraitQTL typeQTL classBTAObserved number of QTLExpected number of QTLP-valueFDR
      FDR = false discovery rate.
      TMProductionAverage daily feed intake6792.75 × 10−96.48 × 10−8
      ReproductionCalving ease (direct)628372.83 × 10−111.33 × 10−9
      Meat and carcassCarcass weight6241373.01 × 10−93.53 × 10−4
      Fat thickness at the 12th rib610215.15 × 10−78.07 × 10−8
      HealthAntibody-mediated immune response234103.76 × 10−43.53 × 10−3
      Cell-mediated immune response23234.41 × 10−43.45 × 10−3
      PMMilkMilk yield, and percentage of protein and fat in milk (EBV)1170.006740.00549
      ReproductionReproductive efficiency23170.010970.00974
      1 FDR = false discovery rate.

      DISCUSSION

      Genetic Parameters

      The heritability for TM and PM estimated in this study was low (TM: 0.14 ± 0.0062 and PM: 0.10 ± 0.0045) but comparable to estimates obtained by others using the pedigree relationship matrix A (0.13,
      • Druet T.
      • Fritz S.
      • Sellem E.
      • Basso B.
      • Gérard O.
      • Salas-Cortes L.
      • Humblot P.
      • Druart X.
      • Eggen A.
      Estimation of genetic parameters and genome scan for 15 semen characteristics traits of Holstein bulls.
      ; 0.12,
      • Yin H.
      • Fang L.
      • Qin C.
      • Zhang S.
      Estimation of the genetic parameters for semen traits in Chinese Holstein bulls.
      ) for PM. High heritability estimates, of up to 0.43 have been reported for motility in Holstein bulls (
      • Al-Kanaan A.
      • König S.
      • Brügemann K.
      Effects of heat stress on semen characteristics of Holstein bulls estimated on a continuous phenotypic and genetic scale.
      ;
      • Butler M.L.
      • Bormann J.M.
      • Weaber R.L.
      • Grieger D.M.
      • Rolf M.M.
      Selection for bull fertility: A review.
      ); however, estimates for PM in this breed are low (0.13–0.15,
      • Butler M.L.
      • Bormann J.M.
      • Weaber R.L.
      • Grieger D.M.
      • Rolf M.M.
      Selection for bull fertility: A review.
      ). The differences in heritability estimates may be the result of environmental and management differences that we found to have a significant effect on these traits. In our study estimated repeatability was moderate (0.31 and 0.33 for TM and PM, respectively) which suggests an effect of the permanent environment and an even higher effect of the residual variance.

      Candidate Genes and QTL Annotation

      In the present study, no loci with significant P-values were found using a single SNP analysis. The SNPs with the smallest P-values would be expected to explain a relatively higher proportion of genetic variation. However, this is not always the case, as P-values are affected by sample size, the magnitude of the effect detected (
      • Greenland S.
      • Senn S.J.
      • Rothman K.J.
      • Carlin J.B.
      • Poole C.
      • Goodman S.N.
      • Altman D.G.
      Statistical tests, p values, confidence intervals, and power: A guide to misinterpretations.
      ), and in particular, by allele frequencies (
      • Aguilar I.
      • Legarra A.
      • Cardoso F.
      • Masuda Y.
      • Lourenco D.
      • Misztal I.
      Frequentist p-values for large-scale-single step genome-wide association, with an application to birth weight in American Angus cattle.
      ).
      • Liu H.
      • Song H.
      • Jiang Y.
      • Jiang Y.
      • Zhang F.
      • Liu Y.
      • Shi Y.
      • Ding X.
      • Wang C.
      A single-step genome wide association study on body size traits using imputation-based whole-genome sequence data in Yorkshire pigs.
      , in a single-step GWAS in swine, found no overlap between the first 20 SNPs with lowest P-values and those explaining the largest variance. These authors attribute this discrepancy to linkage disequilibrium between SNPs and the causal variants and interaction between SNPs. Therefore, we used the variance explained by sliding window intervals to identify genomic regions involved in sperm motility.
      The ten 1-Mb windows with highest additive genetic variance for TM and PM explained, respectively, 6.91% and 6.07% of the additive genetic variance, indicating the highly polygenic nature of these traits. Other studies have also found that several regions across the genome each make small contributions to the total genetic variation for SM. For example,
      • Yin H.
      • Zhou C.
      • Shi S.
      • Fang L.
      • Liu J.
      • Sun D.
      • Jiang L.
      • Zhang S.
      Weighted single-step genome-wide association study of semen traits in Holstein bulls of China.
      found that the most significant regions, located on BTA 1, 3, 4, 5, 11, 17, and 19, explained less than 0.1% of the genetic variance for SM in Holstein bulls. In beef cattle,
      • Sweett H.
      • Fonseca P.A.S.
      • Suárez-Vega A.
      • Livernois A.
      • Miglior F.
      • Cánovas A.
      Genome-wide association study to identify genomic regions and positional candidate genes associated with male fertility in beef cattle.
      identified 10 windows on BTA 9, 13, 20, and 24 that explained 7.17% of the genetic variance of this trait. There is little consistency in the regions identified across studies, even for signals on the same chromosomes;
      • Rezende F.M.
      • Dietsch G.O.
      • Peñagaricano F.
      Genetic dissection of bull fertility in US Jersey dairy cattle.
      identified genes associated with SM on BTA 5, 11, 22, and 25 in Jersey cattle, with no regions shared with other studies.
      The window with highest additive genetic variance we identified was on BTA1 (102,482,169–103,484,977), which explained 1.29% of the genetic variance for TM and 1.62% for PM. This window included the gene E1BHJ0, which codes for profilin. Profilins are ubiquitous proteins present in mammals that regulate actin polymerization and are essential for cell viability (
      • Witke W.
      • Sutherland J.D.
      • Sharpe A.
      • Arai M.
      • Kwiatkowski D.J.
      Profilin I is essential for cell survival and cell division in early mouse development.
      ), fertilization, and normal embryo development in cattle (
      • Rawe V.Y.
      • Payne C.
      • Schatten G.
      Profilin and actin-related proteins regulate microfilament dynamics during early mammalian embryogenesis.
      ). Four profilin isoforms are expressed in mammalian testis, including spermatogenic cells (
      • Behnen M.
      • Murk K.
      • Kursula P.
      • Cappallo-Obermann H.
      • Rothkegel M.
      • Kierszenbaum A.L.
      • Kirchhoff C.
      Testis-expressed profilins 3 and 4 show distinct functional characteristics and localize in the acroplaxome-manchette complex in spermatids.
      ). Profilin-1 and profilin-3 play an important role in mouse sperm metabolism and motility (
      • Vicens A.
      • Lüke L.
      • Roldan E.R.S.
      Proteins involved in motility and sperm-egg interaction evolve more rapidly in mouse spermatozoa.
      ) and the loss of profilin-3 has been associated with impairment of spermiogenesis (
      • Umer N.
      • Arevalo L.
      • Kosterin I.
      • Phadke S.
      • Lohanadan K.
      • Kirfel G.
      • Sons D.
      • Witke W.
      • Schorle H.
      Loss of profilin3 impairs spermiogenesis by affecting acrosome biogenesis, autophagy, manchette development, and mitochondrial organization.
      ). The second most important region (0.98% variance) for PM was on BTA26 (46,055,585–47,048,342) and included cytokinesis dedicator 1 (DOCK1) which encodes a member of the migration and invasion dedicator of the cytokinesis protein superfamily which is involved in cell motility (
      • Valente T.S.
      • Baldi F.
      • Sant’Anna A.C.
      • Albuquerque L.G.
      • Paranhos da Costa M.J.R.
      Genome-wide association study between single nucleotide polymorphisms and flight speed in Nellore cattle.
      ). With regard to TM, the second most important region (0.83% variance) was on BTA 23 (24,376,529–25,372,486) and included IL17A and IL17F which have an immunoregulatory role in the male reproductive tract (
      • Bagri P.
      • Anipindi V.C.
      • Kaushic C.
      The role of IL-17 during infections in the female reproductive tract.
      ). Many of the genes in the other significant windows that fell into BP, CC, and MF categories (Table 3) are associated with spermatogenesis and male fertility [e.g., the peptidylarginine deiminase family (PAD), on BTA 2]. The PAD are a class of Ca+2 dependent enzymes that are widely expressed in mammalian tissues including ovules, ovaries, testes, and in early embryos (
      • Wang S.
      • Wang Y.
      Peptidylarginine deiminases in citrullination, gene regulation, health, and pathogenesis.
      ). The PAD enzymes perform the citrullination of arginine residues, which affects the stability and degradation of proteins, and play a crucial role during cell development, differentiation, and fertility (
      • Fichtner T.
      • Kotarski F.
      • Gärtner U.
      • Conejeros I.
      • Hermosilla C.
      • Wrenzycki C.
      • Taubert A.
      Bovine sperm samples induce different NET phenotypes in a NADPH oxidase-, PAD4-, and Ca++-dependent process.
      ). The enzyme PADI6 is involved in decondensation of sperm chromatin (
      • Wright P.W.
      • Bolling L.C.
      • Calvert M.E.
      • Sarmento O.F.
      • Berkeley E.V.
      • Shea M.C.
      • Hao Z.
      • Jayes F.C.
      • Bush L.A.
      • Shetty J.
      • Shore A.N.
      • Reddi P.P.
      • Tung K.S.
      • Samy E.
      • Allietta M.M.
      • Sherman N.E.
      • Herr J.C.
      • Coonrod S.A.
      ePAD, an oocyte and early embryo-abundant peptidylarginine deiminase-like protein that localizes to egg cytoplasmic sheets.
      ;
      • Alghamdi M.
      • Alasmari D.
      • Assiri A.
      • Mattar E.
      • Aljaddawi A.A.
      • Alattas S.G.
      • Redwan E.M.
      An overview of the intrinsic role of citrullination in autoimmune disorders.
      ). In humans, sperm chromatin organization has been shown to affect sperm function during fertilization and early embryo development (
      • Miller D.
      • Brinkworth M.
      • Iles D.
      Paternal DNA packaging in spermatozoa: More than the sum of its parts? DNA, histones, protamines and epigenetics.
      ) and has been implicated in subfertility (
      • Souza E.T.
      • Silva C.V.
      • Travençolo B.A.N.
      • Alves B.G.
      • Beletti M.E.
      Sperm chromatin alterations in fertile and subfertile bulls.
      ). The BoLAclass 1 genes and JSP.1 are associated with both PM and TM. These genes map on BTA 23 and are part of the major histocompatibility complex involved in antigen presentation (
      • Tijjani A.
      • Utsunomiya Y.T.
      • Ezekwe A.G.
      • Nashiru O.
      • Hanotte O.
      Genome sequence analysis reveals selection signatures in endangered trypanotolerant West African Muturu Cattle.
      ). The major histocompatibility complex antigens play a key role in immune response and affect production, diseases, and fertility traits, in particular the sperm-ovum interaction (
      • Ellis S.A.
      • Hammond J.A.
      The functional significance of cattle major histocompatibility complex class I genetic diversity.
      ). BoLA-NC1 is essential for the establishment and maintenance of pregnancy and in the maternal immune response to the fetus (
      • Shu L.
      • Peng X.
      • Zhang S.
      • Deng G.
      • Wu Y.
      • He M.
      • Li B.
      • Li C.
      • Zhang K.
      Non-classical major histocompatibility complex class makes a crucial contribution to reproduction in the dairy cow.
      ).
      Five genes in the most significant GO biological processes term (GO:BP), were associated exclusively with TM. Glutathione S-transferase alpha (GSTA), which is within one of the regions with the highest effect on BTA 23, belongs to the GST superfamily of multigenic and multifunctional isoenzymes that catalyze glutathione-dependent reactions (GO:BP Glutathione metabolic process and GO:MF Glutathione transferase activity). Expression of the GST has been reported in male reproductive tissues in cattle (
      • Juyena N.S.
      • Stelletta C.
      Seminal plasma: An essential attribute to spermatozoa.
      ), pigs (
      • Pérez-Patiño C.
      • Parrilla I.
      • Li J.
      • Barranco I.
      • Martínez E.A.
      • Rodriguez-Martínez H.
      • Roca J.
      The proteome of pig spermatozoa is remodeled during ejaculation.
      ), and humans (
      • Dutta S.
      • Sengupta P.
      • Slama P.
      • Roychoudhury S.
      Oxidative stress, testicular inflammatory pathways, and male reproduction.
      ) including sperm, seminal plasma, and in testicular somatic cells, particularly in Leydig and Sertoli cells (
      • Benbrahim-Tallaa L.
      • Tabone E.
      • Tosser-Klopp G.
      • Hatey F.
      • Benahmed M.
      Glutathione S-Transferase Alpha expressed in porcine sertoli cells is under the control of follicle-stimulating hormone and testosterone.
      ). GSTA also plays a role in spermatogenesis, sperm capacitation and sperm-oocyte binding, and is involved in the nuclear decondensation of the sperm chromatin (
      • Kumar R.
      • Singh V.K.
      • Atreja S.K.
      Glutathione-S-transferase: Role in buffalo (Bubalus bubalis) sperm capacitation and cryopreservation.
      ;
      • Kwon W.-S.
      • Rahman M.S.
      • Lee J.-S.
      • Kim J.
      • Yoon S.-J.
      • Park Y.-J.
      • You Y.-A.
      • Hwang S.
      • Pang M.-G.
      A comprehensive proteomic approach to identifying capacitation related proteins in boar spermatozoa.
      ;
      • Llavanera M.
      • Mateo-Otero Y.
      • Bonet S.
      • Barranco I.
      • Fernández-Fuertes B.
      • Yeste M.
      The triple role of glutathione S-transferases in mammalian male fertility.
      ).
      Genes on BTA 23 associated with TM included tripartite motifs (TRIM31 and TRIM10) and interleukins (IL17A, IL17F) that code for proteins involved in immune processes in mammals (GO:BP: antigen processing and presentation and immune system process). TRIM31 is a member of the E3 ubiquitin ligases protein family that is involved with innate immunity (
      • Ozato K.
      • Shin D.-M.
      • Chang T.-H.
      • Morse III, H.C.
      TRIM family proteins and their emerging roles in innate immunity.
      ;
      • Tijjani A.
      • Utsunomiya Y.T.
      • Ezekwe A.G.
      • Nashiru O.
      • Hanotte O.
      Genome sequence analysis reveals selection signatures in endangered trypanotolerant West African Muturu Cattle.
      ) and IL17A is a member of the interleukin family that initiate a potent inflammatory response (
      • Zenobia C.
      • Hajishengallis G.
      Basic biology and role of interleukin-17 in immunity and inflammation.
      ). In humans the level of IL17 in seminal fluid has been shown to inversely correlate with SM and mortality (
      • Qian L.
      The relationship between IL-17 and male infertility: Semen analysis.
      ) and to play an important role in successful embryo implantation. In addition, IL-17 regulates the secretion of cytokines (IL-1β and TNF-α), chemokines, and promotes the recruitment of neutrophils (
      • Song X.
      • He X.
      • Li X.
      • Qian Y.
      The roles and functional mechanisms of interleukin-17 family cytokines in mucosal immunity.
      ;
      • Wattegedera S.R.
      • Corripio-Miyar Y.
      • Pang Y.
      • Frew D.
      • McNeilly T.N.
      • Palarea-Albaladejo J.
      • McInnes C.J.
      • Hope J.C.
      • Glass E.J.
      • Entrican G.
      Enhancing the toolbox to study IL-17A in cattle and sheep.
      ) which are related to the modulation of the immune response to sperm antigens and possibly the protection of sperm from the response of the female immune system (
      • Vera O.
      • Vásqucz L.A.
      • Gladys Muñoz M.
      Semen quality and presence of cytokines in seminal fluid of bull ejaculates.
      ).
      The GO:CC terms contained the α-tubulin (ATAT1), tubulin beta class I (TUBB), and flotillin 1 (FLOT1) genes, which code for proteins that have been proposed as biomarkers of male fertility (
      • Gilbert I.
      • Bissonnette N.
      • Boissonneault G.
      • Vallée M.
      • Robert C.
      A molecular analysis of the population of mRNA in bovine spermatozoa.
      ;
      • Fu Q.
      • Pan L.
      • Huang D.
      • Wang Z.
      • Hou Z.
      • Zhang M.
      Proteomic profiles of buffalo spermatozoa and seminal plasma.
      ;
      • Patil S.K.
      • Somashekar L.
      • Selvaraju S.
      • Jamuna K.V.
      • Parthipan S.
      • Binsila B.K.
      • Prasad R.V.
      • Ravindra J.P.
      Immuno-histological mapping and functional association of seminal proteins in testis and excurrent ducts with sperm function in buffalo.
      ). Alpha-tubulin is a cytoskeletal protein that contributes to microtubule organization, sperm morphology, and flagellar function, and is present in the cytoplasm of Sertoli cells, Leydig cells, spermatogonia, primary spermatocytes, secondary spermatocytes, and spermatids of seminiferous tubules (
      • Kalebic N.
      • Sorrentino S.
      • Perlas E.
      • Bolasco G.
      • Martinez C.
      • Heppenstall P.A.
      αTAT1 is the major α-tubulin acetyltransferase in mice.
      ). Tubulin plays a role in microtubule formation in sperm cells of buffaloes (
      • Fu Q.
      • Pan L.
      • Huang D.
      • Wang Z.
      • Hou Z.
      • Zhang M.
      Proteomic profiles of buffalo spermatozoa and seminal plasma.
      ) and is expressed in spermatogonia, Sertoli cells, spermatids, and spermatozoa of pig testes (
      • Luo J.
      • Rodriguez-Sosa J.R.
      • Tang L.
      • Bondareva A.
      • Megee S.
      • Dobrinski I.
      Expression pattern of acetylated α-tubulin in porcine spermatogonia.
      ). In buffalo sperm, α-tubulin has been shown to provide structural support, facilitate cell motility, and to participate in oocyte binding and activation (
      • Patil S.K.
      • Somashekar L.
      • Selvaraju S.
      • Jamuna K.V.
      • Parthipan S.
      • Binsila B.K.
      • Prasad R.V.
      • Ravindra J.P.
      Immuno-histological mapping and functional association of seminal proteins in testis and excurrent ducts with sperm function in buffalo.
      ). JAKMIP1 and PKDH1 are related to sperm function and cell polarization (
      • Hatzirodos N.
      • Irving-Rodgers H.F.
      • Hummitzsch K.
      • Harland M.L.
      • Morris S.E.
      • Rodgers R.J.
      Transcriptome profiling of granulosa cells of bovine ovarian follicles during growth from small to large antral sizes.
      ), embryo development and survival (
      • Gigarel N.
      • Frydman N.
      • Burlet P.
      • Kerbrat V.
      • Tachdjian G.
      • Fanchin R.
      • Antignac C.
      • Frydman R.
      • Munnich A.
      • Steffann J.
      Preimplantation genetic diagnosis for autosomal recessive polycystic kidney disease.
      ), respectively. MDC1, EPHAH2, ATAT1, and FLOT1 code for proteins that are associated with focal adhesion, a process involved in fertilization, ovum transport, and sperm-ovum interaction (
      • Talbot P.
      • Shur B.D.
      • Myles D.G.
      Cell adhesion and fertilization: Steps in oocyte transport, sperm-zona pellucida interactions, and sperm-egg fusion.
      ). Additionally, PPP1R10, PPP1R11, PPP2R2C, and FSCN3 genes were identified and code for proteins related to spermatogenesis, SM, and fertility (
      • Han Y.B.
      • Feng H.L.
      • Cheung C.K.
      • Lam P.M.
      • Wang C.C.
      • Haines C.J.
      Expression of a novel T-complex testis expressed 5 (Tctex5) in mouse testis, epididymis, and spermatozoa.
      ,
      • Han Y.
      • Song X.-X.
      • Feng H.-L.
      • Cheung C.-K.
      • Lam P.-M.
      • Wang C.-C.
      • Haines C.J.
      Mutations of t-complex testis expressed gene 5 transcripts in the testis of sterile t-haplotype mutant mouse.
      ;
      • Cheng L.
      • Pilder S.
      • Nairn A.C.
      • Ramdas S.
      • Vijayaraghavan S.
      PP1γ2 and PPP1R11 are parts of a multimeric complex in developing testicular germ cells in which their steady state levels are reciprocally related.
      ;
      • Gao Y.
      • Li S.
      • Lai Z.
      • Zhou Z.
      • Wu F.
      • Huang Y.
      • Lan X.
      • Lei C.
      • Chen H.
      • Dang R.
      Analysis of long non-coding RNA and mRNA expression profiling in immature and mature bovine (Bos taurus).
      ), PPP2R5A codes for protein phosphatase 2 regulatory subunit B alpha and is correlated with spermatogonial mitosis and it is related to sphingolipid signaling pathway which is associated with SM (
      • Sarakul M.
      • Elzo M.A.
      • Koonawootrittriron S.
      • Suwanasopee T.
      • Jattawa D.
      • Laodim T.
      Characterization of biological pathways associated with semen traits in the Thai multibreed dairy population.
      ;
      • Cai X.
      • Wu S.
      • Mipam T.D.
      • Luo H.
      • Yi C.
      • Xu C.
      • Zhao W.
      • Wang H.
      • Zhong J.
      Testis transcriptome profiling identified lncRNAs involved in spermatogenic arrest of cattleyak.
      ). Genes within the region of BTA 23, which showed an association with PM in the current study, belong to the OR family and KHDRBS2. The OR are involved in chemotaxis and determine speed, strength, and direction of movement of the spermatozoa toward the ovum (
      • Ali M.A.
      • Wang Y.
      • Qin Z.
      • Yuan X.
      • Zhang Y.
      • Zeng C.
      Odorant and Taste Receptors in sperm chemotaxis and cryopreservation: Roles and implications in sperm capacitation, motility, and fertility.
      ). Of the 20 to 66 testicular OR known in mammals, at least 8 are expressed in sperm cells, testis, and epididymus (
      • Milardi D.
      • Colussi C.
      • Grande G.
      • Vincenzoni F.
      • Pierconti F.
      • Mancini F.
      • Baroni S.
      • Castagnola M.
      • Marana R.
      • Pontecorvi A.
      Olfactory receptors in semen and in the male tract: From proteome to proteins.
      ). In humans, the OR influence spermatogenesis, epididymal growth, acrosome reaction, initiate flagella movement and facilitate cell-to-cell communication (
      • Olaniyan O.T.
      • Dare A.
      • Okotie G.E.
      • Adetunji C.O.
      • Ibitoye B.O.
      • Eweoya O.
      • Dare J.B.
      • Okoli B.J.
      Ovarian odorant-like biomolecules in promoting chemotaxis behavior of spermatozoa olfactory receptors during migration, maturation, and fertilization.
      ). The OR in the sperm cell activate proteins that increase intracellular calcium ions and hence motility (
      • Vanderhaeghen P.
      • Schurmans S.
      • Vassart G.
      • Parmentier M.
      Specific repertoire of olfactory receptor genes in the male germ cells of several mammalian species.
      ). KHDRBS2, has been associated with fertility in San Martinero cattle (
      • De León C.
      • Manrique C.
      • Martínez R.
      • Rocha J.F.
      Research article genomic association study for adaptability traits in four Colombian cattle breeds.
      ) and pregnancy status in Brahman beef cattle (
      • Reverter A.
      • Porto-Neto L.R.
      • Fortes M.R.S.
      • McCulloch R.
      • Lyons R.E.
      • Moore S.
      • Nicol D.
      • Henshall J.
      • Lehnert S.A.
      Genomic analyses of tropical beef cattle fertility based on genotyping pools of Brahman cows with unknown pedigree.
      )
      The QTL enrichment analysis for TM revealed QTL associated with feed intake, calving ease, carcass weight, and fat thickness on BTA 6, and antibody-mediated immune response and cell-mediated immune response on BTA 23. Previous studies have identified other QTL on BTA 23 related to genes such as IL4 (interleukin-4), antibody-mediated immune response, and DMI that may be associated with SM (
      • Wang H.
      • Misztal I.
      • Aguilar I.
      • Legarra A.
      • Muir W.M.
      Genome-wide association mapping including phenotypes from relatives without genotypes.
      ).

      CONCLUSIONS

      In this study, we identified a large number of genes associated with SM and involved with sperm and sperm-related traits in Italian Holsteins bulls. These findings suggest that many genes controlling SM could have pleiotropic effects and be functionally correlated with other traits related to fertility. Our results show the highly polygenic nature of SM, with each of numerous candidate genes responsible for a small proportion of the genetic variance. SM influences male fertility, a trait that has probably experienced strong selective pressure even before cattle domestication, about 10,000 years ago. As such, beneficial alleles of genes with major effects on this trait are likely to have reached fixation, leaving only variants with small effect in the domesticated populations, which were identified by GWAS as presented here. The characterization of selection signatures across different timescales, as well as the analysis of ancient genomes from domestic and wild animals, may help in identifying variants related to fertility that have experienced selection before and after domestication events.

      ACKNOWLEDGMENTS

      J.L.W. was supported by the project LEO Livestock Environment Open-data, 16.2 – PSRN 2014-2020 funded through Fondo Europeo Agricolo per lo Sviluppo Rurale (FEASR). A.A. and P.A.M. received support from the Italian Ministry of Education, University and Research (MUR; Rome, Italy) for Progetti PRIN2017 20174BTC4R. The authors have not stated any conflicts of interest.