Breeding for improved protein fractions and free amino acids concentration in bovine milk.

Considerable resources are required to routinely measure detailed milk compositional traits. Hence, an insufficient volume of phenotypic data can hinder genetic progress in these traits within dairy cow breeding programmes. The objective of the present study was to quantify the opportunities for breeding for improved milk protein and free amino acid (FAA) composition by exploiting mid-infrared spectroscopy (MIRS) predictions routinely recorded from milk samples. Genetic parameters for protein fractions and FAA composition were estimated using 134,546 test-day records from 16,166 lactations on 9,572 cows using linear mixed models. Heritability of MIRS-predicted protein fractions ranged from 0.19 (α-lactalbumin) to 0.55 (β-lactoglobulin A), while heritability of MIRS-predicted FAA ranged from 0.08 for glycine to 0.29 for glutamic acid. Genetic correlations among the MIRS-predicted FAA were moderate to strong ranging from -0.44 (aspartic acid and lysine) to 0.97 (glutamic acid and total FAA). Adjustment of the genetic correlations for the genetic merit of 24-h milk yield did not greatly affect the correlations. Results from the current study highlight the presence of exploitable genetic variation for both protein fractions and FAA in dairy cow milk. Besides, the direction of genetic correlations reveals that breeding programmes directly selecting for greater milk protein concentration carry with them favourable improvement in casein and whey fractions.

INTRODUCTION

Although most dairy cow breeding objectives include milk protein and fat concentration and yield at the macro‐level (Cole & VanRaden, 2018; Miglior et al., 2005), few consider the detailed milk composition. However, detailed milk composition affects the quality of a variety of dairy products that can be manufactured (Bonfatti, Di Martino, et al. 2011; Wedholm et al., 2006). For example, the concentration of casein (CN) in milk has a favourable effect on the quantity of protein transferred from milk into cheese curd. High concentrations of αs1‐casein (αs1‐CN), β‐casein (β‐CN), κ‐casein (κ‐CN) and β‐lactoglobulin B (β‐LG B) are known to increase cheese yield (Wedholm et al., 2006).

Milk rich in FAA tends to present unfavourable cheese‐making characteristics mainly due to protein hydrolysis (Auldist et al., 1995; Urech et al., 1999). Nevertheless, presence of FAA is important for some products like infant formulas even though the milk from mammals differs in FAA composition and concentration. In fact, cow milk presents lower quantity of FAA than human milk (McDermott et al., 2016; van Sadelhoff et al., 2020). Some FAA are known to have a key role in the development of the immune system, particularly in infants (van Sadelhoff et al., 2020). However, bovine milk FAA have generally been overlooked thus far, even if some of them are essential or conditionally essential, such as leucine, valine and glycine. For nutritional reasons, FAA‐enriched infant formulas and amino acid‐based products may be valuable when breastfeeding is not possible (van Sadelhoff et al., 2020). Despite the high cost, amino acid‐based supplements provide proteins in non‐allergenic forms and are thus suitable for people with bovine milk protein allergies (Verduci et al., 2019). Given these considerations, it is worth investigating whether and how bovine milk FAA profile could be manipulated through breeding.

Interbreed differences in the concentration of protein fractions and FAA have been documented (Auldist et al., 2004; Lopez‐Villalobos, 2012; McDermott et al., 2017) and individual protein fractions, including α‐lactalbumin (α‐LA) and γ‐casein (γ‐CN), are known to be heritable (Graml & Pirchner, 2003). For example, Schopen et al. (2009) reported a heritability of between 0.25 (β‐CN) and 0.80 (β‐LG) for protein fractions quantified using capillary zone electrophoresis. Similarly, Huang et al. (2012) reported that the proportion of total phenotypic variance attributable to additive genetic variation varied between 0.33 (αs1‐CN, β‐CN and α‐LA) and 0.68 (β‐LG). Differences in the gold standard methods used to quantify milk protein fractions, as well as the characteristics of the studied population such as the breeds investigated, parities and stages of lactation represented, could have contributed to the difference in estimates (Schopen et al., 2009). Heritability estimates calculated by Bonfatti, Cecchinato, et al. (2011) in a population of Italian Simmental dairy cows, for the relative proportions of protein fractions expressed as a percentage of total CN (0.18 for γ‐CN to 0.69 for β‐CN) were greater than those for protein fraction contents expressed as g/L of milk (0.11 g/L for α‐LA to 0.53 g/L for κ‐CN). At present, little is known about the genetic parameters of FAA.

The usefulness of mid‐infrared spectroscopy (MIRS) analysis of milk to predict milk composition and animal health traits on a large scale is now well established (De Marchi et al., 2014; Niero et al., 2016; Rutten et al., 2011; Soyeurt et al., 2011). Despite this, genetic parameters of MIRS‐predicted protein fractions have been reported in only a limited number of studies (Bonfatti et al., 2017; Sanchez et al., 2017), while no estimates are available in the literature for FAA.

Therefore, the objective of the present study was to quantify the potential of breeding for milk protein fractions and FAA in dairy cows, but doing so by exploiting the MIRS predictions routinely available for individual cows.

MATERIALS AND METHODS

Milk samples and datasets

Gold standard dataset and development of MIRS prediction models

Seven hundred and fifteen milk samples were collected from seven research dairy farms operated by the Teagasc Animal and Grassland Research and Innovation Centre (Moorepark, Fermoy, Co. Cork) between August 2013 and August 2014, inclusive. Individual milk proteins (αs1‐CN, αs2‐CN, β‐CN, κ‐CN, α‐LA, and β‐LG variant A and variant B) were determined in 557 of the 715 milk samples using reverse‐phase high‐performance liquid chromatography (HPLC). The FAA‐related traits, namely total FAA, glutamic acid (or glutamate, Glu), glycine (Gly), lysine (Lys), arginine (Arg), aspartic acid (Asp), serine (Ser) and valine (Val), were quantified in all 715 milk samples using cation‐exchange HPLC. Further details regarding gold standard analyses for both protein fractions and FAA are described in detail in McDermott et al. (2016). Each milk sample was also analysed using a mid‐infrared spectrometer MilkoScan FT6000 located at the Teagasc Animal and Grassland Research and Innovation Centre (Moorepark, Fermoy Co. Cork, Ireland) for milk chemical composition quantification and spectra determination. Gold standard data were matched to their respective milk spectra, which were subsequently used to develop equations to predict individual and groups of milk proteins and FAA using partial least squares regression analysis (PROC PLS) in the SAS software v. 9.4 (SAS Institute Inc., Cary, NC). More details on the methodology used to develop the MIRS prediction models are outlined in McDermott et al. (2016). Overall, the accuracy of prediction, defined as the correlation coefficient between the reference and predicted values using (leave‐one‐out) cross‐validation of protein fractions, ranged from 0.43 (β‐LG variant A) to 0.76 (total whey and total LG). Cross‐validation prediction accuracy for FAA varied from 0.51 (Ser) to 0.75 (Gly).

Large dataset with MIRS‐predicted traits

Between the years 2013 and 2015, inclusive, 171,279 milk samples from 17,353 lactations on 10,162 cows were collected from 7 research herds and 69 commercial dairy herds (Visentin et al., 2017). Animals in the research herds were subjected to different experimental treatments, based mainly on different feeding and management strategies. All animals were milked twice daily. The commercial dairy farms (n = 69) were located in the south‐west of Ireland (Munster Region) and represented typical Irish dairy farms.

Individual milk samples were sporadically taken and delivered to the laboratory of the Teagasc Animal and Grassland Research and Innovation Centre. All milk samples collected in both research and commercial herds were analysed using the same Fourier transform infrared spectrometer MilkoScan FT6000 (Foss Electronic A/S, Hillerød, Denmark) based at the Animal and Grassland Research Centre, Teagasc. Somatic cell count (SCC) was determined using a Fossomatic (Foss Electronic A/S, Hillerød, Denmark). Prediction equations developed by McDermott et al. (2016) based on the gold standard dataset were applied to these milk spectra in order to predict (a) groups of proteins (total CN, total whey, total LG), (b) individual proteins (αs1‐CN, αs2‐CN, β‐CN, κ‐CN, α‐LA, β‐LG A and β‐LG B) and (c) FAA traits (total FAA, Glu, Gly, Lys, Arg, Asp, Ser and Val). Spectral variables and information on 24‐h milk yield along with cow calving date and parity number were available for all milk samples. The average number of milk samples available per cow was equal to 72 for cows in the experimental farms (median = 57) and 5.15 for those in the commercial herds (median = 4).

In the final dataset, spectra with a Mahalanobis distance relative to the centroid of the gold standard dataset greater than three were discarded prior the application of MIRS prediction models (Williams, 2007). Furthermore, MIRS‐predicted values of proteins and FAA, as well as 24‐h milk yield and composition, greater than three standard deviations from the mean were set as missing. Moreover, the skewed‐distributed values of Glu, Gly, Lys, Arg, Asp, Ser and Val were loge‐transformed. Breed composition was available for each cow and only data from Holstein, Friesian, Jersey, Norwegian Red, Montbeliarde cows as well as their crosses sampled between 5 and 305 days in milk (DIM) and from parities ≤10 were retained. Parities greater than 6 were grouped together for analysis and only contemporary groups with a minimum of ten records were used. In particular, as both research and commercial farms were present, contemporary groups were generated either by combining the experimental treatment and the test date or by combining the commercial herd ID and the test date. Somatic cell count (SCC, n/μl) was converted to somatic cell score (SCS), as conventionally: SCC = 3 + log2(SCC/100). Following all edits, 134,546 records from 16,166 lactations on 9,572 cows were available; of these, 40,260 were from commercial herds (i.e. 13,337 lactations and 8,260 cows).

Variance components

Pedigree information was retrieved from the Irish Cattle Breeding Federation database and each animal was traced back (where available) up to at least four generations, resulting in a pedigree file of 33,949 and 8,210 individuals for the large dataset and gold standard dataset, respectively. Coefficients of general heterosis and recombination loss were quantified based on the equations presented by both VanRaden and Sanders (2003) and Visentin et al. (2017). Genetic, permanent environmental and residual (co)variances for protein fractions and FAA composition were quantified using linear mixed models in ASReml v. 3.0 (Gilmour et al., 2011).

Large dataset

Fixed effects included in the linear mixed models were the contemporary group, the interaction between parity (1, 2, 3, 4, 5 and ≥6) and stage of lactation (10 classes with length of 30 d, from 5 to 305 DIM), milking time (AM or PM), proportion of cow breed (Friesian, Jersey, Norwegian Red, Montbeliarde and other; Holstein was not included to avoid linear dependencies), and general heterosis and recombination loss coefficients of the cow. The latter two effects were fitted as linear covariates. Random effects included in the models were the direct additive genetic effect of the animal and both a within‐ and an across‐lactation cow permanent environmental effect.

Genetic and phenotypic (co)variances among the MIRS‐predicted protein fractions were estimated using a series of trivariate (3  3) analyses that included two MIRS‐predicted protein fractions as dependant variables, along with the total protein concentration (%). Covariances between MIRS‐predicted FAA were estimated in a similar manner. Thus, two MIRS‐predicted FAA together with the 24‐h milk yield were the dependent variables. In matrix notation, the model was where, for traits 1, 2 and 3, is the vector of phenotypic observations of the dependent variable; is the vector of fixed effects; is the vector of random additive genetic effect; is the vector of permanent environmental effect of the cow within a given lactation; is the vector of permanent environmental effect of the cow across lactations; is the vector of random residuals. Incidence matrices relating the corresponding effects to the dependent variable were X, Z, W and S, respectively. The following structure was assumed for the (co)variance between the random effects: where for the traits 1, 2 and 3: σa12, σa13 and σa23 are the additive genetic covariances, σp12, σp13 and σp23 are the within‐lactation permanent environmental covariances; σc12, σc13 and σc23 are the across‐lactation permanent environmental covariances; σe12, σe13 and σe23 are the residual covariances; σ2 a, σ2 p, σ2 c and σ2 e indicate the additive genetic variance, the within‐ and across‐lactation permanent environmental variance and the residual variance; A was the pedigree‐based relationship matrix of size equal to number of individuals (n = 33,949), while the I denoted identity matrixes of appropriate order, that is equal to number of lactations, number of cows and number of records for σ2 p, σ2 c and σ2 e, respectively.

Genetic correlations among the MIRS‐predicted protein fractions were adjusted for their respective genetic correlation with total protein concentration as follows: where is the partial correlation coefficient between traits x and y independent of trait z, represents the correlation coefficient between traits x and y, is the correlation coefficient between traits x and z, and is the correlation coefficient between traits y and z. In a similar manner, genetic correlations among the MIRS‐predicted FAA were adjusted for their respective genetic correlation with 24‐h milk yield.

The phenotypic and genetic variance adjusted for their respective correlation (phenotypic or genetic) with protein concentration (protein fractions) or 24‐h milk yield (FAA) were calculated for each trait afterwards, as: where is the variance (phenotypic or genetic) estimated by the REML algorithm for each trait, and is the correlation (phenotypic or genetic) between trait x (protein fraction or FAA) and trait y (protein concentration or 24‐h milk yield). Heritability was the ratio of the additive genetic variance (adjusted or unadjusted) to the phenotypic variance calculated by summing up of the additive genetic (adjusted or unadjusted), the within‐ and the across‐lactation permanent environmental and the residual variances.

Gold standard dataset

The model used for genetic parameter estimation in the gold standard dataset differed from that used in the large dataset. In particular, the contemporary group was replaced by the fixed effect of herd (n = 7) and the two random permanent environmental effects were not considered, since only one record per cow (in a single lactation) was present. Heritability was calculated as the ratio between the additive genetic variance to the sum of the additive genetic and residual variance. The genetic correlation between the gold standard trait and its MIRS prediction was estimated using (co)variance components obtained from the bivariate analyses.

RESULTS

Mean and variation

The mean of the gold standard measurement and the mean of the respective MIRS‐predicted milk protein fractions were comparable. In fact, the mean value of the gold standard total CN and the mean of MIRS‐predicted total CN were 36.92 g/L and 36.16 g/L, respectively (Tables 1 and 2, respectively). Similarly, the mean of the gold standard FAA was similar to the mean of their relative MIRS predictions (Tables 1 and 2). The genetic standard deviation of protein fractions ranged from 0.44 g/L (β‐LG A) to 2.47 g/L (total CN; Table 1) for the gold standard data. The MIRS predictions ranged from 0.04 g/L (α‐LA) to 1.85 g/L (total CN) for the unadjusted analyses (Table 2). The genetic standard deviation of the gold standard and the respective MIRS‐predicted FAA variables were also similar, with the exception of total FAA. In fact, gold standard total FAA had a genetic standard deviation of 10.48 μg/mL, which was approximately twice that for the MIRS‐predicted values (5.05 μg/ml). The coefficient of genetic variation differed among traits ranging from 3.63 (β‐CN) to 47.08 (β‐LG B) for the gold standard protein fractions and from 4.02 (Gly) to 14.30 (total FAA) for the gold standard FAA (Table 1).

TABLE 1

Number of records (n), mean, genetic standard deviation (σg), heritability (standard error) and coefficient of genetic variation (CVg) for milk yield and composition, gold standard protein fractions and free amino acids (FAA) and the genetic correlation (r g) between gold standard versus MIRS‐predicted FFA and protein fractions

Trait	N	Mean	σg	Heritability	CVg, %	r g
24‐h Milk yield, kg/d	497	22.08	2.64	0.16 (0.10)	11.96	–
Protein, %	585	3.74	0.22	0.45 (0.12)	5.88	–
Casein, %	585	2.82	0.19	0.51 (0.12)	6.74	–
Fat, %	584	4.71	0.53	0.30 (0.10)	11.25	–
Lactose, %	584	4.61	0.07	0.15 (0.10)	1.52	–
SCS, log2 units	252	2.48	0.23	0.03 (0.12)	9.27	–
Protein fractions, g/L
Total casein	441	36.92	2.47	0.20 (0.12)	6.69	0.69 (0.18)
αs1‐casein	442	14.18	1.19	0.22 (0.12)	8.39	0.71 (0.15)
αs2‐casein	442	3.67	0.51	0.38 (0.14)	13.90	0.82 (0.15)
β‐casein	441	12.95	0.47	0.05 (0.09)	3.63	0.63 (0.15)
κ‐casein	440	6.01	0.60	0.21 (0.14)	9.98	0.64 (0.24)
Total whey	437	6.18	0.97	0.41 (0.15)	15.70	0.73 (0.15)
α‐Lactalbumin	441	1.13	0.27	0.24 (0.16)	23.89	0.45 (0.33)
Total lactoglobulin	438	5.04	0.84	0.38 (0.15)	16.67	0.76 (0.13)
β‐Lactoglobulin A	441	2.59	0.44	0.15 (0.12)	16.99	0.68 (0.17)
β‐Lactoglobulin B	400	2.74	1.29	0.69 (0.16)	47.08	0.77 (0.11)
FAA, μg/mL a
Total FAA	568	73.28	10.48	0.19 (0.11)	14.30	0.53 (0.12)
Glutamic acid	568	27.41	0.29	0.41 (0.13)	5.46	0.83 (0.11)
Glycine	565	5.64	0.15	0.11 (0.09)	4.02	0.78 (0.32)
Lysine	565	3.41	0.29	0.23 (0.12)	8.98	0.38 (0.35)
Arginine	507	3.15	0.20	0.27 (0.13)	6.36	0.51 (0.26)
Aspartic acid	483	2.12	0.29	0.24 (0.13)	10.54	0.52 (0.28)
Serine	488	1.22	0.30	0.37 (0.15)	13.64	0.38 (0.11)
Valine	524	1.39	0.24	0.15 (0.11)	10.32	0.46 (0.17)

Traits were log‐transformed before analysis.

TABLE 2

Number of records (n), mean, genetic standard deviation (σg) with the unadjusted heritability (standard error) and repeatability (standard error) for milk yield and composition, MIRS‐predicted protein fractions and MIRS‐predicted free amino acids (FAA). Heritability and genetic standard deviation after adjustment are also reported

Trait	n	Mean	Unadjusted			Adjusted
			σg	Heritability	Repeatability	σg	Heritability
24‐h Milk yield, kg/day	133,999	20.54	1.53	0.17 (0.02)	0.54 (0.01)	–	–
Protein, %	125,045	3.72	0.18	0.44 (0.02)	0.56 (0.01)	–	–
Casein, %	125,484	2.82	0.15	0.39 (0.02)	0.52 (0.01)	–	–
Fat, %	111,055	4.57	0.36	0.30 (0.01)	0.32 (0.01)	–	–
Lactose, %	110,954	4.77	0.09	0.36 (0.02)	0.49 (0.01)	–	–
SCS, units	63,505	2.50	0.26	0.05 (0.01)	0.47 (0.01)	–	–
Protein, g/L
Total casein	134,395	36.16	1.85	0.42 (0.02)	0.55 (0.01)	0.27	0.11
αs1‐casein	134,409	13.70	0.79	0.44 (0.02)	0.54 (0.01)	0.59	0.36
αs2‐casein	134,482	3.66	0.18	0.36 (0.02)	0.44 (0.01)	0.14	0.28
β‐casein	134,235	12.96	0.54	0.45 (0.02)	0.48 (0.01)	0.41	0.30
κ‐casein	134,463	6.03	0.34	0.36 (0.02)	0.46 (0.01)	0.25	0.27
Total whey	134,383	6.13	0.52	0.37 (0.02)	0.45 (0.01)	0.50	0.36
α‐lactalbumin	134,106	1.11	0.04	0.19 (0.01)	0.24 (0.01)	0.03	0.17
Total lactoglobulin	134,236	5.09	0.51	0.39 (0.02)	0.48 (0.01)	0.48	0.37
β‐lactoglobulin A	134,489	2.39	0.19	0.46 (0.02)	0.54 (0.01)	0.16	0.42
β‐lactoglobulin B	128,391	2.61	0.56	0.42 (0.02)	0.52 (0.01)	0.55	0.42
FAA, μg/ml b
Total FAA	134,546	53.60	5.05	0.24 (0.01)	0.36 (0.01)	5.05	0.24
Glutamic acid	134,426	30.93	4.40	0.29 (0.02)	0.40 (0.01)	4.38	0.29
Glycine	133,650	8.09	0.66	0.08 (0.01)	0.15 (0.01)	0.66	0.08
Lysine	134,105	4.75	0.53	0.17 (0.01)	0.23 (0.01)	0.52	0.16
Arginine	134,425	3.38	0.23	0.15 (0.01)	0.24 (0.01)	0.23	0.15
Aspartic acid	134,433	2.73	0.39	0.29 (0.02)	0.38 (0.01)	0.39	0.28
Serine	112,918	2.74	0.22	0.15 (0.01)	0.20 (0.01)	0.22	0.15
Valine	133,957	1.52	0.16	0.13 (0.01)	0.25 (0.01)	0.16	0.12

Protein factions and FAA were adjusted for protein concentration and 24‐h milk yield, respectively. Heritability considered the genetic variance calculated as per Equation 2.

Traits were log‐transformed before analysis.

Heritability and repeatability estimates

Heritability estimates of milk yield and composition from the gold standard dataset (Table 1), that is the one containing just the reference values of milk protein fractions and FAA, were similar to those from the large dataset with research and commercial cow data (Table 2). The heritability estimates did, however, differ in terms of standard errors, which was due to the fewer records in the gold standard dataset. Heritability estimates for the gold standard protein fractions, albeit characterized by large standard errors, ranged from 0.05 (β‐CN) to 0.69 (β‐LG B; Table 1). Heritability estimates for the MIRS‐predicted protein fractions (Table 2) were higher than those estimated for the gold standard protein fractions (Table 1), although the range of overall heritability estimates for the MIRS‐predicted protein fractions was smaller (0.19 for α‐LA to 0.46 for β‐LG A; Table 2) than the range for the respective gold standard measures (0.05 for β‐CN to 0.69 for β‐LG B; Table 1). Similar to protein fractions, heritability estimates for MIRS‐predicted FAA had a narrower range (0.08 for Gly to 0.29 for Asp and Glu; Table 2) than the respective gold standard range (0.11 for Gly to 0.41 for Glu; Table 1). Repeatability estimates for MIRS‐predicted protein fractions and MIRS‐predicted FAA were moderate (0.24 for α‐LA to 0.55 for total CN; 0.15 for Gly to 0.40 for Glu; Table 2). The estimated genetic standard deviation for each protein fraction trait, genetically independent of protein concentration, was less than the respective unadjusted measure and this was also reflected in lower heritability estimates; only a small reduction in genetic standard deviation and heritability was observed for the whey protein fractions (i.e. α‐LA, β‐LG A and β‐LG B; Table 2). There was also little impact on the heritability estimates for FAA when adjusted for differences in the genetic merit of 24‐h milk yield (Table 2).

Genetic correlations

Moderate‐to‐strong genetic correlations existed between the gold standard protein fractions and their respective MIRS‐predicted protein fractions, ranging from 0.45 (α‐LA) to 0.82 (αs1‐CN; Table 1). The genetic correlation between the gold standard Glu and MIRS‐predicted Glu was strong (0.83), similar to the correlation between the gold standard and MIRS‐predicted Gly (0.78; Table 1). Protein fractions were negatively correlated with 24‐h milk yield, but they were positively associated with protein and CN concentration (Table 3). Protein fractions were weakly correlated with both lactose concentration and SCS. With the exception of both β‐CN and α‐LA, all protein fractions were negatively correlated with lactose concentration. Individual FAA were weakly to moderately genetically correlated with 24‐h milk yield, protein, CN, fat concentration, lactose concentration and SCS (Table 4). Correlations between MIRS‐predicted FAA and SCS in general did not differ from zero with the exception of the genetic correlation between MIRS‐predicted Val and SCS (0.39 ± 0.09; Table 4).

TABLE 3

Genetic correlations (standard errors) of MIRS‐predicted protein fractions with 24‐h milk yield, concentration of protein, casein, fat and lactose, as well as somatic cell score (SCS)

Trait	24‐h milk yield	Protein	Casein	Fat	Lactose	SCS
Total casein	−0.57 (0.04)	0.99 (0.001)	0.99 (0.002)	0.72 (0.02)	−0.09 (0.04)	0.11 (0.08)
αs1‐casein	−0.58 (0.04)	0.99 (0.001)	0.98 (0.002)	0.74 (0.03)	−0.12 (0.04)	0.10 (0.08)
αs2‐casein	−0.59 (0.04)	0.95 (0.01)	0.96 (0.01)	0.79 (0.03)	−0.05 (0.04)	0.10 (0.08)
β‐casein	−0.57 (0.04)	0.94 (0.01)	0.81 (0.03)	0.69 (0.02)	0.06 (0.04)	0.07 (0.08)
κ‐casein	−0.50 (0.04)	0.97 (0.01)	0.94 (0.01)	0.62 (0.03)	−0.20 (0.04)	0.10 (0.08)
Total whey	−0.30 (0.05)	0.49 (0.03)	0.46 (0.03)	0.39 (0.03)	−0.07 (0.04)	0.09 (0.08)
α‐lactalbumin	−0.45 (0.05)	0.55 (0.03)	0.59 (0.03)	0.68 (0.03)	0.35 (0.04)	0.11 (0.08)
Total lactoglobulin	−0.29 (0.04)	0.49 (0.03)	0.46 (0.03)	0.38 (0.03)	−0.12 (0.04)	0.08 (0.08)
β‐lactoglobulin A	−0.56 (0.04)	0.86 (0.01)	0.87 (0.01)	0.75 (0.02)	−0.09 (0.04)	0.15 (0.08)
β‐lactoglobulin B	−0.08 (0.05)	0.18 (0.03)	0.15 (0.04)	0.11 (0.04)	−0.11 (0.04)	0.07 (0.08)

TABLE 4

Genetic correlations (standard errors) between log‐transformed MIRS‐predicted free amino acids (FAA) with 24‐h milk yield, concentration of protein, casein, fat and lactose, as well as somatic cell score (SCS)

Trait	24‐h Milk yield	Protein	Casein	Fat	Lactose	SCS
Total FAA	0.01 (0.06)	−0.10 (0.04)	−0.09 (0.04)	−0.07 (0.04)	0.29 (0.04)	0.07 (0.09)
Glutamic acid	0.10 (0.05)	−0.18 (0.04)	−0.17 (0.04)	−0.15 (0.04)	0.33 (0.04)	0.08 (0.09)
Glycine	0.05 (0.06)	−0.19 (0.04)	−0.17 (0.04)	−0.24 (0.04)	0.28 (0.05)	−0.15 (0.10)
Lysine	−0.33 (0.05)	0.52 (0.03)	0.51 (0.03)	0.51 (0.03)	−0.16 (0.04)	0.07 (0.09)
Arginine	−0.18 (0.06)	0.21 (0.04)	0.17 (0.04)	0.40 (0.04)	−0.33 (0.04)	0.13 (0.12)
Aspartic acid	0.14 (0.05)	−0.19 (0.04)	−0.20 (0.04)	−0.16 (0.04)	0.25 (0.04)	−0.01 (0.09)
Serine	−0.13 (0.06)	−0.15 (0.04)	−0.11 (0.04)	0.24 (0.04)	0.26 (0.04)	0.18 (0.09)
Valine	−0.21 (0.06)	0.32 (0.04)	0.33 (0.04)	0.25 (0.04)	−0.16 (0.04)	0.39 (0.09)

Genetic correlations among the MIRS‐predicted protein fractions were weak to strong (Table 5). An almost unity genetic correlation (0.99) existed between protein concentration and total CN, between total CN and αs1‐CN, and between total LG and total whey. The genetic correlation between total LG and total whey did not change after adjusting for their respective genetic correlation with protein concentration, albeit most of the correlations among the protein fractions weakened when calculated genetically independent of protein concentration. For example, the strong positive genetic correlation that existed between β‐CN and αs2‐CN (0.89) became negative (−0.05) once adjusted for genetic merit for protein concentration (Table 5). Genetic correlations among the MIRS‐predicted FAA were weak to strong and ranged from −0.44 (Asp and Lys) to 0.97 (Glu and Total FAA; Table 6). Adjusting the correlations for the genetic merit of 24‐h milk yield did not greatly affect the correlation estimates (Table 6). For example, the unadjusted correlation between Arg and Val was 0.63 and changed to 0.62 when adjusted for the genetic merit of 24‐h milk yield.

TABLE 5

Genetic correlations (standard errors) among the concentration of MIRS‐predicted protein fractions (above diagonal); estimates adjusted for milk protein concentration are reported below diagonal

Trait	Total CN	αs1‐CN	αs2‐CN	β‐CN	κ‐CN	Total‐W	α‐LA	Total LG	β‐LG A	β‐Lg B
Total CN	–	0.99 (0.00)	0.92 (0.01)	0.97 (0.002)	0.96 (0.002)	0.47 (0.03)	0.60 (0.03)	0.47 (0.03)	0.85 (0.01)	0.17 (0.03)
αs1‐CN	0.91	–	0.92 (0.01)	0.96 (0.002)	0.95 (0.002)	0.51 (0.03)	0.61 (0.03)	0.50 (0.03)	0.81 (0.01)	0.18 (0.03)
αs2‐CN	−0.24	−0.24	–	0.89 (0.01)	0.90 (0.01)	0.42 (0.03)	0.56 (0.03)	0.41 (0.03)	0.88 (0.01)	0.10 (0.03)
β‐CN	0.79	0.63	−0.05	–	0.80 (0.02)	0.48 (0.03)	0.70 (0.02)	0.46 (0.03)	0.86 (0.01)	0.13 (0.04)
κ‐CN	0.07	−0.09	−0.12	−0.22	–	0.47 (0.03)	0.38 (0.04)	0.47 (0.03)	0.72 (0.02)	0.22 (0.03)
Total‐W	0.09	0.16	−0.17	0.05	−0.01	–	0.47 (0.03)	0.99 (0.00)	0.35 (0.03)	0.91 (0.01)
α‐LA	0.53	0.56	0.14	0.63	−0.37	0.25	–	0.44 (0.03)	0.72 (0.02)	0.17 (0.04)
Total LG	0.04	0.12	−0.17	−0.02	0.01	0.99	0.24	–	0.34 (0.03)	0.92 (0.01)
β‐LG A	0.11	0.25	0.44	0.29	−0.79	−0.01	0.58	−0.18	–	−0.03 (0.03)
β‐LG B	−0.05	−0.01	−0.20	−0.13	0.18	0.96	0.09	0.90	−0.36	–

Total casein (Total CN), αs1‐casein (αs1‐CN), αs2‐casein (αs2‐CN), β‐casein (β‐Casein), κ‐casein (κ‐CN), total whey (Total‐W), α‐lactalbumin (α‐LA), total lactoglobulin (Total LG) and β‐lactoglobulin A (β‐LG A) and β‐lactoglobulin B (β‐LG B).

Adjusted genetic correlations calculated as per Equation 1.

TABLE 6

Genetic correlations (standard errors) among the concentration of MIRS‐predicted free amino acids (FAA; above diagonal); estimates adjusted for 24‐h milk yield are reported below diagonal

Trait	Total FAA	Glutamic acid	Glycine	Lysine	Arginine	Aspartic acid	Serine	Valine
Total FAA	–	0.97 (0.002)	0.56 (0.04)	−0.20 (0.04)	0.31 (0.04)	0.88 (0.01)	0.60 (0.03)	0.49 (0.04)
Glutamic acid	0.99	–	0.53 (0.04)	−0.35 (0.04)	0.19 (0.04)	0.91 (0.01)	0.57 (0.03)	0.42 (0.04)
Glycine	0.56	0.53	–	−0.36 (0.04)	−0.27 (0.05)	0.57 (0.03)	0.37 (0.04)	−0.09 (0.05)
Lysine	−0.20	−0.31	−0.36	–	0.56 (0.03)	−0.44 (0.04)	0.06 (0.04)	0.49 (0.04)
Arginine	0.24	0.21	−0.22	0.53	–	0.13 (0.04)	0.32 (0.04)	0.63 (0.03)
Aspartic acid	0.90	0.90	0.57	−0.42	0.16	–	0.32 (0.04)	0.18 (0.04)
Serine	0.38	0.59	0.39	0.02	0.30	0.30	–	0.49 (0.04)
Valine	0.49	0.45	−0.09	0.45	0.62	0.19	0.27	–

Adjusted genetic correlations calculated as per Equation 1.

DISCUSSION

The present study aimed to quantify the extent of genetic variability in detailed milk protein and FAA composition predicted by MIRS in individual cow milk samples. Nowadays, global milk recording programmes rely on infrared tools for determination of fat and protein content of samples intended to herd testing and genetic evaluations (De Marchi et al., 2014; Franzoi et al., 2021). Currently, protein fractions and FAA are not either determined by MIRS on a routine basis or used for within‐country genetic evaluations. For both protein fractions and FAA, the genetic correlations estimated in the present study between both the gold standard and MIRS prediction were moderate to strong, demonstrating that the MIRS predictions are genetically very similar to their corresponding gold standard measures. All traits were heritable and exhibited considerable and exploitable genetic variation; therefore, MIRS is a viable method to collect a large quantity of accurate data for use in genetic evaluations with the goal of improving the added value of granular milk.

Response to genetic selection

Response to genetic selection is a function of the extent of genetic variability, the accuracy of selection, selection intensity and generation interval (Rendel & Robertson, 1950). Therefore, to achieve genetic gain, the trait must exhibit genetic variation and sufficient phenotypic or genomic information on the trait must be available to ensure a high accuracy of selection (Costa et al., 2020; Simm et al., 2020). Considerable genetic variation clearly exists for all milk quality traits examined in the present study with the traits being low to moderately heritable. Therefore, this suggests that there is potential to alter both the protein composition and FAA composition in bovine milk using selective breeding. The greater the heritability, the fewer progeny records required to achieve a high accuracy of selection. In other words, the greater the heritability, the smaller the reference population required to develop genomic predictions (Daetwyler et al., 2008).

Heritability estimates of all the MIRS‐predicted protein fractions in the present study (0.19 to 0.46; Table 2) were higher than the heritability estimate for 24‐h milk yield (0.17; Table 2). Moreover, the heritability of 24‐h milk yield and protein concentration (0.44; Table 2) in the present study were in the range of the heritability estimates for the gold standard protein fractions (0.05 to 0.69; Table 1), but only the heritability of 24‐h milk yield was in the range of the heritability of gold standard FAA (0.11 to 0.41; Table 1) and of the MIRS‐predicted FAA (0.08 to 0.29; Table 2). Previous heritability estimates of protein fractions in bovine milk have been documented to be moderate to high ranging from 0.25 for β‐CN to 0.80 for total LG (Huang et al., 2012; Schopen et al., 2009), although they differed across studies. Recent heritability estimates for gold standard protein fractions in milk from Simmental cows ranged from 0.18 (κ‐CN) to 0.68 (αs1‐CN; Bonfatti, Cecchinato, et al. (2011) and in milk from both Holstein and Holstein‐Jersey crosses the range was between 0.33 (αs1‐CN, β‐CN and α‐LA) and 0.68 (β‐LG; Huang et al., 2012).

In the present study, all protein fractions were positively genetically correlated with protein concentration, contradicting previous findings by Schopen et al. (2009) who estimated genetic parameters for the same traits but determined with capillary zone electrophoresis in milk from 1,940 primiparous cows. By way of example, αs1‐CN, β‐CN, α‐LA and total whey were negatively genetically correlated with the protein concentration (Schopen et al., 2009), but those estimates were characterized by quite large standard errors and sometimes did not differ from zero. Genetic correlations among the protein fractions obtained in the present study were also different from those of Schopen et al. (2009) and Bonfatti,Cecchinato, et al. (2011) in dairy cows. Again, different methods were used to determine milk protein fractions across studies, as well as the characteristics of the studied population such as breeds, parities and seasons represented all of which could have contributed to the difference in correlations. Indeed, Bonfatti, Cecchinato, et al. (2011) used data from protein fractions determined by high‐performance liquid chromatography from 2,167 Simmental cows. In the present study, protein fractions were predicted by MIRS on up to 134,100 milk samples from Holstein, Friesian, Jersey, Norwegian Red and Montbeliarde cows as well as their crosses across a wide range of parities and lactation stages.

Many current national breeding programmes indirectly select for protein concentration through a negative weighting on milk yield concurrent with a positive weighting on protein yield (Cole & VanRaden, 2018). The coefficients of genetic variation for milk yield and protein concentration in the present study were 7.45% and 4.84%, respectively. Such values were similar or even smaller than coefficient of genetic variation for the gold standard total casein and total whey (6.69% and 15.70%, respectively), as well as lower than the coefficient of genetic variation for FAA (14.30%). Favourable trends in genetic gain for milk yield in dairy cows are well documented (Berry, 2008; Berry et al., 2014; Norman & Powell, 1999), with lactation milk yield per cow doubling over the past 40 years (Oltenacu & Broom, 2010). Given the coefficient of variation for milk yield, protein fractions and FAA, a similar genetic gain may be reached for protein fractions and FAA if these traits were included in a selection index with a high accuracy of selection. Satisfactory selection accuracy could be achieved by generating these phenotypes on a representative pool of cows for evaluation purpose as could be the case for a reference population for genomic selection.

However, a relevant question is the expected response to selection when including these traits, detailed proteins or FAA, in the official selection index and among the “true” breeding objectives. The official selection index in Ireland already includes a variety of features such as milk yield and protein concentration. One approach to estimate such a benefit in response is to derive the coefficient of genetic variation of the traits after adjustment for the genetic merit of (1) protein concentration in the case of protein fractions (2) or milk yield in the case of FAA.

In the present study, the near unity correlation between protein concentration and total CN (0.99) and the resulting small coefficient of genetic variation for total CN after adjustment for the genetic merit of protein concentration (0.75%) suggest that directly including CN in the selection index may result in little additional benefit. Similarly, the genetic correlations of protein concentration with αs1‐CN, β‐CN and κ‐CN were strong (0.99, 0.94 and 0.97, respectively) and the coefficient of genetic variation for αs1‐CN, β‐CN and κ‐CN after adjustment was just 4.30%, 3.16% and 4.14%, respectively. Hence, a very high selection pressure would be required for αs1‐CN, β‐CN and κ‐CN, that means reducing both the selection pressure on the other traits present in the index and their genetic gain. The correlation between protein concentration and α‐LA was only 0.55 and the corresponding coefficient of genetic variation after adjustment for the genetic merit was 2.70%. Therefore, even if genetic gain would be slow, it could be advantageous to include α‐LA per se in the selection index. High concentrations of α‐LA are indeed desirable for infant formula and dairy powders manufacturing (Lien, 2003). In some circumstances, the addition of valuable proteins like α‐LA would make sense in dairy cattle selection index, for example in those countries with a wide dairy product portfolio.

The weak genetic correlations between FAA and 24‐h milk yield as well as between FAA and protein concentration signify that existing selection indexes that consider milk yield and protein concentration are not fully exploiting the potential to genetically improve the desirable protein fractions and FAA. The genetic correlation close to 1 between total FAA and Glu (0.99) suggests that one should expect little benefit when including both of them in the selection index. The lack of very strong genetic correlations among the other FAA (i.e. Gly, Lys, Arg, Asp, Ser and Val) plus the existence of a moderate to large coefficient of genetic variation for these traits (Table 1) suggests that direct selection can be pursued simultaneously on all these traits together with traditional traits, like milk yield and protein concentration.

As proposed by Henchion et al. (2016), a potential approach would be to define the emphasis to be placed to a milk quality sub‐index, for example to the milk α‐LA in this study, that finds the consensus of different stakeholders (i.e. breed associations and companies, milk producers and processors, farm advisors and researchers). In Henchion et al. (2016), all the actors of the dairy chain were questioned to gauge their opinion about detailed milk quality traits within an overall main breeding goal. Based on their findings, the relative importance on a milk quality sub‐index was suggested to range from 6% to 10%. Nevertheless, the biggest concern raised by the interviewed stakeholders was their limited knowledge on the fine quality characteristics, followed by the scarce consumer demand for enhanced‐quality milk and dairy foodstuffs. In addition, some categories of consumers are not familiar with nutritional concepts and health benefits associated with consumption of animal products, including milk and dairy. Opinion and perception of different actors must be taken into account to properly quantify the selection index economic weights going in favour of desirable milk traits.

However, given a scenario where certain specific milk traits (e.g. protein fractions) are considered in the payment system, sub‐indexes can be useful to farmers and facilitate more bespoke breeding decisions based on the payment system. A good example is the Trentingrana cheese Consortium (Italy) that incentivizes farmers to produce milk with more favourable cheese‐making characteristics by applying specific rewards (Penasa et al., 2016).

Practical implications

Routine access to vast quantities of phenotypic data for protein fractions and FAA is imperative to achieve a high accuracy of selection and thus genetic gain. In addition, the actual response achievable is a function of how closely the selection index trait (i.e. MIR‐predicted phenotype) reflects the goal trait (i.e the true measure of the trait). This is a function of the genetic correlation between the MIR‐predicted value and its true value. In cattle, the genetic correlation between the reference trait and its MIRS prediction is often less than unity. Therefore, while selection using proxies may still achieve progress in a given direction, the actual response is likely to be less than optimum relative to if the gold standard was actually measured (assuming no difference in heritability).

Using the deterministic equation provided by Goddard (2009) to evaluate the accuracy of genomic evaluations with different sized reference populations, the number of genotyped individuals in the reference population necessary to achieve an animal EBV accuracy of 0.70 is approximately equal to 6,000 for predicted αs1‐CN. On the contrary, a larger reference population size is needed for predicted α‐LA (≈ 14,000), Glu (≈9,000) and Gly (≈33,000) to achieve the same accuracy of genomic evaluations. For this calculation, the effective population size was set at 75 as reported by McParland et al. (2007) for Irish Holsteins. Mid‐infrared spectroscopy is an efficient method commonly used by milk recording organizations worldwide to measure milk gross composiiton. Other than this, the ability of MIRS to predict individual protein fractions and FAA with a reasonable accuracy has been previously documented (McDermott et al., 2016). Thus, MIRS is the most popular, rapid and cost‐effective method for generating substantial quantities of phenotypes related to milk detailed composition, health traits and more (Benedet et al., 2020; Costa et al., 2021; Grelet et al., 2018).

Of potential interest to many processors are the similarities, or lack thereof, between human and bovine milk and the potential strategies to make them more similar, ideally at low cost and preferably even within the farm gate. Human and bovine milk differ in their protein profile; human milk has a whey to casein ratio of approximately 60:40, while bovine milk has a whey to casein ratio of approximately 20:80. While human milk does not contain β‐LG, it is present in large amount in bovine milk; on the contrary, α‐LA is the dominant protein in human milk but is relatively low in concentration in bovine milk (Lien, 2003). Therefore, it may be advantageous for infant formula producers to select bovine milk with a higher concentration of α‐LA and a lower concentration of β‐LG. The amino acid profile in human and bovine milk also differs (Chuang et al., 2005): human milk contains more Glu, Lys, Arg, Asp, Ser and Val than bovine milk (Ghadimi & Pecora, 1963).

Processors also aim to maximize the efficiency of transformation of the milk they purchase into saleable high‐quality and nutritious products. A high concentration of CN to total protein is essential for the transfer of proteins from milk to cheese with high concentrations of αs1‐CN, β‐CN, κ‐CN and β‐LG B being associated with greater cheese yield (Wedholm et al., 2006). The heritability estimates of the MIRS‐predicted CN fractions ranged from 0.36 (κ‐CN and αs2‐CN) to 0.45 (β‐CN) and the genetic standard deviation calculated for gold standard CN fractions in the present study (0.79 g/L for αs1‐CN) were comparable to results obtained by Graml and Pirchner (2003), i.e. 0.94 g/L for αs1‐CN). Variability in predicted CN fractions also existed in the present study, as the genetic standard deviation ranged from 0.18 g/L for αs2‐CN to 0.79 g/L for αs1‐CN. This was true eve after adjusting for differences in protein concentration, indicating an opportunity to improve milk‐specific protein fractions through selective breeding.

CONCLUSIONS

Gold standard and MIRS‐predicted traits were heritable and exhibited considerable exploitable genetic variation in their own right. Nonetheless, variability in some cases was eroded once expressed as relative to total protein concentration, which is present in most dairy cattle breeding goal worldwide. Traditional and genomic evaluations of detailed protein composition make sense if MIRS prediction equations will be implemented for phenotyping and whenever dairy chain stakeholders will conclude that certain protein fractions deserve inclusion in the breeding objective.

CONFLICT OF INTEREST

The authors declare no conflict of interest.