Protein profile of cow milk from multibreed herds and its relationship with milk coagulation properties

Abstract The aims of this study were to estimate the effect of cow breed on milk protein profile, and to assess the relationship between protein profile and milk coagulation properties (MCP) predicted with mid-infrared spectroscopy (MIRS). The dataset included 173,841 test-day observations of 12,533 specialised (Holstein-Friesian [HF] and Brown Swiss [BS]) and dual-purpose cows (Simmental [SI] and Alpine Grey) sampled in 488 multi-breed herds. Fixed effects considered in the mixed model were breed, month of sampling, parity class, stage of lactation, and first-order interactions and random effects were cow nested within breed, herd and the residual. BS had the greatest casein (CN) fractions content and SI the greatest whey fractions content. A lesser clear pattern was observed when considering the proportion of protein fractions. Longer rennet coagulation, longer curd-firming time and firmer curd were observed for Alpine Grey, HF and BS, respectively. Regarding the index of milk aptitude to coagulate (IAC), BS and SI had the greatest value and HF and Alpine Grey (AG) the lowest one. Weak relationships were estimated between protein fractions and MCP. For whey proteins, all the analysed traits were optimised at β-lactoglobulin (β-LG) concentrations next to the maximum of the physiological range and at α-lactalbumin (α-LA) concentrations at the minimum level. For CNs, all the analysed traits were optimised at α-CN and β-CN concentrations closer to the minimum level and at κ-CN concentration at intermediate levels. Results highlighted the lower efficiency of HF milk for cheese manufacturing compared to the other breeds. Highlights Breed, month of sampling, parity order, stage of lactation and their interactions were important to explain the variability of milk protein fractions and coagulation traits. Brown Swiss (BS) had the greatest casein (CN) fractions content and Simmental (SI) the greatest whey fractions content. Longer rennet coagulation was observed for Alpine Grey, whereas firmer curd was observed for BS. Coagulation properties were optimised in milk with low concentration of α-CN and β-CN and for intermediate levels of κ-casein.


Introduction
Milk proteins are important for human nutrition, both for their high biological value and their role as precursors of bioactive proteins and peptides. Milk proteins are known for their antimicrobial and immunomodulatory role and for facilitating absorption of nutrients (Haug et al. 2007;Lisson et al. 2013). Moreover, milk protein profile influences milk coagulation properties (MCP). Several authors have reported that an increase of some casein (CN) fractions content reduces rennet individuals are milk protein polymorphisms (Ketto et al. 2017). The effect of breed on milk protein fraction content has been reported in Estonian Red and Estonian Holstein cattle, with the former having greater b-LG, a S2 -CN and j-CN content than the latter (Jõudu et al. 2008). Nevertheless, the reference method to determine milk protein fractions (high-performance liquid chromatography), is expensive and time consuming, and thus it does not allow large-scale phenotyping.
In the last years, mid-infrared spectroscopy (MIRS) has shown its potential as rapid and cost-effective tool to analyse milk quality traits at population level (De Marchi et al. 2014). Although the development of accurate prediction models for protein fractions is quite challenging (Rutten et al. 2011;De Marchi et al. 2014), they can be enhanced through a novel chemometric approach using uninformative variable elimination applied to milk mid-infrared spectral wavelengths, as proposed by . The application of those improved models to historical test-day spectra for a posteriori prediction in multi-breed herds opens new opportunities to evaluate breed differences in terms of milk protein fractions, which is of interest for the dairy industry to increase processing efficiency and develop new functional products.
In the multifaceted scenario of cattle farming, multibreed dairy herds take on a particular interest, especially thanks to their potential performance in terms of milk production, animal reproduction and concentrate-conversion efficiency compared to single breed herds (Magne et al. 2016). Moreover, a complementarity occurs between lactation performances of breeds with diverse dairy aptitude (Magne et al. 2016) and such differences may also play a key role in the view of a sustainable livestock systems (Dumont et al. 2014).
Therefore, the aim of this study was to estimate sources of variation of milk protein profile of Holstein-Friesian (HF), Brown Swiss (BS), Simmental (SI) and Alpine Grey (AG) cows farmed in multi-breed herds, and to assess its relationship with MCP using data predicted with MIRS.

Data origin
Information on individual cow milk samples from routine cow milk testing was provided by the Breeders Association of Bolzano province (Bolzano, Italy) and the South Tyrolean Dairy Association (Bolzano, Italy).
Due to the orography of this mountainous province of the Italian Alps, farms in Bolzano area are characterised by small herd size (15 cows per herd on average) and cows fed with forage or hay and concentrates, and part of individuals are moved to highland pastures during summer season. The study involved two specialised dairy breeds, HF and BS and two dual-purpose breeds, SI and AG, farmed in multi-breed herds from 2011 to 2019. Data included milk yield (MY, kg/ d), somatic cell count (SCC, cells/mL), fat, crude protein (CP), CN and lactose percentage, milk urea content (MU, mg/dL) pH and spectral information. Each cow was evaluated every 4 or 5 weeks, and individual milk samples (50 mL) added with 200 mL of preservative (Bronysolv; ANA.LI.TIK Austria, Vienna, Austria) were collected and analysed for milk quality traits in the laboratory of the South Tyrolean Dairy Association following the guidelines of the International Committee for Animal Recording. The SCC was measured with Cell Fossomatic (FOSS Electric A/S, Hillerød, Denmark) and transformed to somatic cell score (SCS) as SCS ¼ log 2 (SCC/100) þ 3 (Wiggans and Shook 1987) to normalise the distribution of the data. Until February 2017, milk fat, CP, CN and lactose percentage and MU content were analysed with MilkoScan FT6000 (Foss, Hillerød, Denmark), and from March 2017 onwards with MilkoScan FT7 (Foss, Hillerød, Denmark). According to the manufacturer instructions, MilkoScan instruments were standardised to ensure the comparability of the collected spectra (Feudale et al. 2002;Juhl 2017).
Prediction models (g/100 mL) for a-CN, i.e. the sum of a S1 -and a S2 -CN, b-CN, j-CN, b-LG and a-LA developed by  were used. Briefly, prediction models were built using partial least squares regression analysis after an uninformative variable elimination procedure and validated through leave-one-out cross-validation, matching the spectra with reference chromatographic analysis performed on 114 individual milk samples which were collected in the same region of this study. Coefficient of determination in cross-validation (root mean square error in cross-validation) were 0.88 (1.05 mg/ mL) for a-CN, 0.60 (0.53 mg/mL) for b-CN, 0.74 (0.88 mg/mL) for j-CN, 0.47 (1.10 mg/mL) for b-LG and 0.37 (0.10 mg/mL) for a-LA. In this study, predicted protein fractions were expressed also as percentage of CP (%).
MCP, including RCT (min), curd-firming time (k 20 , min) and curd firmness 30 min after rennet addition (a 30 , mm) were predicted using models developed by Visentin et al. (2016). As for milk proteins, partial least squares regression analysis after an uninformative variable elimination procedure was used. The models were built in a calibration set (80% of the samples) and externally validated by applying the model to the remaining 20% of samples. The reference analysis to determine MCP was carried out using the Formagraph (Foss, Hillerød, Denmark). Coefficients of determination in cross-validation (root mean square error in crossvalidation) were 0.55 (2.86 min), 0.59 (1.00 min) and 0.56 (8.43 mm) for RCT, k 20 and a 30 , respectively . Moreover, the index of milk aptitude to coagulate (IAC) was estimated according to Penasa et al. (2015).

Data editing
Data were edited following the approach of Manuelian, Penasa, Visentin, Zidi, et al. (2018) for mineral composition of cow milk from multi-breed herds. The original dataset (n ¼ 2,760,766) was edited to retain records of purebred BS, HF, SI and AG cows of parity 1-12, between 5 and 305 days in milk (DIM), from cows that did not change herd during their productive life-cycle. Age at calving was checked to ensure consistency within parity order of the animal and cows whose age at calving deviated more than three standard deviations from the respective breed parity mean were discarded from the dataset. Records for the other investigated traits that deviated more than three standard deviations from the respective mean within each breed were treated as missing values. Contemporary groups were defined as cows tested in the same herd and date (herd-test-date, HTD) and HTD with less than 3 animals and cows with less than three observations within lactation were removed. Only records from herds with 2, 3 or 4 breeds and breed combinations present in more than 2 herds were kept, with each breed representing at least 10% of the herd. After editing, 173,841 observations of 12,533 cows in 489 multi-breed herds were available for further analyses. Parity, DIM and herd size averaged 2.80 ± 1.73, 152.16 ± 82.81 d and 25.68 ± 16.74 cows (range: 6-157), respectively. Detailed number of herds and cows involved in the study and relative number of records included in the dataset is reported in Table 1.

Statistical analysis
MY, composition, SCS, MU, protein fractions, MCP and IAC were analysed through a linear mixed model using HPMIXED procedure of SAS version 9.4 (SAS Institute Inc., Cary, NC). The model was: where y ijklmn is the dependent variable; m is the overall intercept of the model; B i is the fixed effect of the ith breed (i ¼ HF, BS, AG, SI); M j is the fixed effect of the jth month of sampling (j ¼ 1-12); S k is the fixed effect of the kth class of stage of lactation (k ¼ 1-30, each class being 10 DIM wide); P l is the fixed effect of the lth parity class (l ¼ 1-5, with the last including parities ! 5); (B Â M) ij is the fixed interaction effect between breed and month of sampling; (B Â S) ik is the fixed interaction effect between breed and stage of lactation; (B Â P) il is the fixed interaction effect between breed and parity; (S Â P) kl is the fixed interaction effect between class of stage of lactation and parity class; H m is the random effect of the mth herd (m ¼ 1-488); C n (B i ) is the random effect of the nth cow nested within the ith breed (n ¼ 1-12,533)$N(0,r 2 C(B) ), where r 2 C(B) is the cow nested within breed variance; and e ijklmn is the random residual $N(0,r 2 e ), where r 2 e is the residual variance. Multiple comparisons of least square means were performed for the main effects of breed, parity class and class of stage of lactation using Tukey-Kramer adjustment. Pearson correlations between the investigated traits and studentised residuals were assessed using the CORR procedure of SAS. Significance was set at p < .01, unless otherwise stated. In order to investigate the effect of protein fractions on MCP, a second linear mixed model was run including protein fractions, i.e. a-CN, b-CN, j-CN, a-LA and b-LG expressed as % of CP, as first, second and third order covariates. Obtained regression coefficients for protein fractions were used to perform a Lagrange optimisation using MCP and IAC as dependant variables, regression coefficients from previous linear mixed model, and protein fractions relative content as independent variables. Optimisation was constrained for single protein content to be in the physiological range and the sum of all protein fractions to be equal to 1 (Madoumier et al. 2019).

Results and discussion
Descriptive statistics and significance of fixed effects Descriptive statistics for MY, chemical composition, SCS, MU, detailed protein composition and MCP are summarised in Table 2. Protein fractions predicted by MIRS were characterised by low to moderate coefficient of variation (CV), similarly to the variation reported by  in single-breed herds located in the same region. Means of RCT, k 20 and a 30 were 21.64, 5.71 min and 16.75 mm, respectively, and the latter exhibited the largest variability among MCP (CV ¼ 59.48%). Average values of RCT and k 20 were greater than those reported by Visentin et al. (2016) from single-breed herds sampled in the same region and measured using the Formagraph. As a consequence, a 30 was higher in the aforementioned study compared to present result which indeed is focussed on multi-breed herds. On the contrary, IAC showed the lowest CV (4.48%) among all the studied traits.
The statistical analysis revealed that fixed effects considered in the model were all significant in explaining the variation of MY, milk composition, pH, SCS, MU, protein fractions, MCP and IAC (Table 3). However, this was somewhat expected due to the large number of observations included in the dataset and confirmed previous results obtained in singlebreed herds of the same area (Franzoi et al. 2019). Focussing on the proportion of CN fractions expressed as % CP, month of sampling showed the highest Fvalue for a-CN, while parity had the highest influence on b-CN and j-CN. Among whey proteins, b-LG was mostly influenced by cow breed and month of sampling, while a-LA was mostly influenced by stage of lactation. Considering protein fractions expressed in g/ 100 mL of milk, a large proportion of a-CN and b-CN variability was explained by stage of lactation, similarly to CP. On the other hand, parity had the strongest effect on j-CN. These results suggest different determinants in the variability among CNs. Finally, month of sampling had the strongest effect on whey protein fractions. In Burlina cattle (Niero, Visentin, et al. 2016), stage of lactation explained a larger proportion of b-CN and j-CN variability than parity; instead, variability of a-CN was explained in a higher proportion by parity effect compared to the stage of lactation. Taking into account MCP, the most significant effect for RCT and IAC was stage of lactation, while breed was the most significant effect for a 30 and k 20 . These results partially agreed with those reported in Rendena and HF cows (Varotto et al. 2015), where breed was the most significant effect for RCT and a 30 . In general, for most of the analysed traits, the random effect of cow explained greater proportion of variability of the investigated traits compared to the random effect of herd, with values that ranged from 14.80% for MU to 44.59% for b-CN. The only exceptions were MY and MU, for which herd effect explained 34.79 and 22.63% of total variance, respectively.

Breed effect
In line with Penasa et al. (2014), HF breed was endowed with the greatest MY (29.37 kg/d), followed Table 2. Number of test day records, mean, standard deviation (SD), range and coefficient of variation (CV) for milk yield, composition, somatic cell score, urea, detailed protein composition and coagulation properties.  (Table 4). BS had greater fat, CP and CN percentages and MU content than other breeds. The lowest SCS was estimated for SI cows. Overall, protein fractions varied significantly across breeds. BS registered the greatest values for all CN groups. Aside from being the highest yielding animals, HF cows had the lowest content of CN fractions, probably due to a Table 3. F-values of fixed effects of breed (B), month of sampling (M), parity order (P), stage of lactation (S) and their first-order interactions and percentage of total variance explained by random effects of herd, cow nested within breed and residual for milk yield, composition, somatic cell score, urea, detailed protein composition and coagulation properties.  dilution effect as reported by Franzoi et al. (2019). Such a lower CP and CN content in HF milk is also in agreement with the lower amount of Ca, Mg and P observed by Manuelian, Penasa, Visentin, Zidi, et al. (2018), due to the positive correlation between minerals and CN micelles (0.23-0.53; p < .05; Toffanin et al. 2015). As regard protein fractions expressed as proportion of CP, SI breed had the greatest content of a-CN, whereas milk of BS and AG breed were the richest in b-CN and j-CN, respectively (p<.05). Milk of HF cows had the greatest content of b-LG and a-LA (9.74 and 2.33% of CP, respectively), whereas BS had the lowest one (7.91 and 2.07% of CP, respectively; p<.05). Breed differences in milk protein profile have been reported also by Gustavsson et al. (2014), who considered Swedish Red, Danish Holstein and Danish Jersey cows, and by Franzoi et al. (2019), who considered HF, BS and SI cows farmed in single-breed herds.
To the authors' knowledge, no studies have attempted to estimate breed differences for detailed protein composition predicted by MIRS in multibreed herds, especially on a large scale. Previous literature focussed on single breeds and on a relatively low number of animals, which were characterised for detailed protein profile through capillary zone electrophoresis technique (Schopen et al. 2009) or high pressure liquid chromatography technique (Cipolat-Gotet et al. 2018). Bonfatti et al. (2017) investigated milk protein composition using phenotypes predicted from a large spectra database of Italian SI cows, and Franzoi et al. (2019) studied the phenotypic variation of predicted protein profile of BS, HF, SI, AG and Pinzgauer cows using data from single-breed herds.
Our results indicated that the four breeds were rather similar in terms of predicted RCT which was the shortest for SI cows (21.32 min) and the longest for AG cows (21.82 min), predicted k 20 which ranged from 5.28 min (BS) to 6.49 min (HF) and predicted IAC which was slightly higher in SI and BS than HF and AG. Compared to the other breeds, BS and HF cows showed the greatest (17.76 mm) and the lowest a 30 (14.15 mm). Moreover, the lower MY and more favourable MCP of BS than in HF is in agreement with findings of Penasa et al. (2014). However, it is worthwhile to note that prediction models used in this study were characterised by moderate accuracy, as evidenced by the coefficient of determination reported by  for detailed milk protein composition, and by Visentin et al. (2016) for MCP. Figure 1 depicts the least squares means of predicted protein composition across stage of lactation for BS, HF, AG and SI. Protein fractions expressed as g/100 mL of milk had the same trend for all the studied breeds. In particular, a-CN, j-CN and b-LG had greater concentration at the beginning of lactation, a swift decline close to the peak of lactation, and a constant increase towards late lactation stages. Such a variation is in agreement with previous authors who reported that the trend of milk protein fractions over lactation is opposite to the trend of MY (McDermott et al. 2017;Franzoi et al. 2019), again confirming the hypothesis of a dilution effect. This was observed also in the case of milk constituents other than protein fractions, with particular regard to major milk minerals (Manuelian, Penasa, Visentin, Zidi, et al. 2018;Visentin et al. 2018). The specific pattern (Figure 1) of milk protein profile along with the variation of major milk minerals across lactation ) may explain the observed trend of MCP, which corroborates with Penasa et al. (2014), i.e. shorter RCT and greater a 30 at the onset of milk production. The trend of a-CN, j-CN and b-LG was only partially mimicked by b-CN, which showed similar concentration at the beginning and at the peak of lactation, and by a-LA, which maintained constant concentration after the peak MY.

Lactation and parity effects
Protein fractions expressed as percentage of CP showed different patterns, with a-CN declining constantly across lactation, while b-CN had a sharp increase up to mid-lactation to subsequently plateau. The j-CN and a-LA fractions peaked between the first and second month of lactation and declined thereafter. Finally, b-LG had a nadir in the early lactation stages and slightly recovered towards late stages. The specific composition in terms of milk protein fractions across lactation likely reflects calf's nutritional needs and cow's physiological requirements. On one hand, a greater proportion of b-LG in early lactation could be linked with the biological function of this protein fraction in calves, due to its capacity to enhance the absorption of hydrophobic ligands such as retinol and fatty acids (Zhang et al. 2015). On the other hand, decreased proteins related to milk component synthesis, such as a-LA in late lactation can be associated with the protection of the mammary gland during the involution process (Zhang et al. 2015).
Trends of predicted protein fractions across parity for BS, HF, AG and SI are depicted in Figure 2. Among CN fractions, both j-CN expressed as g/100 mL of milk and j-CN measured as proportion of CP showed decreasing trend in later parities. An opposite trend was observed for b-CN fraction. As for whey protein fractions, b-LG expressed as g/100 mL increased in older cows, and this trend was confirmed also when considering the same protein measured on CP basis.

Associations between protein fractions and MCP
Pearson correlation coefficients estimated from raw data are presented in Table 5 (above the diagonal). Within protein fractions, correlations were moderate to low. The strongest associations were assessed between j-CN and b-LG (À0.56) and between b-LG and a-LA (0.48), whereas almost null correlation was estimated between b-CN and j-CN. These results are not fully in agreement with McDermott et al. (2016) who reported positive correlations between milk protein fractions in different cow breeds (HF, Jersey, Norwegian Red and crossbreds) that ranged from 0.32 (b-LG-B and a s1 -CN) to 0.92 (b-LG-A and j-CN).
Correlations of protein fractions with MCP and IAC were weak, from À0.24 (j-CN and IAC) to 0.23 (b-LG and IAC). Close to null or null correlations were observed between a-CN and a 30 , a-CN and IAC, a-LA and a 30 and a-LA and IAC. Positive and low correlations were estimated between RCT and CN fractions content; this was somehow expected, since greater CN content is likely to increase RCT due to the higher ratio of substrate available for rennet proteolysis. Overall, such rather weak relationships might be partly argued considering that models used for the prediction of MCP are based on 30 min of lactodynamographic assays, which were taken as reference    Jõudu et al. (2008). In addition to those protein fractions, total amount of j-CN also affects a 30 (Jõudu et al. 2008). Mild to strong correlations were observed between MCP, and between MCP and IAC, with positive associations ranging from 0.39 (RCT and k 20 ) to 0.94 (a 30 and IAC), and negative associations ranging from À0.54 (k 20 and a 30 ) to À0.93 (RCT and IAC). Those correlations are in agreement with the studies conducted in HF, BS, SI, AG  and Pinzgauer  single-breed herds in Bolzano province.
Pearson correlation coefficients calculated from the residuals (Table 5, below the diagonal) estimated by the model were generally stronger than those estimated from raw data. Within protein fractions, the strongest associations were estimated between j-CN and b-LG (À0.61), in a way that greater content of CN fraction hampered that of whey protein fraction, and between b-LG and a-LA (0.28), in a way that greater content of a fraction mutually enhanced the other. Correlations of the residuals between protein fractions and MCP indicate that CN fractions were positively associated with RCT, whereas whey fractions were negatively associated, despite a close to zero correlation between a-LA and RCT (À0.03). The b-CN and j-CN fractions were negatively associated to a 30 (À0.22 and À0.16, respectively) and IAC (À0.20 and À0.25, respectively). Within milk coagulation traits, a strong and positive association between a 30 and IAC (0.96) and a strong negative correlation between RCT and IAC (À0.90) were estimated.

Covariances between protein fractions and MCP
In order to further investigate the effect of protein fractions on MCP and IAC, a linear mixed model was fitted with protein fractions as covariates up to third order. Obtained coefficients were used to perform a constrained Lagrange optimisation, with the aim to identify protein composition optimising MCP and IAC ( Table 6). The proportion of a-CN optimising MCP and IAC was similar for all the analysed traits ranging from 35.96% for k 20 to 39.97% for IAC. On the other hand, remarkable differences were found for b-CN concentrations, with lower values optimising k 20 (16.20%) and higher values optimising RCT (24.30%). On the opposite, k 20 was optimised at higher j-CN concentrations (27.67%), while RCT was minimised at lower j-CN concentrations (16.68%). In general, a 30 was mostly maximised at intermediate values compared to concentration optimising RCT and k 20 . This was expected considering that a 30 is both linked to RCT and k 20 .
Taking into account whey proteins, all the analysed traits were optimised at b-LG concentrations next to the maximum of the physiological range. This could be surprising, considering b-LG is a whey protein, usually considered not playing an important role in milk technological characteristics. Nevertheless, it was previously demonstrated that spiking raw milk with 0.2% b-LG had favourable effect on RCT. At the same time, the effect of increasing b-LG concentration was detrimental after heating milk (Kannan and Jennes 1961). This is in accordance with findings reporting b-LG interacts with j-CN, in particular after heat treatments (Imafidon et al. 1991). Similarly, Abeykoon et al. (2016) reported a favourable significant correlation between b-LG relative content and RCT in individual milk samples from different breeds. On the contrary, a-LA optimised MCP and IAC at its lower physiological limit (1.39%). This evidence could be linked to chelating effect of a-LA on divalent cations (mainly Ca and Zn) involved in paracasein reticulum formation (L€ onnerdal and Lien 2003).