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

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 ).

Illumina BovineSNP50 genotypes for 636 animals were imputed to high-density (777,962) with Beagle 4.1 () 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 () bovine assembly version using Plink 1.9 (). Quality control excluded indels, SNPs with minor allele frequency lower than 0.05 (), 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 (). Low quality imputed variants, with DR2 (dosage R², which is the estimated squared correlation between the estimated allele dose and the true allele dose), lower than 0.3 were excluded (; ), 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.

The genomic BLUP using individuals with repeated measurements were analyzed with a repeatability single trait animal model:

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:

where ${\sigma }_u^2,$ ${\sigma }_{pe}^2,$ and ${\sigma }_e^2$ 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 (). The heritability was estimated as the ratio of additive genetic variance to the total phenotypic variance $({\sigma }_u^2$ + ${\sigma }_{pe}^2$ + ${\sigma }_e^2)$ and the repeatability was calculated as sum of individual variance $({\sigma }_u^2$ + ${\sigma }_{pe}^2)$ divided by phenotypic variance.

The percentage of additive genetic variance explained by the ith SNP window was then calculated as

where ui is the genetic value of the ith genomic region under consideration, ${\sigma }_u^2$ 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 . All procedures were performed using BLUPF90 family programs (). To control false positives, we applied the false discovery rate (FDR) method for multiple testing according to as follows:

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).

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 () 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 () 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 ).

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).

As previously noted, our results indicated that the SNPs with the smallest P-values did not explain a large percentage of genetic variance (Figure 2–3), contrary to expectation. P-values are influenced by sample size () and, in particular, by allele frequencies ().

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.

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.

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).

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, ; 0.12, ) for PM. High heritability estimates, of up to 0.43 have been reported for motility in Holstein bulls (; ); however, estimates for PM in this breed are low (0.13–0.15, ). 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.

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 (), and in particular, by allele frequencies (). , 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, 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, 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; 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 (), fertilization, and normal embryo development in cattle (). Four profilin isoforms are expressed in mammalian testis, including spermatogenic cells (). Profilin-1 and profilin-3 play an important role in mouse sperm metabolism and motility () and the loss of profilin-3 has been associated with impairment of spermiogenesis (). 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 (). 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 (). 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⁺² dependent enzymes that are widely expressed in mammalian tissues including ovules, ovaries, testes, and in early embryos (). 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 (). The enzyme PADI6 is involved in decondensation of sperm chromatin (; ). In humans, sperm chromatin organization has been shown to affect sperm function during fertilization and early embryo development () and has been implicated in subfertility (). 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 (). 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 (). BoLA-NC1 is essential for the establishment and maintenance of pregnancy and in the maternal immune response to the fetus ().

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 (), pigs (), and humans () including sperm, seminal plasma, and in testicular somatic cells, particularly in Leydig and Sertoli cells (). GSTA also plays a role in spermatogenesis, sperm capacitation and sperm-oocyte binding, and is involved in the nuclear decondensation of the sperm chromatin (; ; ).

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 (; ) and IL17A is a member of the interleukin family that initiate a potent inflammatory response (). In humans the level of IL17 in seminal fluid has been shown to inversely correlate with SM and mortality () 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 (; ) 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 ().

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 (; ; ). 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 (). Tubulin plays a role in microtubule formation in sperm cells of buffaloes () and is expressed in spermatogonia, Sertoli cells, spermatids, and spermatozoa of pig testes (). In buffalo sperm, α-tubulin has been shown to provide structural support, facilitate cell motility, and to participate in oocyte binding and activation (). JAKMIP1 and PKDH1 are related to sperm function and cell polarization (), embryo development and survival (), 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 (). Additionally, PPP1R10, PPP1R11, PPP2R2C, and FSCN3 genes were identified and code for proteins related to spermatogenesis, SM, and fertility (, ; ; ), 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 (; ). 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 (). Of the 20 to 66 testicular OR known in mammals, at least 8 are expressed in sperm cells, testis, and epididymus (). In humans, the OR influence spermatogenesis, epididymal growth, acrosome reaction, initiate flagella movement and facilitate cell-to-cell communication (). The OR in the sperm cell activate proteins that increase intracellular calcium ions and hence motility (). KHDRBS2, has been associated with fertility in San Martinero cattle () and pregnancy status in Brahman beef cattle ()

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 ().