Selection strategies focusing on milk yield and associated measures in dairy cattle have led to significant increases in milk production traits. For example, during the past 20 years, increases of 3500 kg of milk, 130 kg of fat and 100 kg of protein per cow per lactation have been achieved in dairy cattle (Shook, Reference Shook2006). However, this success was accompanied by a decline in fertility of the modern high-producing dairy cow, such that infertility has become a major concern of farmers and the dairy industry (Pryce et al. Reference Pryce, Royal, Garnsworthy and Mao2004). As an indication of suboptimal reproductive performance of dairy cows, the first-service pregnancy rate has declined from 70% to 40% (Dobson et al. Reference Dobson, Walker, Morris, Routly and Smith2008). Many reasons account for this decrease in reproductive efficiency including changes in physiology, nutritional management and genetics (Lucy, Reference Lucy2001; Veerkamp & Beerda, Reference Veerkamp and Beerda2007). The antagonistic relationship between production traits and reproductive performance of dairy cows has been shown via the estimation of genetic correlations in quantitative genetic studies (e.g. Dematawewa & Berger, Reference Dematawewa and Berger1998; Kadarmideen et al. Reference Kadarmideen, Thompson and Simm2000; Royal et al. Reference Royal, Flint and Woolliams2002). Different mechanisms may underlie the negative association between yield and fertility including pleiotropic gene effects and linkage or complex physiological associations (Veerkamp et al. Reference Veerkamp, Beerda and van der Lende2003).
Among known genes that affect milk production, milk proteins have received considerable attention in recent years because of possible associations between milk protein genotypes and economically important traits in dairy cattle. More than 95% of the proteins contained in bovine milk are coded by six genes: the two main whey proteins, α-lactalbumin (α-LA) and β-lactoglobulin (β-LG), coded by the LALBA and LGB genes, and the four caseins αS1-CN, αS2-CN, β-CN, and κ-CN coded by the CSN1S1, CSN1S2, CSN2, and CSN3 genes, respectively (Caroli et al. Reference Caroli, Chessa and Erhardt2009). A large number of studies have consistently documented a very high positive genetic correlation between milk and protein yields (Chauhan & Hayes, Reference Chauhan and Hayes1991; Welper & Freeman, Reference Welper and Freeman1992). Moreover, many studies have reported associations between milk production traits and polymorphisms within the casein region (Martin et al. Reference Martin, Szymanowska, Zwierzchowski and Leroux2002; Nilsen et al. Reference Nilsen, Olsen, Hayes, Sehested, Svendsen, Nome, Meuwissen and Lien2009) and whey protein genes (Bovenhuis et al. Reference Bovenhuis, van Arendonk and Korver1992; Bleck & Bremel, Reference Bleck and Bremel1993a).
Given that the six main milk protein genes are directly involved in milk production and hence have been targeted for strong selection for increased milk yield in dairy cattle, we hypothesized that these genes show selection footprints associated with fertility traits. Indeed, the relationships between protein concentration or protein yield and fertility have been extensively investigated (e.g. Dematawewa & Berger, Reference Dematawewa and Berger1998) showing high and antagonistic genetic correlations. However, knowledge of the relationship between milk protein variants and cow fertility is limited and in some cases is contradictory. Demeter et al. (Reference Demeter, Markiewicz, van Arendonk and Bovenhuis2010) did not detect any significant relationship between β-LG, β-CN and κ-CN protein variants or β-κ-CN haplotypes with several fertility traits in Holstein dairy cattle. In contrast, earlier studies reported significant associations between β-CN variants and probability of conception (Hargrove et al. Reference Hargrove, Kiddy, Young, Hunter, Trimberger and Mather1980) and between β-LG variants and age at first conception and age at first calving (Lin et al. Reference Lin, McAllister, Ng-Kwai-Hang, Hayes, Batra, Lee, Roy, Vesely, Wauthy and Winter1987). Moreover, Ruottinen et al. (Reference Ruottinen, Ikonen and Ojala2004) detected a significant association between β-κ-CN haplotypes and the number of days from calving to first service in Finnish Ayrshire cows.
It is important to note that no studies have addressed the association between milk protein genes and fertilization rate and early embryonic survival, the main factors contributing to infertility in dairy cattle (Santos et al. Reference Santos, Thatcher, Chebel, Cerri and Galvao2004; Morris & Diskin, Reference Morris and Diskin2008). Recently, we have developed an in-vitro fertilization (IVF) experimental system in cattle which allowed the identification of several genetic factors involved in embryonic loss, and therefore infertility, at both genomic and expression levels (Khatib et al. Reference Khatib, Maltecca, Monson, Schutzkus, Wang and Rutledge2008a,Reference Khatib, Monson, Schutzkus, Kohl, Rosa and Rutledgeb; Khatib et al. Reference Khatib, Huang, Wang, Tran, Bindrim, Schutzkus, Monson and Yandell2009; Huang et al. Reference Huang, Kirkpatrick, Rosa and Khatib2010a,Reference Huang, Yandell and Khatibb; Huang & Khatib, Reference Huang and Khatib2010). Thus, the main objective of this study was to investigate the association between genetic variations in LALBA, LGB and casein genes and fertilization success and early embryonic development in Holstein cattle.
Materials and Methods
Fertility data collection
Ovaries from mature cows were collected from a local abattoir and immediately used in the IVF experiments as described in Khatib et al. (Reference Khatib, Maltecca, Monson, Schutzkus, Wang and Rutledge2008a,Reference Khatib, Monson, Schutzkus, Kohl, Rosa and Rutledgeb). A total of 6893 in-vitro fertilizations were performed and a total of 4661 embryos were produced using oocytes from 399 ovaries collected from 399 Holstein cows and semen samples from 12 Holstein bulls. For 92 ovaries, oocytes were fertilized by two different bulls each. Fertilization rate was calculated as the number of cleaved embryos 48 h post fertilization out of the total number of oocytes exposed to sperm. Blastocyst rate was calculated as the number of embryos that reached blastocyst stage out of the total number of embryos cultured.
Genotyping
DNA was extracted from ovaries using standard phenol/chloroform protocols. Sixty-six single nucleotide polymorphisms (SNPs), located in both coding and non-coding regions, were genotyped for the six main bovine milk protein genes: α-lactalbumin (LALBA), β-lactoglobulin (LGB) and the four caseins (CSN1S1, CSN2, CSN1S2 and CSN3). Forty-eight of the 66 SNPs are located in the casein region and were reported previously by Nilsen et al. (Reference Nilsen, Olsen, Hayes, Sehested, Svendsen, Nome, Meuwissen and Lien2009). Thirteen SNPs previously reported by Ganai et al. (Reference Ganai, Bovenhuis, van Arendonk and Visker2008) were genotyped in the LGB gene, and five SNPs were genotyped in LALBA gene (these five were previously reported by Bleck & Bremel, Reference Bleck and Bremel1993b; Kaminski et al. Reference Kaminski, Ahman, Rusc, Wójcik and Malewski2005; Martins et al. Reference Martins, Milazzotto, Feitosa, Coutinho, Simoes, Marques, Assumpcao and Visintin2008). A complete description of the SNPs, including the physical position and the marker allele frequencies of each SNP can be found in Supplementary Table 1. Genotyping of SNPs was performed at GeneSeek (Lincoln NE, USA).
Fertilization and blastocyst rates for genotypic classes of α-lactalbumin (LALBA) and β-lactoglobulin (LGB) genes
† NA: these genotypic classes were not considered in the association analysis because of the small number of records
Statistical analysis
Association between each SNP and fertilization and blastocyst rates was analysed using the following mixed linear model:
where y ijk represents in turn, the fertilization or blastocyst rate of oocytes k from ovary i fertilized with semen from bull j; μ represents a general mean for the trait considered, o i represents the random effect of the individual ovary from which oocytes were harvested; b j represents the random effect of the sire used in the fertilization; SNP ijk represents the fixed effect of the genotypic class (3 levels) for the SNP considered; and e ijk represents the residuals, assumed to be normal, independent, and identically distributed with mean 0 an variance Iσe2. Ovaries and bulls were assumed to be uncorrelated with variance structures Iσo2 and Iσb2, respectively. Association between the SNP and fertilization or blastocyst rate was tested using a likelihood ratio test by comparing with a reduced model without the SNP effect. All analyses were performed using the lme4 package of R language/environment (R Development Core Team, 2009).
In order to address multiple testing correction issues, we applied the simpleM method (Gao et al. Reference Gao, Starmer and Martin2008) which is an effective number of independent tests approach based on principal component analysis that takes into account linkage disequilibrium (LD) information among SNPs. First, a correlation matrix for the SNPs was constructed using the composite linkage disequilibrium. Then, eigenvalues were computed through principal component analysis on the composite LD matrix. Finally, the effective number of SNPs (Meff) used in the study were calculated as the number of principal components required to jointly explain 95% of the variance in the SNPs. P values from single SNP association tests were further adjusted for multiple comparisons using the Šidák correction based on Meff (adjusted P value=1−(1−P value)Meff).
Results and Discussion
In this study, we report the association between polymorphisms in LALBA and LGB genes, the two major whey proteins of bovine milk, with fertilization success and blastocyst rate in bovine embryos. To the best of our knowledge, this is the first study that reports associations between variants in LALBA and LGB genes with fertility evaluated in early embryonic development.
Out of 66 SNPs, four were found to be monomorphic in our population, and 8 SNPs had a minor allele frequency less than 0·05 and were not further considered in the study. Overall, 54 SNPs were analysed for association with fertilization and blastocyst rates. Figure 1 shows the linkage disequilibrium, measured with the correlation coefficient, r, among all the SNPs used in the analysis.
A heat map of linkage disequilibrium (r) between the 54 SNPs used in the association analysis. To facilitate the interpretation of the LD map we used only 3 colours where black indicates r⩾0·66, grey for 0·33<r<0·66, and white for r⩽0·33
Two SNPs, located in non-coding regions of the LALBA and LGB genes, showed significant associations with fertilization rate (Table 1). For SNP62, located 1691 bp upstream of the transcription start site in LALBA, oocytes collected from GG cows showed a fertilization rate of 67·0% compared with 58·7% for oocytes collected from GC cows (P=0·0007; adjusted P=0·011). For SNP61, which is located in intron 6 of LGB, oocytes from GG cows showed a fertilization rate of 66·6% compared with 57·4% obtained from oocytes from GA cows (P=0·0026; adjusted P=0·041). Interestingly, this SNP showed moderate linkage disequilibrium (r between 0·33 and 0·36) with other SNPs located in LGB included SNP58 and SNP59, non-synonymous mutations that differentiated the A and B variants of β-lactoglobulin (Fig. 1).
SNP49, located in the promoter region of LGB, showed evidence of association with embryonic blastocyst rate at a 10% significance level (Table 1). Embryos produced from GA cows showed a blastocyst rate of 42·1% compared with 29·9% for embryos produced from GG cows (P=0·0061; adjusted P=0·093). Importantly, this SNP did not show remarkable linkage disequilibrium with the other SNPs located on LGB gene (Fig. 1). Although not reaching statistical significance after adjustment for multiple testing, SNP17 located in CSN2 showed evidence of association with blastocyst rate (P=0·02).
Adjusting for multiple testing is a major issue in association studies owing to the correlation between many of the tests. Given the density of SNPs used in our study, many of these SNPs were found to be in high linkage disequilibrium (Fig. 1), therefore tests performed on each SNP are not independent. In this case, Bonferroni or Šidák corrections, popular approaches for controlling the experiment-wise error rate, would be overly conservative resulting in the overlooking of true positive signals. A valid adjustment for multiple testing must account for the correlation between tests. The simpleM method (Gao et al. Reference Gao, Starmer and Martin2008; Gao et al. Reference Gao, Becker, Becker, Starmer and Province2010) was specifically developed to calculate the number of informative SNPs being tested in a genetic association study using principal components analysis. In this context, the number of independent comparisons can be defined as the number of principal components accounting for a large portion of the variance in the SNP data. In this study, we estimated the effective number of tests (Meff=16) as the number of principal components necessary to account for 95% of the variance. Taking into account the small size of our phenotypic database and also that the simpleM approach still exhibits a trend of increasingly conservative adjustment for less significant P values (Pahl & Schäfer, Reference Pahl and Schäfer2010), we suggest that a cut-off of 0·10 for significant association is reasonable to control the experiment-wise error rate and also to minimize the incidence of false negative associations.
No significant associations were found between SNPs located in the casein region and fertilization or blastocyst rates. These results are in accordance with previous studies that did not detect significant associations between β-LG, β-CN and κ-CN protein variants and cow fertility traits such as calving interval, rate of calving after first insemination, and number of services per conception in Holsteins (Tsiaras et al. Reference Tsiaras, Bargouli, Banos and Boscos2005; Demeter et al. Reference Demeter, Markiewicz, van Arendonk and Bovenhuis2010).
Protein α-LA is essential for the biosynthesis of lactose in mammary gland (Bleck & Bremel, Reference Bleck and Bremel1993b). It interacts with UDP-galactosyltransferase to form lactose synthetase, which catalyses the synthesis of lactose from UDP-galactose and glucose. Lactose is the main osmole in milk and causes movement of water into the secretory vesicles of mammary epithelial cells. The concentration of α-LA is hypothesized to regulate milk volume through this mechanism by increasing or decreasing the amount of lactose secreted into the milk (Voelker et al. Reference Voelker, Bleck and Wheeler1997). In our study, SNP62, located in the 5′-untranslated region of LALBA, showed significant associations with fertilization rate. Variation in this region of the gene could potentially modify LALBA expression by altering the binding of RNA polymerase and various transcription factors that control gene expression. In fact, Lunden & Lindersson (Reference Lunden and Linderson1998) reported a relationship between milk lactose concentration and a mutation in the promoter of LALBA in the Swedish Red and White breed. Theoretically, an increase in lactose concentration could increase milk yield and therefore could be negatively associated with cow fertility. The possible relationship between lactose concentration, and hence milk yield, with SNP62 requires further research and could shed light on the biological mechanism underlying the significant association between the LALBA gene and fertilization success that we report here.
Protein β-LG is the major whey protein in bovine milk, although the biological functions of this gene are still not well understood. It could have a role in metabolism of phosphate in the mammary gland and the transport of retinol and fatty acids in the gut (Kontopidis et al. Reference Kontopidis, Holt and Sawyer2004). SNP49 in the promoter region and SNP61 in intron 6 of LGB showed significant associations with IVF fertility traits (Table 1). Polymorphisms in regulatory regions could affect the expression levels of the gene and therefore its biological function. In fact, Ganai et al. (Reference Ganai, Bovenhuis, van Arendonk and Visker2008) reported a significant association between SNP49 and the relative β-LG protein concentration in the milk of Holstein cows. They found that animals with the AA genotype had lower relative β-LG protein concentration compared with GG animals (Ganai et al. Reference Ganai, Bovenhuis, van Arendonk and Visker2008). It is well known that high concentration of β-LG protein is associated with low concentration of total casein proteins (Lunden et al. Reference Lunden, Nilsson and Janson1997); caseins represent about 80% of total protein in Holstein milk. Thus, it is reasonable to expect that higher concentration of β-LG protein caused by the G allele could be associated with lower concentration of total caseins and hence associated with an increase in milk production. Interestingly, in this study, we showed that embryos produced from GG cows showed a blastocyst rate of 29·9% compared with 42·1% for embryos produced from GA cows. Overall, these results could suggest an antagonistic relationship between milk yield and fertility for SNP49 in the LGB gene. Additional studies are needed to confirm and refine this relationship.
So far, there are no reports on direct involvement of the LALBA and LGB genes in pathways that affect the fertilization process or early embryonic development, although the involvement of these genes in milk production and associated measures has been well documented. It is widely accepted that a negative genetic correlation exists between milk yield traits and fertility in dairy cattle (Veerkamp & Beerda, Reference Veerkamp and Beerda2007). Basically, two different mechanisms may underlie this genetic association: pleiotropic gene effects and linkage disequilibrium. Mutations in many genes (e.g. PRLR, UTMP, GHR and STAT 5A) have been simultaneously associated with an increase in milk production in dairy cattle and a decrease in fertilization success and early embryonic survival (Khatib et al. Reference Khatib, Huang, Wang, Tran, Bindrim, Schutzkus, Monson and Yandell2009). On the other hand, it is well-known that recent events of selection, inbreeding and genetic drift may contribute substantially to the genetic relationships between two separate traits (Lewontin, Reference Lewontin1988). Recent studies have shown that many regions of the Holstein cattle genome, the highest-producing dairy animals in world, display signatures of recent selection (Hayes et al. Reference Hayes, Chamberlain, Maceachern, Savin, McPartlan, MacLeod, Sethuraman and Goddard2009; MacEachern et al. Reference MacEachern, Hayes, McEwan and Goddard2009; Qanbari et al. Reference Qanbari, Pimentel, Tetens, Thaller, Lichtner, Sharifi and Simianer2010). Moreover, these studies reported high values of linkage disequilibrium for some candidate regions harbouring major genes related to dairy quality. Although we cannot exclude a pleiotropic effect of LALBA and LGB genes on fertility traits of dairy cattle, the restructuring of the dairy cattle genome over the last decades caused by intense selection for production traits and a decrease in effective population size (i.e. inbreeding) may have resulted in a hitchhiking effect on a large number of loci affecting fertility and therefore our results could be explained also by the hitchhiking phenomenon.
Conclusion
Polymorphisms in LALBA and LGB genes showed significant associations with fertilization success and early embryonic development in Holstein cattle. No significant associations were detected for mutations located in the casein region. Although the molecular mechanisms underlying the association between whey protein genes and fertility have not yet been characterized, this study provides the first evidence of association between these genes and fertility traits evaluated in early embryonic development. Furthermore, these results may provide further evidence of the antagonistic relationship between milk production traits and reproductive performance.
This study was supported by USDA Hatch grant No. WIS-142-PRJ16JH from the University of Wisconsin-Madison. The authors thank Wen Huang for assistance with DNA preparation.
