Introduction
Excess phosphorus (P) excretion into the environment contributes to eutrophication of surface and groundwater, and is considered a major environmental concern (Smith and Alexander, Reference Smith and Alexander2000). Therefore, several countries have implemented, or are currently implementing, regulations aimed at reducing the amount of manure P applied onto farmland (e.g., EU Water Framework Directive; EC, 2000). In dairy cattle systems, P excretion of cattle and milk P efficiency (output of P in milk as a fraction of P intake) varies widely (Harrison et al., Reference Harrison, Dorigo, Reynolds, Sinclair, Dijkstra and Ray2021). Milk P content is not affected by dietary P content (Keanthao et al., Reference Keanthao, Goselink, Dijkstra, Bannink and Schonewille2021), but milk P efficiency is negatively related to dietary P intake (Klop et al., Reference Klop, Ellis, Bannink, Kebreab, France and Dijkstra2013). Besides, milk P efficiency is related to the P content of milk, but in farm level analyses, a fixed milk P content is usually assumed (e.g., 90 and 100 mg/100 g of milk; NRC, 2001 and CVB, 2005, respectively). A fixed milk P content can be considered convenient for the purpose of calculating P requirements of dairy cattle. However, for both environmental and economic reasons, it is of great practical relevance to accurately assess the amount of P secreted in milk and therefore not ending up in manure. The use of a fixed milk P content to calculate the amount of P secreted in milk is not always accurate (Pfeffer et al., Reference Pfeffer, Beede, Valk, Pfeffer and Hristov2005; Bannink et al., Reference Bannink, Šebek, Dijkstra, Vitti and Kebreab2010; Klop et al., Reference Klop, Ellis, Blok, Brandsma, Bannink and Dijkstra2014). The concentrations of milk protein and milk fat, and to a lesser extent milk lactose, vary greatly in practice and this variation is associated with variation in the milk P content. Indeed, the milk P content was found to be related to milk protein content (Wu et al., Reference Wu and Satter2000; Bijl et al., Reference Bijl, Van Valenberg, Huppertz and Van Hooijdonk2013; Klop et al., Reference Klop, Ellis, Blok, Brandsma, Bannink and Dijkstra2014), milk fat content (NRC, 2001) and milk lactose content (Klop et al., Reference Klop, Ellis, Blok, Brandsma, Bannink and Dijkstra2014). Previously, Klop et al. (Reference Klop, Ellis, Blok, Brandsma, Bannink and Dijkstra2014) attempted to predict milk P content and found that around 30% of the variation in milk P content remained unexplained by their best regression models. However, they used only milk protein, lactose and fat contents as predictor variables. Thus, it is possible that the relatively high unexplained variation in milk P content reflected their somewhat limited set of predictors in their regression models. Furthermore, Klop et al. (Reference Klop, Ellis, Blok, Brandsma, Bannink and Dijkstra2014) used milk production data from individual cows from three separate experiments, which may limit the generalizability of their results to the Dutch dairy cow population. The aim of the current study was therefore twofold; 1) to improve the prediction of the P content in bovine milk, and 2) to obtain predicted P contents of bovine milk representative for the Dutch dairy cow population. To this end, collected bulk tank milk samples were collected from 14 dairy plants located across the Netherlands over one year and used a broader spectrum of potential predictor variables compared to Klop et al. (Reference Klop, Ellis, Blok, Brandsma, Bannink and Dijkstra2014). In the process, information on the seasonal variation in the various measured milk constituents was obtained, analogous to a previous analysis done in 2005/2006 (Heck et al., Reference Heck, Van Valenberg, Dijkstra and Van Hooijdonk2009).
Materials and methods
Collection of samples and chemical analyses
Over one year, from week 14 in 2017 up until week 13 in 2018, bulk milk samples were collected weekly from 14 dairy plants across the Netherlands. After collection, the samples were pooled across the 14 dairy milk plants to obtain one composite sample per week, resulting in a total of 52 pooled samples. Each dairy plant represented a region. Therefore, the collected milk samples, predominantly originating from Holstein-Friesian cows, were considered representative of the average Dutch bovine milk (described previously by Pacheco-Pappenheim et al., Reference Pacheco-Pappenheim, Yener, Heck, Dijkstra and van Valenberg2021). Directly after pooling, milk samples were preserved with 0.03% sodium azide pending chemical analysis. The 52 weekly milk samples were analysed for their selected milk constituents at Qlip laboratory (Qlip, Zutphen, the Netherlands). Milk fat was analysed using the Röse-Gottlieb method (ISO 1211, ISO, 2010). Both the total N content of milk and the N originating from non-protein N compounds (NPN-compounds) were analysed using the Kjeldahl method according to ISO 8968-1 (ISO, 2014). The N associated with the NPN-compounds in supernatants was measured in milk samples treated with trichloroacetic acid to precipitate the milk proteins (i.e. true milk protein). The N associated with true milk protein was subsequently calculated as total milk-N minus N from NPN-compounds. Finally, true milk protein was calculated as 6.38 × N, while NPN-compounds were calculated as 6.25 × N. Milk urea was analysed using an enzymatic method according to ISO 14637 (ISO, 2004b). Milk lactose and casein content were analysed using high-performance liquid chromatography (ISO 5548, ISO, 2004a). The content of non-casein protein (NCP) in milk was calculated as the difference between the true milk protein content and the casein content of milk. The milk contents of Na, Mg, P, K, Ca, Mn and Cu were determined using an inductively coupled plasma mass spectrometer based on ISO 21424 (ISO, 2018).
Statistical analyses
Monthly variation in milk constituents
For each month, the weekly values (n = 52) on milk composition were used, i.e. the 4 or 5 samples that present the month. All data on the monthly milk composition were subjected to ANOVA (IBM SPSS Statistics 27th edition, 2020) using the following model:
$${Y_{\rm{ij}}} = \mu + {M_i} + {e_{ij}},$$
where Yij = the response variable, µ = the overall mean, Mi = month (i = 1 to 12) and eij = residual error. The level of statistical significance was preset at P < 0.050.
Prediction of the P content of milk
First, Pearson’s correlation coefficients and linear regression equations were calculated to detect individual milk variables that were related to the milk P content using SPSS Statistics 24.0 software (IBM SPSS Statistics 27th edition, 2020). Then, multiple regression analyses were performed using SPSS Statistics 24.0 software (IBM SPSS Statistics 27th edition, 2020) with milk P content as dependent variable and the various milk constituents as potential predictor variables. Stepwise regression was performed by incorporating or removing milk constituents showing the highest or lowest partial correlation coefficients for their relation to the residual variance in the milk P content. Both the explained variance (adjusted R 2; R 2adj) and Mallows’ Cp (IBM SPSS Statistics 27th edition, 2020) were used as principal determinants to identify the most suitable models to predict the milk P content. R 2adj assesses the fit of the model to the data, adjusting for the number of predictors, while Mallows’ Cp evaluates the predictive ability of the model by considering both fit and complexity.
The various equations were evaluated using the mean square prediction error (MSPE) which was calculated as: MSPE =
$\mathop \sum \nolimits_{i \,=\, 1}^n {\left( {Oi - Pi} \right)^2}/n$
, where n is the total number of observations and Oi and Pi are the observed and predicted values, respectively. The square root of MSPE (RMSPE) is considered an estimate of the overall prediction error and is expressed as a percentage of the observed mean. The variance inflation factor (VIF) was used to check for multicollinearity between predictor variables in multiple regression models (IBM SPSS Statistics 27th edition, 2020) and was calculated as: VIFi = 1/(1 – R 2i), where R 2i is the square of the i th Pearson correlation coefficient between predictor variables in the model. The Bayesian information criterion (BIC) was used to indicate the goodness of fit, where smaller values indicate a better fit (IBM SPSS Statistics 27th edition, 2020). The residual slopes were tested for significance from zero values to verify for heteroscedasticity. Throughout, the level of statistical significance was preset at P < 0.050.
Results
Variation in milk composition throughout the year
The lactose content of milk differed between months (Table 1). However, the variation in lactose content of milk between months was quantitatively small (Figure 1). In contrast, the fat content of milk varied considerably between months (P < 0.001) and the lowest fat contents of milk were found from the end of June up until half August (i.e. summer). Likewise, the lowest protein contents of milk were also observed during summer but compared to the fat content of milk, the variation in protein contents between months was smaller (P < 0.001). The seasonal variation in the casein contents of milk mirrored the variation in milk protein contents (Figure 1). The content of non-casein protein (Table 1) was found to be significantly affected by month (P < 0.001) but the absolute difference between the minimum and maximum values was small. The contents of NPN-compounds (Table 1) were found to be similar between months (P = 0.237), whereas the urea contents of milk varied between months (P = 0.023), with the highest values being observed in August while the lowest values were found in March and May (Table 1).
Table 1. Descriptive statistics on the pooled composition of bovine raw milk, collected weekly in the Netherlands from April 2017 up until March 2018 from 14 dairy plants located across the country (n = 52). The mean, minimum and maximum values are the monthly means of the 4 or 5 samples that represented the month in question. All units are expressed per 100 g of milk and the P-values were generated with the use of ANOVA
a P-value: the value was obtained by analysis of variance using the following model: Yij = µ + Mi + eij, where Yij = the response variable, µ = the overall mean, Mi = month (i = 1 to 12) and eij = the residual error. The level of statistical significance was preset at P < 0.050.
b Equal values in February and June.
c Equal values in July and August.
d Equal values in March and May.
NCP = Non casein protein, NPN = Non protein nitrogen compounds, SD = Standard deviation, CV = Coefficient of variation, expressed as % of the mean.
Figure 1. Weekly variation in the contents of fat- (▲), protein- (●), casein- (◆) and lactose (■) in Dutch raw bovine milk collected in the Netherlands from April 2017 up until March 2018. Week 14 to 52 indicate April to December 2017; week 1 to 13 indicate January to March 2018 (n = 52).
Except for the Na and Cu content of milk, all other measured minerals in milk were affected by month (Table 1; P < 0.050). With respect to the Ca and P contents of milk, the lowest values were found during summer (Figure 2). Likewise, the lowest Mg content of milk was found in July (Table 1) while the highest values were observed in December. For K, the time window between the lowest and highest content in milk was small, i.e. in March and April, respectively. The Mn contents of milk varied considerably between months (Table 1; P < 0.001) with lowest values in milk during August and greatest values during March.
Figure 2. Weekly variation in the contents of Ca (▲), P (■) and K (●) in Dutch raw bovine milk collected in the Netherlands from April 2017 up until March 2018. Week 14 to 52 indicate April to December 2017; week 1 to 13 indicate January to March 2018 (n = 52).
Linear relationships between individual milk constituents and the milk P content
Except for the lactose- and the Cu content, the contents of all other individual milk constituents were found to be related (P < 0.050) with the P content of milk (Table 2). Except for the NPN and urea content, these milk constituents were positively related to milk P content. The variation in either the fat-, protein- or the casein content of milk each explained ∼55% of the variation in the milk P content (P < 0.001), and in the case of the protein- and casein content of milk, the intercepts of the linear regressions were not different from zero (P ≥ 0.130). The variation in milk proteins other than casein also accounted for variation in milk P (P < 0.010) but to a lesser extent compared to fat, protein or casein (Table 2). In contrast to the aforementioned milk constituents, the urea content of milk was negatively related to the milk P content. The variation in milk urea contents explained only 17% of the variation in milk P (P < 0.001). Likewise, the content of NPN-compounds of milk was also negatively related to the milk P content (P < 0.001) and variation in the contents of NPN-compounds of milk accounted for only 9% of the variation in milk P content.
Table 2. Pearson’s correlation coefficients* and linear regression equations on the relationship between individual milk constituents and the content of phosphorus in milk (Pmilk, n = 52).a,b
* Indicates a non-significant correlation (P > 0.05).
a The contents of lactose, fat, protein, casein, non-casein protein (NCP) and non-protein nitrogen compounds (NPN) are expressed as g/100 g of milk. The contents of urea, P, Ca, Na, K and Mg in milk are expressed as mg/100 g of milk and the contents of Cu and Mn are expressed as µg/100 g of milk.
b Pearson’s R = Pearson’s correlation coefficient, R 2adj = adjusted R 2 where R 2 is the square of Pearson’s R, P constant and P slope = the probability of type I error of the intercept and slope, respectively, of the regression equations.
Both the Ca and Mg content of milk correlated highly with the P content of milk (Table 2), explaining 67 and 62% of the variation in milk P, respectively (P < 0.001). The variation in the contents of Na and K also contributed (P < 0.001) to the explained variation in the milk P content, but to a lesser extent compared to Ca and Mg. The intercepts of the linear regression equations of Na, K, Ca and Mg content with P content did not differ from zero (P ≥ 0.289). The variation in the Mn content of milk accounted for 49% of the variation in milk P content (P < 0.001).
Multiple regression models to predict the milk P content
Various predictor sets were used to identify the most suitable equation to predict the milk P content (Table 3). The eight sets of predictor variables were designed to explore different approaches to predict milk P. The first set includes all potential predictors, thereby providing a comprehensive baseline. The second set comprises routinely measured variables, following the reasoning of Klop et al. (Reference Klop, Ellis, Blok, Brandsma, Bannink and Dijkstra2014). The following three sets (3 to 5) were informed by the results of preceding models, allowing us to refine variable selection iteratively. The last three sets (6 to 8) began with the three minerals showing the highest adjusted R2 when linearly regressed against milk P, in order to evaluate how the inclusion of a strong predictor as the first variable affects the contribution of remaining candidates. To minimize collinearity, candidate variables in sets 2 to 8 were only included when Pearson correlation coefficients were below 0.7.
Table 3. Overview of the candidate predictor variables used for linear stepwise, multiple regression to predict the P content (mg/100 g) of Dutch raw bulk tank milk
a Within predictor sets 2 to 8, Pearson’s correlation coefficient between candidate predictor variables was <0.7. The contents of fat, protein, lactose, casein, non-casein protein (NCP) and non-protein nitrogen compounds (NPN) are expressed as g/100 g of milk. The contents of urea, P, Ca, Mg, Na and K in milk are expressed as mg/100 g of milk and the contents of Cu and Mn are expressed as µg/100 g of milk.
Upon stepwise regression of all measured milk constituents against the milk P content, the fat-, Ca- and K contents of milk were found to significantly contribute to the explained variance in the milk P content (Table 4, equation 1), with the lowest RMSPE, Mallows’ Cp and BIC of all prediction models evaluated. The VIF value of this model was 3.29, which is just below the VIF value threshold of 3.3 indicative of collinearity (Kock and Lynn, Reference Kock and Lynn2012). Therefore, various candidate models were compiled with the aim to reduce multicollinearity between the candidate predictor variables.
Table 4. Overview of the regression equations to predict milk phosphorus content (Pmilk), expressed as mg/100 g of milk. The multiple regression equations were obtained after stepwise regression using the candidate predictor variablesa indicated in Table 3 (for all models P < 0.001, n = 52)
a The contents of fat, protein, lactose, casein and non-casein protein (NCP) are expressed as g/100 g of milk. The contents of urea, Ca, Mg and K in milk are expressed as mg/100 g milk and the Mn content is expressed as µg/100 g of milk.
b R 2adj = adjusted R2 where R2 is the square of Pearson’s correlation coefficient.
c Mallows’ Cp = Mallow’s criterion for a model with p predictor variables.
d VIF = Variance inflation factor.
e RMSPE = Square root of the mean square prediction error and is expressed as a percentage of the observed mean milk P content, i.e. 101.2 mg/100 g of milk.
f Bayesian information criterion.
g NCP = Non casein protein.
All subsequent multiple regression models (Table 4, equations 2–8) had a VIF value < 3.3. The model containing only the protein and urea content of milk as predictor variables (Table 4, equation 2) accounted for only 59% of the variation in the milk P content. The models described with equations 3 to 7 (Table 4) were all similar to explain the variation in milk P content (76% < R 2adj < 79%) with varying performance when based on Mallows’ CP or BIC. The model containing, amongst others, Mn (equation 8) explained the variation in the milk P content to a somewhat lower extent (R 2adj = 73%) compared to equations 3 to 7. The 4 most promising equations to predict P content of milk, i.e. equations 1, 3, 5 and 7, were assessed for heteroscedasticity by regressing the residuals against the predicted values (Figure 3). Visual inspection of the graphs indicated no evidence of heteroscedasticity for any of the 4 equations (Figure 3). When predicted milk P contents were plotted against observed milk P contents (Figure 3), minor underprediction was apparent for equation 7, whereas equations 1, 3 and 5 slightly overpredicted milk P content.
Figure 3. Predicted vs observed milk P content (left) and residuals (i.e., observed – predicted) vs predicted values for milk P content (mg/100 g) of selected multiple regression models from models 1, 3, 5 and 7 (right) (n = 52). In each residual plot, the solid line represents a fitted regression. R2 adj is the adjusted R2 where R2 is the square of Pearson’s correlation coefficient and the P-value represents the probability of type I error of the slope of the regression line.
Finally, the prediction model reported by Klop et al. (Reference Klop, Ellis, Blok, Brandsma, Bannink and Dijkstra2014), i.e. P in milk (g/kg) = −0.64 + 0.0223 × milk protein (g/kg) + 0.0191 × milk lactose (g/kg), was applied to the current dataset to estimate milk P content (Figure 4). Visual inspection of the graph, however, indicated a poor prediction performance and evidence of heteroscedasticity when the model of Klop et al. (Reference Klop, Ellis, Blok, Brandsma, Bannink and Dijkstra2014) was used to predict the milk P content (Figure 4).
Figure 4. Predicted vs observed milk P content (left) and residuals (i.e., observed – predicted) vs predicted values for milk P content (g/kg) when the regression equation reported by Klop et al. (Reference Klop, Ellis, Blok, Brandsma, Bannink and Dijkstra2014) is used to predict the milk P content, i.e., P in milk (g/kg) = −0.64 + 0.0223 × milk protein (g/kg) + 0.0191 × milk lactose (g/kg), (n = 52). In the residual plot, the dotted line represents a fitted regression. R2adj is the adjusted R2 where R2 is the square of Pearson’s correlation coefficient and the P-value represents the probability of type I error of the slope of the regression line.
Discussion
Variation in milk composition throughout the year
Holstein-Friesian is the predominant dairy breed in the Netherlands (CRV, 2021) and the current data reflects the milk composition of Holstein-Friesian cows. In the Netherlands, where calving occurs generally all year round, changes in milk composition between months are likely related to seasonal changes in the cows’ diet (Heck et al., Reference Heck, Van Valenberg, Dijkstra and Van Hooijdonk2009). Cows may graze outside from ∼ April until September/October and during the grazing season, cows consume fresh grass either or not in combination with ensiled feed (maize silage and grass silage) and, depending on the milk yield, supplemental concentrates (CBS, 2022). In late autumn and winter, the diet consumed by cows in the Netherlands mainly includes ensiled feed and concentrate (CBS, 2022). Pacheco-Pappenheim et al. (Reference Pacheco-Pappenheim, Yener, Heck, Dijkstra and van Valenberg2021) evaluated variation in milk fat content and milk fatty acid profile based on the same milk samples as obtained in the present analysis in relation to seasonal changes in diet composition. They argued that the lower milk fat content in summer periods is likely related to increased concentrations of fatty acid intermediates that arise from biohydrogenation of polyunsaturated fatty acids in grass herbage, leading to decreased de novo synthesis of fatty acids in the mammary gland. The increased milk protein content in winter months compared with summer months may be related to changes in glucogenic and aminogenic nutrient supply (Lemosquet et al., Reference Lemosquet, Guinard-Flament, Raggio, Hurtaud, Van Milgen and Lapierre2010). A higher ratio of concentrate to roughage in the diet in the winter season is associated with lower levels of fibre and higher levels of starch in the diet, which stimulates the production of propionic acid in the rumen (Bannink et al., Reference Bannink, Kogut, Dijkstra, France, Kebreab, Van Vuuren and Tamminga2006) and increases microbial protein supply (Dijkstra et al., Reference Dijkstra, France and Davies1998). Increasing glucogenic energy at various levels of protein supply may improve postabsorptive amino acid transfer efficiency from the gut to milk protein (Rius et al., Reference Rius, McGilliard, Umberger and Hanigan2010). Indeed, both increased propionic acid supply and increased amino acid supply stimulated milk protein content (Raggio et al., Reference Raggio, Lobley, Lemosquet, Rulquin and Lapierre2006). In the current study, the variation in the lactose content of milk was small, which is in line with other studies (Costa et al., Reference Costa, Lopez-Villalobos, Sneddon, Shalloo, Franzoi, De Marchi and Penasa2019). Lactose is the main osmotic regulator between blood and mammary gland alveolar lumen, with lactose determining the amount of absorbed water from blood and thus milk volume produced (Costa et al., Reference Costa, Lopez-Villalobos, Sneddon, Shalloo, Franzoi, De Marchi and Penasa2019). Milk urea concentrations were lowest in March and May and highest in August, likely reflecting seasonal variation in dietary protein and energy balance. It is plausible that the relatively high protein content of summer pasture, combined with lower OM digestibility and thus energy availability, leads to less efficient N utilization and increased excretion as urea (Spek et al., Reference Spek, Dijkstra, van Duinkerken, Hendriks and Bannink2013). By contrast, diets in early spring have a higher digestibility with more energy available both in the rumen to enable efficient microbial protein synthesis as well as post-absorptive to enable efficient synthesis of milk protein, culminating in lower milk urea levels (Van Vuuren et al., Reference Van Vuuren, Van den Pol-Van Dasselaar, Elgersma, Dijkstra and Tamminga2006).
Seasonal variation in milk mineral concentrations may reflect dietary mineral supply, pasture composition and environmental factors. Indeed, previous work shows associations between feeding practices, supplementation and milk composition (Alothman et al., Reference Alothman, Hogan, Hennessy, Dillon, Kilcawley, O’Donovan, Tobin, Fenelon and O’Callaghan2019; Dunshea et al., Reference Dunshea, Walker, Williams and Doyle2019), suggesting these factors could contribute to the observed seasonal patterns. However, an elaboration on the potency of the dietary mineral supply to affect milk mineral concentrations is considered not opportune because our study lacks data on feed mineral content.
The values of P, Ca, Mg, K, Na and Cu contents of bovine milk adopted by the Central Bureau for Livestock Feeding (CVB, 2005), are 100, 120, 12, 150, 46 mg/100 g and 4 µg/100 g, respectively. Except for the Na and Cu content of milk, the overall mean mineral contents of bovine milk in the current study are thus in line with the values adopted by the CVB (2005). It is, however, of interest to note that the mean P content of milk was found to be 1.2% greater than the fixed milk P content (100 mg/100 g of milk; CVB, 2005) or values used in farm level analyses to calculate P excretion with manure (e.g., Ireland, 90 mg P/100 g of milk, Mihailescu et al., Reference Mihailescu, Murphy, Ryan, Casey and Humphreys2015; Australia, 91 mg P/100 g of milk, Gourley et al., Reference Gourley, Dougherty, Weaver, Aarons, Awty, Gibson, Hannah, Smith and Peverill2012). Under the assumption that the whole farm P intake of the cows is accurately assessed, it thus appears that the calculated P excretion with manure is systematically overrated.
The Na contents of milk in the current study were ∼24% lower compared to the value reported by the CVB, whereas the Cu content of milk was ∼48% greater in the current study. The CVB (2005) values on Na and Cu contents of milk are, however, means estimated on the basis of literature studies. The current observations on the Na and Cu contents of bovine milk fall within the ranges reported by Schonewille and Beynen (Reference Schonewille and Beynen2005) and Van den Top (Reference Van den Top2005), respectively. The variation of milk Na, Ca, P and Mg generally resembles the variation in the protein content of milk and to a somewhat lesser extent the variation in milk fat content. This similarity in variation is most likely related to the process of casein and milk fat formation (see also next paragraph) in the mammary gland (Shennan and Peaker, Reference Shennan and Peaker2000; Mather, Reference Mather, Fuquay, Fox and McSweeney2011). The highest and lowest milk K contents were found in March and April, respectively. In April, the majority of dairy cows start grazing grass herbage again, and the high K content of grass herbage (especially in spring) compared with K content in maize silage and concentrate (CBS, 2022) may explain the observed variation in the milk K content. With respect to the Mn content of milk, a reference value is not provided by the CVB (2005). The current mean Mn content is somewhat lower than that reported by Castillo et al. (Reference Castillo, St-Pierre, Silva del Rio and Weiss2013), i.e. 3.0 µg/100 g. There are indications that Mn is involved in the process of protein and fat synthesis in the mammary gland (Lönnerdal et al., Reference Lönnerdal, Keen and Hurley1984) but numerically the difference between the minimum and maximum value of milk Mn content is small and is not in line with the corresponding values of protein and fat content of milk.
Linear relationships between individual milk constituents and the milk P content
Inorganic P is released during the formation of lactose in the Golgi apparatus (Shennan and Peaker, Reference Shennan and Peaker2000), thereby indicating that variation in milk P content can be, at least partly, explained by variation in the lactose content of milk. In the current study, however, no significant relationship was found between the contents of lactose and P in milk which is in contrast with results reported by Klop et al. (Reference Klop, Ellis, Blok, Brandsma, Bannink and Dijkstra2014). This discrepancy in results is probably caused by the difference in variation around the mean lactose and P content of milk between the two studies. In the study reported by Klop et al. (Reference Klop, Ellis, Blok, Brandsma, Bannink and Dijkstra2014), the weighed mean lactose content of milk was calculated to be 45.5 g/kg of milk with a coefficient of variation (CV) of 3.8%. In the current study, the mean lactose content was 45.1 g/kg milk with CV of only 0.49%. Likewise, the weighed mean milk P content and CV in the study reported by Klop et al. (Reference Klop, Ellis, Blok, Brandsma, Bannink and Dijkstra2014) were 0.993 g/kg milk and 8.0%, respectively, while the mean P content of milk and CV were, respectively, 1.012 g/kg milk and 3.5% in the current study.
The variations in the contents of fat, protein or casein in milk explained a significant amount of variation in the milk P content. Milk fat is known to be excreted in the form of droplets surrounded by a monolayer of phospholipids (Shennan and Peaker, Reference Shennan and Peaker2000; Mather, Reference Mather, Fuquay, Fox and McSweeney2011), which might explain the observed positive relationship between the fat and P content of milk. The process of casein synthesis in the Golgi vesicles is associated with the release of inorganic P (Shennan and Peaker, Reference Shennan and Peaker2000), and esterification of calcium phosphate with the protein matrix of casein micelles also occurs (Holt, Reference Holt2004; Bijl et al., Reference Bijl, Van Valenberg, Huppertz and Van Hooijdonk2013). The latter may explain the observed positive relationship between both casein and protein content and the milk P content because the contents of milk protein and casein are highly correlated (Pearson’s r = 0.985). The current data corroborate the notion that the milk P content is primarily related to the casein content of milk because the variation in the non-casein protein content explained only 34% of the variation in milk P, while the variation in the contents of NPN-compounds or that of urea explained ≤17% of the variation in the milk P content. The binding of casein micelles with calcium phosphate (Ahmadi et al., Reference Ahmadi, Markoska, Huppertz and Vasiljevic2024) also explains the high correlation between the Ca and P content of milk. In fact, the variation in Ca content of milk was found to explain the greatest amount of variation in milk P content compared to the other selected milk constituents. Next to Ca and phosphate, Mg also is, at least partly, associated with the casein micelles (Holt, Reference Holt2004). This finding is in line with the current data in that the Ca and Mg contents were highly correlated (Pearson’s r = 0.943). The high correlation between the Ca and Mg content of milk is corroborated by data reported by Bijl et al. (Reference Bijl, Van Valenberg, Huppertz and Van Hooijdonk2013).
Potassium likely contributes to charge balance in milk independent of casein-bound minerals, such as P and Ca, by compensating for negatively charged anions like citrate, sulphate and chloride. In our dataset, milk K content showed only a weak positive correlation with Ca and the relationship with free Ca could not be determined. Nevertheless, the inclusion of K as a predictor of milk P is physiologically plausible and supported by its role in maintaining ionic balance in milk.
Next to lactose, Na and K are involved in the osmotic pressure of milk, which typically resembles the osmotic pressure of blood (Bijl et al., Reference Bijl, Van Valenberg, Huppertz and Van Hooijdonk2013). From the perspective of the osmotic pressure of milk, an inverse relationship between Na or K and lactose is expected. However, in the current study, both the Na content and K content of milk were not significantly correlated with the milk lactose content (P ≥ 0.117), most likely due to the low variation in the lactose content of milk. Visentin et al. (Reference Visentin, Niero, Berry, Costa, Cassandro, De Marchi and Penasa2019) reported a positive genetic correlation between milk fat and K contents. In the current dataset however, variation in milk K content explained only 1.8% of the variation in milk fat content (P = 0.334), indicating no meaningful phenotypic association. Furthermore, in contrast to the well-known role of Na and K in relation to the osmotic pressure of milk, there is no published information on their relationship, if any, with the formation of casein micelles or milk fat droplets. Thus, the observed positive relationship between the milk Na or K content and the milk P content is difficult to explain. Similarly, the current data implies a positive relationship between the Mn content of milk and milk P, but the underlying mechanistic principle, if any, is ambiguous.
Predictions of milk P content by multiple regression
Across the five indices used to evaluate the current multiple regression models, 4 equations emerged as the most promising for predicting the P content of milk, i.e. equations 1, 3, 5 and 7. Although equation 7 produced the greatest range of residuals, overall prediction errors were broadly similar for equations 3, 5 and 7. Despite a slight overprediction of milk P content, equation 1, which includes milk fat, Ca and K content as predictors, performed best overall with respect to R2adj, Mallows’ Cp, RMSPE, and BIC, with VIF (just) below the threshold of 3.3. Milk protein or casein content did not enter this equation, likely because of the high correlation between protein or casein content and Ca content.
In equation 7, three out of the five predictor variables (i.e., K, urea and non-casein protein) do not have a clear physiological relationship with milk P content (Klop et al., Reference Klop, Ellis, Blok, Brandsma, Bannink and Dijkstra2014) and the use of lactose in equation 7 is somewhat peculiar in view of its complete lack of relationship with the milk P content. Therefore, the use of equation 7 is not preferred to predict to milk P content. Likewise, equations 5 and 6 can be considered less opportune because in these equations only the milk fat or Ca content, respectively, can be causally linked to the milk P content. Following this reasoning, equations 3 and 4 are also not preferred (Klop et al., Reference Klop, Ellis, Blok, Brandsma, Bannink and Dijkstra2014). In view of its lowest prediction error and with 2 of the 3 variables, viz. fat and Ca, having a clear physiological relationship with milk P content, equation 1 is preferred to predict the milk P content. Although it is difficult to causally link the K content of milk with that of P, the outcome of the multiple regression analysis indicates that the milk K content is an important predictor variable in the current dataset. Indeed, in 7 out of the 8 models presented the milk K content was found to be significantly contributing to the explained variance in the milk P content.
The Klop et al. (Reference Klop, Ellis, Blok, Brandsma, Bannink and Dijkstra2014) model had lower prediction performance compared with the optimal equation derived in the current study. Since the contents of protein, fat and lactose of milk are routinely monitored in the Netherlands, Klop et al. (Reference Klop, Ellis, Blok, Brandsma, Bannink and Dijkstra2014) limited their candidate predictor variables to these components. However, in view of the results, the use of a broader spectrum of potential predictor variables to predict the milk P content appears justified, at least for scientific reasons. Conversely, using a broader rather than narrower set of predictor variables increases costs of chemical analyses and may amplify variation due to analytical errors associated with each predictor variable. In our preferred multiple regression model (i.e. equation 1), milk Ca and K contents are required to predict the P content of milk. Considering the potential for increased analytical errors and costs, it is suggested that a direct measurement of milk P content is preferable when calculating the mandatory amount of manure P.
Conclusions
Our study demonstrated significant seasonal variation in the P content of bovine milk. The highest milk P content was found during winter and the lowest during summer. The contents of fat, protein, casein, Ca and Mg in milk were all highly and positively correlated with milk P content and have physiological relevance. The preferred multiple regression equation to predict the milk P content included milk fat, Ca and K content as predictor variables and explained 80% of the variation (R 2adj) in the milk P content. Future studies are warranted to assess the performance of the currently preferred multiple regression equation using independent data on milk composition. While the present study enhances our understanding of physiological relationships involved in milk synthesis, direct measurement of milk P content is preferred for calculating the amount of manure P excreted into the environment.
Acknowledgments
This paper builds upon work done for the first author’s PhD thesis, entitled: Reappraisal of the Phosphorus Requirement of Lactating Dairy Cows. Utrecht University, the Netherlands, 2022.
Author contributions
KL and DJ: were responsible for data collection. KP and SJT: performed statistical analyses. KP, SJT and DJ: prepared the original draft of the manuscript. All authors reviewed and edited the draft. All authors have read and agreed to the final version of the manuscript.
Funding statement
Funding of this research by the Ministry of Agriculture, Nature and Food Quality (research theme ‘Manure and Environment’ BO-20-004-136, project KD-2017-121) is gratefully acknowledged. Open access funding provided by Utrecht University
Competing interests
The authors declare that there are no conflicts of interest.
Ethical standards
Not applicable.
Declaration of generative AI in scientific writing
Generative AI tools were not used at all to prepare the manuscript.
Open access
