Genetic background of calcium and phosphorus phases predicted from milk mid-infrared spectra of Holstein cows

Abstract

In bovine milk, Ca and P are partitioned between micellar (MP) and soluble phase (SP), both having important effects on coagulation properties. In particular, greater mineral content in the MP translates into better milk coagulation ability. Nevertheless, the high analytic costs of gold standard methods hamper the possibility to deepen partition of minerals in MP and SP on a large scale. In this study, MP and SP of Ca and P were predicted from 111,653 milk mid-infrared spectra of Holstein cows to investigate genetic parameters. The average coefficient of determination of the prediction models in cross validation was 0.73. Heritability estimates of MP and SP of Ca and P ranged from 0.472 to 0.548 and the two phases of the same mineral were negatively correlated. The MP of Ca was genetically associated with protein yield (0.284) and content (0.658); in the case of MP of P, the latter were equal to 0.262 and 0.808, respectively. The current selection index of Italian Holstein places positive emphasis on protein percentage and yield, thus it is likely that the MP of the investigated minerals is increasing at the expense of the SP. In perspective, it would be important to assess genetic correlations between measured and predicted phenotypes to corroborate the use of such predictions for management and breeding purposes.

Highlights

• Micellar and soluble fractions of Ca and P predicted from milk spectra are heritable.

• Micellar fraction of Ca and P is favourably genetically associated with traits of interest for the dairy industry, such as milk protein and coagulation properties.

• The Italian Holstein selection index emphasises protein percentage and yield, and thus it is indirectly improving micellar phase of Ca and P.

Keywords:

• Milk mineral
• coagulation
• casein micelle
• curd
• genetic correlation

Introduction

Milk and dairy products are important sources of minerals in human nutrition because they provide around 55 and 25% of recommended daily intakes of Ca and P, respectively (Górska-Warsewicz et al. 2019). The amount of certain minerals in milk largely affects technological properties such as enzymatic coagulation and synaeresis. Moreover, soluble Ca acts as molecular bridge between micelles during the formation of paracasein reticulum. Previous studies highlighted differences in micellar content of Ca and P between naturally occurring ‘good’ and ‘poor’ coagulating samples (Malacarne et al. 2014).

Mid-infrared spectroscopy (MIRS) is a cost-effective and non-destructive technique for the analysis of several milk features, including difficult-to-measure traits (Benedet et al. 2020; Grelet et al. 2021) and milk coagulation properties (MCP) (Chessa et al. 2014; De Marchi et al. 2014). Some authors have demonstrated the feasibility of MIRS to predict minerals from milk spectra (Soyeurt et al. 2009; Visentin et al. 2016) but they did not distinguish between micellar (MP) and soluble phase (SP) (Malacarne et al. 2018). On the other hand, a recent study has demonstrated that SP and MP of Ca and P of individual bovine milk can be predicted using backward interval partial least squares analysis (Franzoi et al. 2019).

To the best of our knowledge, no studies have characterised the genetic background of SP and MP of Ca and P in cow milk. Therefore, the aim the present study was to estimate genetic parameters of such new MIRS-predicted traits and to assess their associations with traits of interest for the dairy industry.

Materials and methods

A total of 473,816 test-day records of 23,450 HO cows collected from January 2011 to February 2019 during the monthly recording procedures in Bolzano province were originally available. Each record had information on milk yield (MY, kg/d), percentage of protein (P%) and fat (F%), somatic cell count (cells/μL), parity, days in milk (DIM), birth date, and calving date. Protein yield (PY, kg/d) and fat yield (FY, kg/d) were calculated from MY and P% and F%. Parity ranged from 1 to 15 and DIM from 5 to 305 d. Milk samples were analysed according to International Committee for Animal Recording guidelines using MilkoScan FT6000 until February 2017 and MilkoScan FT7 afterwards (FOSS Electric A/S, Hillerød, Denmark). Standard samples were routinely used to calibrate instruments according to manufacturer’s instructions (Juhl 2017). Somatic cell count was determined using Fossomatic (FOSS Electric A/S, Hillerød, Denmark) and was log₂-transformed to somatic cell score (SCS, units).

Spectral data of all samples from 5000 to 900 cm⁻¹ were stored to allow a posteriori application of the available MIRS models to predict MCP, total mineral composition (Visentin et al. 2016), and mineral phases (Franzoi et al. 2019). Coefficients of determination (root mean square error) in external validation were 0.54 (2.90 min), 0.56 (1.22 min), and 0.52 (9.00 mm) for rennet coagulation time (RCT), curd-firming time (k₂₀), and curd firmness 30 min after enzyme addition to milk (a₃₀), respectively, and 0.67 (122.00 mg/kg) and 0.68 (88.12 mg/kg) for total Ca and total P, respectively.

The models developed by Franzoi et al. (2019) were used for prediction of MP, SP, and the MP to SP ratio (MP/SP) for Ca and P. Briefly, 93 HO samples covering DIM from 6 to 373 were analysed for the content of Ca and P in both SP and MP (Franzoi et al. 2018). Calibration models were developed using backward interval partial least squares algorithm, ending up with coefficients of determination (root mean square error in cross-validation) of 0.77 (2.69) mg/100 mL for SP of Ca, 0.76 (5.30) mg/100 mL for MP of Ca, 0.69 (0.36) for MP/SP of Ca, 0.73 (2.56) mg/100 mL for SP of P, 0.73 (4.76) mg/100 mL for MP of P, and 0.68 (0.22) for MP/SP of P. Belay et al. (2017) reported that such accuracies allow the predictions to be considered good estimates of the actual phenotype.

Spectra with Mahalanobis distance greater than three from spectral centroid were considered as outliers. Moreover, observations with predicted values outside the range of the reference data used for calibrations were set to missing to avoid data extrapolation. Traits deviating more than three standard deviations from the corresponding mean were set to missing. Lactations and contemporary groups (herd-test-date, HTD) with less than 5 test-days and 5 cows, respectively, were discarded from the dataset. The final data consisted of 111,653 test-day records of 9519 HO cows in 338 herds.

Variance components and solutions for the fixed effects were estimated in ASReml v4.1 (Gilmour et al. 2015) through univariate analyses, and covariance components between the traits were estimated through bivariate analyses. The general form of the linear model, in matrix notation, was: $$y\mathrm{ }=\mathbf{X}b\mathrm{ }+{\mathbf{Z}}_{\mathbf{w}}w\mathrm{ }+{\mathbf{Z}}_{\mathbf{a}}a\mathrm{ }+\mathrm{ }e,$$ where y is the vector of dependent variable(s); b is the vector of fixed effects of HTD (n = 10,504), parity (1, 2, 3, 4, and ≥5), and stage of lactation (30 classes); w is the vector of solutions for the random permanent environmental effect of the cow; a is the vector of solutions for the random additive genetic effect of the animal; and e is the vector of random residuals. The incidence matrices X, Z_w, and Z_a linked the corresponding effects to the dependent variable(s) y. Random effects were assumed to be normally distributed with means zero and variance–covariance structures of additive genetic, permanent environmental, and residual effects in the bivariate analyses that were G⊗A, W⊗I, and R⊗I, respectively, where G is the 2 × 2 additive genetic (co)variance matrix, W is the 2 × 2 (co)variance matrix of permanent environmental effects, R is the residual (co)variance matrix, A is the additive genetic relationship matrix among individuals (n = 31,645) which included cows with phenotypic records and 6 generations of ancestors, I is an identity matrix of appropriate order, and ⊗ is the Kronecker product. Heritability (h²), repeatability (t), phenotypic correlation (r_p), and genetic correlation (r_g) were computed as reported by Costa et al. (2019b).

Assuming a given number of offspring per bull, it is possible to estimate the ratio between correlated (ΔG_C) and direct response (ΔG) to selection for a specific trait (y). This is of interest whether two traits (x and y) are correlated, and genetic selection is already in progress for a trait x (Syrstad 1970). For the calculations of this study, the number of offspring per bull was set to 70 in the formula: $$\frac{{\mathbf{\Delta }\mathbf{G}}_{\mathbf{C}}}{\mathbf{\Delta }\mathbf{G}}=\mathrm{ }{\mathit{i}\mathrm{ }·\mathrm{ }\mathbf{r}}_{\mathbf{g}}·\sqrt{\frac{{\mathbf{h}}_{\mathit{x}\mathrm{ }}^2·\mathrm{ }\left[\mathrm{ }4\mathrm{ }+\left(\mathbf{n}-1\right)\mathrm{ }·\mathrm{ }{\mathbf{h}}_{\mathit{y}}^2\right]}{{\mathbf{h}}_{\mathit{y}\mathrm{ }}^2·\mathrm{ }\left[\mathrm{ }4\mathrm{ }+\left(\mathbf{n}-1\right)\mathrm{ }·\mathrm{ }{\mathbf{h}}_{\mathit{x}}^2\right]}},$$ where i is the selection intensity; r_g is the genetic correlation between x and y; n is the number of daughters with information; $h_x^2$ is the heritability of trait x (i.e., MIRS-predicted minerals in the present paper); and $h_y^2$ is the heritability of trait y.

Results and discussion

Average MY and composition (Table 1) agreed with official national statistics of Italian HO breed (AIA 2020). As regards mineral traits, de la Fuente et al. (1996) reported that MP/SP of Ca and P are strongly affected by the method used for their quantification. Previous studies reported values of MP/SP for Ca that ranged from 2.20 to 2.68 w/w, using filtration methods (Holt and Jenness 1984; de la Fuente et al. 1996; Malacarne et al. 2014), and 4.06 w/w based on rennet coagulation procedure (de la Fuente et al. 1996). The MP/SP of P (1.55 w/w) was very close to the value of 1.45 w/w obtained on 12 milk samples by de la Fuente et al. (1996), but higher than values (0.80–1.26 w/w) reported in other studies using different quantification methods (Holt and Jenness 1984; de la Fuente et al. 1996; Malacarne et al. 2014). For the present study, calibrations were developed from reference data using a modified version of the rennet coagulation procedure (Franzoi et al. 2018) and thus results intermediate to those available in the literature were somehow expected.

Table 1. Mean, range, and coefficient of variation (CV) of milk yield and quality traits, mineral fractions, and coagulation properties predicted by mid-infrared spectroscopy.

Trait^a	Mean	Range	CV, %
Milk yield, kg/d	29.71	47.00	25.35
Protein, kg/d	0.96	1.90	0.11
Fat, kg/d	1.16	3.16	9.72
Protein, %	3.31	2.27	10.55
Fat, %	4.02	4.73	17.15
SCS	2.56	13.27	73.82
Minerals
Ca
Total, mg/kg	1351.05	864.52	10.83
Micellar, mg/kg	963.42	871.56	14.25
Soluble, mg/kg	282.54	284.18	15.75
Micellar/soluble, w/w	3.44	4.98	23.01
P
Total, mg/kg	1032.35	774.78	13.37
Micellar, mg/kg	602.94	579.97	15.31
Soluble, mg/kg	365.19	313.66	13.57
Micellar/Soluble, w/w	1.55	2.72	27.55
Milk coagulation properties
RCT, min	22.20	18.31	16.42
k20, min	6.22	7.53	19.32
a30, mm	14.00	42.13	63.24

^aSCS: somatic cell score; RCT: rennet coagulation time; k₂₀: curd-firming time; a₃₀: curd firmness 30 min after enzyme addition to milk.

All the fixed effects included in the statistical model significantly affected the variation of Ca and P in MP and SP (p<.001). The greatest MP/SP was observed in first-parity cows for Ca and in second-parity cows for P (Figure 1). The MP of Ca and P decreased from parity 1 to parity ≥5, whereas the SP of Ca had an opposite trend to that of P (Figure 1). The increase of Ca SP with parity order was likely due to a reduction of casein and thus to a shift of the equilibrium to the SP. The pattern of Ca and P in MP across DIM (Figure 2) resembled the typical lactation curve of F% and P%. Ca and P in SP decreased from early to late lactation; as a result, the MP/SP of both minerals decreased during the first 30 DIM and increased thereafter.

Figure 1. Effect of parity on micellar and soluble phases of Ca and P, and their ratio.

Figure 2. Effect of lactation stage on micellar and soluble phases of Ca and P, and their ratio.

In this study, h² and t of MIRS-predicted MP, SP, and MP/SP of Ca and P from univariate analyses averaged 0.520 and 0.630, respectively, and they were very close to estimates from bivariate models (Table 2). Heritability of mineral phases mirrored that of P% (0.469 ± 0.014) and it was greater than h² of F% (0.356 ± 0.011). The lowest h² was estimated for MP/SP of P (0.472 ± 0.013). Costa et al. (2019a) reported h² of 0.423, 0.430, 0.446, and 0.531 for MIRS-predicted P%, casein percentage, total Ca, and total P, respectively, in Italian HO. Similarly, in the same population, Visentin et al. (2019) estimated average h² of 0.540 and 0.420 for MIRS-predicted Ca and P, respectively. Apart from the MP/SP of P, the t of other Ca and P fractions was greater than 0.600, suggesting that few predicted records are adequate to explain the variability of the phenotype and get reliable estimates of cows’ mean and breeding values for the phenotype itself.

Table 2. Additive genetic variance (σ²_a), permanent environmental variance (σ²_w), heritability (h²), and repeatability (t) of Ca and P phases from univariate analyses^a, and average, minimum and maximum heritability estimates from bivariate analyses.

Trait	Univariate analyses				Bivariate analyses
	σ^2a	σ^2w	h^2	t	Average	Minimum	Maximum
Ca
Micellar, mg/kg	55.637	10.316	0.520	0.616	0.519	0.518	0.521
Soluble, mg/kg	7.149	1.272	0.518	0.610	0.517	0.516	0.518
Micellar/soluble, w/w	0.238	0.044	0.548	0.650	0.541	0.541	0.553
P
Micellar, mg/kg	30.744	6.775	0.513	0.626	0.512	0.511	0.516
Soluble, mg/kg	8.115	1.634	0.533	0.640	0.532	0.531	0.534
Micellar/soluble, w/w	0.047	0.009	0.472	0.561	0.471	0.470	0.475

^aStandard errors of heritability and repeatability ranged from 0.013 to 0.014 and 0.004 to 0.005, respectively.

The strongest and weakest r_p were estimated between the two MP (0.896 ± 0.002) and between the two SP (–0.132 ± 0.008), respectively. As regards r_g, the MP and MP/SP of Ca were positively strongly correlated with MP and MP/SP of P, and for both minerals the SP was negatively genetically associated with the MP and the ratio (Table 3). The strong relationship between MP of Ca and MP of P is likely attributable to their complementary and structural role in casein micelle, whereas the weak association between SP of Ca and SP of P suggests that they have different genetic background. Moreover, within mineral, the MP was negatively genetically correlated with the SP, meaning that it is possible to adopt selection strategies to shift mineral equilibrium according to the desired objective. In particular, from cheesemakers’ point of view, a shift to MP is preferred.

Table 3. Phenotypic (below diagonal) and genetic correlations^a (above diagonal) between Ca and P phases and their ratio.

Trait	Ca			P
	Micellar	Soluble	Micellar/Soluble	Micellar	Soluble	Micellar/Soluble
Ca
Micellar, mg/kg		–0.299	0.735	0.906	–0.317	0.897
Soluble, mg/kg	–0.256		–0.834	–0.530	–0.210	–0.240
Micellar/soluble, w/w	0.720	–0.809		0.871	–0.075	0.642
P
Micellar, mg/kg	0.896	–0.477	0.853		–0.086	0.760
Soluble, mg/kg	–0.396	–0.132	–0.160	–0.177		–0.643
Micellar/soluble, w/w	0.886	–0.243	0.643	0.768	–0.670

^aStandard errors ranged from 0.002 to 0.008 and 0.010 to 0.025 for phenotypic and genetic correlations, respectively.

The r_p and r_g of MIRS-predicted minerals fractions with MY, P%, PY, F%, FY, SCS, and MCP are displayed in Table 4. The MP of Ca and P was negatively genetically associated with MY, whereas the r_g between MY and SP was 0.277 ± 0.036 in the case of Ca and close to zero in the case of P. As a result, the MP/SP of Ca and P was negatively associated with MY. However, the ratio was positively genetically associated with PY and FY for both Ca and P. In general, the negative correlation between MP and MY confirms the lower concentration of micellar Ca and P in high-producing cows, in accordance with studies that dealt with r_g and r_p of MY with P% and casein content (Cassandro et al. 2008). The MP of both minerals was positively associated with P% and F%, both phenotypically and genetically (Table 4), which is favourable for milk destined to produce cheese. Correlations of predicted Ca and P fractions with SCS were weak or close to zero (Table 4). Finally, all fractions of Ca and P were favourably associated with RCT, except for SP of P, while MP and MP/SP of both minerals were favourably correlated with k₂₀ and a₃₀ (Table 4).

Table 4. Phenotypic (r_p) and genetic (r_a) correlations^a of milk Ca and P phases (mg/kg) and their ratio (w/w) with yield and quality traits, and coagulation properties.

Trait^b		Ca			P
		Micellar	Soluble	Micellar/Soluble	Micellar	Soluble	Micellar/soluble
Milk yield, kg/d	rp	–0.256	0.088	–0.233	–0.270	0.059	–0.211
	ra	–0.406	0.277	–0.480	–0.524	0.013	–0.318
Protein, kg/d	rp	0.041	–0.099	0.074	0.089	0.036	0.026
	ra	0.284	0.022	0.098	0.262	–0.020	0.198
Fat, kg/d	rp	0.053	0.053	0.039	0.045	–0.042	–0.10
	ra	0.444	–0.124	0.354	0.407	–0.198	0.370
Protein, %	rp	0.647	–0.417	0.665	0.783	–0.066	0.519
	ra	0.658	–0.449	0.694	0.808	–0.002	0.520
Fat, %	rp	0.399	–0.038	0.332	0.402	–0.1364	0.261
	ra	0.649	–0.358	0.654	0.700	–0.163	0.551
SCS	rp	–0.018	–0.004	–0.023	–0.039	0.064	–0.049
	ra	–0.098	–0.042	–0.046	–0.106	0.113	–0.132
RCT, min	rp	–0.134	–0.030	–0.091	–0.171	0.241	–0.243
	ra	–0.201	–0.020	–0.117	–0.207	0.236	–0.326
k20, min	rp	–0.509	0.221	–0.466	–0.650	0.109	–0.413
	ra	–0.576	0.295	–0.535	–0.732	0.083	–0.524
a30, mm	rp	0.314	–0.108	0.288	0.406	–0.186	0.375
	ra	0.376	–0.134	0.317	0.441	–0.159	0.422

^aStandard errors ranged from 0.003 to 0.008 and 0.009 to 0.044 for phenotypic and genetic correlations, respectively.

^bSCS: somatic cell score; RCT: rennet coagulation time; k₂₀: curd-firming time; a₃₀: curd firmness 30 min after enzyme addition to milk.

An indirect favourable selection for the MP of Ca and P likely exists in the Italian HO, since the official selection index emphasises P% (3%), PY (36%), F% (2%), and FY (8%), and does not place any weight to MY (ANAFIJ 2019). Supporting this, in the case of P% the ratio between ΔG_C and ΔG was 0.65 for Ca and 0.80 for P; in the case of PY, the ratio was 0.23 for Ca and 0.21 for P. Such results based on MIRS-predicted traits, suggested that an indirect selection for MP of both minerals seems to be currently present and favourable. As expected, for both minerals, indirect selection based on P% or PY was less efficient than direct selection on the MP itself.

Conclusions

In the present study, MIRS-predicted micellar and soluble fractions of Ca and P showed exploitable genetic variation and were associated with traits of interest for the dairy industry. The current selection index of Italian Holstein places positive emphasis on protein yield and content thus, according to the findings obtained for the MIRS-predicted traits, this is expected to indirectly increase the micellar Ca and P and to favourably affect milk coagulation ability. Future studies should investigate the feasibility of enhancing the robustness of MIRS models to predict minerals fractions. This could be achieved by increasing the reference dataset and testing different prediction algorithms. Having a larger dataset would allow to assess genetic correlations between measured and predicted phenotypes, and corroborate the use of such predictions for management and breeding purposes.

Acknowledgments

The South Tyrolean Dairy Association (Bolzano, Italy), the Breeders Association of Bolzano Province (Bolzano, Italy), and the Italian Holstein and Jersey Breeders Association (Cremona, Italy) are gratefully acknowledged for providing data used in this study.

Ethical approval

The approval of the Ethical Committee for the Care and Use of Experimental Animals was not required for this study because the data were extrapolated from an existing database.

Disclosure statement

This article has been accepted for publication in Italian Journal of Animal Science, published by Taylor & Francis. The authors are fully responsible for the content of the manuscript and declare that there is no conflict of interest associated with the paper.