Bayesian-genome-wide association study and post-GWAS on reproductive traits of Holstein dairy cattle

Abstract

Bayesian multiple regression models are often used for genomic selection, where all markers are adjusted simultaneously as random effects to reduce the false-discovery rate. The purpose of this research was to identify individual candidate genes for some reproductive traits in Holstein cattle through genome wide association studies (GWAS) using a Bayesian method in combination with post-GWAS analysis. Reproductive traits included: days open (DO), pregnancy rate (PR), calving interval (CI) and age at first calving (AFC). The animals were genotyped using single-nucleotide polymorphism (SNP) panels of different densities imputed to a 50 K SNP density. After quality control, we included 2400 genotyped animals. According to the Bayesian analysis, there were 19 windows with an explained additive genetic variance of >0.1 percent for CI, DO, AFC and PR in Holstein cattle, which were 3, 3, 6 and 7, respectively. Using Bayesian analysis, 79 genes were located within or nearby (250-kb) 19 significant SNPs/windows in the Bos taurus autosomes. Among these genes, we identified 25 candidate genes for reproductive traits, namely CHD7, CLVS1, EVX2, MAT2B, NUDCD2, GPR39, NCKAP5, LYPD1, HOXD13, SEMA5B, CCNG1, SEMA5A, BRF1, PSEN2, CACHD1, SUGTA, ELF1, SNORA70, AKT1, TM2D1, SLF1, MCTPA, PAB2A, MTRF1 and ADCY2. Additionally, another nine candidate genes (CLVS1, GPR39, CENPF, AMOT, ARF1, CCDC186, ADCY2, BDP1 and AMOTL1) were identified in the network cluster analysis as hub genes for reproductive traits. The results of gene set enrichment analysis (GSEA) and pathway analysis, suggest that the most important gene ontology term involving cellular metabolic process was related to the AFC trait. To summarize, Bayesian methods were used to identify SNPs and candidate genes that could be useful in genomic selection to improve reproductive traits of Holstein dairy cattle.

In recent times, genetic progress in production animal breeding has increasingly become the result of selection based on genomic information of important economic traits (Miglior et al., 2017). Genome-wide association studies (GWAS) can be used to identify genomic regions or causative mutations that account for a significant proportion of the genetic variance associated with a trait. The GWAS approach is based on the assessment of correlations between genetic markers and the trait of interest in a large number of samples (Jiang et al., 2019). The ability to identify single-nucleotide polymorphisms (SNPs) that are associated with economically important traits is the main advantage of GWAS in dairy cattle breeding. Various statistical models are used in GWAS, including Bayesian multiple regression and generalized linear mixed models. Since most economic traits are controlled by multiple loci, estimates of marker effects by linear mixed models have low accuracy. Another challenge of GWAS is the multiple testing correction. Many loci may fail to overcome the strict criterion of significance test because the normal Bonferroni correction is too cautious (Wang et al., 2016). When the goal is reduction of experimental power to achieve high stringency, traditional multiple testing corrections such as Bonferroni do not take into account the genotype dependencies among SNPs in linkage disequilibrium (Fernando et al., 2004). An alternative approach is to use Bayesian multiple regression models for genomic selection, where all markers are adjusted simultaneously as random effects to reduce the false-discovery rate (Kemper et al., 2015). A clear advantage of multi-locus models is that they do not require Bonferroni correction, making them more robust approaches for GWAS. Although Bayesian methods use shrinkage approaches to handle a large sample size of markers, they find difficulty when the number of markers increases due to memory allocation or computation time limitations. Furthermore, these models suffer from multicollinearity problems in situations where marker density is unusually high (Segura et al., 2012). Some of these problems can be solved using the Bayesian multiple regression model, which was initially designed for genome prediction. Bayesian methods allow the effects of all genotyped markers to be adjusted simultaneously by accounting for different prior distributions of marker effects. Furthermore, additional sources of information such as omics data or annotation information can be incorporated (Sollero et al., 2017). The purpose of the research described here was to identify individual candidate genes for some reproductive traits in Holstein cattle through a GWAS using a Bayesian method in combination with post-GWAS analysis.

Materials and methods

Phenotypic data

Pedigree and phenotypic records for insemination and calving dates of Iranian Holstein dairy cattle were obtained from the Animal Breeding Center of Iran. The records of Holstein heifers and cows (lactation 1–3) were collected between 2000 and 2021. Reproductive traits included were days open (DO), pregnancy rate (PR), calving interval (CI) and age at first calving (AFC). A summary of the pedigree file that was used in genetic analysis is presented in online Supplementary Table S1. Descriptive statistics for AFC and reproductive traits are presented in online Supplementary Table S2.

Genotypic data

The animals were genotyped using SNP panels of different densities (Illumina and GeenSeek Genomic companies) imputed to a 50 K SNP density. The following quality control (QC) criteria were used to evaluate the processes: deviations from Hardy–Weinberg equilibrium (P-value < 10⁻⁶), animal call rate (> 90%), SNP call rate (> 95%) and minor allele frequency (> 0.01). We performed quality control using the PLINK software. Additionally, duplicate, twin and animals who failed parentage tests were eliminated from the analysis. QC was initially carried out on the 54 001 SNP array and SNPs on chromosome X or with an unclear physical location were eliminated. Thus, 51 185 markers were found on the 54 001 SNP chip. The density of each SNP chip after imputation was 51 185. The 2400 animals genotyped with 41 099 SNP markers demonstrated success with QC filtering for association analysis. The SNP number and animal number of the final analysis are presented in online Supplementary Table S3.

Statistical model

To estimate genetic parameters for the AFC trait, the following sire model was fitted:

(1)$$y_{{\rm ik}} = {\rm \mu } + HYS_{\rm i} + {\rm \alpha }_k + e_{{\rm ik}}$$

Where: y _ikis the kth observation; μ is the population mean; fixed effects include: HYS _i: herd-years-season; α_k is additive genetic α_k~N (0, ${\rm A \alpha }_a^2 $), and e _ik: is the random residual error e _ik~N (0, I).

Using the following repeatability single-trait model, the genetic variance components, breeding values, and corresponding reliability reproductive traits were estimated:

(2)$$\eqalign{y_{{\rm ijlmk}} = {\rm \mu } + HYS_{\rm i} + b( Age_{\rm j}-\overline {Age_ \ldots \;) } + P_l + Pe_{\rm m} + {\rm \alpha }_k + e_{{\rm ijlmk}}}$$

Where Age _j: is the kth covariate of age at calving; Pe _m: is the random effect of the permanent environment; P _l: is the fixed effect of parity (3 levels); other terms are defined above. Variance components and heritability were estimated by Bayesian method and using Gibbs2f90 software from the BLUPF90 program family. We used the deregressed proofs (DRP) as pseudo phenotypes for DGV and GWAS prediction in the next multiple-step method. To calculate DRP, DRP packages from the R programming language were utilized.

Bayesian genome-wide association analysis

Using the Bayes-CP _i method, GWAS were performed with GenSel software, according to the presented model:

(3)$$y_i = {\rm \mu } + \mathop \sum \limits_{\,j = 1}^k Z_jU_{\,j\;}{\rm \delta }_j + e_i$$

where y _i is DRP on animal i for the respective trait; μ is the population mean; k is the number of markers; Z _j is the vector includes genotypes for the SNPj (coded as: −10 homozygous dominant, −10 homozygous recessive and 0 heterozygous); U _j is the random substitution effect for marker j; δ_j is a random 0/1 variable showing the absence or presence of marker j, and e _i is the random residual effects assumed normally distributed N (0,$\alpha _e^2 $). We used the Markov Chain Monte Carlo (MCMC) method with the GenSel implementation to obtain the posterior distributions of the SNP effects. After a burn-in period of 1000 iterations in which data were discarded, 51 000 iterations were performed to determine the posterior mean effect of each SNP. Multiple regression models with π = 0.9 were fitted simultaneously in each MCMC iteration of the Bayesian variable selection process. Finally, we displayed the Manhattan plot of genetic variance using the R program.

Genome-wide association study and functional enrichment analyses

Following GWAS analysis, gene set enrichment analysis (GSEA) was used to identify molecular pathways and putative functional categories associated with reproductive traits in Holstein cattle. Information from the NCBI platform (National Science for Biotechnology Information) was utilized in this context. The Ensembl (https://www.ensembl. org; NCBI, https://www.ncbi.nlm.nih.gov/genome/; genome-browser https://genome.ucsc.edu) platform was used to assess the biological functions of the genes and genetic ontology. STRING (https://string-db.org) databases were employed for GSEA. The BiomaRt package was used for positional genes were located within or nearby 250 kb of the significant SNPs.

Results and discussion

Genetic parameter estimation

Since reproductive traits are known to have low heritability, achieving increased genetic progress compared to higher heritable variables is more challenging. Estimates of heritability for reproductive traits (AFC, CI, DO, and PR) in Iranian Holstein cattle are presented in Table 1. Compared to other reports on Holstein cows the heritability value for the AFC trait was higher, however, other traits showed comparable heritability. Previous studies have estimated heritability similarly for the AFC trait (0.19) in dairy cattle (Ali et al., 2019), and for reproductive traits in Iranian Holstein cows (Faraji et al., 2011). The heritability of the AFC trait varies between 0.05 and 0.74 in different reports (Pirlo et al., 2000). Moreover, heritability is a specific property that is unique for population, environmental conditions and imposed management decisions.

Table 1. Variance components, and heritability estimates for fertility traits in Iranian Holstein cattle

Genome-wide association analysis

Using GWAS based on Bayesian methods, we were able to determine which SNP marker regions were most significant in the Illumina BovineSNP50 panel. The relationship between SNPs and reproductive traits in Holstein cattle is represented by a Manhattan plot in Fig. 1, and Table 2 shows the results of the GWAS for reproductive traits. GWAS results include chromosomal location and proportion of variance of 1-Mb genome windows by informative regions (more than 1.0%). According to the Bayesian analysis, there were 19 windows with an explained additive genetic variance of >0.1 percent for CI, DO, AFC and PR in Holstein cattle, which were 3, 3, 6, and 7, respectively. Quantitative trait loci (QTL) associated with the target traits were defined as windows that account for ≥ 1.0% of the total variance.

Figure 1. Manhattan plots of the genome-wide association studies for reproductive traits. The x-axis and the y-axis represent the chromosome number and the observed windows, respectively. The SNPs were plotted against their genomic positions. (a) DO, (b) AFC, (c) CI, (d) PR.

Table 2. Significant SNPs identified for reproductive traits in Iranian Holstein cows using Bayes-GWAS approach

Using Bayesian analysis for reproductive traits in Holstein cows

Using Bayesian analysis, 79 genes were located within or nearby (250-kb) to 19 significant SNPs/windows in the Bos taurus autosomes. The windows with the most informative located within or nearby (250-kb) detected CI, DO, AFC and PR traits on (BTA7) at 75 Mb 0.18%, (BTA3) at 84 Mb 0.23%, (BTA3) at 17 Mb 0.21% and (BTA21) at 69 Mb 0.28%, respectively. Significant SNPs reproductive traits in Iranian Holstein cattle are shown in Table 2.

Due to the high frequency of QTL controlling reproductive traits in this region, we identified 132 important SNPs through Bayesian analysis. There is particular interest in the AFC trait on BTA3, as shown in Fig. 1b. The region of BTA3 at 1-Mb was the most significant window region. However, significant SNPs were also associated with heifer reproductive trait within the SPRR4, CLVS1, SNORA70, CCNG1, NUDCD2 and MAT2B genes. A total of 23 genes were found within or near the 250 kb of these SNPs.

The association analysis detected significant SNPs related to the PR trait, and we identified 152 significant SNP, whislt BTA21 at 1-Mb was the most significant window region as shown in Fig. 1d. Among the significant SNPs associated with the PR trait, SNP/windows on chromosome 21 were located within the BRF1 gene. BRF1 gene is one of the important genes affecting fertility traits (Butler et al., 2020).

BTA7 at 1-Mb was the most significant window region, with 70 significant SNPs identified for the CI trait, as shown in Fig. 1c. We identified 57 significant SNP for the DO trait, and BTA3 at 1-Mb was the most significant window region (Fig. 1a). Our results for the 1-Mb and significant regions are in agreement with previously found regions (BTA3, BTA14, BTA1, BTA21, BTA17 and BTA12) in the animal QTL database that may have an impact on reproductive traits. With Bayesian-GWAS, a substantially higher level of genetic variance allowed for the identification of the more significant essential regions.

Finding the genomic regions in GWAS that influence reproductive traits is interesting because dairy cow fertility is a complex trait with low heritability that involves multiple biological pathways. However, this study provided new information on many genes and molecular networks affecting reproductive traits in Holstein cattle, as shown in Table 2. The CHD7 gene, which encodes the protein known as chromodomain helicase DNA binding protein 7, was one of the important positional candidate genes within 1-Mb upstream and downstream of this overlapping window. In addition, the genes on BTA14 (AFC) have previously been reported by Mohammadi et al. (2020). In our research, the CHD7 gene was identified in the AFC (BTA14) trait, consistent with the findings of Mota et al. (2017), which plays an important role in reproductive traits. In addition, this gene influences the growth trait and fertility of Simmental cattle (Fonseca et al., 2018), meat quality traits in Brahman cattle (Braz et al., 2021) and lactation persistency in Canadian Holsteins (Do et al., 2017). In Nellore cattle, the CLVS1 gene plays an essential role in determining the age at first calving trait (Mota et al., 2017). The SEMA5A and MCTP1 genes play a significant antibacterial role in mammary gland infection by controlling neutrophil migration to the infection site and also affect carcass traits, beef tenderness and oviduct–embryo interactions, once again in Nellore cattle (Santana et al., 2015). The ELF1 gene has been reported as a candidate gene related to feed efficiency traits that help to reduce environmental consequences from livestock production, enhancing beef cattle profitability (De Lima et al., 2020). The small proline-rich proteins (SPRR4) gene affects reproductive traits, and 175-kb flanking regions of BTB-00661215 SNP were detected in Iranian Holstein cows as a potential candidate gene for the AFC trait (Mohammadi et al., 2022). Additionally, in Brahman cattle, the PAB2A gene is associated with scrotal growth and male fertility (Soares et al., 2017). The genes HOXD, EVX2 and KIAA1715 were found to play important roles in limb function, digit number and lumbar vertebrae in the potential genome region of sheep (Ren et al., 2016). In the peri-implantation phase of pregnancy, the MAT2B gene regulates the expression of critical genes required for cell-specific metabolism and transport of polyamines in the sheep endometrium, affecting fertility and offspring health, growth rates and changes in the embryonic stage of bovine fibroblast cells with DNA methylation, but also has relationships with milk production characteristics of Chinese Holstein cows (Ye et al., 2022). The expression of cow genes related to the cell cycle, membrane stability and apoptosis, as well as other embryonic genes, is influenced by the NudC domain containing 2 (PDCD2) gene (Hoelker et al., 2012). The G protein-coupled receptor 39 gene (GPR39), as in other mammals, is crucial for digestive and metabolic processes in cows (Yamamoto et al., 2008). The NCK associated protein 5 (NCKAP5) gene was studied in various cattle breeds and crosses to assess its influence on temperamental traits. This gene encodes a transcription factor essential for osteoblast development and chondrocyte maturation, processes that may influence crucial phenotypes of the Limousin breed (Mariadassou et al., 2020). In addition, HOXD13 and SLF1 genes were found to affect sperm morphometry through inbreeding in a cattle breed and related to sperm DNA methylation and reproductive traits in cattle (Liu et al., 2019). We found a significant SNP on BTA20 for PR traits in Holstein cows. This SNP is located in the semaphorin5A (SEMA5A) gene. The SEMA5B gene is an essential paralog of SEMA5A and is crucial for fertility pathways in cows (Cai et al., 2019) as well as improvements in mastitis resistance and there are indications that the FEZL-SEMA5A pathway may regulate innate immunity and neuronal development in cattle. According to Han and Peñagaricano (2016), the BRF1 gene has been associated with bull fertility in Holstein cattle. The CACHD1 and CDC42BPA genes reportedly affect the quality of meat produced beef cattle (Ryu and Lee, 2016). According to Li et al. (2020), the PSEN2 and SUGT1 genes are linked to a haplotype that in dairy cows affects fertility and leads to altered conception, early embryonic losses, pregnancy losses and mutations. This haplotype has also been linked to heritable abnormalities called BMS (bovine male subfertility) in the Russian Simmental population. The biological mechanisms underlying a number of economically important traits in cattle breeds are linked to the SNORA70 gene (Taye et al., 2017). Furthermore, activation of the AKT1 gene in the epithelial cells of the mammary gland plays a role in the transmission of lactation signals (Wan et al., 2011). Research has shown that the TM2D1 gene is associated with adaptation and fertility traits in Colombian creole cattle breeds, as well as with apoptosis of ovarian follicles in cattle and buffalo (Zoheir et al., 2021). The KIAA0825 gene was found to have an impact on endoparasite infection, general health and fertility in German Holsteins (Wolf et al., 2021). The adenylate cyclase 2 (ADCY2) gene is one of the primary candidates for reproductive traits, having effects on age at first calving and first calving interval in water dairy buffaloes and connected to several biological processes, and the Yanbian cattle have previously shown an association with meat tenderness (Silva et al., 2020). Mutations in the MTRF1 gene have related to a significant effect on body weight and other quantitative traits, as reported by (Reynolds et al., 2022).

Post-GWAS (GSEA and network cluster analysis)

Gene enrichment analysis was used in GWAS to identify potential functional categories and molecular pathways associated with reproductive traits in Iranian Holstein cattle, and these data are presented in online Supplementary Table S4. A total of 43 biological terms were significantly associated with reproductive traits (FDR < 0.05). Given that many gene ontology (GO) terms are correlated, the data need to be interpreted with caution. However, after the Bonferroni correction for post-GWAS, the top 43 terms were highly enriched in the overexpressed-up gene set. As concerns significant GO terms for biological processes for reproductive traits in our Iranian Holstein cows, the most significant terms were ‘GO:0006412’ Translation involved the cellular metabolic process related to the AFC trait. Network clusters, and interactions between potential genes for reproductive traits are shown in Fig. 2. The genes were classified into three clusters by the network analysis using the K-means algorithm technique (red, green and blue). In the first (red) cluster CLVS1 and GPR39 (in relation to AFC), AMOT and AMOTL1 (in relation to DO), ARF1 and CCDC186 (in relation to CI), ADCY2 and BDP1 (in relation to PR) genes were the hub genes of the network. Moreover, the hub genes of the second cluster (green) were CHD7 and ICT1 (in relation to AFC), ANKRD32 and ATM (in relation to DO), GHRL and GPR39 (in relation to CI), C14orf79 and CEPT1 (in relation to PR). Also, BHMT and BHMT2 hub genes of the third cluster were related to AFC, FAM81B and PITPNB (in relation to DO), CA8 and ITPR1 (in relation to CI) and FASTKD3 and JAG2 (in relation to PR). The first cluster's annotated genes were significantly involved in the regulation of body weight, gastrointestinal mobility, hormone secretion and cell death (in relation to AFC), and may play a role in the polarity, proliferation and migration of endothelial cells (in relation to DO), affect the cell cycle and differentiation of various tissues (in relation to CI), and be involved for instance in cell cycle control and signal transduction (in relation to PR).

Figure 2. Network clusters between the identified genes for (a): age at first calving (AFC), (b) days open (DO), (c) calving interval (CI), and (d) pregnancy rate (PR) traits in Iranian Holstein cattle.

Gene ontology (GO) terms for biological processes

The biological process (‘GO:0006412’) associated with the regulation of gene expression, protein translation and indeed any process related to the AFC trait are presented in online Supplementary Table S4 and the important role of this terms in feed efficiency in beef cattle and bovine oocyte quality and female fertility (Macaulay et al., 2016) has already been reported. The association of embryonic limb morphogenesis (‘GO:0030326’) term for the DO trait is a biological process that occurs in the embryo by which the anatomical structures of the limb are generated and organized and the biosynthetic process (‘GO:0044272’) terms related with mastitis induced by E. coli was demonstrated by Darang (Darang et al., 2022). Additionally, the role of this GO:0051649 term in meat performance in pigs has been reported by Wimmers et al. (2010). The genes included in (‘GO:1901566’) affect organonitrogen compound biosynthetic processes, the decline in follicle-stimulating hormone, the bovine granulosa cell pathway and the mitochondrial translational termination pathway, whilst (‘GO:0070126’) affects pre-implantation embryo development in cattle (Clark, 2021).

Network cluster analysis

Network analysis of the first cluster CLVS1 and GPR39 hub genes associated with the AFC trait is shown in Fig. 2a. The GPR39 gene is known for its important role in gastrointestinal and metabolic functions in cattle and other mammals. Additionally, the CLVS1 gene is a potential candidate gene for carcass traits due to its biological functions. Regarding the gene hub-related PR trait (Fig. 2d), the AMOT gene has a significant impact on regulating gene expression in mouse, cow and human morula embryos. The role of the hub gene, ADP ribosylation factor 1 (ARF1) in cell proliferation and lactation of bovine mammary epithelial cells has been investigated. Shown in Fig. 2c, the CCDC186 hub gene (first cluster involved in CI) helps in milk production traits in water buffalo, as well as affecting of feed efficiency in mid-lactation Holstein dairy cows. In addition, the role of the BDP1 gene associated with white blood cell count in Chinese Holstein cows has been reported by Chen et al. (2022). The CHD7 gene (second cluster, associated with AFC) is involved in the organization of chromatin and has associated GO terms related to biological processes including in-utero embryonic development and cow genitalia development and may play an important role in regulating puberty and reproduction as well as having effects on in vitro fertilization, embryo transfer and superovulation in Holstein cattle. Additionally, in beef cattle, the CEPT1 gene may have an important regulatory function in lipid metabolism. Finally, the ITPR1 gene, a member of the third cluster, was associated with both male and female fertility and also with embryo size and morphology in both mice and humans (Cole et al., 2018).

In conclusion, using Bayesian analysis, 79 genes were located within or near to 19 significant SNPs/windows in the Bos taurus autosomes. Among these genes, 25 candidate genes for reproductive traits were identified, namely CHD7, CLVS1, EVX2, MAT2B, NUDCD2, GPR39, NCKAP5, LYPD1, HOXD13, SEMA5B, CCNG1, SEMA5A, BRF1, PSEN2, CACHD1, SUGTA, ELF1, SNORA70, AKT1, TM2D1, SLF1, MCTPA, PAB2A, MTRF1 and ADCY2. Additionally, the candidate genes CLVS1, GPR39, CENPF, AMOT, ARF1, CCDC186, ADCY2, BDP1 and AMOTL1 were identified in the network cluster analysis as hub genes for reproductive traits. The results of gene set enrichment analysis and pathway analysis suggest that the most important gene ontology term involving cellular metabolic process was related to the AFC trait. A combination of Bayesian GWAS and post-GWAS analysis can be effectively used to identify potential candidate genes and causative mutations that could improve reproductive traits in dairy cattle.