Causal relationships between milk quality and coagulation properties in Italian Holstein-Friesian dairy cattle

Background Recently, selection for milk technological traits was initiated in the Italian dairy cattle industry based on direct measures of milk coagulation properties (MCP) such as rennet coagulation time (RCT) and curd firmness 30 min after rennet addition (a30) and on some traditional milk quality traits that are used as predictors, such as somatic cell score (SCS) and casein percentage (CAS). The aim of this study was to shed light on the causal relationships between traditional milk quality traits and MCP. Different structural equation models that included causal effects of SCS and CAS on RCT and a30 and of RCT on a30 were implemented in a Bayesian framework. Results Our results indicate a non-zero magnitude of the causal relationships between the traits studied. Causal effects of SCS and CAS on RCT and a30 were observed, which suggests that the relationship between milk coagulation ability and traditional milk quality traits depends more on phenotypic causal pathways than directly on common genetic influence. While RCT does not seem to be largely controlled by SCS and CAS, some of the variation in a30 depends on the phenotypes of these traits. However, a30 depends heavily on coagulation time. Our results also indicate that, when direct effects of SCS, CAS and RCT are considered simultaneously, most of the overall genetic variability of a30 is mediated by other traits. Conclusions This study suggests that selection for RCT and a30 should not be performed on correlated traits such as SCS or CAS but on direct measures because the ability of milk to coagulate is improved through the causal effect that the former play on the latter, rather than from a common source of genetic variation. Breaking the causal link (e.g. standardizing SCS or CAS before the milk is processed into cheese) would reduce the impact of the improvement due to selective breeding. Since a30 depends heavily on RCT, the relative emphasis that is put on this trait should be reconsidered and weighted for the fact that the pure measure of a30 almost double-counts RCT. Electronic supplementary material The online version of this article (doi:10.1186/s12711-015-0123-7) contains supplementary material, which is available to authorized users.

Conclusions: This study suggests that selection for RCT and a 30 should not be performed on correlated traits such as SCS or CAS but on direct measures because the ability of milk to coagulate is improved through the causal effect that the former play on the latter, rather than from a common source of genetic variation. Breaking the causal link (e.g. standardizing SCS or CAS before the milk is processed into cheese) would reduce the impact of the improvement due to selective breeding. Since a 30 depends heavily on RCT, the relative emphasis that is put on this trait should be reconsidered and weighted for the fact that the pure measure of a 30 almost double-counts RCT.

Background
In recent years, increasing efforts have been made to enhance efficiency in the Italian dairy industry and dairy cattle breeding organizations have started selecting for a wide range of novel traits. Milk coagulation properties (MCP) have been included in the data recording system and breeding values are routinely produced for Italian Holstein bulls [1]. Milk coagulation properties, namely rennet coagulation time (RCT) and curd firmness after 30 min from rennet addition (a 30 ), have been shown to be good predictors of milk technological quality and cheese yield [2][3][4], which are key factors in dairy industries where most of the milk produced is processed into cheese. In particular, a 30 is the trait that has the strongest impact on Grana Padano cheese processing [4].
Generally, selection for RCT and a 30 is based either on correlated traits such as somatic cell score (SCS), fat, protein and casein percentages [5][6][7][8] or on direct measures of RCT and a 30 [1,9]. These two traits can vary in terms of curd firmness at different time points in the milk coagulation process and depend heavily on each other. This is inherent to the test used (see Bittante [10] and Bittante et al. [11] for a review of current knowledge), i.e. RCT (in min) measures the amount of time between rennet addition and the beginning of the coagulation process, whereas a 30 measures curd firmness 30 min after rennet addition. The longer the milk takes to start coagulating, the softer the curd will be at the end of the test, and vice versa. Somatic cell score and milk casein percentage (CAS) are considered to affect RCT and a 30 [12][13][14] and are correlated at the genetic level [5,15]. Pretto et al. [7] suggested that the genetic correlation that exists between SCS and CAS could be used as a predictor in breeding programs that focus on improving MCP.
The overall genetic effects that influence MCP are probably distributed into multiple causal paths: on the one hand, some genes may affect MCP directly, while, on the other hand, some genes may affect other milk quality parameters, which in turn affect the ability of milk to coagulate. Alternatively, a causal path that involves MCP may exist. For instance, a strong association between a 30 and RCT could support a causal hypothesis that variability in a 30 is mostly explained by the influence of RCT, while there is no strong direct genetic effect on a 30 (i.e. they are independent). In other words, some genes may not strongly and directly affect both RCT and a 30 , but only RCT.
As discussed by [16], in the classical genetic evaluation scenario, breeding values of candidate individuals are predicted by fitting multiple trait models (MTM), which neglect the causal network that influences phenotypic traits. Structural equation models (SEM) [17][18][19][20][21][22][23][24][25][26][27][28][29][30][31] can help "dissect" the overall genetic effects expressed by MTM into distinct sources of genetic variation, by separating the common sources of variation that affect directly two or more traits in the system (e.g. the genetic correlation between SCS and RCT) from the causal effect that one phenotypic trait plays on the other (e.g. the causal effect of SCS on RCT). In addition, using nonintervened data (such as field data routinely collected for genetic evaluations), SEM would be able to predict genetic effects for scenarios for which interventions on the phenotype are performed [16]. For example, let us consider a scenario in which the goal is to predict the individuals' genetic effects on RCT or a 30 when milk quality traits are physically controlled (e.g., by filtering somatic cell load [32,33] or standardizing casein percentage [34,35]). Such a scenario will take only the individuals' genetic effects on RCT and a 30 into account, since the genetic influence mediated by SCS and CAS is blocked. Alternative approaches to assess the impact of SCS and CAS on MCP require that additional experimental records, this time under the given intervention, and then genetic effects based on these data can be predicted. Therefore, the genetic parameters and breeding values for a 30 that are estimated with SEM result in interpretations that differ from those obtained with MTM. The breeding value of a given cow that is estimated with SEM for a 30 indicates its genetic merit in terms of firmness of the milk produced that is not mediated by SCS and CAS and that cannot be obtained with standard MTM. Moreover, if indirect effects play a major role in the variation of a 30 and if the traits through which the effect is mediated are in the selection index, it might not be necessary to consider a 30 since that would represent redundant information.
The aim of this study was to: (1) infer the magnitude of the causal effects of two traditional milk quality traits (SCS and CAS) on MCP, (2) estimate the causal effect of RCT on a 30 and the genetic variation in the latter when the causal effect of the former is removed, and (3) estimate genetic and phenotypic variation of a 30 by taking the causal effect of all other traits into account, i.e. by assessing the magnitude of the variance components that would hold if some milk quality parameters are controlled.

Data collection and editing procedure
Routine assessment of MCP via mid-infrared spectroscopy began in September 2011 in the Veneto region of Italy [1,36]. Approximately 25 000 cows are currently under monthly routine control. A panel of traits is routinely assessed with Milko-Scan FT6000 (Foss Electric A/S, Hillerød, Denmark), including fat, protein and casein percentages, milk coagulation properties (such as RCT, curd firming time and a 30 ) and fatty acid profile of milk. Milk somatic cell count (SCC) is determined with Cell Fossomatic 250.
For this study, we extracted data from the same dataset as in Tiezzi et al. [36]. Traditional milk quality parameters were chosen, i.e. SCS (as log-transformation of SCC) and CAS, while we used RCT and a 30 as measures of MCP. We retained only the records from early-lactation (5 to 125 days in milk) on first-lactation cows in order to avoid accumulation of carry-over effects of deteriorated milk quality on coagulation properties (e.g. an identical decrease in SCS may have a different impact on RCT in early and late lactation because of the accumulation of the effect of SCS over lactation, such that late lactation RCT and a 30 may be affected by early lactation SCS, late lactation SCS and their interaction). Therefore, data editing was similar to that in [36], except that only the records from earlylactation (5 to 125 days in milk) first-parity cows were considered. For statistical analysis, 8783 records collected on 3266 first-lactation Italian Holstein cows across the period from January to December 2012 were used for statistical analyses. Cows were sired by 128 AI bulls and reared in 309 herds.

Statistical analysis
We fitted three SEM and a single MTM. The baseline MTM (M0) was as follows: where the index 'i' indicates correspondence to the i th trait, y 1, , y 2 , y 3 and y 4 are the vectors reporting the four traits (SCS, CAS, RCT and a 30 , considered in this order), X and b i are the incidence matrix and the corresponding vector of fixed effects (intercept and four classes of stage of lactation, namely 5 to 34, 35 to 64, 65 to 94 and 95 to 125 days in milk), Z h and h i are the incidence matrix and corresponding vector of herd random effect (309 levels), Z p and p i are the incidence matrix and vector of cow permanent environmental random effect (3266 levels), Z s and s i are the incidence matrix and vector of sire random additive genetic effect (128 sires, 1254 total individuals in the sire-MGS pedigree), and e i are random residuals.
Different causal structures with varying complexity in terms of number of causal connections were assigned to each model. Causal connections are represented in the same manner as in Wu et al. [26]. For instance, λ yx indicates a causal effect of x on y, traits are coded as follows: 1 for SCS, 2 for CAS, 3 for RCT and 4 for a 30 . Model 1 (M1, Figure 1) takes the effect of SCS on RCT and a30 (λ 31 and λ 41 , respectively) and the effect of CAS on RCT and a 30 (λ 32 and λ 42 , respectively) into account. This model is represented as follows: In model 2 (M2, Figure 2), only the effect of RCT on a 30 (λ 43 ) was considered, In model 3 (M3, Figure 3), effects of SCS, CAS and RCT on a 30 were considered (λ 41, λ 42 and λ 43 , respectively).
where y 1 , y 2 , y 3 , y 4 , X, b i , Z h , h i , Z p , p i , Z s , s i and e i are defined as for the MTM. Furthermore, two additional models derived from model M1 were used to avoid possible confounding between the effects of SCS and CAS on RCT and a 30 : one took only the effect of SCS on RCT and a 30 into account, while the other took only the effect of CAS on RCT and a 30 into account. However, since the estimated values of the causal effects were similar with model M1 and the two derived models, results from these models are not presented.
Analyses were implemented in a Bayesian framework using the software SIRBAYES [24,26]. For structural coefficients, a multivariate normal prior distribution was assumed as N(1λ 0 , Iτ 2 ), where hyperparameters were  λ 0 = 0 and τ 2 = 10 000. For fixed effects, the prior distribution was normal, with mean 0 and variance 10 000. Prior distributions for sire, cow permanent environmental and herd effects were multivariate normal with the following covariance structures: for sire effect s~N(0, G ⊗ A) where A was the numerator relationship matrix and G is the sire effect covariance matrix, and for cow and herd effects and residuals, it was assumed that p~N(0, P ⊗ I), h~N(0, H ⊗ I), e~N(0, R ⊗ I) where I is an identity matrix and H, P and R are the respective covariance matrices. Prior distributions for covariance matrices G, H, and P were independent inverse-Wishart invWish(ν, S), where ν are the number of degrees of freedom and S is the scale. Prior distribution for R was an independent inverse-Wishart invWish(ν, S) only for the MTM, while an inverted chi-square invChisq(ν, S) was used for the SEM since R was forced to be diagonal. For all priors, the number of degrees of freedom (ν) was set to 6. In this study, R was assumed as diagonal for SEM, i.e., all residual covariances were constrained to 0. This assumption is required to identify structural coefficients. The meaning of this commonly adopted parametric constraint conflicts with the quantitative genetics that underlie sire MTM, i.e. alleles that are inherited from the dam may have associated effects in more than one trait, and they are expected to be absorbed by the residual covariance. In fact, this theoretical contradiction in SEM that are based on sire models was largely neglected in previous studies [19][20][21][22][23][24][25][26][27][28][29][30][31]. However, the possible confounding is not expected to be important here because of the structure of the data and the model: a model that takes into consideration dam effect would not take this effect into account well since most cows do not share the same dam in the dataset.
Structural coefficients were sampled using the Metropolis-Hastings algorithm, and the remaining parameters were sampled using Gibbs sampling [28]. For each model, 120 000 iterations were run, discarding the first 20 000 as burn-in and retaining one every 10 samples for inferences. Posterior means and 95% highest probability density intervals were calculated on the remaining 10 000 samples. Convergence was assessed by visual inspection of the trace and running mean plots and estimates of autocorrelation and effective samples size were obtained using the 'coda' package [37] in R (http://cran.r-project.org).
Since the interpretation of the parameters that are estimated with the SEM (contained in G, P, H, and R) differs from that of the analogous parameters with a MTM [16], further transformation is required to be able to compare (co)dispersion of overall random effects between the four models fitted. For each model, transformation for the estimated covariance matrices to the MTM scale was performed as: where G, P, H, R and Λ are defined as above. Genetic and phenotypic correlations were calculated in the usual way from the (co)variance components in G*, P*, H* and R*. Heritability (h 2 ) was computed as: where σ 2 s is the sire additive genetic variance, σ 2 p is the cow permanent environmental variance, σ 2 h is the herd environmental variance and σ 2 e is residual variance. These variance components are obtained from G*, P*, H* and R*. For easier interpretation, posterior means of causal effects were transformed to standard deviation units by applying the formula yx is the transformed value, λ yx is the posterior mean of the causal effect of x on y, sd (x) is the standard deviation of the independent variable and sd (y) is the standard deviation of the dependent variable.
For a 30 , we computed the difference in sire additive genetic and phenotypic variance (sum of sire additive genetic, cow permanent environmental, herd and residual variances) between model M0 and each of the considered models under the causal effect (RCT under model M1 and a 30 from all models), expressed as the relative difference with model M0 (difference between the variance components divided by the respective variance component of model M0 and scaled to 100). In addition, we computed the heritability from the SEM variance components (i.e., from G, P, H and R) according to the formula above.

Descriptive statistics and observed correlations
Descriptive statistics and observed correlations are in Table 1. Means (SD) for SCS, CAS, RCT and a 30 were equal to 2.35 (1.66), 2.46 (0.23), 18.9 (3.80) and 23.0 (8.53), respectively. The data originated from the same dataset as in Tiezzi et al. [36] but observations were restricted to early-lactation first-parity cows. Descriptive statistics were in partial agreement with the previous study. Summer et al. [38] also found a significant increase in casein percentage from early to late lactation in Italian Friesian cows. Correlations of RCT with SCS and CAS were equal to 0.087 and -0.021, respectively, while correlations of a 30 with the other traits were -0.849 for RCT, -0.107 for SCS and 0.346 for CAS.

Heritabilities, genetic and phenotypic correlations
Heritabilities (on the diagonal) and genetic and phenotypic correlations (above and below the diagonal, respectively) estimated with each model are in Table 2. These parameter estimators should be interpreted as the standard parameters obtained with a MTM.
Heritabilities for all traits considered were consistent across models: SCS ranged from 0.021 in model M1 to 0.030 in model M0, CAS ranged from 0.141 in M1 to 0.157 in M2, RCT from 0.112 in M3 to 0.167 in M0 and a 30 ranged from 0.139 in M2 to 0.187 in M0. Also, genetic and phenotypic correlations did not vary significantly across models, since in most cases, the posterior mean of one model fell within the 95% HPD (highest posterior density interval) intervals of the other models. Genetic correlations that involved RCT were almost null with SCS (posterior means ranging from -0.081 to 0.081), very moderate and negative with CAS (-0.241 to -0.155), and strong and negative with a 30 (-0.933 to -0.883). The trait a 30 presented weak negative correlations with SCS (-0.149 to -0.072) and moderate positive correlations with CAS (0.374 to 0.572). Genetic correlation between SCS and CAS was null, ranging from -0.180 to 0.059 across models. Phenotypic and genetic correlations were similar in direction and magnitude: RCT was moderately and positively correlated with SCS (0.147 to 0.182) and strongly and negatively correlated with a 30 (-0.851 to -0.58), while correlation with CAS was null (-0.125 to -0.021). a 30 was weakly and negatively correlated with SCS (-0.192 to -0.151) and moderately and positively correlated with CAS (0.211 to 0.494), while correlation between SCS and CAS was null (0.021 to 0.042).
Estimates of heritabilities obtained with the MTM were lower than those reported by Tiezzi et al. [36], who considered whole lactations up to the ninth parity. In our study, restricting the dataset to early-lactation firstparity cows led to decreased heritabilities for all traits: from 0.093 to 0.030 for SCS, from 0.283 to 0.157 for CAS, from 0.210 to 0.167 for RCT and from 0.238 to 0.187 for a 30 . To the best of our knowledge, there are no studies on the heterogeneity of variance components across lactation and parities for milk coagulation properties; however, Muir et al. [39] found a lower heritability for SCS in the first lactation than in later lactations for Italian Holsteins (0.165, 0.211 and 0.252 for first, second and third lactations, respectively), while Odegard et al. [40] reported a heritability for SCS less than 0.08 at the beginning of lactation and a value of 0.10 in late lactation, although they used Norwegian Red cattle data. Estimates of genetic and phenotypic correlations are in agreement with Tiezzi et al. [36], therefore stage and number of lactation do not appear to affect correlations.
As mentioned above, we found no significant differences in estimates of genetic and phenotypic correlations between MTM and SEM. However, exceptions were observed between CAS and a 30 , i.e., models M0 and M2 led to lower values (0.291 and 0.211, respectively) while models M1 and M3 led to the highest values (0.494 and 0.463, respectively). Including causal effects between CAS and a 30 increased the correlation estimates. This was observed in several other studies: (1) Konig et al. [29] who analyzed correlations between milk yield and claw disorders in German Holstein, found that SEM resulted in lower genetic correlations compared to MTM; (2) de los Campos et al. [21] reported that, in dairy goats, genetic correlations between milk yield and SCS differed between MTM and SEM; and (3) Wu et al. [26] showed that including causal effects modified the genetic correlations between milk yield and SCS, while heritabilities varied little across models.

Causal effects on RCT and a 30
The causal effects estimated with the three SEM and their transformations to the scale of standard deviation units for both traits are in Table 3. Model M1 took the effects of traditional milk quality parameters (SCS and CAS) on milk coagulation measures (RCT and a 30 ) into account. According to the posterior mean of this parameter, an increase of 1 unit in SCS (e.g. from 2.00 to 3.00)  The effects inferred from model M3 were weaker than those from other models, although they agreed in sign. The posterior means of the effect of SCS, CAS and RCT on a 30 were equal to -0.267 (95% HPD intervals: -0.327 to 0.207, -0.0520 SD units), 12.845 (95% HPD intervals: 12.443 to 13.232, 0.3465 SD units) and -1.792 (95% HPD intervals: -1.819 to -1.764, -0.7983 SD units), respectively.
Overall, increasing milk quality (i.e. lower SCS and higher CAS) led to better milk coagulation ability (i.e.   lower RCT and greater a 30 ), although the traits used and the strength of the relationship varied across studies.
Reviewing the effects of SCS on cheese process and quality, Le Marechal et al. [14] found that, in most of the studies, a high SCS was associated with extended rennet clotting time and lower curd firmness, and Mazal et al. [41] showed that curd firmness of milk decreased as SCS decreased from~800 000 to 170 000 cells per mL). Politis and Ng-Kwai-Hang [13] reported a regression coefficient of a 30 on CAS of 12.92, which is very close to the value found here (12.845) although the methodologies used differ and the required parametric interpretation does not allow straightforward comparisons. However, in a study conducted by Grandison and Ford [12], a correlation of 0.807 was estimated between SCC and coagulum strength. Causal relationships between somatic cell score, casein percentage and milk coagulation properties

Impact of causal effects on curd firmness
Our results suggest that even under a hypothetical scenario in which SCS and CAS are maintained at constant values by external intervention and their influence on variability is nullified, the firming process is expected to vary, with part of the variability being attributable to the additive genetic component. In an experimental scenario in which CAS and SCS are standardized across samples, we would still find some additive genetic variation in RCT and a 30 . However, the causal pattern appears to differ between RCT and a 30 . For RCT, both direct sire additive genetic and phenotypic variances obtained with model M1 were essentially equal to those obtained with the MTM scenario (results reported in Table S2 [See Additional file 1: Table S2]), which is what would be expected if variation in RCT was weakly mediated by traditional milk quality parameters. In fact, RCT showed low phenotypic correlations with SCS and CAS, which was translated as weak inferred causal effects, especially from CAS. Weak statistical dependences generally suggest the absence of strong causal effects.
Curd firmness can be considered as the most pertinent coagulation measure that can account for cheese yield under certain processing conditions [7]. Assuming model M1, 18.1% of the additive genetic variance and 14.4% of the phenotypic variance can be explained by causal effects of SCS and CAS, which leaves a large proportion of the variation explained by other sources. The moderate correlations between a 30 and CAS inferred with MTM (an observed phenotypic correlation of 0.346 (Table 1) and a genetic correlation of 0.374 (Table 2)) can be considered to result from a causal effect rather than from a common source of variation i.e., the difference in variability for a 30 between real and causality-free scenarios was noticeable for all variance components but relatively larger for the additive genetic variance. More evidence is provided by the fact that the square of the phenotypic and genetic correlations between a 30 and CAS is close to the drop in variance for a 30 under model M2 (Table 4). This scenario suggests that the overall genetic association, which is observed between CAS and a 30 with MTM, is mostly due to a phenotypic causal effect of the former on the latter, rather than a pleiotropic effect on both traits directly. This hypothesis is supported by the null values estimated for genetic covariances under model M1 [See Additional file 1: Table S1].
Here, we found causal dependencies between RCT and a 30 . Our results confirm that curd firmness is intrinsically connected with RTC and depends causally on it [10], which suggests that it virtually cannot vary when the latter is held constant. In this study, model M2 accounted for this dependency as an effect of RCT on a 30 . This type of model makes it possible to investigate how the system would react if the samples were physically set to have the same coagulation time, i.e. removing variability due to the effect of starting coagulation time. Variance components and heritabilities for scenarios that involve such interventions can be inferred from information provided by Λ, G, P, H and R pertaining to model M2 before transformation to the standard MTM, and are in Table 4. It should be noted that the heritability of a 30 (0.097) is low, but non-null, and the posterior mean (95% HDP intervals) of its genetic covariance with CAS is 0.030 (0.011; 0.052), similar to the corresponding genetic covariance of 0.052 (0.001; 0.105) obtained with model M0 [See Additional file 1: Table S1]. These results indicate that the association between the genetic effect of CAS and the genetic effect of a 30 that is not mediated by RCT is still present. Thus, under model M2, curd firmness can be considered as a trait by itself and should not be ignored in selection indices, despite its low heritability (0.097), which will constrain genetic progress.
The almost complete loss of variance for a 30 with model M3 suggests that the causal effects of SCS, CAS and RCT absorb a large part of the phenotypic variation of a 30 . Subsequently, independent sources of environmental variation are scarce and direct additive genetic variance is negligible. Under the causal assumptions applied here, we show that the observed genetic associations between a 30 and the other traits are due to phenotypic causal effects rather than to common sources of genetic variation, i.e. genes with pleiotropic effects. Under these circumstances, the importance of a 30 in selection indices should be downplayed if SCS, CAS and RCT are already taken into account.

Conclusions
This study inferred causal relationships between two traditional milk quality measures (somatic cell score and casein percentage) and milk coagulation properties (rennet coagulation time and curd firmness). Results from this study suggest that the additive genetic variance of milk coagulation properties does not depend on traits such as somatic cell score or casein percentage. If selection is performed on these traits, coagulation properties will be improved only indirectly and to a small extent. This means that selection for milk coagulation properties should be performed on their direct measures, and cannot rely entirely on correlated traits since external interventions on the correlated traits may break down the causal path. In addition, including both rennet coagulation time and curd firmness in genetic evaluations appears redundant, considering that the latter depends largely on the former. If rennet coagulation time is the only selection objective, specific models that include right-censoring would probably better suit this purpose, because 5 to 10% samples of milk do not start coagulation within the first 30 min after rennet addition. Otherwise, it is necessary to demonstrate that curd firmness can have an impact on cheese yield and quality of products even if coagulation time is physically set to a constant value. This trait could represent an additional selection objective, but further research is needed.

Additional file
Additional file 1: Table S1. Estimates 1 of co-variance components for sire additive genetic effect, cow permanent environmental effect, herd effect and residual as estimated with different models 2 . 1 Estimates are the means (lower and upper bound of the 95% HPD interval) of the marginal posterior distributions. 2 The models differ in the structural coefficients considered: M0 is the standard multiple trait model; in M1 are considered the causal effects of both SCS and CAS on RCT and a 30 ; in M2 is considered the causal effects of RCT on a 30 ; in M3 the causal effects of SCS, CAS and RCT on a 30 are considered.