BDNF Contributes to the Genetic Variance of Milk Fat Yield in German Holstein Cattle

The gene encoding the brain-derived neurotrophic factor (BDNF) has been repeatedly associated with human obesity. As such, it could also contribute to the regulation of energy partitioning and the amount of secreted milk fat during lactation, which plays an important role in milk production in dairy cattle. Therefore, we performed an association study using estimated breeding values (EBVs) of bulls and yield deviations of German Holstein dairy cattle to test the effect of BDNF on milk fat yield (FY). A highly significant effect (corrected p-value = 3.362 × 10−4) was identified for an SNP 168 kb up-stream of the BDNF transcription start. The association tests provided evidence for an additive allele effect of 5.13 kg of fat per lactation on the EBV for milk FY in bulls and 6.80 kg of fat of the own production performance in cows explaining 1.72 and 0.60% of the phenotypic variance in the analyzed populations, respectively. The analyses of bulls and cows consistently showed three haplotype groups that differed significantly from each other, suggesting at least two different mutations in the BDNF region affecting the milk FY. The FY increasing alleles also had low but significant positive effects on protein and total milk yield which suggests a general role of the BDNF region in energy partitioning, rather than a specific regulation of fat synthesis. The results obtained in dairy cattle suggest similar effects of BDNF on milk composition in other species, including man.


INTRODUCTION
The gene encoding the brain-derived neurotrophic factor (BDNF ) is among those genes that have been repeatedly associated with human obesity. Since a significant association between the V66M BDNF variant and body mass index was reported in healthy adults (Gunstad et al., 2006), an additional six SNPs located downstream of the BDNF gene were identified to be associated with body mass index in a genome-wide association study (Thorleifsson et al., 2009). Recently, the BDNF effect on body mass index was confirmed in the GIANT population with 249,796 individuals of European ancestry (Speliotes et al., 2010). The allele substitution effect of BDNF on the body mass index was small, 0.19 kg/m 2 , explaining 0.07% of the phenotypic variance in the population (Speliotes et al., 2010).
In both humans and knockout mice it was shown that BDNF can influence eating behavior and energy expenditure by the regulation of appetite and energy balance in the hypothalamus (Kernie et al., 2000;Unger et al., 2007). Higher BDNF levels are more likely to lower feed intake and increase physical activity (Monteleone et al., 2005). Stimulation of neuronal cells has shown that BDNF may change synaptic plasticity (Soliman et al., 2010). BDNF is well investigated for its role in neuronal differentiation in humans (Binder and Scharfman, 2004), however, its role in energy metabolism is not yet clear.
If the genetic determinants not only lead to variation in deposition of fat in adipose tissues but also to higher secretion of fat during lactation, we would also expect effects on fat yields (FYs) in milk. Therefore, we used estimated breeding values (EBVs) of bulls and own production performance data of cows of German Holstein dairy cattle to test the BDNF effect on milk FY.
The effect of BDNF observed in humans makes it a promising candidate gene to test in cows. However, restricting analyses to the genomic region of BDNF as a candidate gene potentially affecting milk FY, reduces the complexity of multiple testing, and even offers the opportunity to detect small genetic effects.
Two independent association tests were performed. In a bull population, we tested the association between SNPs in the genomic region of BDNF with the EBVs for milk FY. The EBV of a bull represents the mean value of own production performances of his daughters. In a cow population, the most significant results obtained in bulls were verified using own production performance data. They were adjusted for environmental effects and used as deviations of the mean of the population (yield deviation). The results provide evidence that the BDNF region contributes to the variation in milk FY, while also having an effect on protein yield (PY) and milk yield (MY). This leads to the conclusion that BDNF not only favors fat synthesis but generally has a role in the regulation of energy balance in cattle. Furthermore, genetic effects of different haplotype groups that were identified in the BDNF region provide evidence for at least two mutations segregating in the BDNF region in German Holstein cattle. Thus, the combined analyses of EBVs and own production performance data can increase the power to reliably detect effects, even with a low effect size.

ANIMALS
The association study was carried out with 2,402 breeding bulls and 1,476 cows of the German Holstein population. In bulls, family patterns were identified with 40 full sibs and 563 half sibs. Cows originated from 296 bulls. The average number of offspring per sire was 22.8; with a minimum of 1 (126 cases) to a maximum of 79 offspring (1 case). Fifty-six of those bulls were breeding bulls of the analyzed bull population. The cows were managed in three herds in the Northeast of Germany (Strucken et al., 2010).
For cows, own production performance data for MY, P%, PY, F%, and FY was available for lactations 1-3. For the association study yield deviations were used. Yield deviations for all milk composition traits were calculated using a restricted maximum likelihood (REML) approach with the following model: Where Y = milk production record; m = overall mean; h = fixed effect of herd (three classes); c = fixed effect of the calving season (46 classes); h × c = interaction between herd and calving season; f = linear regression on age at first calving; e = random residual.

GENOTYPIC DATA
The bulls were genotyped with the Illumina Bovine SNP50K Bead-Chip (Matukumalli et al., 2009) containing 54,001 SNPs. SNP data from this chip were subjected to rigorous validation by a remapping procedure suggested by (Schmitt et al., 2010) against the current Btau 4.0 assembly. In total 2,017 ambiguous SNP positions were defined as missing due to substantial deviations between the mapping strategy of the manufacturer and our own. A quality check of obtained genotypes, revealed that 8,748 SNPs were removed either because they failed genotyping in more than 10% of the animals (749 SNPs) or due to a minor allele frequency below 1% (7,998 SNPs). Finally, 45,334 SNPs passed the quality check for the subsequent association study. Out of 2,402 bulls, 48 were removed for low genotyping (>10% missing SNPs).
Cows were genotyped for five selected SNPs, where four SNPs compose the haplotype block with the lowest p-value in the BDNF region in the bull population. The fifth SNP is located in the DGAT1 region. Genotyping was performed using KBioscience assays. Primers for the genotyping of the selected SNPs are listed in Table 1.

HAPLOTYPE INFERENCE AND BLOCK COMPUTATION
Haplotype analyses for the bulls were carried out on a population of 2,354 bulls and 672 German Holstein bull dams, genotyped on the same Illumina platform. Based on more stringent filter criteria (<3% missing genotypes, minor allele frequency >5%, <5% missing SNP calls), haplotypes were derived using the software fastPHASE (Scheet and Stephens, 2006). The program was run with 10 random starts (parameter -T) and 25 iterations (parameter -C) per chromosome. Phased genotyping data was partitioned into haplotype blocks using the solid spine algorithm implemented in the software HAPLOVIEW v4.1 (Barrett et al., 2005). A block was defined if all markers within a region were in linkage disequilibrium (D > 0.8) with the first and last marker of that region but not necessarily with each other. According to the solid spine approach, nine haplotype blocks were found in the BDNF region.
The programme SimWalk (Weeks et al., 1995) was used to infer the phase of the haplotypes in the cow population. Haplotypes were inferred for the HTB4, which gave significant results in the bull population. For the generation of haplotype phases, we used the same four SNPs that contributed to HTB4 in the bull population.

ASSOCIATION ANALYSIS
In a first step, association analyses were performed for the EBVs for milk FY in bulls and milk FY deviations in cows. In a second step, pleiotropic effects of the indentified genomic region (haplotype and SNPs) on MY and other milk composition traits were tested in bulls and cows. BDNF is located on BTA15 at 58,220,662-58,221,703 bp. For the analysis of the bull population, the SNP (ARS-BFGL-NGS-112654) from the Illumina Bovine SNP50K BeadChip that was most adjacent to the genomic position of BDNF was defined as the center of the candidate region. Thirty-six SNPs that were located within 1 Mb up-and downstream of this center-SNP (57,251,933 and 59,189,700 bp) were used for association analyses. SNP association was tested using the software PLINK 1.06 (Purcell et al., 2007). A linear regression model estimated the additive effects of the SNP alleles. Since DGAT1 is already known to have a large effect on milk FY in dairy cattle (Grisart et al., 2002), the most significant DGAT1 SNP was fitted in the model to account for the allelic dosage of the DGAT1-effect. To find the most significant DGAT1-SNP on the chip, association analysis was performed with all SNPs within the 2 Mb region around the center-SNP of DGAT1.
To adjust for population stratification, the pairwise population concordance test (PPC) in PLINK was performed, based on an identity-by-state-similarity matrix. A total of 124 significant clusters (p < 10 −4 ) were identified. Those were fitted in the model using the multi dimensional scaling (MDS) approach. Association results were adjusted for multiple testing using the Bonferroni correction (p < 0.05). For significant SNPs, effects between genotype groups were tested for significance using ANOVA and subsequent Tukey-Kramer tests as implemented in the statistic programme SAS version 9.2 (SAS 2008).
Haplotype association analyses were performed with all blocks that were indentified in the 2 Mb candidate region. All haplotypes having a population frequency >1% were included in the analyses. The same adjustments as used in the single SNP association analyses were fitted in the model. For generating haplotype effect plots, the phased data obtained by fastPHASE was used. A Tukey-Kramer test was performed to test for differences between the haplotype groups. Diplotypes were generated and those occurring more than 10 times in the population were included to estimate the residual variance.
For the validation of the significant effects found in the bull population, an association test was performed in the cow population. For testing SNP and haplotype effects, a linear model (PROC MIXED) was applied using the SAS programme version 9.2 (SAS 2008). The model accounted for the effects of DGAT1 and the population substructure. The same DGAT1-SNP that was used in the bull population was genotyped in the cow population and fitted as a fixed effect to correct for DGAT1. Due to the high number of 296 different sires in the cow population, sires were added as a random effect in the model to account for the observed population structure. Haplotype association tests were performed as described for bulls. Genetic variances that can be explained by alleles, genotypes, haplotypes, or diplotypes were calculated. The residual variance explained by the genetic effect of single SNPs or haplotype blocks was estimated as a reduction of residual sum of squares in the full model (including all genetic effects under consideration) and the reduced model.

SNP ASSOCIATION WITH MILK FAT YIELD IN THE BULL POPULATION
Out of 36 SNPs in the 2 Mb BDNF candidate region, six SNPs were significantly associated with the average EBV for FY over the first three lactations (Figure 1). The most significant SNP for milk FY was Hapmap53286-rs29015961 (p = 3.362 × 10 −4 ). The difference between the two homozygous classes AA and GG of this SNP was 10.26 kg FY over the first three lactations. The genotype GG was associated with the lowest mean value for EBV of FY (12.57) and had the lowest frequency (5.81%) in the bull population (Figure 2A). All genotypes at this locus were significantly different from each other (p < 0.05). This genotype explained 1.72% of the total genetic variance in the bull population.

HAPLOTYPE ASSOCIATION WITH MILK FAT YIELD IN THE BULL POPULATION
The most significant SNP belonged to haplotype block 4, which was also most significantly associated (p < 10 −5 ) with the EBV of FY over the first three lactations and their average. This block, consisting of four SNPs, covered 187.5 Kb. The average D was 0.95 (Figure 1). The regression coefficient β for the haplotype GGAG with the lowest mean for the average EBV FY was −3.97. This means that each extra GGAG haplotype reduces the mean phenotype by 3.97 kg fat. The first SNP (BTB-00608390) of this block was also significantly associated (p < 0.047) with the average EBV for FY over the first three lactations and had similar genetic effects as the most significant SNP at the second position of this haplotype block. The two SNPs Hapmap43566-BTA-37197 and ARS-BFGL-NGS-36995 at the end of this haplotype block were not significantly associated with the EBV of FY when analyzed separately. Furthermore, their pattern of inheritance was different from the first two SNPs in the block (distance between the second and third SNP was 86.8 kb; Figure 2A). Only 25 bulls out of 2,402 (1.04%) had a haplotype (AGGA) that harbors the low fat alleles at the most significant SNP and the SNP at the fourth position. This implies that these two low effect alleles were derived independently from each other and are rarely inherited together.
In the most significant haplotype block 4, eight different haplotypes occurred in the bull population, of which five had a frequency above 1% in the population and were included in the association analyses. The effects of the five haplotypes on EBV FY over lactations 1-3 clustered in three distinct groups: group I consisted of haplotype GGAG, group II of haplotypes AAGA and AAAG and group III of haplotypes GAAG and AAAG ( Figure 2C). The differences in the mean EBV for milk FY were significant between all groups (p < 10 −3 ). The haplotype effects within a group did not differ. The haplotype AAAG occurring most frequently (33.34%) in the population belonged to group II, which showed a medium EBV for FY over the first three lactations of 20.98 kg (Figure 2C).

SNP ASSOCIATION WITH MILK FAT YIELD IN THE COW POPULATION
The significant association between the SNP Hapmap53286-rs29015961 and milk FY in the bull population could be verified in the cow population for the first lactation and the average over the www.frontiersin.org first two lactations. The direction and magnitude of the genotype effects in the cow population was consistent with those observed in the bull population ( Figure 2B). All genotype classes were significantly different from each other (p < 0.05). The difference between the homozygous GG and AA genotypes was 13.55 kg in the cow population. The FY of heterozygous cows was between the homozygous classes, and thus the mode of inheritance was additive. The genotype GG with the low milk FY occurred with a frequency of 4.68% in the population which is very similar to the bull population. This SNP explained 0.60% of the genetic variance in the cow population which is 1.12% less than in the bull population. Furthermore, the direction and magnitude of effects of the first and fourth SNP at the haplotype block 4 were consistent with those obtained in the bull population. The differences between the genotype classes for the third SNP in the haplotype block (Hapmap43566-BTA-37197 ) were rather small and inconsistent over all lactations (Figure 2B).

HAPLOTYPE ASSOCIATION WITH MILK FAT YIELD IN THE COW POPULATION
Using the same four SNPs that contributed to haplotype block 4 in the bull population, 16 different haplotypes were inferred in the cow population. Out of these 16 haplotypes, 9 occurred with a frequency below 1% in the population and therefore were excluded from further analyses. Among those were all five haplotypes present in the bull population. Those two haplotypes that did not occur in the bull population (GGGG and AAAA) had low frequencies (1.82 and 1.37%, respectively) in the cow population. After fitting the sire as a random and DGAT1 as a fixed effect, no significant associations of the haplotype with milk FY was found in any lactation. However, without adding those effects, significant associations were found for all lactations and their average, except for lactation three. Consistent with the bull population, the haplotype with the lowest milk FY was GGAG. In the cow population, a second haplotype (GGGG) occurred with low mean (Figure 2D). The haplotypes with the highest means were GAAG (13.20 kg) and AAGG (8.10 kg). As in the bull population, the cow haplotypes could be divided into three phenotypic effect groups (p < 0.01; Figure 2D).

ASSOCIATION WITH OTHER MILK TRAITS
The SNPs belonging to haplotype block 4, which was significantly associated with milk fat, were separately tested for association with MY and other milk composition traits. In bulls, Hapmap53286-rs29015961, the second SNP of haplotype block 4, was also significantly associated with PY in lactation two (p = 0.005) and the Frontiers in Genetics | Genetic Architecture average of all lactations (p = 0.033) and with MY for all lactations and the average of them, except for the second lactation ( Table 2). The minor allele of the SNP had a negative effect on yield traits and fat content, but positive effects on protein content in lactation one and for the average of all three lactations. In the second and third lactations, the effect of the minor allele of this SNP was nearly 0.
In cows, consistent with bulls, Hapmap53286-rs29015961 had also significant effects on PY (p = 0.042) and MY which were significant in lactation two (p = 0.020). The direction and magnitude of effects were similar with those observed in the bull population.
In both, the bull and cow population, haplotype effects were tested on additional milk composition traits. In bulls, haplotype block 4 had similar effects on the EBV for PY for the second lactation and the average of the first three lactations and on MY for all lactations and the average of them. For all yield traits, the haplotype GGAG had a negative and haplotype AAGG a positive effect. No significant association was found with milk content traits. In cows, significant associations between the haplotypes and milk composition traits disappeared after fitting the whole model with DGAT1 as a fixed effect and the sire as a random effect. However, when the diplotypes were analyzed, significant association with PY, P%, and MY for lactation 2 were found. Furthermore, if only the sire was fitted as a random effect, significant associations between diplotypes and F% for lactations 2 and 3 and the average of the first two lactations were evident.

VARIANCE CONTRIBUTION OF SIGNIFICANT GENETIC EFFECTS
In the bull population, the subsequent fitting of alleles of the four SNPs of the haplotype block 4 in the model decreased the residual variance and thus increased the genetic variance explained by the four SNPs (Table 3). Fitting haplotypes into the model resulted in a smaller effect than fitting the alleles of the four SNPs. Genotypes explained more of the variance than the haplotypes or alleles. The diplotypes explained the highest proportion (2.22%) of the genetic variance (Table 3). This pattern was found throughout all lactations.
The proportion of genetic variance explained by genotypes, haplotypes, and diplotypes was also calculated within the cow population. Except for the yield deviation of milk FY in lactation 2, the diplotypes explained most of the genetic variance for all milk composition traits. The maximal genetic variance explained by www.frontiersin.org  diplotypes was 1.50% for milk FY in the average of the first two lactations. Just as diplotypes most accurately captured the genetic structure in bulls and cows at the most significant locus, they also seem to capture the genetic variance most effectively.

DISCUSSION
The results provide evidence that the BDNF gene region harbors SNPs and haplotypes which are significantly associated with FY in German Holstein cattle. SNPs harbored in the most significant haplotype block were significantly associated with yield traits, but not content traits. For all of these four SNPs the minor alleles were associated with the lower milk FY. This is probably a result of selection for high milk FY over many generations. In addition to variation in milk FY the analyzed region also showed effects on PY and MY. As MY is highly dependent on the amount of osmotically active lactose in the milk, it can be assumed that BDNF also regulates the level of lactose. However, EBVs were not available for lactose in milk. This leads to the suggestion that BDNF not only favors fat synthesis but generally might have a role in the regulation of energy balance. The variance that can be explained by the most significant SNP in bulls was 1.7%, which was low but higher than in the human obesity study (0.9% phenotypic variance for BMI; Speliotes et al., 2010). This shows that the size of the bull population and the accuracy of breeding values are large enough to detect loci with very small genetic effects. The genetic variance explained by the most significant SNP in the cow population was 0.6%. The yield deviations in the cow data were also adjusted for environmental effects, but the number of observations per individual is much lower than in the bull population, where records of many daughters for one individual are available.
Adding DGAT1, which accounts for 30% of the variation in milk F% (Grisart et al., 2002) as a fixed effect in both populations, increased the p-values. In the cow population, no significant associations were found after fitting DGAT1 into the model. However, in the bull population, significant associations with other milk traits than milk FY disappeared when DGAT1 was not fitted in the model. By the inclusion of DGAT1 the small proportion of variance explained by BDNF might appear or disappear depending on the trait. Therefore, it is debatable whether it is necessary to add DGAT1 as a fixed effect in the model with the risk of losing some putative significant SNPs, or to add it as an effect to account for the major effect of DGAT1. An epistatic interaction between the DGAT1 and BDNF loci was not found in our populations (data not shown).
Our data suggest the presence of more than one genetic variation in the BDNF region responsible for the observed variation in milk FY. In both the bull and the cow population, we could distinguish three haplotype groups having different effects on milk FY. For the existence of three haplotype groups at least two mutations with different effects must occur. In the bull population eight haplotypes were identified of which five occurred with frequencies above 1%. In the cow population, 16 haplotypes were identified, of which seven occurred with a frequency above 1%. This indicates that a large genetic variability exists in the cow population.
The candidate gene approach that we used for the analysis of effects of the BDNF region reduced the stringency of multiple testing by including prior knowledge and thus supported the identification of significant associations between few SNPs and milk traits. According to the infinitesimal model (Visscher, 2008), many genes with small effects are expected to contribute to milk FY and the other milk composition traits. Here, we demonstrate the identification of such a locus with a rather small effect in a bull population and the successful validation of the effects in a cow population.
The analyses of SNP and haplotype effects provide evidence that haplotypes captured more of the underlying genetic variation than the sum of all single SNPs in the haplotype block. The diplotype level is the genetic structure in which the biology of the organism functions whereas the haplotypes is the genetic level which is inherited. Partitioning into diplotypes showed that more genetic variance can be explained by reflecting the genetic architecture as precisely as possible.

CONCLUSION
The study has shown that bull and cow populations having a large number of records for different milk traits are powerful to test candidate genes with low effects. The additional information gained by including haplotypes and diplotypes into the analysis may contribute to a higher accuracy in determining positional genetic effects. This in turn may enhance breeding programs aiming on genomic selection, which currently use additive effects of single SNPs only. Finally, the detailed information on the genetic structure and effects of the BDNF region can be used for fine mapping studies and resequencing of related haplotypes to facilitate the search for the causal variants.