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.

Keywords:

• Cattle breed
• mid-infrared spectroscopy
• phenomics
• milk coagulation properties
• proteins

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 coagulation time (RCT), and increase curd firmness (Vallas et al. 2012; Poulsen et al. 2013) and acid coagulation properties (Ketto et al. 2017).

About 80% of milk proteins are CNs, classified as α_S1-CN, α_S2-CN, β-CN and κ-CN, and the remaining 20% are whey proteins, mainly β-lactoglobulin (β-LG) and α-lactalbumin (α-LA). However, milk protein composition differs among individuals (Ketto et al. 2017; Franzoi et al. 2019) and breeds (Haug et al. 2007), and one of the factors associated with differences among 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 β-LG, α_S2-CN and κ-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 Niero, Penasa, et al. (2016). 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.

Materials and methods

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/µL), 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 µL 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₂(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 α-CN, i.e. the sum of α_S1- and α_S2-CN, β-CN, κ-CN, β-LG and α-LA developed by Niero, Penasa, et al. (2016) 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 α-CN, 0.60 (0.53 mg/mL) for β-CN, 0.74 (0.88 mg/mL) for κ-CN, 0.47 (1.10 mg/mL) for β-LG and 0.37 (0.10 mg/mL) for α-LA. In this study, predicted protein fractions were expressed also as percentage of CP (%).

MCP, including RCT (min), curd-firming time (k₂₀, min) and curd firmness 30 min after rennet addition (a₃₀, 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 cross-validation) were 0.55 (2.86 min), 0.59 (1.00 min) and 0.56 (8.43 mm) for RCT, k₂₀ and a₃₀, respectively (Visentin et al. 2016). 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.

Table 1. Number (frequency between brackets) of herds and cows involved in the study, and number of records (frequency between brackets) within herd and breeds.

Item^a	Number (frequency, %)	Number of records (frequency, %)
Herds
BS + HF	148 (30.27)	57,358 (32.99)
HF + SI	113 (23.11)	40,388 (23.23)
BS + SI	84 (17.18)	23,247 (13.37)
BS + HF + SI	68 (13.91)	28,794 (16.56)
SI + AG	28 (5.73)	10,513 (6.05)
HF + AG	22 (4.50)	5470 (3.15)
HF + SI + AG	11 (2.25)	3968 (2.28)
BS + AG	8 (1.64)	2501 (1.44)
BS + HF + AG + SI	4 (0.82)	711 (0.41)
BS + HF + AG	3 (0.61)	883 (0.51)
Overall	489 (100.00)	173,841 (100.00)
Cows
BS	4483 (35.77%)	64,429 (37.14%)
HF	4013 (32.02%)	52,127 (30.05%)
SI	3319 (26.48%)	46,387 (26.74%)
AG	718 (5.73%)	10,898 (6.28%)
Overall	12,533 (100.00%)	173,481 (100.00%)

^aBS: Brown Swiss; HF: Holstein-Friesian; SI: Simmental; AG: Alpine Grey.

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: $$y_{\mathit{ijklmn}}=\mathrm{\mu }+{\mathrm{B}}_i+{\mathrm{M}}_j+{\mathrm{S}}_k+{\mathrm{P}}_l+{\left(\text{B × M}\right)}_{ij}+{\left(\text{B × S}\right)}_{ik}+{\left(\text{B × P}\right)}_{il}+{\left(\text{S × P}\right)}_{kl}+{\mathrm{H}}_m+{\mathrm{C}}_n\left({\mathrm{B}}_i\right)+{\mathrm{e}}_{\mathit{ijklmn}},$$ where y_ijklmn is the dependent variable; µ 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,σ²_C(B)), where σ²_C(B) is the cow nested within breed variance; and e_ijklmn is the random residual ∼N(0,σ²_e), where σ²_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. α-CN, β-CN, κ-CN, α-LA and β-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 Niero, Penasa, et al. (2016) in single-breed herds located in the same region. Means of RCT, k₂₀ and a₃₀ 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₂₀ 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₃₀ 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.

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.

Trait	No of records	Mean	SD	Range	CV, %
Milk yield, kg/d	173,841	26.51	7.45	51.20	28.10
Milk composition, %
Fat	171,742	4.14	0.68	4.68	16.46
Crude protein	173,484	3.50	0.38	2.48	10.94
Casein	173,397	2.75	0.31	2.01	11.24
Somatic cell score, units	173,841	2.40	1.87	9.59	78.06
Milk urea, mg/dL	173,841	21.36	7.07	44.9	33.09
Protein fractions, g/100 mL
α-Casein	173,678	1.54	0.19	2.40	12.50
β-Casein	169,187	1.04	0.18	1.57	16.79
κ-Casein	169,188	0.58	0.13	1.03	21.78
β-Lactoglobulin	170,557	0.30	0.11	1.07	36.61
α-Lactalbumin	173,803	0.08	0.01	0.11	8.01
Protein fractions, % of crude protein
α-Casein	172,911	43.91	2.81	19.20	6.39
β-Casein	168,652	29.84	3.93	26.66	13.19
κ-Casein	168,707	16.57	3.35	20.64	20.18
β-Lactoglobulin	169,977	8.52	2.97	17.40	34.90
α-Lactalbumin	173,187	2.18	0.26	1.77	12.01
Milk coagulation properties^a
RCT, min	159,206	21.64	3.86	19.62	17.83
k20, min	168,563	5.71	1.29	8.05	22.55
a30, mm	127,123	16.75	9.97	48.96	59.48
IAC	125,814	100.58	4.50	24.41	4.48

^aRCT: rennet coagulation time; k₂₀: curd-firming time; a₃₀: curd firmness 30 min after rennet addition. IAC: index of milk aptitude to coagulate

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 single-breed 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 F-value for α-CN, while parity had the highest influence on β-CN and κ-CN. Among whey proteins, β-LG was mostly influenced by cow breed and month of sampling, while α-LA was mostly influenced by stage of lactation. Considering protein fractions expressed in g/100 mL of milk, a large proportion of α-CN and β-CN variability was explained by stage of lactation, similarly to CP. On the other hand, parity had the strongest effect on κ-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 β-CN and κ-CN variability than parity; instead, variability of α-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₃₀ and k₂₀. 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₃₀. 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 β-CN. The only exceptions were MY and MU, for which herd effect explained 34.79 and 22.63% of total variance, respectively.

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.

	Fixed effects								Random effects
Trait	B	M	P	S	B × M	B × S	B × P	S × P	σ^2herd, %	σ^2cow (breed), %	σ^2residual, %
Milk yield, kg/d	1078.04	269.55	3456.84	3600.22	5.78	10.48	21.00	81.10	34.79	24.42	40.78
Milk composition, %
Fat	188.71	316.65	55.93	451.25	4.78	8.97	14.00	5.51	12.46	25.97	61.57
Crude protein	1107.83	958.20	396.88	5147.91	8.51	32.98	9.72	10.70	13.37	43.07	43.56
Casein	1079.55	962.77	465.25	4822.38	7.09	31.41	14.06	11.87	13.48	44.22	42.30
Somatic cell score, units	38.19	7.64	433.46	254.79	2.34	3.37	10.36	3.96	6.68	19.72	73.60
Milk urea, mg/dL	349.60	41.97	76.93	105.65	10.58	4.07	2.56	3.88	22.63	14.80	62.56
Protein fractions, g/100 mL
α-Casein	817.82	1447.12	462.33	1784.73	9.62	15.19	15.31	4.83	15.51	35.24	49.25
β-Casein	372.77	171.97	719.31	2577.06	11.70	8.74	4.34	6.15	12.72	44.59	42.69
κ-Casein	636.01	316.92	2373.15	195.22	11.53	3.26	8.57	2.47	7.45	39.55	53.00
β-Lactoglobulin	155.13	796.64	193.68	544.63	4.76	12.35	5.04	4.82	5.45	35.99	58.56
α-Lactalbumin	193.71	268.27	52.53	104.87	2.76	9.53	15.23	3.29	11.19	33.30	55.50
Protein fractions, % of crude protein
α-Casein	122.73	778.81	382.37	580.35	5.91	8.47	11.20	2.59	12.17	22.43	65.40
β-Casein	140.34	125.45	960.77	639.77	11.49	4.60	3.09	3.99	10.96	43.73	45.31
κ-Casein	219.83	709.81	2156.74	108.32	11.73	1.65	5.79	2.53	8.43	37.12	54.46
β-Lactoglobulin	452.70	582.23	262.84	205.62	6.24	10.36	7.68	3.98	5.51	35.13	59.36
α-Lactalbumin	1327.26	275.62	241.77	3135.61	7.57	24.43	12.61	5.62	8.64	42.39	48.97
Milk coagulation properties^a
RCT, min	5.03	67.93	176.08	325.09	9.03	5.75	20.17	3.31	7.66	28.98	63.37
k20, min	138.61	94.36	86.54	138.1	2.17	3.11	14.47	2.72	5.78	24.59	69.62
a30, mm	979.99	168.45	564.51	717.95	4.73	9.25	6.35	2.45	8.10	38.94	52.95
IAC	15.79	59.52	156.01	164.13	2.87	4.11	16.42	3.09	5.31	28.51	66.18

^aRCT: rennet coagulation time; k₂₀: curd-firming time; a₃₀: curd firmness 30 min after rennet addition. IAC: index of milk aptitude to coagulate.

All fixed effects were highly significant (p<.001).

Breed effect

In line with Penasa et al. (2014), HF breed was endowed with the greatest MY (29.37 kg/d), followed by SI (25.67 kg/d), BS (25.16 kg/d) and AG (20.79 kg/d) breeds (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 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 α-CN, whereas milk of BS and AG breed were the richest in β-CN and κ-CN, respectively (p<.05). Milk of HF cows had the greatest content of β-LG and α-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.

Table 4. Least squares means (standard error) of milk yield, composition, somatic cell score, urea, detailed protein composition and coagulation properties.

Trait	Brown Swiss	Holstein-Friesian	Alpine Grey	Simmental
Milk yield, kg/d	25.16 (0.17)^c	29.37 (0.17)^a	20.79 (0.25)^d	25.67 (0.17)^b
Milk composition, %
Fat	4.27 (0.013)^a	4.04 (0.013)^c	3.86 (0.024)^d	4.12 (0.014)^b
Crude protein	3.62 (0.006)^a	3.30 (0.006)^c	3.53 (0.012)^b	3.51 (0.007)^b
Casein	2.84 (0.005)^a	2.58 (0.005)^c	2.76 (0.01)^b	2.75 (0.005)^b
Somatic cell score, units	2.50 (0.029)^ab	2.47 (0.029)^b	2.52 (0.057)^a	2.19 (0.031)^c
Milk urea, mg/dL	23.33 (0.168)^a	20.15 (0.168)^d	21.87 (0.252)^b	21.09 (0.172)^c
Protein fractions, g/100 mL
α-Casein	1.58 (0.004)^a	1.44 (0.004)^d	1.54 (0.007)^c	1.56 (0.004)^b
β-Casein	1.12 (0.003)^a	1.02 (0.003)^d	1.08 (0.007)^b	1.03 (0.004)^c
κ-Casein	0.60 (0.002)^a	0.47 (0.002)^d	0.58 (0.005)^b	0.54 (0.003)^c
β-Lactoglobulin	0.29 (0.002)^b	0.32 (0.002)^a	0.29 (0.004)^b	0.33 (0.002)^a
α-Lactalbumin	0.074 (0.0001)^c	0.076 (0.0001)^b	0.077 (0.0002)^a	0.077 (0.0001)^a
Protein fractions, % of crude protein
α-Casein	43.67 (0.051)^b	43.59 (0.051)^b	43.62 (0.093)^b	44.45 (0.054)^a
β-Casein	31.05 (0.081)^a	30.93 (0.081)^ab	30.62 (0.163)^b	29.33 (0.087)^c
κ-Casein	16.63 (0.066)^a	15.05 (0.067)^c	16.68 (0.135)^a	15.37 (0.071)^b
β-Lactoglobulin	7.91 (0.047)^c	9.74 (0.048)^a	8.08 (0.102)^c	9.23 (0.052)^b
α-Lactalbumin	2.07 (0.004)^b	2.33 (0.004)^c	2.21 (0.008)^a	2.22 (0.004)^a
Milk coagulation properties^A
RCT, min	21.50 (0.068)^b	21.52 (0.069)^b	21.82 (0.138)^a	21.32 (0.073)^c
k20, min	5.28 (0.022)^d	6.49 (0.022)^a	5.56 (0.046)^c	5.86 (0.024)^b
a30, mm	17.76 (0.162)^a	14.15 (0.17)^c	16.27 (0.337)^b	16.41 (0.179)^b
IAC	100.93 (0.074)^a	100.49 (0.078)^b	100.32 (0.157)^b	100.97 (0.082)^a

^ARCT: rennet coagulation time; k₂₀: curd-firming time; a₃₀: curd firmness 30 min after rennet addition. IAC: index of milk aptitude to coagulate.

Least squares mean with different superscript letters within a row are significantly different (p<.01).

To the authors’ knowledge, no studies have attempted to estimate breed differences for detailed protein composition predicted by MIRS in multi-breed 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₂₀ 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₃₀ (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 Niero, Penasa, et al. (2016) for detailed milk protein composition, and by Visentin et al. (2016) for MCP.

Lactation and parity effects

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, α-CN, κ-CN and β-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 (Visentin et al. 2018) may explain the observed trend of MCP, which corroborates with Penasa et al. (2014), i.e. shorter RCT and greater a₃₀ at the onset of milk production. The trend of α-CN, κ-CN and β-LG was only partially mimicked by β-CN, which showed similar concentration at the beginning and at the peak of lactation, and by α-LA, which maintained constant concentration after the peak MY.

Figure 1. Least squares mean of (A) α-casein (α-CN), (B) β-casein (β-CN), (C) κ-casein (κ-CN), (D) β-lactoglobulin (β-LG) and (E) α-lactalbumin (α-LA) across stage of lactation for Brown Swiss (

), Holstein-Friesian (

), Alpine Grey (

), and Simmental (

) cows, expressed as g/100 mL of milk (left side) or percentage of crude protein (right side) of the figure.

Protein fractions expressed as percentage of CP showed different patterns, with α-CN declining constantly across lactation, while β-CN had a sharp increase up to mid-lactation to subsequently plateau. The κ-CN and α-LA fractions peaked between the first and second month of lactation and declined thereafter. Finally, β-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 β-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 α-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 κ-CN expressed as g/100 mL of milk and κ-CN measured as proportion of CP showed decreasing trend in later parities. An opposite trend was observed for β-CN fraction. As for whey protein fractions, β-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.

Figure 2. Least squares mean of (A) α-casein (α-CN), (B) β-casein (β-CN), (C) κ-casein (κ-CN), (D) β-lactoglobulin (β-LG) and (E) α-lactalbumin (α-LA) across parity class for Brown Swiss (

), Holstein-Friesian (

), Alpine Grey (

), and Simmental (

) cows, expressed as g/100 mL of milk (left side) or percentage of crude protein (right side) of the figure.

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 κ-CN and β-LG (−0.56) and between β-LG and α-LA (0.48), whereas almost null correlation was estimated between β-CN and κ-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 (β-LG-B and α_s1-CN) to 0.92 (β-LG-A and κ-CN).

Table 5. Pearson correlation coefficients for protein fractions (measured as % of crude protein) and milk coagulation properties.

Trait	α-Casein	β-Casein	κ-Casein	β-Lactoglobulin	α-Lactalbumin	RCT	k20	a30	IAC
Protein fractions
α-Casein	–	−0.49	−0.22	0.19	0.33	0.06	0.13	0.01	0.01
β-Casein	−0.52	–	0.04	−0.32	−0.15	0.17	0.10	−0.18	−0.14
κ-Casein	−0.32	0.18	–	−0.56	−0.46	0.21	−0.12	−0.10	−0.24
β-Lactoglobulin	0.21	−0.33	−0.61	–	0.48	−0.22	−0.22	0.22	0.23
α-Lactalbumin	0.24	−0.05	−0.34	0.28	–	0.06	0.11	−0.01	0.00
Milk coagulation properties^a
RCT	0.05	0.19	0.27	−0.38	−0.03	–	0.39	−0.74	−0.93
k20	0.12	0.20	−0.05	−0.26	0.39	0.53	–	−0.54	−0.44
a30	0.07	−0.22	−0.16	0.30	−0.04	−0.77	−0.49	–	0.94
IAC	0.08	−0.20	−0.25	0.32	0.03	−0.90	−0.47	0.96	–

Above the diagonal, Pearson correlation coefficients calculated from raw data; under the diagonal, Pearson correlation coefficients calculated from studentised residuals estimated by the model.

^aRCT: rennet coagulation time (min); k₂₀: curd-firming time (min); a₃₀: curd firmness 30 min after rennet addition (mm).

IAC: index of milk aptitude to coagulate

All correlations differed from zero (p<.001).

Correlations of protein fractions with MCP and IAC were weak, from −0.24 (κ-CN and IAC) to 0.23 (β-LG and IAC). Close to null or null correlations were observed between α-CN and a₃₀, α-CN and IAC, α-LA and a₃₀ and α-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 measurements. In this view, curd of the late-coagulating samples may have had not enough time to develop, thus providing only a limited description of the potential gel firmness (Niero et al. 2018). Weak to moderate correlations between MCP and protein fractions were also reported by McDermott et al. (2016). Moreover, a significant effect of total α_s1-CN, β-CN and β-LG content on RCT milk from Estonian Red, Red-and-White Holstein and Estonian Holstein was observed by Jõudu et al. (2008). In addition to those protein fractions, total amount of κ-CN also affects a₃₀ (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₂₀) to 0.94 (a₃₀ and IAC), and negative associations ranging from −0.54 (k₂₀ and a₃₀) to −0.93 (RCT and IAC). Those correlations are in agreement with the studies conducted in HF, BS, SI, AG (Visentin et al. 2018) and Pinzgauer (Manuelian, Penasa, Visentin, Cassandro, et al. 2018) 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 κ-CN and β-LG (−0.61), in a way that greater content of CN fraction hampered that of whey protein fraction, and between β-LG and α-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 α-LA and RCT (−0.03). The β-CN and κ-CN fractions were negatively associated to a₃₀ (−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₃₀ 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 α-CN optimising MCP and IAC was similar for all the analysed traits ranging from 35.96% for k₂₀ to 39.97% for IAC. On the other hand, remarkable differences were found for β-CN concentrations, with lower values optimising k₂₀ (16.20%) and higher values optimising RCT (24.30%). On the opposite, k₂₀ was optimised at higher κ-CN concentrations (27.67%), while RCT was minimised at lower κ-CN concentrations (16.68%). In general, a₃₀ was mostly maximised at intermediate values compared to concentration optimising RCT and k₂₀. This was expected considering that a₃₀ is both linked to RCT and k₂₀.

Table 6. Results from Lagrange constrained optimisation of protein relative concentration for milk coagulation properties and index of milk aptitude to coagulate.

Protein fractions, % of crude protein	Optimisation range	Optimising concentrations
		RCT^a	k20	a30	IAC
α-Casein	34.67–53.87	38.85	35.96	38.00	39.97
β-Casein	16.20–42.86	24.30	16.20	20.65	22.38
κ-Casein	7.03–27.67	16.68	27.67	21.34	17.91
β-Lactoglobulin	1.38–18.78	18.78	18.78	18.62	18.35
α-Lactalbumin	1.39–3.16	1.39	1.39	1.39	1.39

^aRCT: rennet coagulation time (min); k₂₀: curd-firming time (min); a₃₀: curd firmness 30 min after rennet addition (mm).

IAC: index of milk aptitude to coagulate.

Taking into account whey proteins, all the analysed traits were optimised at β-LG concentrations next to the maximum of the physiological range. This could be surprising, considering β-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% β-LG had favourable effect on RCT. At the same time, the effect of increasing β-LG concentration was detrimental after heating milk (Kannan and Jennes 1961). This is in accordance with findings reporting β-LG interacts with κ-CN, in particular after heat treatments (Imafidon et al. 1991). Similarly, Abeykoon et al. (2016) reported a favourable significant correlation between β-LG relative content and RCT in individual milk samples from different breeds. On the contrary, α-LA optimised MCP and IAC at its lower physiological limit (1.39%). This evidence could be linked to chelating effect of α-LA on divalent cations (mainly Ca and Zn) involved in paracasein reticulum formation (Lönnerdal and Lien 2003).

Conclusions

Our results demonstrated that breed, month, parity, stage of lactation and first-order interactions were important contributors to milk protein fractions and coagulation traits. BS and SI had the greatest total amount of CN and whey fractions, respectively. However, a lesser clear pattern was observed when considering the distribution of the protein fractions within total CP. The greatest RCT, k₂₀ and a₃₀ was obtained for AG, HF and BS, respectively. Moreover, two groups of breeds were observed when looking at the IAC; BS and SI had the greatest value, whereas HF and AG the lowest one. Weak relationships were assessed between protein fractions and MCP. For whey proteins, all the analysed traits were optimised at β-LG concentrations next to the maximum of the physiological range and at α-LA concentrations at the minimum level. For CN proteins, all the analysed traits were optimised at α-CN and β-CN concentrations closer to the minim level, and at κ-CN concentration at intermediate levels. Results highlighted the lower efficiency of HF milk for cheese manufacturing compared to other breeds, even if reared in multi-breed herds. Nevertheless, it is well established that HF cows represent a fundamental resource for farmers and herds profitability, mainly due to their high genetic merit and great MY; in this light, future research is needed to evaluate how and to which extent specialised dairy breeds should be included in respect to dual purpose, and to define optimal breed herd composition within the context of cheese manufacturing.

Ethical approval

Procedures adopted in this study do not fall into the scope of an animal ethics evaluation, as they do not reach the thresholds established in the ‘Directive 2010/63/EU (Art. 1) of the European Parliament and of the council of 22 September 2010 on the protection of animals used for scientific purposes’.

Acknowledgements

The authors gratefully thank the South Tyrolean Dairy Association (Bolzano, Italy) and the Breeders Association of Bolzano province (Bolzano, Italy) for providing the data used in the study.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

None of the data were deposited in an official repository. The data that support the findings presented in this study are available from the first author or corresponding author upon reasonable request.

Additional information

Funding

This research was funded by AGER GA No. 2017-1153-INNOVAMILK and by Italian Ministry of Economic Development CD00628 “LatTraZione 2021 SABELLI”