An R Shiny WebApp was developed to display results from this research and it is freely available at
https://alcantara.shinyapps.io/prodollar. In addition to the tables and figures presented herein, this WebApp features extra information, such as the full data set, density plots (number of daughters, RBV, and reliabilities), result summaries from regression analysis, AIC values from SBR, full results from PCA (all eigenvalues and eigenvectors), and extra plots from PCA (explained variance, variable contribution, quality of representation of traits, and interactive biplot).
3.1 Multiple polynomial regression analysis
Estimates of changes in Pro$ were expressed as net estimate values (
Table 3), which are the regression coefficients for the linear terms or the addition of the linear and quadratic terms for intermediate-optimum traits. The full model explained 72.5% of the variance in Pro$, with heel depth (HD), dairy capacity (DC), and body depth (BD) having the largest effects, whereas bone quality (BQ), udder depth (UD), and rump angle (RA) were associated with the smallest amount of variation in Pro$ (
Table 3).
Positive effects of DC on Pro$ indicate that a more angular and capacious cow would have higher profit compared to the average dairy cow in Canada. An increase in one unit of RBV for DC is estimated to increase Pro$ by $59.00. A possible explanation would be due to the contribution that DC has on the longevity and production of dairy cows, since accumulated profit is the core of Pro$. A strong positive relationship between DC and functional longevity has been reported in Canadian Holsteins (
Sewalem et al. 2004), where nonangular cows (score 1) were 2.47 times more likely to be culled than intermediate angular animals (score 5). Similar positive correlations are also observed in Iranian (
Dadpasand et al. 2008), Czech (
Zavadilová et al. 2011), and American Holsteins (
Caraviello et al. 2004), showing lower relative risk of culling for more angular cows.
Previous research found high genetic correlations between DC and milk yield to 240 days of lactation (0.48) in an Irish Holstein-Friesian population (
Berry et al. 2004). In a study with Czech Holsteins (
Zink et al. 2014), genetic correlations between DC and yield traits were moderate to high (0.32, 0.34, and 0.42 for milk (MY), protein (PY), and fat yield (FY), respectively). In an Italian Brown Swiss population study (
Samoré et al. 2010), moderate genetic correlations for DC and production traits were found (0.36, 0.39, and 0.23 for MY, FY, and PY, respectively).
Cows with well-sprung ribs tend to stay longer in the herd since such conformation is thought to facilitate the cow's ability to process large volumes of roughage while sustaining high production levels (
Atkins et al. 2008). The positive correlation between DC and Pro$ seems to play an important role in balancing the negative effects that deeper bodies have on Pro$, whereby high feed-consuming cows with deeper bodies might be able to offset some of the feed costs by being more efficient in processing such large amounts of feed towards milk production.
Body depth appeared to contribute the most to lowering Pro$ index value, decreasing Pro$ by −$61.95 for every increase in one unit of RBV for BD, meaning that deeper animals are associated with a lower Pro$ index value. The contribution that BD has on other traits such as feed effciency (
Manafiazar et al. 2016), longevity (
Zavadilová and Stipkova 2012), and fertility (
Jagusiak et al. 2014) in dairy cows could be a possible explanation to this effect. A recent study on the feed efficiency of dairy cows (
Manafiazar et al. 2016) showed that inefficient animals with high residual feed intake and high dry matter intake (DMI) tend to have greater body size, including a deeper body.
Manzanilla-Pech et al. (2016) also reported a moderate genetic correlation between BD and DMI in Holsteins from both the Netherlands (0.26) and the United States (0.49). Different countries, such as Australia, New Zealand, and the Netherlands, have incorporated indirect measures of feed efficiency into their selection indices, including live or predicted live body weight, body condition score, and size (
Brito et al. 2020;
Pryce et al. 2014,
2015). It is known that feed accounts for as high as 60% of dairy farm expenses (
Connor 2015), therefore, one could expect that BD would have noteworthy impact on profit of dairy farms.
Zavadilová and Stipkova (2012) studied correlations between BD and length of productive life (LPD), i.e., number of days between first calving and culling, and number of lactations (NL). They showed that BD had unfavorable genetic correlations with both LPD and NL, either when correcting LPD and NL for milk production in first lactation (−0.22 and −0.21, respectively) or not (−0.23 and −0.26, respectively). Thus, regardless of milk yield in first lactation, cows with a deeper body would have fewer days between first calving and culling, and have less lactations than an average dairy cow.
Zink et al. (2014) showed that deeper animals tend to perform poorly for fertility traits, with unfavorable genetic correlations increasing from first to second lactation for days open (0.14 to 0.37) and first service to conception (0.23 to 0.43), suggesting that a deeper animal would stay open for more days and would take longer to conceive.
Jagusiak et al. (2014) reported corroborating results, finding that daughters of bulls with a higher breeding value for BD had worse nonreturn rates in primiparous cows, with an unfavorable genetic correlation of –0.41, implying that deeper daughters would need to be re-bred more times than shallower daughters. While conception is delayed and calving interval lengthens due to unsuccessful inseminations, cows could spend a significant portion of their lactation at suboptimal production levels, typically affecting milk-based income. The period a cow stays open is also often extended, while maintenance costs are not offset by any income from milk production.
3.2 Stepwise backward regression analysis
Stepwise backward regression was used to identify the most important predictive traits associated with Pro$ according to the AIC, for which lower values indicate better fitting models. In practice, if a trait does not significantly contribute to the variation observed in Pro$, its removal from the model will either have no impact in AIC value or will decrease it.
Results in
Table 3 show that there was virtually no difference between adjusted
R2 obtained in the full (0.7238) and best reduced model (0.7239). However, the best reduced model showed a lower AIC (122 066) compared to the full model (122 075) and the full model relative likelihood (RL) was only 0.01, indicating that the full model was 0.01 times as probable as the best reduced model to minimize the information loss. Therefore, the reduced, more parsimonious, model showed better goodness of fit after penalizing for the number of estimated parameters to discourage overfitting. All the covariates in the best reduced regression model were retained with a
P < 0.01, except for TL
2 (
P = 0.136), RA
2 (
P = 0.028), and UD (
P = 0.065).
As expected from the estimated effect on Pro$, removing either HD, BD, or DC would negatively impact the model the most, leading to the highest loss of information (AIC = 122 986, 122 833, and 122 833, respectively). The quadratic terms for RLSV, CW, RTP, and HFE were removed from the best reduced regression model, but this was not true for their respective linear terms, suggesting that these traits might have a linear association with Pro$ when all other conformation traits are taken into consideration.
Furthermore, in addition to BQ having the smallest effect on Pro$ (−$1.02), SBR results showed that BQ was removed from the best reduced model. Recent genetic parameters estimated for conformation traits in Canadian Holsteins (
Oliveira et al. 2021) showed moderate genetic and phenotypic correlations between BQ and DC (0.44 and 0.36, respectively) and UT (0.43 and 0.35, respectively), indicating that the contribution from BQ to Pro$ might have been captured indirectly through other traits in the model. Only the quadratic term of UD was retained in the best reduced model, suggesting that UD has a perfect quadratic correlation with Pro$, and an intermediate-sized udder is most desirable.
3.3 Principal component analysis
Eigenvalues, proportion of the explained variance, and their cumulative sum over the 26 principal components are shown in
Table 4. The first seven PCs explained 67.7% of the total variance among traits and all of them had eigenvalues greater than one, which suggests that they could be considered in a reduced model to predict Pro$ without significant loss predictive ability (
Kaiser 1960). Nonetheless, only components one and two will be discussed here.
The first principal component (PC1) explained 28.45% of the total variance among traits, opposing variables characterized by a strongly positive or negative coordinate on the
x-axis (
Fig. 2). Given the large number of animals in the data set, only a subset of 600 animals comprising the highest and lowest Pro$ ranked animals was plotted (
Fig. 2). It was noticeable that PC1 separated the data in two clusters of bulls based on their shared characteristics among trait RBVs (
Fig. 2). The right cluster contains 300 bulls with the lowest Pro$ index value, whereas the left one clusters the top 300 bulls ranked according to their Pro$ index value. The same clustering pattern persisted even when all bulls were added to the biplot, as it can be visualized on the ShinyApp.
In a PCA biplot, the cosine of the angle between pairs of vectors is inversely proportional to the correlation between the corresponding variables. Therefore, we observed a high correlation between MSL, FA, FTP, RTP, RLRV, HFE, HD, and FAN (
Fig. 2), meaning that animals with high Pro$ index values would have similar RBV for such traits. The clustering of BQ among traits contributing towards highly profitable bulls is noteworthy, and it corroborates with the results from the SBR analysis that the contribution of BQ to the variance in Pro$ was indirectly captured by other traits.
Variable contributions to each PC are directly proportional to the size of their eigenvector coefficients, whereby coefficients with different signs contribute in opposite directions to their PC. Traits such as PW, LS, THP, RA, UD, and ST had very low eigenvector coefficients (−0.05 to 0.03) for PC1 (
Table 4), therefore suggesting they contribute very little to the differentiation between high and low Pro$ bulls. Most traits from the MS related to udder shape and capacity, e.g., UF, MSL, FA, RAH, and RAW, showed the highest contribution to PC1 (
Table 4), and are expected to have similar RBV among high Pro$ bulls (
Fig. 2).
Conversely, positioning of TL among low Pro$ bulls (
Fig. 2) indicate that these bulls have similar RBVs for TL. Even though TL is an intermediate-optimum trait, the negative net estimate from the regression analysis (−5.07) (
Table 3) shows that long teats are expected to have a greater negative effect on Pro$ compared to short teats. In addition to longer teats being attributed to a higher risk of injury from housing and handling (
Berry et al. 2004), milking teats that are excessively longer than the teatcup liner may suffer deterioration of blood circulation, form edema and cyanosis, and have an overall reduction in immunity of the tissue (
Gašparík et al. 2019). Even though short teats are not expected to have such a high negative economic impact, they are not desirable from a management perspective since short teats are associated with increased kick offs, incomplete milkings, and attachment failures in automated milking systems (
Carlström et al. 2016;
Dechow et al. 2020). Therefore, the concerning genetic trend towards smaller TL in Canadian Holsteins (
Oliveira et al. 2021) prompted changes to the weighting given to TL in the genetic evaluation to promote selection for longer teats, while selecting for improved MS (
Fleming and van Doormaal 2021).
The second PC (PC2) accounted for 8.64% of the total variance in Pro$ and contrasting coefficients can be seen between traits from the MS and other scorecard traits (
Table 4). Variables with opposite signs for eigenvector coefficients have a different direction of variation, which suggests a subclustering of animals among high Pro$ bulls that would have a strong performance for MS versus DS and FL (
Fig. 2). Even though HD had the highest effect on the variance in Pro$, it seemed to have nearly zero contribution to PC2 (−0.002), indicating that this trait might not add to the differentiation among highly profitable bulls.