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.


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 orska-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.  (Juhl 2017). Somatic cell count was determined using Fossomatic (FOSS Electric A/S, Hillerød, Denmark) and was log 2 -transformed to somatic cell score (SCS, units).

Materials and methods
Spectral data of all samples from 5000 to 900 cm À1 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 20 ), and curd firmness 30 min after enzyme addition to milk (a 30 ), 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 (herdtest-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: 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 GA, WI, and RI, 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 2 ), 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 (DG C ) and direct response (DG) 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: 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 2 x is the heritability of trait x (i.e., MIRS-predicted minerals in the present paper); and h 2 y is the heritability of trait y.

Results and discussion
Average MY and composition ( 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.
In this study, h 2 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 2 of F% (0.356 ± 0.011). The lowest h 2 was estimated for MP/SP of P (0.472 ± 0.013). Costa et al. (2019a) reported h 2 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 2 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.
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 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 2. Additive genetic variance (r 2 a ), permanent environmental variance (r 2 w ), heritability (h 2 ), and repeatability (t) of Ca and P phases from univariate analyses a , and average, minimum and maximum heritability estimates from bivariate analyses.  Table 3. Phenotypic (below diagonal) and genetic correlations a (above diagonal) between Ca and P phases and their ratio.  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.  (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 20 and a 30 (Table 4).
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 DG C and DG 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.

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.