Skip to main content
  • Research article
  • Open access
  • Published:

Genome-wide association study on Fourier transform infrared milk spectra for two Danish dairy cattle breeds

Abstract

Background

Infrared spectral analysis of milk is cheap, fast, and accurate. Infrared light interacts with chemical bonds present inside the milk, which means that Fourier transform infrared milk spectra are a reflection of the chemical composition of milk. Heritability of Fourier transform infrared milk spectra has been analysed previously. Further genetic analysis of Fourier transform infrared milk spectra could give us a better insight in the genes underlying milk composition. Breed influences milk composition, yet not much is known about the effect of breed on Fourier transform infrared milk spectra. Improved understanding of the effect of breed on Fourier transform infrared milk spectra could enhance efficient application of Fourier transform infrared milk spectra. The aim of this study is to perform a genome wide association study on a selection of wavenumbers for Danish Holstein and Danish Jersey. This will improve our understanding of the genetics underlying milk composition in these two dairy cattle breeds.

Results

For each breed separately, fifteen wavenumbers were analysed. Overall, more quantitative trait loci were observed for Danish Jersey compared to Danish Holstein. For both breeds, the majority of the wavenumbers was most strongly associated to a genomic region on BTA 14 harbouring DGAT1. Furthermore, for both breeds most quantitative trait loci were observed for wavenumbers that interact with the chemical bond C-O. For Danish Jersey, wavenumbers that interact with C-H were associated to genes that are involved in fatty acid synthesis, such as AGPAT3, AGPAT6, PPARGC1A, SREBF1, and FADS1. For wavenumbers which interact with –OH, associations were observed to genomic regions that have been linked to alpha-lactalbumin.

Conclusions

The current study identified many quantitative trait loci that underlie Fourier transform infrared milk spectra, and thus milk composition. Differences were observed between groups of wavenumbers that interact with different chemical bonds. Both overlapping and different QTL were observed for Danish Holstein and Danish Jersey.

Background

There is a large number of applications for trait predictions utilizing Fourier transform infrared (FT-IR) milk spectra from the mid-infrared range. Fourier-transform infrared spectroscopy determines light absorbance across the infrared spectrum. Light is absorbed when it interacts with chemical bonds. Wavenumber 1690 cm− 1, for example, interacts with C=O of amide I, and 1600 cm− 1 is involved in N-H bending of amide II [1, 2]. These chemical bonds are typical for protein molecules. Wavenumbers from the lower energy region that ranges from 1150 to 1040 cm− 1 interact with C-OH, which is abundantly present in sugar molecules [1, 2]. This chemical bond, however, is also present in fat and protein molecules, but more scarcely. Wavelengths from the infrared region that ranges from 2950 to 2850 cm− 1 induce C-H stretching [1, 2]. Triglyceride molecules are rich in C-H bonds, but C-H bonds are also present in many other molecules.

Mid-infrared light is commonly used in combination with the principal least square regression method to analyse chemical composition of milk [3]. The major milk components fat, protein, and lactose have been successfully predicted with FT-IR milk spectra [3]. In addition, minor milk components have been predicted with FT-IR milk spectra, such as fatty acids [4,5,6], protein fractions [7, 8], and ketone bodies [9,10,11]. Concentration of ketone bodies in milk can be used as an indicator for subclinical ketosis [9,10,11], or energy balance [12, 13].

Associations to genomic regions have been observed for both milk composition, and infrared milk spectra. Fatty acid composition, for example, has been associated to many different genomic regions [14, 15]. FT-IR milk spectra have been linked to genes that have been associated to milk composition previously, such as Diacylglycerol O-acyltransferase 1 (DGAT1) or beta-lactoglobulin (PAEP) [16,17,18]. A genome wide association study (GWAS) on a subset of wavenumbers revealed associations for individual wavenumbers to a variety of genes [18]. Examples are the gene for the growth hormone receptor, or the gene UMPS [18]. FT-IR milk spectra are also moderately to strongly heritable [17, 19,20,21]. To get a better understanding of the genetic background of FT-IR milk spectra, it is necessary to further study the association between milk spectra and the genome.

Cattle breed influences milk composition [19, 22,23,24], and the genetic architecture of milk composition [25,26,27]. These breed differences in milk composition are reflected in the FT-IR milk spectra. Heritability of FT-IR milk spectra varied across breeds [17, 19,20,21]. Not much is known about the breed differences in the genes that indirectly underlie FT-IR milk spectra. Enhanced knowledge on breed differences in the genetic architecture of FT-IR milk spectra could provide us with a better understanding of differences in milk composition across breeds. Finally, it could facilitate future application of FT-IR milk spectra in across breed prediction of novel phenotypes.

The aim of this study is to perform a GWAS on a selection of wavenumbers in two dairy cattle breeds, Danish Holstein and Danish Jersey.

Results

Selection of wavenumbers

After removal of wavenumbers which interact with water molecules, 530 waveunumbers were left. For these 530 wavenumbers, correlations were calculated. The correlation matrices were nearly identical for Danish Holstein and Danish Jersey. The heatmap of the correlation matrix for Danish Holstein is presented in Fig. 1. For both breeds 17 blocks of highly correlated neighbouring wavenumbers were observed, where the correlation between wavenumbers was > 0.95. From each block, the wavenumber with the highest correlation sum was selected. For four blocks, the selected wavenumber was different for Danish Holstein and Danish Jersey (Table 1). For both Danish Holstein and Danish Jersey, 15 out of the 17 selected wavenumbers had a heritability > 0.05. For Danish Holstein and Danish Jersey separately, Table 1 presents an overview of the selected wavenumber per block, the chemical bond with which the wavenumber interacts, heritability of the selected wavenumber, number of quantitative trait loci (QTL) identified for the selected wavenumber, number of QTL unique for the selected wavenumber, and chromosomes on which QTL were located. A QTL was defined as one, or several overlapping groups of 100 neighbouring SNPs (SNP group), for which each individual SNP group explained > 0.35% of the total additive genetic variation. A peak was defined as the SNP group within a QTL, which explained most of the total additive genetic variation.

Fig. 1
figure 1

Heatmap of the phenotypic correlation matrix for wavenumbers in Danish Holstein. The upper and lower triangle are identical. Seventeen blocks of strongly positively correlated neighbouring wavenumbers are indicated with black dashed square boxes. The upper left corner represents wavenumber 3008 cm− 1 and wavenumber group 1, and the lower right represents wavenumber 925 cm− 1 and wavenumber group 17

Table 1 Fifteen selected wavenumbers (wvn) from fifteen positively correlated wavenumber blocks (see Fig. 1), the chemical bond (CB) with which the selected wavenumber interacts, heritability of the selected wavenumber, total number of QTL for the selected wavenumber, number of QTL unique for the selected wavenumber, and chromosomes on which QTL were located

Peak regions

Table 2 shows an overview of genomic regions, which were associated to groups of wavenumbers, which interact with different chemical bonds. For each group of wavenumbers, genomic regions of 100 consecutive SNPs which explained > 0.35% of the total additive genetic variation are listed. This genomic region is referred to as the “peak region”. There can be more peak regions on one chromosome. Table 2 gives an overview of the highest peak region for each chromosome, meaning that only one peak region per chromosome is described. A peak region is not necessarily associated to all wavenumbers of a group. For each peak region, those wavenumbers are presented for which the proportion of explained additive genetic variation by the peak region > 0,35%. Candidate genes located within the peak region are named in the final column.

Table 2 Top SNP groups explaining most total additive genetic variation for wavenumbers, which interact with different chemical bonds

An overview of the number of QTL per chromosome, for Danish Holstein and Danish Jersey separately, and the number of overlapping QTL between the two breeds are shown in Table 3. Results are presented for all wavenumbers combined, and for groups of wavenumbers based on the chemical bond with which they interact (Table 1).

Table 3 Number of QTL per chromosome observed for Danish Holstein (DH) and Danish Jersey (DJ), for wavenumbers which interact with different chemical bonds. Overlap (ol) indicates number of overlapping QTL between the two breeds

QTL and wavenumbers interacting with different chemical bonds

Wavenumbers interacting with alkanes

For Danish Holstein, the three peak regions explaining most additive genetic variation for wavenumbers interacting with alkanes were positioned on BTA 6 (0.54%) harbouring the casein (CSN) cluster, on BTA 14 (2.04%), and on BTA 29 (0.48%). For Danish Jersey, the three peak regions were positioned on BTA 6 (2.25%) harbouring the CSN cluster, on BTA 14 (2.10%), and on BTA 20 (0.67%) harbouring GHR, and MRPS30. The CSN cluster is a genomic region on BTA 6 containing genes, which code for the milk protein casein.

Wavenumbers interacting with C=O

For Danish Holstein, the three peak regions explaining most additive genetic variation for wavenumbers interacting with C=O were positioned on BTA 5 (0.43%) harbouring MGST1, on BTA 14 (9.76%), and on BTA 17 (0.39%). For Danish Jersey, the four peak regions were positioned on BTA 12 (0.54%), on BTA 14 (5.11%), on BTA 19 (0.54%) harbouring SREBF1, and on BTA 20 (0.60%) harbouring GHR, and MRPS30.

Wavenumbers interacting with C-H

For Danish Holstein, the two peak regions explaining most additive genetic variation for wavenumbers interacting with C-H were positioned on BTA 5 (0.39%) harbouring MGST1, and on BTA 14 (7.88%). For Danish Jersey, the three peak regions were positioned on BTA 12 (0.63%), on BTA 14 (5.38%), and on BTA 20 (0.64%) harbouring GHR, and MRPS30.

Wavenumbers interacting with C-O

For Danish Holstein, the three peak regions explaining most additive genetic variation for wavenumbers interacting with C-O were positioned on BTA 14 (8.55%), on BTA 19 (0.52%), and on BTA 29 (0.54%). For Danish Jersey, the four peak regions were positioned on BTA 1 (0.71%), and AGPAT3, on BTA 6 (0.88) harbouring the CSN cluster, on BTA 14 (5.27%), and on BTA 19 (0.68%) harbouring SREBF1.

Wavenumbers interacting with CO-N

For Danish Holstein, the three peak regions explaining most additive genetic variation for wavenumbers interacting with C-ON were positioned on BTA 5 (0.36%) harbouring MGST1, on BTA 6 (0.35%) harbouring the CSN cluster, and on BTA 14 (5.19%). For Danish Jersey, the three peak regions were positioned on BTA 6 (1.04%) harbouring the CSN cluster, on BTA 14 (1.90%), and on BTA 29 (0.49%) harbouring FADS1.

Wavenumbers interacting with N-H

For Danish Holstein, the three peak regions explaining most additive genetic variation for wavenumbers interacting with N-H were positioned on BTA 6 (0.55%) harbouring the CSN cluster, on BTA 14 (9.73%), and on BTA 20 (0.53%) harbouring ANKH. For Danish Jersey, the three peak regions were positioned on BTA 3 (0.73%) harbouring GBA, on BTA 6 (2.32%) harbouring the CSN cluster, and on BTA 14 (4.89%).

Wavenumbers interacting with –OH

For Danish Holstein, the three peak regions explaining most additive genetic variation for wavenumbers interacting with -OH were positioned on BTA 14 (1.12%), on BTA 20 (0.40%) harbouring ANKH, and on BTA 29 (0.45%). For Danish Jersey, the three peak regions were positioned on BTA 6 (1.09%) harbouring the CSN cluster, on BTA 14 (0.52%), and BTA 16 (0.46%).

Breed differences

Breed differences are clearly visible in Tables 1, 2 and 3, and the Manhattan plots in Additional files 1 and 2. Overall, more QTL were observed for Danish Jersey compared to Danish Holstein (Tables 2 and 3). For Danish Holstein, most QTL were located on BTA 19 and BTA 20, and for Danish Jersey on BTA 5 and BTA 16. For both breeds, most QTL were observed for wavenumbers interacting with C-O. Heritability of wavenumbers was slightly lower for Danish Jersey than for Danish Holstein (Table 1). The proportion of explained variation by the peak region of DGAT1 was higher in Danish Holstein compared to Danish Jersey.

Overlapping QTL

Overlapping peak regions were observed on BTA 5 (91.1–94.9 Mbp) harbouring MGST1, on BTA 6 (81.1–84.6 Mbp) harbouring the CSN cluster, on BTA 19 (32.1–37.6), on BTA 20 (56.0–60.8 Mbp) harbouring ANKH, on BTA 21 (6.2–11.0 Mbp) harbouring IGFIR, and on BTA 25 (0.1–4.7 Mbp). Most overlapping QTL were observed for wavenumbers interacting with C-O and N-H. No overlapping QTL between the two breeds were observed for wavenumbers interacting with C=O, C-H, or –OH.

Discussion

To get a better understanding of the genetics of milk composition, this study aimed at performing a GWAS on a selection of wavenumbers interacting with different chemical bonds in two dairy cattle breeds, Danish Holstein and Danish Jersey.

For each breed separately, fifteen wavenumbers were selected from blocks of strongly positively correlated neighbouring wavenumbers based on the maximum correlation sum within block, and a minimum heritability of 0.05. The correlation between wavenumbers within one block were close to one, and analysis of all wavenumbers within one block would most probably result in similar findings. For four blocks, different wavenumbers were selected for Danish Holstein and Danish Jersey (Table 1). The selected wavenumbers were within the same infrared region. Therefore, we assumed that results of e.g. 1295 cm− 1 for Danish Holstein are comparable to results of 1299 cm− 1 for Danish Jersey.

BTA 14

A QTL on Bos Taurus autosome (BTA) 14 in the genomic region of DGAT1 was associated to all wavenumbers, with the exception of 1557 cm− 1 in Danish Jersey. The QTL in DGAT1 explained most additive genetic variation for 14 out of 15 wavenumbers in Danish Holstein, and for 9 out of 15 wavenumbers in Danish Jersey (Table 2). Because DGAT1 is a well-known major milk gene, the genomic region of BTA 14 will not be thoroughly discussed.

Wavenumbers interacting with alkanes

The wavenumber 1449 cm− 1 is known to interact with alkanes [1, 2]. The chemical bonds present in alkanes resemble those of saturated fatty acids. For both breeds, a QTL on BTA 19 (19_a) was identified. This genomic region harbours the gene SREBF1. The gene SREBF1 is known as a key player in fatty acid synthesis [14]. In line with this, Bouwman et al. [28] observed a QTL for saturated fatty acids in milk in the same genomic region. For both Danish Holstein and Danish Jersey, a QTL on BTA 6 (6_b) harbouring the CSN cluster was observed. This QTL has previously been associated to protein percentage [27, 29,30,31], caseins, whey proteins [30], and cheese yield [32]. In a GWAS on wavenumbers in Dutch Friesian Holstein, Wang and Bovenhuis [18] also observed an association between the QTL 6_b and wavenumber 1469 cm− 1.

Wavenumbers interacting with C=O

The chemical bond C=O typically appears in fat molecules and protein molecules, and interacts with 1735 and 1696 cm− 1 [1, 2]. The QTL associated to wavenumbers interacting with C=O have been associated to a variety of milk production traits. In Danish Holstein, a QTL on BTA 5 (5_a) harbouring MGST1 has previously been associated to milk composition [29, 30, 33], and fatty acid composition [28, 34]. In Danish Jersey, the QTL on BTA 6 (6_a) harbouring PPARGC1A, and BTA 19 (19_a) harbouring SREBF1 have both been associated to fat percentage, and fatty acid composition in milk [14, 28, 34, 35]. Furthermore, for Danish Jersey, several QTL were identified that were linked to protein in milk previously. A QTL on BTA 11 (11_b) harbouring PAEP has been strongly associated to beta-lactoglobulin in milk, and protein composition [27, 30].

Wavenumbers interacting with C-H

The chemical bond C-H is present in many molecules, such as fat, protein, and lactose. The C-H bond strongly interacts with 2988 and 2872 cm− 1 [1, 2]. The C-H bond is most abundantly present in the fatty acid tails of fat molecules. This is why wavenumbers in the region of 2988 and 2872 cm− 1 are used for prediction of fat percentage in milk [1, 36]. In Danish Holstein, the QTL on BTA 5 (5_b) harbouring MGST1, and the QTL on BTA 17 (17_a) have been associated to fatty acid composition in milk previously [28, 34]. For Danish Jersey, many QTL were located in genomic regions of genes, which have been associated to milk fatty acid synthesis [14, 37]. Examples of these genes are AGPAT3 on BTA 1 (1_b), PPARGC1A on BTA 6 (6_a), SREBF1 on BTA 19 (19_a), AGPAT6 on BTA 27 (27_b), and FADS1 on BTA 29 (29_b) [14, 37]. The gene AGPAT6 on BTA 27 is described as one of the key links in milk fatty acid synthesis [37]. Interestingly, the genomic region of AGPAT6 was only associated to wavenumbers that interact with C-H (Table 2). An additional QTL on BTA 20 (20_b) harbouring ANKH was observed for Danish Holstein. This QTL has been strongly associated to alpha-lactalbumin [27], and lactose percentage in milk [18]. For Danish Jersey, two QTL (11_b and 21_b) were found. Within this genomic region, two genes were located that have been linked to proteins in milk [27, 30, 38].

Wavenumbers interacting with C-O

The chemical bond C-O is abundantly present in sugar molecules, and it interacts with wavenumbers in the infrared region from 1250 to 950 cm− 1 [1, 2]. This infrared region and the infrared region that ranges from 1400 to 1250 cm− 1 (see next section) are used for prediction of lactose in milk [1, 36]. For Danish Holstein, the observed QTL did not reveal a strong link between this infrared region and lactose in milk. The QTL on BTA 5 (5_b) harbouring MGST1, however, has been associated to milk composition, including lactose percentage [29]. Most of the QTL observed for Danish Holstein, however, have been associated to fatty acids or groups of fatty acids, such as the QTL on BTA 17 (17_a), BTA 19 (19_c), BTA 21 (21_a and 21_c), and BTA 28 (28_a) [34, 39]. For Danish Jersey, on the other hand, many of the currently observed QTL have been linked to lactose in milk. Four QTL (19_a, 19_b, 22_b, and 28_b) have been associated to lactose percentage in milk [18, 38]. In addition, the QTL on BTA 1 (1_b) harbouring SLC37A1 and AGPAT3, and the QTL on BTA 5 (5_a) were both associated to alpha-lactalbumin in milk [27, 30]. Alpha-lactalbumin is a milk protein that plays a critical role in converting glucose into lactose [40]. Finally, the QTL 22_b was associated to wavenumbers, which interact with C-O exclusively. The QTL 22_b is harbours the gene lactotransferrin (LTF). The protein lactotransferrin is a selective antibacterial milk protein that is involved in the mucosal protection of the mammary gland, and possibly protects against mastitis [41, 42].

Wavenumbers interacting with CO-N and N-H

The chemical bonds CO-N, and N-H are present in protein molecules. These chemical bonds interact with the infrared region that ranges from 1550 to 1500 cm− 1, and infrared region around 1600 cm− 1, respectively [1, 2]. These infrared regions are used for prediction of protein percentage in milk [1, 36]. The two groups of wavenumbers interacting with CO-N and N-H have many overlapping QTL, and therefore will be discussed together. Firstly, for both breeds, a strong association was observed between CO-N and N-H interacting wavenumbers and the CSN cluster on BTA 6 (6_b; Table 2). The CSN cluster has been associated to many traits related to protein in milk, such as protein percentage [18, 27, 30, 33], and protein composition [27, 30]. Secondly, a QTL on BTA 20 (20_b) was observed for both breeds. This QTL harbours the gene ANKH, which is strongly associated to alpha-lactalbumin, and it is expressed in mammary tissue [27]. Finally, one more QTL was observed for both breeds, which was located on BTA 17 (17_a). This QTL has been associated to alpha-S2-casein in milk [30].

For Danish Jersey, additional QTL were identified that have been associated to milk protein previously. Firstly, three QTL on BTA 3 (3_a), BTA 10 (10_b), and BTA 20 (20_a) have been associated to protein percentage in milk [30, 33]. Secondly, the QTL on BTA 11 (11_b) is located within a genomic region, which harbours several genes that control beta-lactoglobulin in milk [30, 43, 44]. Thirdly, the QTL no BTA 24 has been associated to beta-lactoglobulin previously as well [30]. Finally, like both the QTL on BTA 20 (20_b), and the QTL on BTA 5 (5_a) have been linked to alpha-lactalbumin [30].

Wavenumbers interacting with –OH

Like C-O, the chemical bond -OH is abundantly present in sugar molecules. The chemical bond –OH interacts with wavenumbers in the infrared region that ranges from 1500 to 1250 cm− 1 [1, 2]. For Danish Holstein, wavenumbers from this infrared region were associated to the QTL on BTA 20 (20_b). This QTL harbours ANKH, which has been strongly associated to alpha-lactalbumin in milk [27]. Alpha-lactalbumin, as discussed earlier, has been described as a key player in lactose synthesis [40]. In Danish Jersey, the QTL explaining most variation was located on BTA 6 (6_b), which harbours the CSN cluster. Another QTL was located on BTA 1 (1_b) harbouring SLC37A1 and AGPAT3, and has been associated to alpha-lactalbumin [27]. Two other QTL, which were positioned on BTA 19 (19_b) harbouring FASN and CCDC57 and on BTA 28 (28_b), have been linked to lactose percentage in milk [18]. These two QTL have also been associated to wavenumbers surrounding wavenumber 1299 cm− 1 [18].

Breed differences

Breed has an effect on milk composition [24, 45], FT-IR milk spectra [5, 23], and the heritability of FT-IR milk spectra [17, 20, 21, 46]. In the current study more QTL were observed for Danish Jersey than for Danish Holstein (Table 1). A reason for this observation could be that DGAT1 explained more additive genetic variation in Danish Holstein than in Danish Jersey. The less dominant role of DGAT1 for Danish Jersey could have allowed for smaller effects to be visible. This could have resulted in the seemingly more polygenic character of milk spectra in Danish Jersey. Differences in allele frequency for the DGAT1 gene have been described before [20, 25]. The fact that not the same QTL were observed for both breeds could have been caused by differences in allele frequencies for SNPs between the two breeds, or even the complete absence of SNPs in one breed [25, 47]. When applying milk spectra directly for estimating breeding values of milk components, these breed differences in allele frequencies may cause reduced prediction accuracy, when predicting across breeds.

Conclusion

The current study observed multiple QTL for FT-IR milk spectra. Different QTL were observed for wavenumbers interacting with different chemical bonds. Wavenumbers that interact with the same chemical bond were often associated to the same QTL, yet some QTL were observed for small subsets of wavenumbers. Different QTL were observed for Danish Holstein and Danish Jersey.

Methods

Study population

The study population consisted of 3274 Danish Holstein cows from 354 herds, and 3408 Danish Jersey cows from 175 herds. For Danish Holstein, 3001 cows were in their first parity, and 273 in their second. For Danish Jersey, 3125 cows were in their first parity, and 283 in their second. For Danish Holstein, 19,656 morning-milk records were provided. For Danish Jersey, 20,228 morning milk records were provided. Cows had between one and twenty milk records with on average 32 days between records. Milk records were collected from October 1st 2015 to September 30th 2016. The year was split into summer, from April 1st 2016 through September 30th 2016, and winter, from October 1st 2015 through March 31st 2016. Milk records were collected from 1 through 400 days in milk (DIM). Obvious outlying milk records with a fat% > 8.0, or a protein% < 2.5 or > 5.0 in Danish Holstein, and a protein% > 5.5 in Danish Jersey were removed from the dataset.

Morning milk records were collected and provided by RYK (Aarhus N, Denmark), the Danish milk recording organization. Infrared spectral analysis was performed by Eurofins-Steins laboratory (Vejen, Denmark) with the MilkoScan FT+ (Foss, Hillerød, Denmark). Transmittance values for 1060 wavenumbers in the infrared region of 5008–925 cm− 1 were provided.

Genotypes

The study population was genotyped with the EuroG10K custom SNP chip. The EuroG10k SNP chip is composed of two parts: (1) SNP from the BovineLD Genotyping BeadChip v.2 [48], and (2) a custom part of selected SNP from sequence data as part of 1000 Bull Genomes Project Run 4 [49] based on their functional annotation or based on GWAS results [50]. Genotypes were imputed from the EuroG10K custom SNP chip to the 50 K using BEAGLE 4 [51]. Reference populations for imputation consisted of 4000 cows for Danish Holstein, and 4576 cows for Danish Jersey. Reference cows were genotyped on the Illumina 50 K BovineSNP50 v.2 BeadChip (Illumina Inc., San Diego, CA). Only autosomal SNPs which were present in both the Danish Holstein reference population and the Danish Jersey reference population were selected. During quality control, SNPs with more than 40% missing genotypes or with a minor allele frequency (MAF) of < 0.01 were excluded. After quality control, genotypes of Danish Holstein cows were imputed from 10,353 to 43,807 SNPs, and genotypes of Danish Jersey cows from 9749 to 39,235 SNPs. Median distance between SNPs was 41 kb for Danish Holstein, and 43 kb for Danish Jersey. All SNPs used for analysis are present on the Illumina BovineSNP50 v.2 BeadChip (Illumina Inc., San Diego, CA).

Phenotypes

FT-IR Milk spectra

Transmittance values were provided for 1060 wavenumbers in the mid-infrared region of 5008–925 cm− 1. Wavenumbers in the infrared regions from 5008 to 3008 cm− 1, and from 1669 to 1623 cm− 1 interact with water molecules, and were excluded from the analysis. A total of 530 wavenumbers were left for further analysis.

Selection of wavenumbers

Selection of wavenumbers was done for each breed separately. Correlations between 530 wavenumbers corrected for season, parity, days in milk, and herd were calculated in R software [52]. The correlation matrix was used to make a heatmap, where axes were sorted in order of wavenumber from 3008 cm− 1 through 925 cm− 1 (Fig. 1). Blocks of strongly positively correlated neighbouring wavenumbers were defined by visual inspection of the heatmap. Within each block, correlation sums were calculated for each wavenumber individually, and the wavenumber with the highest correlation sum was selected for further analysis.

Genetic analysis

Model description

Analysis of selected wavenumbers was done with the Bayz software package (http://www.bayz.biz/) [53]. We used the model:

$$ {y}_{ijkl}=\mu +{Parity}_i+{Season}_j+{\beta}_1{DIM}_{ijkl}+{\beta}_2{e}^{-0.05{DIM}_{ijkl}}+{Herd}_k+{CowA}_l+{CowPE}_l+{E}_{ijkl.} $$
(1)

Where yijkl is the transmittance value for one selected wavenumber; μ is mean transmittance value; Parityi corrects for the fixed effect of parity (i = 1 or 2); Seasonj corrects for the fixed effect of season during which the milk sample was collected (j = summer or winter); β1DIMijkl and \( {\beta}_2{e}^{-0.05{DIM}_{ijkl}} \) correct for lactation stage (Wilmink function) [54], where DIMijkl is dimijkl /365 (dimijkl = 1…365). For all fixed effects and regressors, a uniform prior distribution was assumed, where ~ UNI(0,+∞); Herdk is a random herd effect, for which a normal prior distribution was assumed, where Herd ~ N(0, \( {\sigma}_{Herd}^2 \)); CowPEl is a permanent environmental effect of cow l, for which a normal prior distribution was assumed, where CowPE ~ N(0, \( {\sigma}_{PE}^2 \)); and Eijkl is the residual variance, for which a normal prior distribution was assumed, where E ~ N(0, \( {\sigma}_E^2 \)). CowAl is the additive genetic effect of cow l, and was modeled using a hierarchical model to depend on SNP effects:

$$ {CowA}_l=\sum m amglm $$
(2)

Where am is the additive effect of SNP m; glm is the allele dosage for SNP m of cow l. Allele dosages were centred. For the additive genetic value, a normal prior distribution was assumed, where CowA ~ N(0, \( {\sigma}_A^2 \)), and all SNP variance parameters had a uniform prior distribution ~ UNI(0,+∞).

A Metropolis-Hastings sampler was used, with 70,000 iterations, including a burn-in of 30,000 iterations.

For all selected wavenumbers, heritability was calculated as:

$$ {h}^2=\frac{\sigma_A^2}{\sigma_{Herd}^2+{\sigma}_A^2+{\sigma}_{PE}^2+{\sigma}_E^2} $$
(3)

Where \( {\sigma}_A^2 \) is the additive genetic variance; \( {\sigma}_{Herd}^2 \) is the variance explained by herd; \( {\sigma}_{PE}^2 \) is the permanent environmental variance; \( {\sigma}_E^2 \) is the residual variance. Wavenumbers with a heritability < 0.05 were excluded from further analyses.

Grouping SNPs

Within each chromosome, SNPs were divided into groups of 100 consecutive SNPs [55]. The grouping procedure was repeated five times for each chromosome, starting with counting at SNP 1, 21, 41, 61, or 81 on the chromosome. Between the five repeated procedures, SNP groups overlapped, yet SNP groups were never identical. Groups with < 80 SNPs were excluded from analysis.

For each group of 100 SNPs, variance of the genomic estimated breeding value was calculated with the gbayz function of Bayz software (http://www.bayz.biz/) [53]. Proportion of total additive genetic variance explained per SNP group was calculated as:

$$ \%{\sigma}_{A, ij}^2=\frac{\sigma_{gEBV, ij}^2}{\sigma_{A,i}^2}\ast 100\% $$
(4)

Where \( \%{\sigma}_{A, ij}^2 \) is the percentage of total additive genetic variance of selected wavenumber i explained by SNP group j; \( {\sigma}_{gEBV, ij}^2 \) is the variance of the genomic estimated breeding value for selected wavenumber i of SNP group j; and \( {\sigma}_{A,i}^2 \) is the total additive genetic variance for selected wavenumber i.

Visual inspection was done on Manhattan plots of \( \%{\sigma}_{A, ij}^2 \), where \( \%{\sigma}_{A, ij}^2 \) of a group was represented by the middle SNP as orientation point (Additional files 1 and 2). For each selected wavenumber, QTL were collected.

Availability of data and materials

The raw datasets that were used in the current study are not available to the public, but can be requested for on reasonable grounds from the responsible co-author (AJB, bart.buitenhuis@mbg.au.dk).

SNPs that were used in the current study were present on the 50 k Bovine SNP array from Illumina. SNP names, and the position of SNPs can be found on:

http://support.illumina.com/downloads.html.

Gene annotation was performed for all SNPs in Ensemble (92), using the UMD3.1 assembly, and the Variant Effect Predictor function [56].

Abbreviations

BTA :

Bos Taurus autosome

DIM :

Days in milk

FA(s):

Fatty acid(s)

FT-IR :

Fourier transform infrared

GWAS :

Genome wide association study

QTL :

Quantitative trait locus/loci

References

  1. Sun D-W. Fourier transform infrared (FT-IR) spectroscopy. In: Infrared spectroscopy for food quality analysis and control: Academic Press/Elsevier; London, UK; 2009. p. 145–78.

  2. Williams DH, Fleming I. Infrared Spectra. In: Spectroscopic methods in organic chemistry: McGraw-Hill; London, New York; 1980.

  3. Luinge HJ, Hop E, Lutz ETG. Determination of the fat, protein and lactose content of milk using Fourier transform infrared spectrometry. Anal Chim Acta. 1993;284:419–33.

    Article  CAS  Google Scholar 

  4. Soyeurt H, Dardenne P, Dehareng F, Lognay G, Veselko D, Marlier M, et al. Estimating fatty acid content in cow milk using mid-infrared spectrometry. J Dairy Sci. 2006;89:3690–5.

    Article  CAS  PubMed  Google Scholar 

  5. Eskildsen CE, Rasmussen MA, Engelsen SB, Larsen LB, Poulsen NA, Skov T. Quantification of individual fatty acids in bovine milk by infrared spectroscopy and chemometrics: understanding predictions of highly collinear reference variables. J Dairy Sci. 2014;97:7940–51.

    Article  CAS  PubMed  Google Scholar 

  6. Fleming A, Schenkel FS, Chen J, Malchiodi F, Bonfatti V, Ali RA, et al. Prediction of milk fatty acid content with mid-infrared spectroscopy in Canadian dairy cattle using differently distributed model development sets. J Dairy Sci. 2017;100:5073–81.

    Article  CAS  PubMed  Google Scholar 

  7. Rutten MJM, Bovenhuis H, Heck JML, Arendonk JAM. Predicting bovine milk protein composition based on Fourier transform infrared spectra. J Dairy Sci. 2011;94:5683–90.

    Article  CAS  PubMed  Google Scholar 

  8. Ferrand M, Miranda G, Larroque H, Leray O, Guisnel S, Lahalle F, et al. Determination of protein composition in milk by mid-infrared spectrometry. In: Proceedings of the international strategies and new developments in milk analysis VI ICAR Reference Laboratory Network Meeting: 28 May 2012; Cork, vol. 16; 2012. p. 41–5.

    Google Scholar 

  9. Heuer C, Luinge HJ, Lutz ETG, Schukken YH, Maas JH, Wilmink H, et al. Determination of acetone in cow milk by Fourier transform infrared spectroscopy for the detection of subclinical ketosis. J Dairy Sci. 2001;84:575–82.

    Article  CAS  PubMed  Google Scholar 

  10. De Roos APW, Bijgaart HJCM, Hørlyk J, Jong G. Screening for subclinical ketosis in dairy cattle by Fourier transform infrared spectrometry. J Dairy Sci. 2007;90:1761–6.

    Article  PubMed  CAS  Google Scholar 

  11. Grelet C, Bastin C, Gelé M, Davière JB, Johan M, Werner A, et al. Development of Fourier transform mid-infrared calibrations to predict acetone, beta-hydroxybutyrate, and citrate contents in bovine milk through a European dairy network. J Dairy Sci. 2016;99:4816–25.

    Article  CAS  PubMed  Google Scholar 

  12. McParland S, Banos G, Wall E, Coffey MP, Soyeurt HH, et al. The use of mid-infrared spectrometry to predict body energy status of Holstein cows. J Dairy Sci. 2011;94:3651–61.

    Article  CAS  PubMed  Google Scholar 

  13. McParland S, Kennedy E, Lewis E, Moore SG, McCarthy B, O'Donovan MO, et al. Genetic parameters of dairy cow energy intake and body energy status predicted using mid-infrared spectrometry of milk. J Dairy Sci. 2015;98:1310–20.

    Article  CAS  PubMed  Google Scholar 

  14. Bionaz M, Loor JJ. Gene networks driving bovine milk fat synthesis during the lactation cycle. BMC Genomics. 2008;9:366.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  15. Schennink A, Bovenhuis H, Leon-Kloosterziel KM, Arendonk JAM, Visker MHPW. Effect of polymorphisms in the FASN, OLR1, PPARGC1A, PRL and STAT5A genes on bovine milk-fat composition. Anim Genet. 2009;40:909–16.

    Article  CAS  PubMed  Google Scholar 

  16. Rutten MJM, Bovenhuis H, Heck JML, Arendonk JAM. Prediction of beta-lactoglobulin genotypes based on milk Fourier transform infrared spectra. J Dairy Sci. 2011;94:4183–8.

    Article  CAS  PubMed  Google Scholar 

  17. Wang Q, Hulzebosch A, Bovenhuis H. Genetic and environmental variation in bovine milk infrared spectra. J Dairy Sci. 2016;99:6793–803.

    Article  CAS  PubMed  Google Scholar 

  18. Wang Q, Bovenhuis H. Genome wide association study for milk infrared wavenumbers. J Dairy Sci. 2018;101:2260–72.

    Article  CAS  PubMed  Google Scholar 

  19. Soyeurt H, Misztal I, Gengler N. Genetic variability of milk components based on mid-infrared spectral data. J Dairy Sci. 2011;93:1722–8.

    Article  CAS  Google Scholar 

  20. Bittante G, Cecchinato A. Genetic analysis of the Fourier-transform infrared spectra of bovine milk with emphasis on individual wavelengths related to specific chemical bonds. J Dairy Sci. 2013;96:5991–6006.

    Article  CAS  PubMed  Google Scholar 

  21. Zaalberg RM, Shetty N, Janss L, Buitenhuis AJ. Genetic analysis of Fourier transform infrared milk spectra in Danish Holstein and Danish Jersey. J Dairy Sci. 2019;102:1–8.

    Article  CAS  Google Scholar 

  22. De Marchi M, Toffanin V, Cassandro M, Penasa M. Invited review: mid-infrared spectroscopy as phenotyping tool for milk traits. J Dairy Sci. 2014;97:1171–86.

    Article  PubMed  CAS  Google Scholar 

  23. Eskildsen CE, Skov T, Hansen MS, Larsen LB, Poulsen NA. Quantification of bovine milk protein composition and coagulation properties using infrared spectroscopy and chemometrics: a result of collinearity among reference variables. J Dairy Sci. 2016;99:8178–86.

    Article  CAS  PubMed  Google Scholar 

  24. Maurice-Van Eijndhoven MHT, Hiemstra SJ, Calus MPL. Short communication: Milk fat composition of 4 cattle breeds in the Netherlands. J Dairy Sci. 2011;94:1021–5.

    Article  CAS  PubMed  Google Scholar 

  25. Maurice-Van Eijndhoven MHT, Veerkamp RF, Soyeurt H, Calus MPL. Heritability of milk fat composition is considerably lower for Meuse-Rhine-Yssel compared to Holstein Friesian cattle. Livest Sci. 2015;180:58–64.

    Article  Google Scholar 

  26. Soyeurt H, Dardenne P, Gillon A, Croquet C, Vanderick S, Mayeres P, et al. Variation in fatty acid contents of milk and milk fat within and across breeds. J Dairy Sci. 2006;89:4858–65.

    Article  CAS  PubMed  Google Scholar 

  27. Sanchez MP, Govignon-Gion A, Croiseau P, Fritz S, Hozé C, Miranda G, et al. Within-breed and multi-breed GWAS on imputed whole-genome sequence variants reveal candidate mutations affecting milk protein composition in dairy cattle. Genet Sel Evol. 2017;49:68.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  28. Bouwman AC, Visker MHPW, Arendonk JAM, Bovenhuis H. Genomic regions associated with bovine milk fatty acids in both summer and winter milk samples. BMC Genet. 2012;13:93.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  29. Littlejohn MD, Tiplady K, Fink TA, Lehnert K, Lopdell T, Johnson T, et al. Sequence-based association analysis reveals an MGST1 eQTL with pleiotropic effects on bovine milk composition. Sci Rep. 2016;6:25376.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  30. Schopen GCB, Visker MHPW, Mullaart PDE, Arendonk JAM, Bovenhuis H. Whole-genome association study for milk protein composition in dairy cattle. J Dairy Sci. 2011;94:3148–58.

    Article  CAS  PubMed  Google Scholar 

  31. Heck JML, Valenberg HJF, Dijkstra J, Hooijdonk ACM. Seasonal variation in the Dutch bovine raw milk composition. J Dairy Sci. 2009;92:4745–55.

    Article  CAS  PubMed  Google Scholar 

  32. Dadousis C, Biffani S, Cipolat-Gotet C, Nicolazzi EL, Rosa GJM, Gianola D, et al. Genome-wide association study for cheese yield and curd nutrient recovery in dairy cows. J Dairy Sci. 2017;100:1259–71.

    Article  CAS  PubMed  Google Scholar 

  33. Raven L-A, Cocks BG, Hayes BJ. Multibreed genome wide association can improve precision of mapping causative variants underlying milk production in dairy cattle. BMC Genomics. 2014;15:62.

    Article  PubMed  PubMed Central  Google Scholar 

  34. Bouwman AC, Bovenhuis, Visker MHPW, Arendonk JAM. Genome-wide association of milk fatty acids in Dutch dairy cattle. BMC Genet. 2011;12:43.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  35. Schennink A, Heck JML, Bovenhuis H, Visker MHPW, Valenberg HJF, Arendonk JAM. Milk fatty acid unsaturation: genetic parameters and effects of Stearoyl-CoA Desaturase (SCD1) and acyl CoA: Diacylglycerol Acyltransferase 1 (Dgat1). J Dairy Sci. 2008;91:2135–43.

    Article  CAS  PubMed  Google Scholar 

  36. Andersen SK. Vibrational spectroscopy in the analysis of dairy products and wine. In: Chalmers JM, Griffiths PR, editors. Handbook of Vibrational Spectroscopy. Chichester: Wiley; 2002.

    Google Scholar 

  37. 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. J Nutr. 2008;138:1019–24.

    Article  CAS  PubMed  Google Scholar 

  38. Cecchinato A, Ribeca C, Chessa S, Cipolat-Gotet C, Maretto F, Casellas J, et al. Candidate gene association analysis for milk yield, composition, urea nitrogen and somatic cell scores in Brown Swiss cows. Animal. 2014;8:1062–70.

    Article  CAS  PubMed  Google Scholar 

  39. Li C, Sun D, Zhang S, Wang S, Wu X, Zhang Q, et al. Genome wide association study identifies 20 novel promising genes associated with milk fatty acid traits in Chinese Holstein. PLoS One. 2014;9:e96186.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  40. Fitzgerald DK, Brodbeck U, Kiyosawa I, Mawal R, Colvin B, Ebner KE. Alpha-lactalbumin and the lactose synthetase reaction. J Biol Chem. 1970;245:2103–8.

    CAS  PubMed  Google Scholar 

  41. Fuquay JF. Milk Proteins. In: Encyclopedia of Dairy Sciences: Academic Press/Elsevier; London, UK; 2011. p. 758.

    Chapter  Google Scholar 

  42. Hagiwara S-I, Kawai K, Anri A, Nagahata H. Lactoferrin concentration in milk from normal and subclinical mastitic cows. J Vet Med Sci. 2003;65:319–22.

    Article  CAS  PubMed  Google Scholar 

  43. Ganai NA, Bovenhuis H, Arendonk JAM, Visker MHPW. Novel polymorphisms in the bovine ?-lactoglobulin gene and their effects on ?-lactoglobulin protein concentration in milk. Anim Genet. 2009;40:127–33.

    Article  CAS  PubMed  Google Scholar 

  44. Bedere N, Bovenhuis H. Characterizing a region on BTA11 affecting ?-lactoglobulin content of milk using h igh-density genotypingand haplotype grouping. BMC Genet. 2017;18:17.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  45. Poulsen NA, Gustavsson F, Glantz M, Paulsson M, Larsen LB, Larsen MK. The influence of feed and herd on fatty acid composition in 3 dairy breeds (Danish Holstein, Danish Jersey, and Swedish red). J Dairy Sci. 2012;95:6362–71.

    Article  CAS  PubMed  Google Scholar 

  46. Soyeurt H, Gillon A, Vanderick S, Mayeres P, Bertozzi C, Gengler N. Estimation of heritability and genetic correlations for major fatty acids in bovine milk. J Dairy Sci. 2007;90:4435–42.

    Article  CAS  PubMed  Google Scholar 

  47. Hayes BJ, Bowman PJ, Chamberlain AJ, Goddard ME. Invited review: genomic selection in dairy cattle: Progress and challenges. J Dairy Sci. 2009;92:433–43.

    Article  CAS  PubMed  Google Scholar 

  48. Boichard D, Chung H, Dassonneville R, David X, Eggen A, Fritz S, et al. Design of a Bovine low-density SNP array optimized for imputation. PLoS One. 2012;7:e34130.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  49. Daetwyler HD, Capitan A, Pausch H, Stothard P, Binsbergen R, Brøndum RF, et al. Whole-genome sequencing of 234 bulls facilitates mapping of monogenic and complex traits in cattle. Nat Genet. 2014;46:858–65.

    Article  CAS  PubMed  Google Scholar 

  50. Boichard D, Boussaha M, Capitan A, Rocha D, Hozé C, Sanchez MP, et al. Experience from large scale use of the EuroGenomics custom SNP chip in cattle. In: 11th World Congress of Genetics Applied to Livestock Production; 2018. p. 675.

    Google Scholar 

  51. Browning BL, Browning SR. Genotype imputation with millions of reference samples. Am J Hum Genet. 2016;98:116–26.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  52. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. 2018. https://www.R-project.org/. Accessed on 20 Nov 2018.

  53. McArt JAA, Poulsen NA, Larsen MK, Larsen LB, Janss LL, Buitenhuis B. Genetic parameters for milk fatty acids in Danish Holstein cattle based on SNP markers using a Bayesian approach. BMC Genet. 2013;14:79.

    Google Scholar 

  54. Wilmink JBM. Adjustment of test-day milk, fat and protein yield for age, season and stage of lactation. Livest Prod Sci. 1987;16:335–48.

    Article  Google Scholar 

  55. Gebreyesus G, Lund MS, Buitenhuis B, Bovenhuis H, Poulsen NA, Janss LG. Modeling heterogeneous (co)variances from adjacent?SNP groups improves genomic prediction for milk protein composition traits. Genet Sel Evol. 2017;49:89.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  56. McLaren W, Gil L, Hunt SE, Riat H, Ritchie GRS, Thormann A, et al. The Ensembl variant effect predictor. Genome Biol. 2016;17:122.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

Download references

Acknowledgements

We would like to acknowledge Goutam Sahana (QGG, Aarhus University) for his ideas and constructive feedback.

Funding

This study was financially supported by the projects: “FT-IR spektre i mælk: Genetisk variation og effekt på sundhed, frugtbarhed og energibalance” of the Milk Levy Fund, Denmark (2015–2017), and Robust and Efficient Dairy Cows- REFFICO, project number: 34009–14-0848 of Grønt Udviklings- og Demonstrations Program (GUDP) (2015–2018). RMZ was financially supported by The Graduate School of Science and Technology (GSST) of Aarhus University. The funders played no part in the decision making regarding study design, collection of data, data analysis, decision to publish, or preparation of the manuscript.

Author information

Authors and Affiliations

Authors

Contributions

AJB contributed in arranging funding, shaping the idea of the study, and editing the paper. RMZ was responsible for imputation of genotypes, data analysis, and preparation of the manuscript. LJ was responsible for the design of the statistical analysis. All authors read and approved the final manuscript.

Corresponding author

Correspondence to R. M. Zaalberg.

Ethics declarations

Ethics approval and consent to participate

All procedures to collect the Danish samples followed the protocols according to the National Guidelines for Animal Experimentation, and the Danish Animal Experimental Ethics Committee. Infrared spectra were available from milk samples collected from routine records of dairy farms. Genotypes were available from farms, which were part of the routine genotyping for genomic selection, and hence, no specific permission was required.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary information

Additional file 1.

Manhattan plots of % of explained additive genetic variation for Danish Holstein. Scale of y-axis runs from 0 to 1%. The horizontal line indicates the cut-off at 0.35%, which was used to define and select QTL.

Additional file 2.

Manhattan plots of % of explained additive genetic variation for Danish Jersey. Scale of y-axis runs from 0 to 1%. The horizontal line indicates the cut-off at 0.35%, which was used to define and select QTL.

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Zaalberg, R.M., Janss, L. & Buitenhuis, A.J. Genome-wide association study on Fourier transform infrared milk spectra for two Danish dairy cattle breeds. BMC Genet 21, 9 (2020). https://doi.org/10.1186/s12863-020-0810-4

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12863-020-0810-4

Keywords