In a previous study, we performed a genome-wide association study for bull fertility using a total of 38,650 SNPs spanning the entire bovine genome (Peñagaricano et al., 2012
). Associations between SNPs and SCR were analyzed using a mixed linear model that included a random polygenic effect and each SNP genotype, one at a time, as a linear covariate or a categorical variable. The first analysis had more power to detect significant SNPs with additive effects, whereas the second approach allowed the detection of significant markers that show some degree of dominance. Out of the 38,650 genetic markers, 3,353 and 3,096 showed significant associations with SCR (nominal P
0.05) using SNP genotypes as linear covariates or categorical variables, respectively. In total, 4,450 different SNPs showed significant associations with bull fertility in one or both models used for the analysis.
The first step in this pathway-based analysis of a whole-genome association study was to assign SNPs to genes. A total of 21,325 of the 38,650 examined SNPs were located within annotated genes or within 20
kb upstream or downstream from annotated genes. This set of SNPs defined a total of 16,819 genes annotated in the UMD 3.1 bovine genome sequence assembly, which in turn were evaluated for the pathway analysis. In addition, 2,767 of the 16,819 genes had one or more SNPs with nominal P
0.05; hence, these genes were defined as significantly associated with bull fertility.
The next step in the analysis was to assign the 16,819 genes into pathways. We tested GO categories and IP entries with more than 30 genes and, in the case of GO terms, located between levels 5 and 9 in the GO hierarchy. A total of 662 GO terms and 248 IP entries met these requirements and hence were tested by a hypergeometric test for enrichment of significant genes associated with bull fertility.
Twenty GO terms showed significant overrepresentation of genes statistically associated with SCR. Table shows the 11 GO terms that belong to the biological process domain, and Table displays the nine GO terms classified into the molecular function (4) and cellular component (5) domains, respectively. These tables show the number of genes in each functional category (out of 16,819 genes), the expected number of significant genes under the null hypothesis (i.e., with no significant overrepresentation), and the actual number of significant genes per category. Only GO terms with q
0.05 were considered significant. In addition, one entry of the IP database, Haloacid dehalogenase-like hydrolase
(IPR005834) showed significant enrichment (q
0.027) of genes associated with bull fertility. This entry, with a total of 37 genes, showed 16 significant genes, which is 10 more than expected under the null hypothesis. Genes statistically associated with SCR (nominal P
0.05) within each significant category (i.e., the 20 GO terms and IPR005834) are listed in File S1 in Supplementary Material.
Biological process terms significantly overrepresented with genes statistically associated with bull fertility.
Molecular function and cellular component terms significantly overrepresented with genes statistically associated with bull fertility.
Six of the 11 significant GO terms classified into the biological process domain show a very close relationship in the GO hierarchy (Table ). Indeed, Rho protein signal transduction
(GO:0007266) is part of Ras protein signal transduction
(GO:0007265) and Regulation of Ras protein signal transduction
(GO:0046578) is part of Regulation of small GTPase mediated signal transduction
(GO:0051056). Furthermore, Ras protein signal transduction
is part of Small GTPase mediated signal transduction
(GO:0007264), which in turn is part of Intracellular signaling cascade
(GO:0007242). Interestingly, GO:0046578 and GO:0051056 negatively regulate GO:0007265 and GO:0007264, respectively. Overall, these results point out that an intracellular signaling pathway involving small GTPases and their negative regulators was significantly associated with bull fertility. Interestingly, a recent microarray study showed that Regulation of small GTPase mediated signal transduction
term was significantly enriched with genes differentially expressed between human sperm that achieved successful fertilization and sperm that failed the intracytoplasmic sperm injection (García-Herrero et al., 2011
Small GTPases coordinate diverse cellular functions, including actin reorganization, junction dynamics, cell movement, cell cycle, cell transformation, and gene transcription (Bustelo et al., 2007
). Several recent studies have shown that small GTPases are relevant molecules for spermatogenesis. For instance, some studies have reported that the translocation of elongating/elongate spermatids across the seminiferous epithelium is assisted by small GTPases present at the site of the ectoplasmic specialization, a unique anchoring junction type located in the testes (Lui et al., 2003
). In addition, other studies have shown that small GTPases are involved in regulation of the mammalian acrosome reaction by controlling the membrane fusion system in the sperm (Iida et al., 1999
). The acrosome reaction allows spermatozoa to penetrate the zona pellucida and fuse with the oocyte membrane and, hence, this is a crucial step during the fertilization (Brucker and Lipford, 1995
). Overall, these results provide further evidence for the possible role of small GTPases and their regulators in spermatogenesis and, hence, in male fertility.
Three significant GO terms classified into the biological process domain are related to neurongenesis (Table ). In fact, Neuron development
(GO:0048666) is part of Neuron differentiation
(GO:0030182), which in turn is part of Neurogenesis
(GO:0022008). One well-known connection between neurons and male fertility is provided by the hypothalamic–pituitary–gonadal axis. Indeed, the maintenance of spermatogenesis depends upon stimulation of the testes by the gonadotropic hormones, follicle-stimulating hormone, and luteinizing hormone, both produced by the pituitary gland in response to gonadotropin-releasing hormone from the hypothalamus (Zirkin, 1998
). It is important to note that recent studies have reported de novo
neurogenesis in postnatal and adult mammalian hypothalamus (Pierce and Xu, 2010
). Another well documented connection between neurogenesis and spermatogenesis is the fact that these two biological processes share several key molecular regulators. Heat shock factors and proteins, bone morphogenetic proteins, and inhibitor of DNA-binding/differentiation (Id) proteins, among others, play relevant roles in both processes (Yokota, 2001
; Abane and Mezger, 2010
; Björk and Sistonen, 2010
; Bragdon et al., 2011
). Heat shock factors, a family of transcriptional regulators, trigger the expression of genes encoding proteins (i.e., heat shock proteins) that function as molecular chaperones, contributing to establish a cytoprotective state in response to different stress conditions (Björk and Sistonen, 2010
). In addition to promoting cell survival under stressful conditions, heat shock factors and proteins are important for developmental processes such as spermatogenesis and neurogenesis (Abane and Mezger, 2010
). Bone morphogenetic proteins, potent growth factors belonging to the transforming growth factor beta superfamily, are responsible for many biological processes including embryogenesis, neurogenesis, and spermatogenesis (Bragdon et al., 2011
). Id proteins comprise a family of proteins that are positive and negative regulators of proliferation and differentiation, respectively (Yokota, 2001
). Id proteins are key regulators of several developmental processes such as neurogenesis, lymphoid organogenesis, mammary gland development, and spermatogenesis (Yokota, 2001
). Altogether, these results present further evidence of the tight connection between neurogenesis and male fertility.
Four GO terms classified into the molecular function domain showed significant association with bull fertility (Table ). Three of these GO terms showed a very close relationship in the GO hierarchy: Nucleoside-triphosphatase activity
(GO:0017111) is part of Pyrophosphatase activity
(GO:0016462), which in turn is part of Hydrolase activity, acting on acid anhydrides, in phosphorus-containing anhydrides
(GO:0016818). Nucleoside-triphosphatases are enzymes that catalyze the hydrolysis of a nucleoside triphosphate in a nucleoside diphosphate and phosphate. Pyrophosphatases are enzymes that catalyze the hydrolysis of a pyrophosphate bond between two phosphate groups. Both types of enzymes, which belong to the family of hydrolases, participate in purine and thiamine metabolism (Zöllner, 1982
). Interestingly, Sertoli cells, which are located in the seminiferous tubules and play an essential role in the maintenance and control of spermatogenesis, have shown remarkable activity of these two enzymes (Casali et al., 2001
). In fact, these enzymes play an important role in the physiological control of extracellular levels of nucleotides and nucleosides inside the seminiferous tubules, which in turn can modulate Sertoli cell responses (Casali et al., 2001
). Thus, findings of this study provide further support for the relevant role of these enzymes (i.e., nucleoside-triphosphatases and pyrophosphatases) in spermatogenesis.
The calcium ion binding
(GO:0005509) term showed a significant overrepresentation of genes associated with bull fertility (Table ). The relationship between calcium and sperm physiology is well documented. Mammalian spermatozoa possess multiple voltage-gated calcium channels and use calcium ion signals to control several physiological responses (Darszon et al., 1999
). Indeed, calcium is considered a regulator of sperm motility, a participant in capacitation, and an essential second messenger for the acrosome reaction (Wennemuth et al., 2003
). Therefore, our findings provide more evidence of the strong association between calcium and sperm biology and, hence, male fertility.
Three of the five significant GO terms grouped within the cellular component domain are related to the cytoskeleton biology. Indeed, Actin cytoskeleton
(GO:0015629) and Cytoskeletal part
are part of Cytoskeleton
(GO:0005856). The close relationship between cytoskeleton and spermatogenesis is well documented. Sertoli cells are characterized by their well-developed cytoskeleton, which is responsible for the collective organization of the seminiferous epithelium (Mruk and Cheng, 2004
). Morphological studies have shown that the Sertoli cell cytoskeleton, among other functions, positions, anchors, and aids in the movement of developing germ cells and participates in the release of mature spermatids from the seminiferous epithelium at spermiation (Mruk and Cheng, 2004
). In fact, the different cellular events that occur during spermatogenesis are associated with extensive changes in cell shape and size and germ cell movement, in which the cytoskeleton has a key role (Lie et al., 2010
). Interestingly, García-Herrero et al. (2010
) have reported that the Cytoskeleton
term was significantly enriched with genes differentially expressed between sperm of infertile and fertile human males. Therefore, our study provides further evidence of the close relationship between cytoskeletal dynamics and spermatogenesis.
The terms Secretion by cell
(GO:0046903), (Table ) and Extracellular matrix part
(GO:0044420), and Plasma membrane part
(GO:0044459; Table ) also showed significant enrichment of genes statistically associated with bull fertility. Despite the wide range of functions of these categories, it is worth noting that spermatogenesis and fertilization involve biological processes related to these terms. For instance, many of the functions of the Sertoli cells are carried out by their protein secretions: these cells synthesize and secrete a number of proteins that have functions such as hormone-like or growth factor-like activity, transport functions, enzymatic activities, as well as proteins that contribute to the basement membrane (Griswold, 1998
). In addition, it seems that the plasma membrane plays an important role in Sertoli-Sertoli and Sertoli-germ cell interactions (Mruk and Cheng, 2004
). Also, it has been reported that the extracellular matrix regulates Sertoli cell differentiation and germ cell development (Hadley et al., 1985
InterPro is a database that integrates diverse information about families, domains, and functional sites of proteins. One entry of the IP database, Haloacid dehalogenase-like hydrolase
(IPR:005834), showed significant overrepresentation of genes associated with SCR. IPR:005834 entry contains a group of hydrolase enzymes including l
-2-haloacid dehalogenase, epoxide hydrolases, and phosphatases. These enzymes are involved in diverse cellular functions such as membrane transport, signal transduction, metabolism, and DNA repair (Arand et al., 2003
). Interestingly, Iguchi et al. (2006
), studying gene expression profile in meiotic male germ cells, reported that one haloacid dehalogenase-like hydrolase is highly expressed in the testes. Thus, our finding provides further evidence of the possible involvement of these enzymes in male fertility.
Previous studies have derived gene networks and pathways related to complex traits in livestock species. For instance, Fortes et al. (2012b
) have used genome and trait associations, hypothalamic-transcriptome data and transcription factors in order to find candidate genes and gene–gene interactions related to first service conception in Brangus heifers. Moreover, Lewandowska-Sabat et al. (2012
) have integrated genomic and gene expression data aimed at identifying significant pathways involved in the immune response to mastitis in dairy cattle. Here, we have extended our previous genome-wide association study to identify biologically relevant pathways involved in bull fertility. Previously, we have reported eight SNPs associated with SCR with a genome-wise adjusted P
-value below 0.10 (Peñagaricano et al., 2012
). Thus, beyond the enumeration of possible candidate genes (i.e., genes harboring the most significant genetic markers) done in our previous report, in this study we have performed a gene-set enrichment analysis to characterize quantitative trait pathways related to bull fertility. The identified pathways could contribute to a better understanding of the genetic architecture of male fertility in dairy cattle. Moreover, as discussed by Snelling et al. (2012
), genetic markers within and surrounding genes initially identified by gene-set enrichment analysis may facilitate the incorporation of marker-assisted selection for complex traits in commercial breeding schemes.
One potential limitation of this study is the limited number of GO categories and IP entries that were analyzed. To avoid testing overly narrow or broad functional categories, while at the same time address the multiple-testing issue, we have analyzed only functional categories with more than 30 genes and, in the case of GO terms, those that were located between levels 5 and 9 in the GO hierarchy. As a result, a total of 662 GO categories and 254 IP entries were interrogated. It is important to note that a key factor in these analyses is the total number of genes, which determines the total number of functional categories that can be inferred. Here, using data from 38,650 SNPs spanning the entire bovine genome, we could analyze a total of 16,819 annotated genes. Future studies, however, could analyze a larger numbers of genes using a high-density SNP chips (e.g., Illumina Bovine HD 770K), which in turn will increase the number of functional categories and provide even further insight into the genetic architecture of bull fertility.