Genetic parameters for milk urea concentration and milk traits in Polish Holstein-Friesian cows

Milk urea concentration (MU) used by dairy producers for management purposes can be affected by selection for milk traits. To assess this problem, genetic parameters for MU in Polish Holstein-Friesian cattle were estimated for the first three lactations. The genetic correlation of MU with milk production traits, lactose percentage, fat to protein ratio (FPR) and somatic cell score (SCS) were computed with two 5-trait random regression test-day models, separately for each lactation. Data used for estimation (159,044 daily observations) came from 50 randomly sampled herds. (Co)variance components were estimated with the Bayesian Gibbs sampling method. The coefficient of variation for MU in all three parities was high (40–41 %). Average daily heritabilities of MU were 0.22 for the first parity and 0.21 for the second and third lactations. Average genetic correlations for different days in milk in the first three lactations between MU and other traits varied. They were small and negative for protein percentage (from −0.24 to −0.11) and for SCS (from −0.14 to −0.09). The weakest genetic correlation between MU and fat percentage, and between MU and lactose percentage were observed (from −0.10 to 0.10). Negative average genetic correlation with the fat to protein ratio was observed only in the first lactation (−0.14). Genetic correlations with yield traits were positive and ranged from low to moderate for protein (from 0.09 to 0.33), fat (from 0.16 to 0.35) and milk yield (from 0.20 to 0.42). These results suggest that the selection on yield traits and SCS tends to increase MU slightly.


Introduction
Milk urea concentration (MU) reflects the efficiency of protein synthesis and provides information for dairy producers about the balance between crude protein and energy in the diet (Oltner and Wiktorsson 1983). In Poland, since 2000, MU is being recorded on a national scale for management purposes. A surplus of nitrogen in relation to energy in the rumen is a source of ammonia. Ammonia, because of its toxicity, is rapidly converted to urea in the liver. Urea molecules dissolve in body fluids, so milk urea nitrogen (MUN) is highly correlated with blood urea nitrogen and urinary nitrogen excretion (Butler et al. 1996;Broderick and Clayton 1997). High MUN level indicates inefficient protein conversion, which causes higher feeding costs and environmental pollution (Jonker et al. 1998). Furthermore, a negative energy balance causes health problems and decline in fertility (Butler and Smith 1989;Collard et al. 2000;Oikonomou et al. 2008).
The phenotypic variability of MUN is high, and the coefficient of variation exceeds 30 % (Stoop et al. 2007;Bastin et al. 2009;Yazgan et al. 2010; Hossein-Zadeh and Ardalan 2011a). MUN depends not only on nutrition but also on other factors, such as year and calendar month when the milk samples were collected, sample type (morning vs. evening), stage of lactation, age at calving, parity group and level of milk production (Godden et al. 2001;Wattiaux et al. 2005;Bastin et al. 2009). In preliminary studies, we found that the calendar month when the milk samples were collected and air temperature during that month were the major factors affecting MU (Rzewuska and Strabel 2013). MU was high from May to October, and the highest values were observed in summer months. However, MU increased with increasing temperature, regardless of the season. This indicates that air temperature should also be taken into account during winter months.
The variability of MUN is partly due to genetic factors and its heritability ranges from 0.13 up to 0.59 (Wood et al. 2003;Miglior et al. 2007;Bastin et al. 2009;Hossein-Zadeh and Ardalan 2011a). Hence, the level of MUN depends, to some extent, on the individual's genetic potential. The possibility of using MU as a predictor trait for indirect selection is of great importance (Wood et al. 2003). Some economically important traits are difficult to measure and their heritabilities are low. Therefore, the use of correlated traits may be beneficial for the selection process. Genetic relationships of MUN with fitness traits, such as reproductive performance and health, have been under investigation (Mitchell et al. 2005;König et al. 2008; Hossein-Zadeh and Ardalan 2011b; Mucha and Strandberg 2011). The estimation of genetic correlations with diseases is currently not possible, as this kind of data is not recorded in Poland. However, the somatic cell score (SCS) and the fat to protein ratio (FPR) are easily available from test-day data. SCS is an indicator of mastitis and FPR is an indicator of cow energy balance and describes its metabolic stability (Buttchereit et al. 2012). The value of FPR depends on the cow's diet and adaptation to metabolic stress at the beginning of lactation.
For good fit of the model, environmental factors affecting MU have to be estimated accurately. The most important environmental factor influencing MUN was herd-test-day (htd) effect, which is usually attributed to the high impact of herd management (Arunvipas et al. 2003;Wattiaux et al. 2005;Stoop et al. 2007;Bastin et al. 2009;Hossein-Zadeh and Ardalan 2011a). With a test-day model, modelling of the shape of individual and group lactation curves is also possible. For the genetic evaluation of dairy cattle in Poland, Legendre polynomials have been used to model both types of curves (Strabel et al. 2005). The shape of the lactation curve for MUN differs depending on the population. For example, some authors found that the MUN curve resembled a lactation curve for milk (Jonker et al. 1998;Stoop et al. 2007) and, in other studies, the curve for MUN had a mirror shape (Wood et al. 2003;Miglior et al. 2006;Hossein-Zadeh and Ardalan 2011a). In a preliminary study, we found that, for Polish Holstein-Friesian cows, MU was the lowest during the first month of lactation and its peak occurred in the fifth month of lactation (Rzewuska and Strabel 2013).
The direct effect of selection on MU would be a reduction of nitrogen emission from dairy farms. Possible selection on MU would have an impact on environmental protection but, currently, there is no legal consideration motivating farmers to select animals to minimise the emission of greenhouse gases. Interest in selection on MU is caused by its relationship with functional traits (health and fertility). The possible use of MU in the selection index requires knowledge of its relation with the traits under selection, such as production traits. Hence, the aim of this study was to use a test-day model to estimate genetic parameters of test-day MU, SCS, FPR and milk production traits for the first three parities of Polish Holstein-Friesian cows.

Materials and methods
MU was measured by infrared spectroscopy in test-day samples routinely collected from 2001 to 2009. MU was not converted to MUN because dairy farmers in Poland use phenotypic information about MU (mg/l). The raw data set was edited to eliminate records with unrealistic measures. The threshold of 0.1 % removed extreme values was applied for lactose percentage and protein percentage, while a threshold of 0.5 % was used for fat percentage and MU. For SCS, we used the threshold of 0.1 % for observations with the highest values and 0.5 % for observations with the lowest values, as the distribution of this trait is positively skewed. Means with standard deviations for all traits from the data set used in the analysis are shown in Table 1. Only records from 5 to 305 days in milk (DIM) in the first three parities, with data for all production traits, were subject to analysis. Records with missing birth data or non-physiological length of calving interval were removed. Cows with records from the second or third parities were required to have records from all previous lactations. Data were restricted to lactation with at least five records per lactation and herds with at least 50 first lactations. Although for cows in the third parity the effect of age at calving was not statistically significant, for cows in the first and second lactations, two classes of age at calving were created: 18-24 and over 24 months for cows in the first parity and 28-36 and over 36 months for cows in the second parity. The average air temperature in a given month was calculated with measurements from seven weather stations and divided into ten classes with 3-°C intervals, starting from −7.9°C. To account for interaction between calendar month and average air temperature in that month, classes of month-temperature were created. In each of the months, only 2-4 classes of temperature were observed. In total, 34 classes of month-temperature were defined. For variance components estimation, 50 herds were randomly sampled from all the available herds in the data set. The description of the data for consecutive lactations is presented in Table 1.
The test-day MU records were analysed using the following multiple-lactation model: where y ijklmnop is the milk urea concentration on test day o of cow m in parity p within herd-test-day effect i belonging to herd-year j and k class of age at calving, milked in l class of month-temperature; htd i is a random herd-test-day effect; HY is a fixed effect of the jth herd-year; AGE is a fixed effect of kth age at calving; b lnp are fixed regression coefficients on DIM specific to the l month-temperature class; z mnop are Legendre polynomials on DIM; a mnp are genetic random regression coefficients specific to the animal effect of cow m; p mnp are random regression coefficients for the permanent environmental effect; and e ijklmnop is the residual effect for each observation. Fixed lactation curves within subclasses of month-temperature were modelled using Legendre polynomials with four parameters. For random regressions, the same function of the second order was used.
Residual variance and variance of the htd effect were assumed to be constant across lactations. The variance structure of the random effects of the model was as follows: where I hp is the identity matrix of order equal to the number of htd levels in parity p; σ hp 2 the variance of the htd effect in parity p; A the additive genetic relationship matrix; G the additive genetic covariance matrix; I the identity matrix of order equal to the number of cows with records; P the permanent environmental covariance matrix; I p the identity matrix of order equal to the number of test-day records in lactation p; σ ep 2 the residual variance in parity p; while ⊗ and ∑ + are the direct matrix product and sum operator, respectively.
Due to a large number of analysed traits and lactations, genetic correlations between MU and other traits were computed with two 5-trait random regression test-day models. In the first model, we analysed MU, milk yield, protein yield, fat yield and SCS, while in the second model, MU was analysed with protein percentage, fat percentage, FPR and lactose percentage.
The model employed for each group of traits was analogous to that used for analysing test-day MU records; however, in those analyses, each lactation was analysed separately.
Covariance components were estimated by Bayesian methods with Gibbs sampling. Gibbs sampling was used to obtain the marginal posterior distribution for variance components of each random effect from the model from 100,000 samples, after discarding 20,000 samples as the burn-in period. Genetic (co)variances on each test day were calculated using the equation described by Jamrozik and Schaeffer (1997). Daily heritability was defined as the ratio of genetic variance to the sum of genetic, permanent environmental and residual variances at a given DIM. The proportion of variance attributable to the htd effect was calculated as the ratio of the htd effect variance and total variance. Results and discussion

Descriptive statistics
The average MU was 230.1, 230.3 and 226.2 mg/l in the first, second and third lactations, respectively (Table 1). Means for other milk components and milk yield in the first parity, presented in Table 1, are in agreement with the literature data (Wood et al. 2003;Miglior et al. 2007;Stoop et al. 2007;Yazgan et al. 2010; Hossein-Zadeh and Ardalan 2011a; Loker et al. 2012). The differences between MUN in subsequent lactations are small, which is consistent with some previous studies (Arunvipas et al. 2003;de Roos and de Jong 2006;Ramírez et al. 2008;Cao et al. 2010). However, contradictory results regarding the influence of parity were also found. In some analyses, the highest MU was observed in the first lactation (Jonker et al. 1998;Cao et al. 2010), in contrast to other reports, where the MUN level was the lowest in the first lactation (Arunvipas et al. 2003;Wood et al. 2003;Hojman et al. 2004;Miglior et al. 2006;Ramírez et al. 2008). The lowest MUN in the first lactation might be due to the animal's growth at this age, which might be associated with protein efficiency. Some other factors might also be important, e.g. a lack of a separate feeding group for primiparous cows (which have special requirements because of their growth) and herd hierarchy (which can be a reason for a lower dry matter intake). Also, lower body weight and lower level of milk production in the first lactation can have an impact on the lower MUN concentration (Jonker et al. 1998). Due to the large number of factors influencing MUN concentration, it is hard to explain why differences among average MU in subsequent lactations observed in our study were negligible. The coefficient of variation (CV) for MU in all three parities was high (40-41 %, Table 1), which confirmed the results from previous studies, where the CV ranged from 33 to 42 % (Stoop et al. 2007;Bastin et al. 2009; Hossein-Zadeh and Ardalan 2011a). By contrast, Yazgan et al. (2010), who used test-day records provided by one of the Polish commercial dairy farms, observed a greater variation (CV about 60 %).

Milk urea variances
Genetic and environmental variances were estimated over DIM during the whole lactation period by each parity. The pattern of variances for MU throughout the first lactation is presented in Fig. 1. On average, genetic variances were quite similar in all three lactations and ranged from 712.8 in the first lactation to 741.6 in the second. The mean daily permanent environmental variance was lower than genetic variance in the first lactation but similar in the second and third lactations. Genetic variances in the first lactation were flattened in the central part of lactation, with higher values in the first 2 months after calving, and increased in the last 2 months of lactation. In the second and third lactations, a slight decrease during lactation was observed, until the seventh month of lactation, when genetic variances suddenly increased. This increase, however, might be due to a lower number of observations, resulting in artefacts related to fitting the polynomials (Misztal et al. 2000). Higher estimates of permanent environmental variation were observed at the beginning of the first parity and at the beginning and end of the second and third lactations. In mid-lactation, the values increased, and in the second part of lactation, permanent environmental variations were greater than genetic variations in all parities. The studies that considered other pop- ulations of Holsteins reported genetic variances higher than permanent environmental variances over all lactations (de Roos and de Jong 2006;Miglior et al. 2007;Bastin et al. 2009), but Yazgan et al. (2010) found lower permanent environmental variation mainly in later lactations. In our study, residual variances were higher than both genetic and permanent environmental variances, throughout all lactations. That may be due to many effects influencing MU, which could not be included in the model because no information on sample type, dry matter intake or feeding groups was available. Measurement errors, usually relatively high for MU, might have also contributed to high residual variance, because infrared spectroscopy, which is used to collected data, is an indirect measure of MU. However, the htd effect explained the greatest part of variance, on average 53, 49 and 47 % in the first, second and third lactations, respectively. Other researchers reported that residual variance was greater than both genetic and permanent environmental variances (de Roos and de Jong 2006; Bastin et al. 2009). Only in the study conducted by Miglior et al. (2007) were residual variances lower than genetic variance, which proves the goodness of fit of the model. Wood et al. (2003) found a statistically significant htd effect, which was confirmed in later studies, where a large influence of the htd effect on the total variance for MUN was reported (de Roos and de Jong 2006; Stoop et al. 2007;Bastin et al. 2009). In the cited studies, the fraction of variance due to the htd effect ranged from 46 % up to 58 %, which is consistent with our results. Those authors suggested that management is an important factor influencing the MUN level, so the MUN concentration can be regulated by management adjustments at the herd level.

Heritabilities
The mean daily MU heritability was 0.22 for the first lactation and 0.21 for the second and third lactations. The heritabilities of daily MU are presented in Fig. 2. The lactation curves for all parities are similar in level and shape. The values for the second lactation were constant over the first 100 days, then decreased and reached the nadir (0.16) at 194 DIM, and, finally, gradually increased to 0.40 at 305 DIM.
The highest values of daily heritabilities were estimated at the end of all three lactations, due to an increase in genetic variance in the last 2 months of lactation. In the first and third lactations, the heritabilities after calving were higher than in the second lactation, but the changes over DIM in later lactations have a lower magnitude. Wood et al. (2003) observed an increased heritability of MUN with higher DIM, but other authors reported a curvilinear pattern, with a rapid increase until 35 DIM, followed by a "U" shape with a minimum at mid-lactation (Yazgan et al. 2010; Hossein-Zadeh and Ardalan 2011a). Hossein-Zadeh and Ardalan (2011a) reported the average heritability for MUN to range from 0.14 to 0.21 over the parities. Those authors estimated heritability after including the htd effect in the sum of variances. When, in our study, heritability was computed with the same method, the averages for all three lactations were equal and amounted to 0.11. The same tendency of lower value of daily heritability after including the htd variance, compared with the results not taking into account this variance (0.13 and 0.29, respectively), was observed by Bastin et al. (2009) in the first lactation. Heritabilities similar to our results were found by Mitchell et al. (2005). Those authors obtained different heritabilities for MUN measured with two methods-infrared (IR) and wet chemistry-and confirmed that they are genetically different traits. The genetic correlation between the two methods was only 0.38 and 0.23 for the first and second lactations, respectively. Higher heritability of IR MUN was explained by the specificity of the MUN measurement laboratory procedures. The levels of other milk components have an influence on MUN estimation because IR is an indirect method of measuring MUN. The lowest heritability was observed in studies conducted on relatively small data sets (Stoop et al. 2007;König et al. 2008;Yazgan et al. 2010) because the elimination of fixed effects is not very effective. MUN heritabilities for the first three parities (0.38 to 0.59), which were higher than those estimated in the present study, were noticed in Canadian investigations, where multiple-trait random regression test-day models were used (Wood et al. 2003;Miglior et al. 2007). However, our results confirm the genetic potential of individuals to produce MU and indicate that it should be taken into consideration when MU is used as an indicator of protein-energy balance in the diet. The mean heritabilities estimated for milk production traits ranged from 0.17 to 0.35, from 0.18 to 0.32 and from 0.20 to 0.30 for the first, second and third lactations, respectively. The highest values were observed for lactose percentage and protein percentage, and the lowest for fat yield ( Table 2). The heritabilities for SCS were lower than the heritabilities for the production traits, and were equal to 0.13, 0.14 and 0.16 for the first three parities, respectively.
The heritabilities for production traits and SCS were within the ranges reported for Polish Holstein-Friesians (Strabel et al. 2004;Strabel and Jamrozik 2006;Yazgan et al. 2010;Jesiołkiewicz et al. 2011;Rzewuska et al. 2011;Ptak et al. 2012) and all-except the heritability for SCS-were lower than those obtained for Canadian and Dutch populations (Wood et al. 2003;Miglior et al. 2007;Stoop et al. 2007). A similar relation between heritabilities estimated for the production traits of the Polish and Canadian populations, with higher values for the latter, was indicated in an earlier study (Strabel et al. 2004). The results of that study implied that differences in heritabilities are a consequence of differences in data structure and environmental conditions. In all lactations, the heritabilities for percentage traits (fat and protein) were higher than the heritabilities for yield traits. This is consistent with the tendency reported by Miglior et al. (2007). The average heritability for FPR in the first lactation in our study was 0.40 and ranged from 0.23 to 0.51. Our estimation was similar to the results for other populations of Holsteins: 0.42 for Canadian (Loker et al. 2012) and 0.30 for German cattle (Buttchereit et al. 2012). Lower heritabilities, ranging from 0.13 to 0.25, were reported for Nordic Red (Negussie et al. 2013), which can be explained by between-breed differences, type of model and effects included in the model.

Genetic correlations
The average genetic correlation between MU at the same DIM presented in Fig. 3 were 0.81 between the first and second lactations and 0.83 between the first and second as well as between the second and third lactations, and were in the range of the correlations estimated in previous studies, which was between 0.74 and 0.97 (Miglior et al. 2007;Hossein-Zadeh and Ardalan 2011b). In a Canadian population, genetic variances were high (Miglior et al. 2007), which can be a reason for the high correlation between lactations.  Hossein-Zadeh and Ardalan (2011b) suggested that MUN at different parities were genetically similar traits. The results of the current study confirm this statement. However, the selection based only on the data from the first lactation could limit the response in traits of interest for later lactations. We found the highest correlation between the most distant lactations, the first and the third, which was in contradiction to the findings of other authors (Miglior et al. 2007; Hossein-Zadeh and Ardalan 2011b). For the first and second parities, the correlation was the lowest after calving and was increasing up to the seventh month of lactation. In contrast, for other parities, the lowest values were observed at the end of lactation, with rather similar correlations in the central part of lactation. The highest correlations were found in different stages of lactation for different pairs of parities (Fig. 3). The maximum values were observed at 218 DIM (0.95), 78 DIM (0.90) and 189 DIM (0.92) between the first and second, the first and third, and the second and third parities, respectively. Genetic correlations between MU at the end of a given lactation and MU at the beginning of the subsequent one were close to zero (results not shown).
Genetic correlations between MU at selected DIM and at the remaining part of lactation in the first parity are shown in Fig. 4. The highest genetic correlations within lactation were observed between MU on adjacent test days at the same stage of lactation. Only in mid-lactation were high values for more distant days found. The days at the end of lactation were weakly correlated with the days in the third month of the first parity. In other parities, correlations between these stages of lactation were negative (e.g. −0.19 between days 70 and 305 of the second lactation). However, for management purposes, the first part of lactation is the most important period due to problems with meeting energy requirements resulting in a frequent occurrence of metabolic disorders. No information about MU genetic correlations within lactation is available in the literature. In this study, the relatively low values of genetic correlations between MU within lactation and between MU in subsequent lactations suggest that possible selection should be based on records from more than one, full lactation.
Daily genetic correlations between MU and production traits ranged from low to moderate (Table 3)  In the first parity, genetic correlation between MU and yield traits increased with lactation trajectory and reached the highest value at the end of lactation. The average genetic correlation in the second parity was lower than those for the first and third lactations. The weakest relationships with fat percentage and lactose percentage were observed (from −0.10 to 0.10). Genetic correlations between MU and protein percentage (from −0.24 to −0.11) were small and negative. The strongest negative correlation was observed in the first parity between MU and protein percentage, with the highest value at the beginning of lactation (−0.48 for day 5). This may be related to the efficiency of protein metabolism and nitrogen utilization, especially during the energy deficit phase of lactation. Also, the average genetic correlations among MU and SCS were negative: −0.14, −0.16 and −0.09 in the first, second and third parities, respectively. In the first lactation, higher negative correlation in the first half of lactation (−0.32 at day 5) and in the last part of lactation (−0.43 at day 305) was observed. In the second parity, the correlation between MU and SCS was stable in mid-lactation, with higher values in the last part of lactation. This indicates that the selection for mastitis resistance based on SCS may lead to an increase in MU. Genetic correlation between MU and FPR in the first lactation is negative and stable over DIM until the sixth month of lactation. This could be explained by changes in the metabolism. In the early lactation, the cow is in a negative energy balance, which results in adipose tissue mobilisation. This leads to an increase in FPR because the fat percentage in milk increased. Prolonged energy deficiency results in hepatic lipidosis. That makes MU synthesis in the liver difficult and results in a decreased MU in milk. However, the genetic correlation between MU and FPR in the next parities showed a completely different pattern. At the beginning of lactation, the correlation is close to zero and increases Genetic correlation DIM Fig. 6 Genetic correlations between milk urea concentration (MU) and fat to protein ratio (♦), protein percentage (■), fat percentage (▲) and lactose percentage (X) within days in milk (DIM) in the first lactation thereafter, reaching the maximum values in mid-lactation (0.24 at day 154 and 0.41 at day 106 of the second and third parities, respectively). The interpretation of correlation with FPR is difficult because that trait has an intermediate optimum. Larger changes in correlations between MU and other milk traits over DIM in the second and third parities were observed. Partially, that may be a result of using fewer data for this part of lactation, but it may also confirm that MU in the subsequent lactation is not the same trait.
In most previous studies, a positive correlation between MUN and milk yield was found, and only Miglior et al. (2007) did not confirm this rule. Contrary to our results, positive genetic correlations of MU with protein percentage, fat percentage (Miglior et al. 2007;Stoop et al. 2007; Hossein-Zadeh and Ardalan 2011a) and FPR (Loker et al. 2012) were reported. At the same time, Miglior et al. (2007) found a negative genetic correlation with SCS (−0.19), which was consistent with our results.

Conclusions
The estimated heritabilities of milk urea concentration (MU) in the milk of Polish Holstein-Friesians in the first three parities were moderate. Genetic correlations between MU within lactation and between MU in the same days of subsequent lactations were low. The positive genetic correlations with protein yield and fat yield and negative genetic correlation with the somatic cell score (SCS) confirm that MU may slightly increase due to selection for production traits and udder health. MU is genetically correlated with the fat to protein ratio (FPR) but the direction of the relationship is inconsistent: negative in the first lactation and positive in the second and third lactations. The ranges of heritabilities of MU suggest that they can be used in a genetic evaluation for that population. Therefore, further studies are needed in order to assess the genetic correlations with fertility traits, which are a low heritable but economically important trait.