Genetic effects and correlations between production and fertility traits and their dependency on the lactation-stage in Holstein Friesians
BMC Genetics volume 13, Article number: 108 (2012)
This study focused on the dynamics of genome-wide effects on five milk production and eight fertility traits as well as genetic correlations between the traits. For 2,405 Holstein Friesian bulls, estimated breeding values (EBVs) were used. The production traits were additionally assessed in 10-day intervals over the first 60 lactation days, as this stage is physiologically the most crucial time in milk production.
SNPs significantly affecting the EBVs of the production traits could be separated into three groups according to the development of the size of allele effects over time: 1) increasing effects for all traits; 2) decreasing effects for all traits; and 3) increasing effects for all traits except fat yield. Most of the significant markers were found within 22 haplotypes spanning on average 135,338 bp. The DGAT1 region showed high density of significant markers, and thus, haplotype blocks. Further functional candidate genes are proposed for haplotype blocks of significant SNPs (KLHL8, SICLEC12, AGPAT6 and NID1). Negative genetic correlations were found between yield and fertility traits, whilst content traits showed positive correlations with some fertility traits. Genetic correlations became stronger with progressing lactation. When correlations were estimated within genotype classes, correlations were on average 0.1 units weaker between production and fertility traits when the yield increasing allele was present in the genotype.
This study provides insight into the expression of genetic effects during early lactation and suggests possible biological explanations for the presented time-dependent effects. Even though only three markers were found with effects on fertility, the direction of genetic correlations within genotype classes between production and fertility traits suggests that alleles increasing the milk production do not affect fertility in a more negative way compared to the decreasing allele.
Time dependency in milk production traits in dairy cattle has been known and analyzed over the last century. Just like the phenotype, underlying genetic effects follow a dynamic expression over time. To account for the dynamic effects of genotypes, functional mapping has been introduced for the detection of QTLs, but has been applied mainly in human and plant genetics [1–4]. In livestock however, time dependency of traits is often accounted for when modeling genetic effects, but reported results are static in the sense of that cumulated 305-day breeding values are made public or that gene effects are given for a whole lactation. Thus, only scarce information about time dependent genetic effects in livestock is known so far. A few studies were carried out for pigs and sheep reporting dynamic QTLs for growth and weight [5–8]. Recently, it has been shown that the known effects of the K232A locus in the DGAT1 (diacylglycerol O-transferase 1) gene on milk yield and the protein production [9–11] are less pronounced or even reversed during the first 40 days of lactation . Furthermore, it was shown that genetic loci significant for the main milk production traits change from early, peak to late lactation .
The first weeks of lactation are a crucial time, especially in high yielding dairy cows. During that time, an energy deficiency manifests due to a drastic increase in milk production and a physiologically restricted energy intake. In those first weeks, the cow needs to draw energy from its adipose stores and in some cases even from its muscles, leading to a loss in body weight. The nutritional status of a cow after calving affects disease resistance and reproductive performance and regions affecting energy status were reported to overlap with regions significant for fertility traits . On one hand, for many populations, an unfavorable correlation between fertility traits and milk production has been reported due to the competition between these traits for the same body resources [15–19]. Fertility traits have a reportedly low heritability and a loss in fertility is so far mainly managed by optimizing environmental conditions . On the other hand, there are also reports questioning the genetic connection between high performance and a decline of fertility [21, 22].
The objectives of this study were to analyze how the genetic influence of genomic regions changes during the most critical interval between early and peak lactation and how changes of genetic effects in milk production could affect fertility. For this purpose, estimated breeding values (EBVs) for the five main milk production traits of 2,405 German Holstein Friesian bulls assessed in 10-day intervals over the first 60 lactation days, cumulated 305-day records and eight fertility traits were used. Association analyses were run between the production and fertility traits and the genome-wide markers of a 50 k SNP chip array. Finally, genetic correlations over the first 60 lactation days and within genotype classes were studied.
For the production traits, we identified 43 genome-wide significant markers for the 10-day intervals and 43 genome-wide significant markers for the 305-day records, of which 34 and 36 markers were located on chromosome 14 (Additional file 1: Table S1 and Additional file 2: Table S2). Seven markers located on chromosomes 6, 18 and 27 were unique for the 10-day intervals, affecting fat yield as well as fat and protein content mainly in the early 10-day intervals (Additional file 1: Table S1). Additionally seven markers located on chromosomes 5 and 14 were unique for the 305-day records, significant for fat content and protein yield. Whilst most markers on chromosome 14 were not significant once the effect of the DGAT1 locus was deducted, an additional ten and seven markers for the 10-day intervals and 305-day records, respectively, were discovered on chromosomes 5 and 27 (Additional file 3: Table S3). The allele effects were doubled compared to the effects when the DGAT1 locus was not deducted, which points to an overestimation for these marker effects. Therefore, given allele effects for those markers on chromosomes 5 and 27 were obtained without the DGAT1 locus in the model to ensure comparability between the sizes of found effects (Additional file 3: Table S3).
For the fertility traits, only three markers on chromosomes 6, 28 and 24 were significantly associated with non-return rate in heifers and overall fertility index, respectively (Additional file 4: Table S4). The DGAT1 locus had no influence on the significance or on effect sizes. None of the three significant fertility markers was significant for the production traits or located close to markers for production traits.
Based on whether allele effects increased or decreased over the first 60 lactation days, we divided the significant markers from the genome-wide association study (GWAS) in three groups:
Group 1 showed increasing effect sizes for all traits (Figure 1A). This group consisted of five markers located in a region between 0.05 and 2.6 Mb on chromosome 14 in which DGAT1 resides (at 50,872; 856,889; 1,307,998; 2,580,414; 2,607,583 bp; Figure 2). None of these five markers clustered together in a haplotype block. Increases in effect sizes of those SNP alleles were highly significant (P < 0.0001) for milk yield between all 10-day intervals as well as for fat content (P < 0.05), except for the first two 10-day intervals. Fat yield and protein content showed only significant increases between time points further apart than 10 days (Additional file 5: Table S5).
Group 2 was associated with decreasing allele effects for all traits across the first 60 days of lactation (Figure 1B). The eight markers of this group reside on chromosomes 6, 18 and 27 (Figure 2). Again, no clustering of markers within a single haplotype block was observed. Due to the low number of significant markers per trait from the GWAS, it was only possible to estimate LSM differences for protein content for markers on chromosome 6 (close to the Casein-Cluster) and 18. Significant differences (P < 0.05) were found between 10-day intervals 1 and 3, 4 and 5 as well as between intervals 2 and 5 (Additional file 5: Table S5). Other changes shown in Figure 1B indicate trends.
Group 3 showed an increase in effect size with later intervals apart from fat yield, where effect sizes became constantly smaller (Figure 1C). This group includes markers on chromosome 14 except the markers from group 1 (Figure 2). The majority of those markers were located in seven haplotype blocks (Table 1). Differences in allele effects for milk yield and fat content were similar to group 1. Fat yield and protein content showed significant differences in allele effects mostly in the first three intervals (Additional file 5: Table S5). For protein yield, not enough markers were significant in the GWAS to estimate LSMs, therefore, only a trend is indicated in Figure 1C.
The differences in allele effects between the 10-day intervals and the 305-day records were highly significant (P < 0.001) for fat and protein content with larger effects for the 305-day records for groups 1 and 3 (Additional file 6: Table S6). Differences between the 10-day intervals and the 305-day records were not significant for group 2. Because the EBVs for the 305-day records were 100x higher for milk yield and 10x higher for fat and protein yield compared to the EBVs for the 10-day intervals, differences for the yield traits were not considered informative.
On the basis of significant SNPs and linked protein coding genes, most significant results from a gene enrichment analysis using GENCODIS[23, 24] were related to chemokines, peptidases, antigens. The genomic regions connected to those gene ontology (GO) terms were located on chromosomes 6 (91 to 92 Mb), 14 (1 and 65 Mb), 18 (56 Mb) and 27 (39 Mb; Table 1).
To refine the genomic regions, haplotype blocks (HTBs) were constructed; 38 out of 53 significant markers (10-day intervals and 305-day records) were located in 22 haplotype blocks spanning on average 135,338 bp. Considering the markers that were significant when accounting for the DGAT1 effect, 45 out of 65 significant markers were found in 30 haplotype blocks (Additional file 7: Table S7). Eight HTBs including more than one significant marker were found on chromosome 14, and one on chromosome 5 (Table 2).
The DGAT1 gene region was the only one for milk production traits with a haplotype block harboring more than one significant marker. Most markers on chromosome 6 were located around the Casein-gene cluster (CSN1S1, CSN2, CSN1S2, and CSN3), however, the four significant markers are not within a single haplotype block and only CSN1S1 is located within the boundaries of a haplotype block. An additional marker on chromosome 6 is located 0.9 Mb upstream of the ABCG2 (ATP-binding cassette, subfamily G, member 2) gene. This marker had a very low MAF (<0.01) and was therefore not included in a haplotype block (Additional file 7: Table S7).
The highest density of protein coding genes (25 genes) was found for the haplotype block on chromosome 14 that included the DGAT1 gene. Another gene rich HTB was located on chromosome 18 containing 17 genes. On chromosome 18, the significant marker ARS-BFGL-NGS-109285 is a polymorphism within the SICLEC12 (sialic acid binding Ig-like lecithin 12) gene which might be responsible for changes in protein content in early lactation. Apart from an additional seven significant markers with polymorphisms in gene regions, a direct functional relation between these genes to the traits under investigation was not obvious.
Taking HTBs around significant SNPs and the next gene up- and down-stream of a haplotype block into account, the following genes with a reported or likely function on the traits under investigation can be proposed: the KLHL8 (Kelch-like-8) on chromosome 6 located 68.91 Kb down-stream of the next significant marker, CYP11B (cytochrome P450, subfamily XI B, polypeptide), which itself harbors a significant marker that does not contribute to a HTB on chromosome 14 at 1.29 Mb, AGPAT6 (1-acylglycerol-3-phosphate O-acyltransferase 6) on chromosome 27 located 79.38 Kb upstream of the next significant marker and NID1 (nidogen 1) on chromosome 28 located 303.97 Kb downstream of the next significant marker (Table 2). Despite the fact that several peaks of P-values were found for markers on chromosome 14 (Figure 2), no single candidate gene can be proposed for those regions on the basis of known gene annotations.
In general, yield traits were positively correlated with each other (0.3 to 0.8; P < 0.0001) and whilst milk yield was always negatively correlated with the content traits (−0.1 to −0.3; P < 0.0001), fat and protein yield were positively correlated with protein content (0.03 to 0.3; P < 0.007). The majority of fertility traits were highly genetically correlated with each other (0.79 to 0.94; P < 0.0001).
Referring to the 305-day records, negative correlations were found between all fertility traits and the yield traits. Fat content, however, showed low but positive correlations (0.03 to 0.08) though not always significant, and protein content showed a mix of positive and negative correlations, but equally low and only significant for non-return rate in heifers, calving to first insemination and days open (Table 3).
Over the first 60 lactation days, EBVs for fertility traits became more negatively correlated with EBVs for milk and protein yield, whilst fat yield showed smaller negative correlations except for calving to first insemination and days open (Figure 3). The content traits showed decreasing negative and increasing positive correlation over the first 60 lactation days. Fat content even showed a change from negative to positive genetic correlation between the third and fourth 10-day interval for all fertility traits except first to successful insemination in heifers (Figure 3).
For the analysis of genetic correlations within genotype classes, the most significant markers were used. On chromosome 14, analyses were restricted to the five SNPs with highest P-values at 0.44, 0.68, 0.89, 1.04 and 1.49 Mb; all except one marker were located in a haplotype block (Table 2, Figure 2). The minor allele frequency (MAF) of those markers ranged between 0.12 and 0.41. On one hand we found five markers, where the major allele was associated with an increase in EBVs for milk and protein yield and that showed a weaker negative correlation (on average −0.23 to −0.36) to most fertility traits when both copies of the allele were present. On the other hand, four markers, where the major allele was associated with a decrease in EBVs for milk and protein yield, showed a stronger negative correlation (on average −0.29 to −0.23) to most fertility traits when both copies of the allele were present (Table 2). If markers were significantly associated with content but not yield traits, allele effects for content traits were treated as opposing effects to milk yield due to the negative correlation between milk yield and content traits. Thus, two markers on chromosomes 5 (ARS-BFGL-NGS-116999, Hapmap49734-BTA-74577) and 6 (ARS-BFGL-NGS-112872, ARS-BFGL-NGS-118182), respectively, were exceptions to the above described pattern. For those markers, both copies of the major allele led to a stronger negative correlation between milk yield and fertility traits (on average −0.27 to −0.17) whilst the major allele was presumed to have an increasing effect on milk yield (Table 2). Correlations between EBVs for fat yield and fertility traits showed in most cases an opposite direction compared to the correlations for milk and protein yield, but fat yield also had an opposing allele effect compared to the other two yield traits.
For the content traits, most correlations within genotype classes were significant for calving to first insemination and days open. Those two fertility traits showed a stronger positive correlation (+0.05) if both copies of that allele were present that was associated with an increase in EBVs for fat or protein content. Exceptions were again the same two markers on chromosome 6 as well as one marker on chromosome 5 (Hapmap60021ss46526426) and 27 (ARS-BFGL-NGS-57448).
Exemplarily, the correlations within genotype classes between the production traits and the overall fertility index are shown for marker ARS-BFGL-NGS-57448 on chromosome 27 (Figure 4). Chromosome 27 was chosen due to the high MAF and the promising candidate gene AGPAT6 proposed for that region.
The maximum or minimum of the lactation curve for the different milk production traits has been stated to be between 28 and 56 days after calving in Holstein Friesian cattle [25, 26], which would be during the last three 10-day intervals for 31–60 days in milk (DIM) in our analysis.
We could identify changes in allele effects over the first 60 lactation days and compared to 305-day records. One reason for the change in allele effects could be that they only follow the change in EBVs due to the strong correlation (>0.99) in the response variable between consecutive intervals. However, the change in allele effects seems to be independent from the change in EBVs which is particularly obvious for milk yield, where the EBVs decreased over the first 60 lactation days but most of the allele effects increased. Furthermore, the variation in EBVs was fairly constant and only marginally increased for all traits over the first 60 lactation days, therefore presented comparisons of absolute allele effects over time is reasonable.
Three groups of markers were identified that showed different directions in allele effect changes. Group 1 consisted of five markers in the DGAT1 region between 0.05 and 2.6 Mb on chromosome 14 that showed increasing effect sizes for all traits during the first 60 days of lactation. The rest of 29 markers on chromosome 14 were regarded as group 3 with increasing effect sizes, apart from fat yield. DGAT1 is a major gene for a QTL on chromosome 14 with opposing effects on milk yield and milk fat [9, 10, 27, 28]. Even though it was shown that the enzyme activity of the DGAT1-protein can depend on the mutation variant [29, 30], it was not yet demonstrated that functional changes or quantity of expression is altered during the course of a lactation [29, 31, 32]. Therefore, other loci in that QTL region may explain the time-dependent changes of effect sizes in early lactation as well as the opposing development of allele effects for fat yield in group 1 and 3, as found in this study. For example, the CYP11B gene, which is located 1 Mb upstream of the DGAT1 gene and on the edge of a haplotype block of three significant markers, showed an opposite effect on fat yield compared to DGAT1. This opposing effect could decrease the overall effect of the DGAT1 region on fat yield, which could explain the herein reported decreasing effect size of markers associated with fat yield in group 3.
Overlapping with a haplotype close to CYP11B was also a region containing eight genes encoding proteins for the lymphocyte complex. Markers in that region were significant for the first 60 lactation days as well as for the 305-day records and strongly affecting most milk production traits. This region could therefore play a role in the immune response of the mammary gland and prevents inflammation during lactation.
Group 2 included markers on chromosomes 6, 18 and 27. The significant markers showed decreasing effects for all traits, which were significant for protein content for markers on chromosomes 6 and 18. The markers on chromosome 6 in the Casein-gene cluster were only significant during the first 60 lactation days with largest effects at the beginning of the lactation. This co-occurrence with changes in gene expression patterns of the casein genes during early lactation . The same pattern of changes of genetic effects was found for a marker on chromosome 18 that is located in the SICLEC12 gene and linked to SICLEC genes 10 and 14 which are all involved in immunoglobulin production. The co-occurrence of changes in allele effects over time of markers in the Casein-gene cluster and the SICLEC genes can be explained by their biological functions. In the first days after calving, when the colostrum is produced, milk is rich in proteins, especially in immunoglobulins, which have bioactive functions and help to activate the calf’s immune system [34, 35]. Furthermore, other genes of the SICLEC family have been reported to affect productive live and fertility in Holstein cattle and are linked to leptin signaling [36, 37]. Leptin plays a major role in regulating foot intake and thus the energy status of a cow which could affect both fertility and productivity [38, 39].
The significant region on chromosome 27, with decreasing effects on fat yield and content during the first 60 days of lactation, was supported by only one marker in the initial model, but five more markers when accounting for the DGAT1 effect. The 1-acylglycerol-3-phosphate O-acyltransferase 6 gene (AGPAT6) is located 79.38 Kb upstream of the initially only significant marker on the edge of a haplotype block. AGPAT6 has a similar function as DGAT1. It is an enzyme in the phospholipid and triglyceride biosynthesis, and thus, contributes to the production of milk fat . The expression of AGPAT6 in the mammary gland was shown to increase drastically over the first 60 lactation days and decreases afterwards . The change in the expression profile is not consistent with the herein reported trend of decreasing allele effects over the first 60 days of lactations but the fact that markers in the AGPAT6 region were only significant at the beginning of lactation showed that the impact on milk fat must be diminishing in late lactation.
Additionally, the GO terms chemokines and peptidases occurred significantly often for our marker set and especially in the region of markers on chromosome 6, 18 and 27 which were only significant in early lactation. Chemokines attract lymphocytes and macrophages and are thus related to immune response, either to help the offspring or prevent infections of the udder itself in a time of high activity and metabolic stress . Something similar can be said about peptidases, especially the kallikrein-related peptidases that were predominant in the region on chromosome 18; connection to immune response and also associations with breast cancer have been reported [42–44]. Chemokines and peptidases however, show no direct link to the traits that were investigated in this study.
Not included in the three described effect groups were markers on chromosomes 5. Only one marker was initially significant for the last 10-day interval, and thus, no change in effect size could be determined. However, five more markers for the 10-day intervals and six markers for the 305-day records became significant when accounting for the DGAT1 effect. The markers had a tendency towards increasing allele effects for fat yield and content over the first 60 lactation days. This, and the fact that initially most of the markers on chromosome 5 were significant for the 305-day records concurs with a decrease in fat content when energy availability is low during early lactation .
Only a few studies have focused on time-dependent genetic associations in livestock to date and, as this study shows, the investigation of association at certain lactation stages seem to be a promising approach to detect loci with only small overall effects [13, 46–48]. Thus, significant associations can sometimes only be found when the phenotype was recorded at the right time. Because the impact of loci changes over time, an investigation at certain developmental or lactation stages might be a contribution to detect parts of the otherwise missed genetic variance and may be better in detecting quantitative trait genes with overall small effect.
Finally, we want to propose two candidate genes for markers detected for the non-return rate in heifers (NRH). The first gene is the kelch-like 8 (KLHL8), located on chromosome 6 within 68.91 Kb of the nearest significant marker. The KLHL8 gene was reported to be preferentially expressed in the female gonads of fish (zebrafish and carbiomedaka) where it may play a role in oogenesis . Even though fish are rather distantly relate to cattle, the gene function might be evolutionary conserved. The second gene is nidogen 1 (NID1) on chromosome 28 and within 303.97 Kb of the nearest significant marker. NID1 is increasingly expressed in the focimatrix (follicular basal laminas) when follicles enlarge before ovulation . Thus, both genes, KLHL8 and NID1, could affect the non-return rate by regulating the ovulation process.
An antagonistic relation between production and fertility traits reflects the competition between these traits for the same body resources [16, 17]. We confirm previous studies reporting an unfavorable genetic correlation between yield and fertility traits [17, 51–53]. The EBVs for content traits were mostly positively correlated with the EBVs for fertility which might be due to a spurious effect resulting from both content and fertility traits, being negatively correlated with milk yield. Additionally, we report an increase of the negative genetic correlation over the first 60 lactation days for the EBVs of milk and protein yield and smaller negative and larger positive correlations for the EBVs of content traits. The stronger negative correlation between EBVs for yield and fertility traits could result from the also increasing milk production in early lactation, depleting important body resources needed for most fertility traits.
Additionally to the time dependency of genetic correlation between production and fertility traits, we reported differences for genetic correlation within genotypes of the most significant markers from the GWAS. With only few exceptions, alleles that significantly increased the EBV for milk yield showed a weaker negative correlation with most of the fertility traits if this allele was present in the genotype. Though for different markers compared to this study, Pimentel et al. (2011) reported SNPs with favorable effects on production and fertility traits and came to the same conclusion that a selection for higher performance based on the genotypes of certain loci does not further influence the fertility negatively .
The content traits showed mainly stronger positive correlations with most of the fertility traits if the yield decreasing allele was present in the genotype. Again, a spurious effect due to the negative correlation of content traits as well as fertility traits with milk yield cannot be ruled out.
This study provided evidence for new loci and confirmed less known loci with effects on milk production traits. Furthermore, significant associations can sometimes only be found when the phenotype was recorded at the right time. Because the impact of loci changes over time, an investigation at certain developmental or lactation stages might be a contribution to detect parts of the otherwise missed genetic variance and may be better in detecting quantitative trait genes with overall small effect. By analyzing the first 60 lactation days more closely, we could show that changes in effect sizes mirror results from gene expression studies and might, therefore, be a less elaborate method in addition to determine the consequences of gene expression. Lactation-stage specific association studies provide information independently of knowing where the genes exert their activity in the body.
Correlations between fertility and production traits, positive or negative, became stronger with progressing lactation. Depending on the genotype, the correlations between traits differed. The correlation was less negative if the yield increasing allele of a candidate marker was present in the genotype. If these results can be confirmed in further studies, there is a good chance that selection for higher performance based on the here presented marker genotypes would not increase the negative influence on fertility.
Estimated breeding values for the five main production traits milk yield (MY), fat and protein yield (FY, PY), and fat and protein content (FC, PC) of 2,405 German Holstein Friesian bulls were averaged over the first three lactations. All records were provided by VIT (Vereinigte Informationssysteme Tierhaltung, Germany) and the production EBVs are raw breeding values not normalized to the average breeding value of a birth cohort.
On average, each bull had 752 daughters contributing to the breeding value estimation, with a median of 114, a maximum of 85,393 and a minimum of 50 (1. Quartile = 93; 3. Quartile = 116). Data for the production traits were given as cumulated 305-day records (Additional file 8: Table S8) and cumulated 10-day intervals (1–5; interval 1 = day 11–20) over the first 60 lactation days. No accuracies for the 10-day interval EBVs were available; instead, the number of records per sire for each 10-day interval is given in Table 4. The heritabilities of the production traits are generally reported to be 0.25 for yield traits and 0.5 for content traits. Even though it was reported that heritabilies increase with DIM, the amount of daughter information and the stable record number between the intervals let us assume that accuracies for the EBVs are high and possible changes in heritabilites did not affect our results. The EBVs for milk yield and fat and protein content were decreasing and for fat and protein yield increasing over the first 60 lactation days (Additional file 9: Figure S1). The standard deviation only marginally increased with later intervals.
For fertility traits, breeding values were available for the overall fertility index (RZR; summarizing all following traits), interval from first to successful insemination classified into cows and heifers (FLc and FLh), non-return rate to 56 days for cows and heifers (NRc and NRh), conception (CON; summarizing first to successful insemination and non-return rate), interval from calving to first insemination for cows (CFc) and days open (DO; Additional file 10: Table S9). Heritabilities for fertility traits range between 0.01-0.1 .
The genotypic information for all bulls was obtained from the bovine 50 k BeadChip (Illumina Inc., San Diego, California, USA). The chip features 54,001 SNPs. The array information was remapped against the Btau 4.0 assembly according to and cleaned of duplicated and not yet allocated markers (leaving 51,983 markers) . 2,339 animals and 43,628 markers passed quality control (minor allele frequency >0.01, call rate >0.9, IBS < =0.95) resulting in an average mean autosomal heterozygosity of 0.30 per SNP, a mean identity by state of 0.72 and a minor allele frequency (MAF) of 0.26 in our data set.
A genome wide association study was conducted with the GenABEL package in R  where we treated the 10-day intervals and the 305-day EBVs as separated traits. Since the strong effect of the DGAT1 region allowed the detection of only a few additional markers outside that QTL, a supplemental analysis was run using the marker ARS-BFGL-NGS-4939, which is closest to the gene, as a fixed effect and thus, to deduct its effects. This marker had no effects on the analysis of fertility traits. A Genome-wide Rapid Analysis using Mixed Models And Score test (GRAMMAS), as implemented in GenABEL (c.f. GenABLE Tutorial, 29.11.2009), was used to account for population stratification based on the genomic kinship matrix [57, 58]. The following polygenic model was used with and without the DGAT1 locus in the model:
where y ijkl are the EBVs for 10-day intervals, 305-day records or fertility traits, respectively; μ is the mean of the population; G i is the polygenic contribution of individual i; s j is the fixed effect of SNP j; DGAT1 j is the fixed effect of the genotype j (AA, KA, and KK) of marker ARS-BFGL-NGS-4939 (only in a supplemental analysis); and e ijkl is a vector of random residuals.
Beside GRAMMAS, the lambda factor was used on the P-values to achieve a non-conservative test statistic. Lambda is defined as the ratio between observed and expected median of the chi2-distribution of the test statistic. Finally, Bonferroni correction was applied to account for multiple testing with α = 0.05, resulting in an adjusted P-value of 1.15-6 (−log10 P-value = 5.94).
Based on the significant markers in the GWAS, least square means differences between allele effects of each 10-day interval as well as between allele effects of the 10-day intervals and 305-day records were estimated with the following model in SAS (version 9.2; SAS Institute Inc., Cary, North Carolina, USA):
with y ijk being the allele effects of significant markers estimated from the GWAS (Additional file 1: Table S1); μ is the mean of the selected marker effects; I i is the 10-day interval i (1–5); m j is a random effect of the marker j to account for the different number of markers and possible relations between the markers; and e ijk is a vector of random residuals. The test of LSM differences between the 10-day intervals was performed for each effect group that was found in the GWAS separately.
The genetic correlations between all traits were estimated as Pearson’s correlation coefficients between the EBVs. Additionally, genetic correlations within genotype classes of markers with most significant P-values within a haplotype block (HTB) were used.
Haplotype blocks over all markers featured on the chip array were derived using the software fastPhase with filtering criteria of <3% missing genotype calls, <5% missing SNP calls, MAF >5%, 10 random starts and 25 iterations , and HAPLOVIEW v4.1  using the solid spine algorithm implemented in the software. A haplotype block was defined if the first and last markers of a region were in linkage disequilibrium (D’ > 0.8) with all intermediate markers whereas the latter had not necessarily to be in LD with each other.
Protein coding genes within the haplotype blocks were taken from the Btau 4.0 assembly. GENCODIS was used for a gene enrichment analysis to suggest possible functional candidate genes within HTBs [23, 24]. GENCODIS integrates different sources of information to search for annotations that frequently co-occur in a set of genes and rank them by their statistical significance.
Yang R, Gao H, Wang X, Zhang J, Zeng ZB, Wu R: A semiparametric approach for composite functional mapping of dynamic quantitative traits. Genetics. 2007, 177 (3): 1859-1870. 10.1534/genetics.107.077321.
Li N, Das K, Wu R: Functional mapping of human growth trajectories. Journal of theoretical biology. 2009, 261 (1): 33-42. 10.1016/j.jtbi.2009.07.020.
Wu R, Lin M: Functional mapping - how to map and study the genetic architecture of dynamic complex traits. Nature reviews Genetics. 2006, 7 (3): 229-237.
Min L, Yang R, Wang X, Wang B: Bayesian analysis for genetic architecture of dynamic traits. Heredity. 2011, 106 (1): 124-133. 10.1038/hdy.2010.20.
Lopez S, France J, Gerrits WJ, Dhanoa MS, Humphries DJ, Dijkstra J: A generalized Michaelis-Menten equation for the analysis of growth. Journal of animal science. 2000, 78 (7): 1816-1828.
Schinckel AP, Ferrel J, Einstein ME, Pearce SM, Boyd RD: Analysis of pig growth from birth to sixty days of age. Professional Animal Scientist. 2004, 20: 79-86.
Lambe NR, Navajas EA, Simm G, Bunger L: A genetic investigation of various growth models to describe growth of lambs of two contrasting breeds. Journal of animal science. 2006, 84 (10): 2642-2654. 10.2527/jas.2006-041.
Hadjipavlou G, Bishop SC: Age-dependent quantitative trait loci affecting growth traits in Scottish Blackface sheep. Animal genetics. 2009, 40 (2): 165-175. 10.1111/j.1365-2052.2008.01814.x.
Winter A, Kramer W, Werner FA, Kollers S, Kata S, Durstewitz G, Buitkamp J, Womack JE, Thaller G, Fries R: Association of a lysine-232/alanine polymorphism in a bovine gene encoding acyl-CoA:diacylglycerol acyltransferase (DGAT1) with variation at a quantitative trait locus for milk fat content. Proceedings of the National Academy of Sciences of the United States of America. 2002, 99 (14): 9300-9305. 10.1073/pnas.142293799.
Thaller G, Kramer W, Winter A, Kaupe B, Erhardt G, Fries R: Effects of DGAT1 variants on milk production traits in German cattle breeds. Journal of animal science. 2003, 81 (8): 1911-1918.
Rahmatalla S, Muller U, Strucken EM, Brockmann GA: Der Effekt von DGAT1-Genvarianten in Deutschen Holstein Kühen unter Produktionsbedingungen. Züchtungskunde. 2008, 80: 473-484.
Strucken EM, de Koning DJ, Rahmatalla SA, Brockmann GA: Lactation curve models for estimating gene effects over a timeline. Journal of dairy science. 2011, 94 (1): 442-449. 10.3168/jds.2009-2932.
Strucken EM, Bortfeldt RH, De Koning DJ, Brockmann GA: Genome wide association for investigating time dependent genetics effects for milk production traits in dairy cattle. Animal genetics. 2011, 10.1111/j.1365-2052.2011.02278.x.
Tetens J, Seidenspinner T, Buttchereit N, Thaller G: Whole genome association study for energy balance and fat protein ratio in German Holstein bull dams. Animal genetics. 2012, in press
Dematawewa CM, Berger PJ: Genetic and phenotypic parameters for 305-day yield, fertility, and survival in Holsteins. Journal of dairy science. 1998, 81 (10): 2700-2709. 10.3168/jds.S0022-0302(98)75827-8.
Rauw WM, Kanis E, Noordhuizen-Strassen EN, Grommers FJ: Undesirable side effects of selection for high production efficiency in farm animals: a review. Livest Prod Sci. 1998, 56: 15-33. 10.1016/S0301-6226(98)00147-X.
Buch LH, Sorensen MK, Lassen J, Berg P, Jakobsen JH, Johansson K, Sorensen AC: Udder health and female fertility traits are favourably correlated and support each other in multi-trait evaluations. Journal of animal breeding and genetics = Zeitschrift fur Tierzuchtung und Zuchtungsbiologie. 2011, 128 (3): 174-182. 10.1111/j.1439-0388.2010.00904.x.
Dunn TG, Ingalls JE, Zimmerman DR, Wiltbank JN: Reproductive performance of 2-year-old Hereford and Angus heifers as influenced by pre- and post-calving energy intake. Journal of animal science. 1969, 29 (5): 719-726.
Van Niekerk A: The effect of body condition as influenced by winter nutrition, on the reproductive performance of the beef cow. South African Journal of Animal Science. 1982, 12: 383-387.
Tsuruta S, Misztal I, Huang C, Lawlor TJ: Bivariate analysis of conception rates and test-day milk yields in Holsteins using a threshold-linear model with random regressions. Journal of dairy science. 2009, 92 (6): 2922-2930. 10.3168/jds.2008-1707.
LeBlanc S: Assessing the association of the level of milk production with reproductive performance in dairy cattle. The Journal of reproduction and development. 2010, 56 (Suppl): S1-S7.
Pimentel EC, Bauersachs S, Tietze M, Simianer H, Tetens J, Thaller G, Reinhardt F, Wolf E, Konig S: Exploration of relationships between production and fertility traits in dairy cattle via association studies of SNPs within candidate genes derived by expression profiling. Animal genetics. 2011, 42 (3): 251-262. 10.1111/j.1365-2052.2010.02148.x.
Carmona-Saez P, Chagoyen M, Tirado F, Carazo JM, Pascual-Montano A: GENECODIS: a web-based tool for finding significant concurrent annotations in gene lists. Genome biology. 2007, 8 (1): R3-10.1186/gb-2007-8-1-r3.
Nogales-Cadenas R, Carmona-Saez P, Vazquez M, Vicente C, Yang X, Tirado F, Carazo JM, Pascual-Montano A: GeneCodis: interpreting gene lists through enrichment analysis and integration of diverse biological information. Nucleic acids research. 2009, 37 (Web Server issue): W317-W322.
Stanton TL, Jones LR, Everett RW, Kachman SD: Estimating milk, fat, and protein lactation curves with a test day model. Journal of dairy science. 1992, 75 (6): 1691-1700. 10.3168/jds.S0022-0302(92)77926-0.
Pollott GE: Deconstructing milk yield and composition during lactation using biologically based lactation models. Journal of dairy science. 2004, 87 (8): 2375-2387. 10.3168/jds.S0022-0302(04)73359-7.
Grisart B, Coppieters W, Farnir F, Karim L, Ford C, Berzi P, Cambisano N, Mni M, Reid S, Simon P, et al: Positional candidate cloning of a QTL in dairy cattle: identification of a missense mutation in the bovine DGAT1 gene with major effect on milk yield and composition. Genome research. 2002, 12 (2): 222-231. 10.1101/gr.224202.
Kuehn C, Edel C, Weikard R, Thaller G: Dominance and parent-of-origin effects of coding and non-coding alleles at the acylCoA-diacylglycerol-acyltransferase (DGAT1) gene on milk production traits in German Holstein cows. BMC genetics. 2007, 8: 62-
Grisart B, Farnir F, Karim L, Cambisano N, Kim JJ, Kvasz A, Mni M, Simon P, Frere JM, Coppieters W, et al: Genetic and functional confirmation of the causality of the DGAT1 K232A quantitative trait nucleotide in affecting milk yield and composition. Proceedings of the National Academy of Sciences of the United States of America. 2004, 101 (8): 2398-2403. 10.1073/pnas.0308518100.
Smith SJ, Cases S, Jensen DR, Chen HC, Sande E, Tow B, Sanan DA, Raber J, Eckel RH, Farese RV: Obesity resistance and multiple mechanisms of triglyceride synthesis in mice lacking Dgat. Nature genetics. 2000, 25 (1): 87-90. 10.1038/75651.
Furbass R, Winter A, Fries R, Kuhn C: Alleles of the bovine DGAT1 variable number of tandem repeat associated with a milk fat QTL at chromosome 14 can stimulate gene expression. Physiological genomics. 2006, 25 (1): 116-120. 10.1152/physiolgenomics.00145.2005.
Ron M, Israeli G, Seroussi E, Weller JI, Gregg JP, Shani M, Medrano JF: Combining mouse mammary gland gene expression and comparative mapping for the identification of candidate genes for QTL of milk production traits in cattle. BMC genomics. 2007, 8: 183-10.1186/1471-2164-8-183.
Kaupe B, Brandt H, Prinzenberg EM, Erhardt G: Joint analysis of the influence of CYP11B1 and DGAT1 genetic variation on milk production, somatic cell score, conformation, reproduction, and productive lifespan in German Holstein cattle. Journal of animal science. 2007, 85 (1): 11-21. 10.2527/jas.2005-753.
Le A, Barton LD, Sanders JT, Zhang Q: Exploration of bovine milk proteome in colostral and mature whey using an ion-exchange approach. Journal of proteome research. 2011, 10 (2): 692-704. 10.1021/pr100884z.
Farrell HM, Jimenez-Flores R, Bleck GT, Brown EM, Butler JE, Creamer LK, Hicks CL, Hollar CM, Ng-Kwai-Hang KF, Swaisgood HE: Nomenclature of the proteins of cows’ milk–sixth revision. Journal of dairy science. 2004, 87 (6): 1641-1674. 10.3168/jds.S0022-0302(04)73319-6.
Cole JB, VanRaden PM, O’Connell JR, Van Tassell CP, Sonstegard TS, Schnabel RD, Taylor JF, Wiggans GR: Distribution and location of genetic effects for dairy traits. Journal of dairy science. 2009, 92 (6): 2931-2946. 10.3168/jds.2008-1762.
Patel N, der Brinkman-Van Linden EC, Altmann SW, Gish K, Balasubramanian S, Timans JC, Peterson D, Bell MP, Bazan JF, Varki A, et al: OB-BP1/Siglec-6. a leptin- and sialic acid-binding protein of the immunoglobulin superfamily. The Journal of biological chemistry. 1999, 274 (32): 22729-22738. 10.1074/jbc.274.32.22729.
Banos G, Woolliams JA, Woodward BW, Forbes AB, Coffey MP: Impact of single nucleotide polymorphisms in leptin, leptin receptor, growth hormone receptor, and diacylglycerol acyltransferase (DGAT1) gene loci on milk production, feed, and body energy traits of UK dairy cows. Journal of dairy science. 2008, 91 (8): 3190-3200. 10.3168/jds.2007-0930.
Liefers SC, Veerkamp RF, te Pas MF, Delavaud C, Chilliard Y, Platje M, van der Lende T: Leptin promoter mutations affect leptin levels and performance traits in dairy cows. Animal genetics. 2005, 36 (2): 111-118. 10.1111/j.1365-2052.2005.01246.x.
Beigneux AP, Vergnes L, Qiao X, Quatela S, Davis R, Watkins SM, Coleman RA, Walzem RL, Philips M, Reue K, et al: Agpat6–a novel lipid biosynthetic gene required for triacylglycerol production in mammary epithelium. Journal of lipid research. 2006, 47 (4): 734-744. 10.1194/jlr.M500556-JLR200.
Bionaz M, Loor JJ: ACSL1, AGPAT6, FABP3, LPIN1, and SLC27A6 are the most abundant isoforms in bovine mammary tissue and their expression is affected by stage of lactation. The Journal of nutrition. 2008, 138 (6): 1019-1024.
Clarkson RW, Wayland MT, Lee J, Freeman T, Watson CJ: Gene expression profiling of mammary gland development reveals putative roles for death receptors and immune mediators in post-lactational regression. Breast cancer research: BCR. 2004, 6 (2): R92-R109. 10.1186/bcr754.
Angelopoulou K, Karagiannis GS: Structural characterization and expression of five novel canine kallikrein-related peptidases in mammary cancer. Mammalian genome: official journal of the International Mammalian Genome Society. 2010, 21 (9–10): 516-524.
Barret AJ, Rawlings ND, Woessner JF: Handbook of proteolytic Enzymes. 2004, Academic, San Diego
de Vries MJ, Veerkamp RF: Energy balance of dairy cattle in relation to milk production variables and fertility. Journal of dairy science. 2000, 83 (1): 62-69. 10.3168/jds.S0022-0302(00)74856-9.
Lund MS, Sorensen P, Madsen P, Jaffrezic F: Detection and modelling of time-dependent QTL in animal populations. Genetics, selection, evolution: GSE. 2008, 40 (2): 177-194. 10.1186/1297-9686-40-2-177.
Forni S, Gianola D, Rosa GJ, de Los Campos G: A dynamic linear model for genetic analysis of longitudinal traits. Journal of animal science. 2009, 87 (12): 3845-3853. 10.2527/jas.2008-1514.
Suchocki T, Szyda J: Statistical modelling of growth using a mixed model with orthogonal polynomials. Journal of applied genetics. 2011, 52 (1): 95-100. 10.1007/s13353-010-0012-2.
Gautier A, Le Gac F, Lareyre JJ: The gsdf gene locus harbors evolutionary conserved and clustered genes preferentially expressed in fish previtellogenic oocytes. Gene. 2011, 472 (1–2): 7-17.
Irving-Rodgers HF, Harland ML, Rodgers RJ: A novel basal lamina matrix of the stratified epithelium of the ovarian follicle. Matrix biology: journal of the International Society for Matrix Biology. 2004, 23 (4): 207-217. 10.1016/j.matbio.2004.05.008.
Windig JJ, Calus MP, Beerda B, Veerkamp RF: Genetic correlations between milk production and health and fertility depending on herd environment. Journal of dairy science. 2006, 89 (5): 1765-1775. 10.3168/jds.S0022-0302(06)72245-7.
Gonzalez-Recio O, Alenda R, Chang YM, Weigel KA, Gianola D: Selection for female fertility using censored fertility traits and investigation of the relationship with milk production. Journal of dairy science. 2006, 89 (11): 4438-4444. 10.3168/jds.S0022-0302(06)72492-4.
Ikonen T, Morri S, Tyriseva AM, Ruottinen O, Ojala M: Genetic and phenotypic correlations between milk coagulation properties, milk production traits, somatic cell count, casein content, and pH of milk. Journal of dairy science. 2004, 87 (2): 458-467. 10.3168/jds.S0022-0302(04)73185-9.
VanRaden PM, Sanders AH, Tooker ME, Miller RH, Norman HD, Kuhn MT, Wiggans GR: Development of a national genetic evaluation for cow fertility. Journal of dairy science. 2004, 87 (7): 2285-2292. 10.3168/jds.S0022-0302(04)70049-1.
Schmitt AO, Bortfeldt RH, Brockmann GA: Tracking chromosomal positions of oligomers - a case study with Illumina’s BovineSNP50 beadchip. BMC genomics. 2010, 11: 80-10.1186/1471-2164-11-80.
Aulchenko YS, Ripke S, Isaacs A, van Duijn CM: GenABEL: an R library for genome-wide association analysis. Bioinformatics. 2007, 23 (10): 1294-1296. 10.1093/bioinformatics/btm108.
Aulchenko YS, de Koning DJ, Haley C: Genomewide rapid association using mixed model and regression: a fast and simple method for genomewide pedigree-based quantitative trait loci association analysis. Genetics. 2007, 177 (1): 577-585. 10.1534/genetics.107.075614.
Amin N, van Duijn CM, Aulchenko YS: A genomic background based method for association analysis in related individuals. PloS one. 2007, 2 (12): e1274-10.1371/journal.pone.0001274.
Scheet P, Stephens M: A fast and flexible statistical model for large-scale population genotype data: applications to inferring missing genotypes and haplotypic phase. American journal of human genetics. 2006, 78 (4): 629-644. 10.1086/502802.
Barrett JC, Fry B, Maller J, Daly MJ: Haploview: analysis and visualization of LD and haplotype maps. Bioinformatics. 2005, 21 (2): 263-265. 10.1093/bioinformatics/bth457.
The project and R.H. Bortfeldt were supported by GenoTrack (FUGATO). We thank A. Bünger and F. Reinhardt (VIT, Verden, Germany) for providing the EBVs.
EMS conceived of the study, carried out statistical analysis of the data and drafted the manuscript. RHB carried out the position correction of SNP data and provided the haplotype blocks. JT and GT provided the genotypes. GAB participated in conceiving of the study, its design and coordination and helped drafting the manuscript. All authors read and approved the final manuscript.
Electronic supplementary material
Additional file 1 Table S1: Complete list of significant markers over the first 60 lactation days separated after 10-day intervals. DIM: days in milk; MY: milk yield; FY: fat yield; PY: protein yield; FC: fat content; PC: protein conten. (DOC 610 KB)
Additional file 2 Table S2: Complete list of markers significant for the 305-day records. MY: milk yield; FY: fat yield; PY: protein yield; FC: fat content; PC: protein content. (DOC 168 KB)
Additional file 3 Table S3: Additional markers detected when accounting for the DGAT1 locus. DIM: days in milk; MY: milk yield; FY: fat yield; PY: protein yield; FC: fat content; PC: protein content. Effects, SD and P-values are given as the results from the analysis without DGAT1 in the model and in brackets with DGAT1 in the model. (DOC 60 KB)
Additional file 4 Table S4: Significant markers and allele effects for fertility traits. RZR: fertility index, summarizing all fertility traits, NRh: non-return rate for heifers. (DOC 26 KB)
Additional file 5 Table S5: LSM-differences between the first 60 lactation days separated after effect groups. P-values *** < 0.0001 ** < 0.001 * < 0.01 † < 0.05. LSM: least square means; DIM: days in milk; MY: milk yield; FY: fat yield; FC: fat content; PC: protein content. (DOC 40 KB)
Additional file 6 Table S6: LSM- differences between the 10-day intervals over the first 60 lactation days and 305-day records after effect groups. P-values *** < 0.0001 ** < 0.001 * < 0.01 † < 0.05. LSM: least square means; DIM: days in milk; FC: fat content; PC: protein content. (DOC 31 KB)
Additional file 7 Table S7: Haplotype blocks, marker and candidate gene information over all approaches. HTB: haplotype block; DIM: days in milk; MAF: minor allele frequency; MY: milk yield; FY: fat yield; PY: protein yield; FC: fat content; PC: protein content; RZR: fertility index, summarizing all fertility traits; NRh: non-return rate for heifers; ns: non-significant. * only significant when DGAT1 region was deducted; allele effects are given as original results due to the bias when a major locus is deducted. (DOC 93 KB)
Additional file 8 Table S8: EBVs for 305-day production traits. EBV: estimated breeding value; MY: milk yield; FY: fat yield; PY: protein yield; FC: fat content; PC: protein content. (DOC 32 KB)
Additional file 9 Figure S1: Average EBVs over first 60 lactation days. EBV: estimated breeding value; filled triangle milk yield; filled circle fat yield; filled square protein yield; empty circle fat content; empty square protein content (axis for milk yield is on the right hand side). (DOC 52 KB)
Additional file 10 Table S9: EBVs for fertility traits. EBV: estimated breeding value; RZR: fertility index, summarizing all following traits, FLc and FLh: 1st to successful insemination separated after cows and heifers, NRc and NRh: non-return rate for cows and heifers, CON: conception (summarizing FLc and FLh, NRc and NRh), CFc: calving to first insemination; DO: and days open. (DOC 32 KB)
About this article
Cite this article
Strucken, E.M., Bortfeldt, R.H., Tetens, J. et al. Genetic effects and correlations between production and fertility traits and their dependency on the lactation-stage in Holstein Friesians. BMC Genet 13, 108 (2012). https://doi.org/10.1186/1471-2156-13-108