2.1 Data

Data for this study were provided by Lactanet (Guelph, Canada) and contained Pro$ values and estimated relative breeding values (RBVs), as well as trait heritabilities, for 26 conformation traits from 9351 proven bulls from the August 2021 genetic evaluation. According to Fleming and van Doormaal (2020a), RBVs are standardized estimated breeding values (mean = 100 and standard deviation = 5) expressed relative to the genetic base of the Holstein breed, which is defined as proven bulls born in the most recent complete 10 year period.

Trait abbreviations, definitions, scoring system, and ideal scores are available in Table 1 (Holstein Canada 2017). Descriptive statistics for RBVs and their reliabilities are presented in Table 2. For further reading about conformation traits and the descriptive statistics of their phenotypic values (i.e., classification scores), please refer to a recent study published by Oliveira et al. (2021).

Table 1.

Scorecard sections, trait names, scoring system and ideal scores, and short trait definitions.

Table 2.

Descriptive statistics of relative breeding values (RBVs) of Canadian Holstein conformation traits and their respective reliabilities.

A minimum reliability of 0.14 for all RBV was imposed to make sure all bulls had a minimum effective daughter contribution of one, as derived from Mrode (2014) (eq. 1).

where rel_min is the calculated minimum reliability (rounded up), EDC_min is the desired minimum effective daughter contribution (one in this case), and h² is the heritability of the trait (0.53 in this case, which is the heritability of stature that was the highest value among all studied conformation traits).

.

2.2 Statistical analysis

2.2.2 Stepwise backward regression analysis

Stepwise backward multiple regression (SBR) model selection was carried out to identify the most important predictive variables associated with Pro$ according to the Akaike information criterion (AIC) using the stepAIC() function implemented in the R MASS package (Ripley et al. 2021), as exemplified in Fig. 1.

Fig. 1.

Diagram of stepwise backward regression carried out to identify the most important predictive variables associated with Pro$ according to the Akaike information criterion (AIC).

Relative likelihood (RL) was used to assess the likelihood that the alternative models would minimize information loss relative to the best model (the one with the lowest AIC value), as defined in eq. 2 (Burnham and Anderson 2002):

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).

Table 3.

Results from the full model (FM) and best reduced model (RM) for the multiple polynomial regression of relative breeding value RBV for 26 conformation traits on Pro$ selection index values. Net estimate considering both linear and quadratic effects on Pro$ are also shown.

Heel depth was the most important trait contributing positively to changes in Pro$, whereby a deeper heel on the outside hoof would lead to higher overall profit in terms of Pro$. An increase in one unit of RBV for HD is estimated to increase Pro$ by $79.13. This is quite a compelling result if one considers the relationship between heel depth, horn lesions, and lameness. A study on Canadian Holstein cows found that extremely shallow heel depth was associated with a significantly higher incidence of horn lesions (e.g., sole and toe ulcer, sole hemorrhage, and white line disease) when compared with the average incidence in the herd (Chapinal et al. 2013). Such lesions are among the main causes of lameness in dairy cattle (Leach et al. 2010; Solano et al. 2015; Cartwright et al. 2017), which is one of the major causes of economic losses in dairy farms (Enting et al. 1997; Ettema and Østergaard 2006). There is a tight association between lameness and higher culling rates (Booth et al. 2004; Hultgren et al. 2004), decreased milk production (Warnick et al. 2001; Green et al. 2002) and poor fertility (Hernandez et al. 2001; Garbarino et al. 2004; Dobson et al. 2008), not to mention the significant negative impact of laminess on animal welfare (Leach et al. 2010).

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 R² 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² (P = 0.136), RA² (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.

Table 4.

Eigenvectors, eigenvalues, proportion of the explained variance (EV), and their cumulative sum (CEV) of principal components (PC) 1 to 7 calculated from RBVs of 26 conformation traits from 9351 proven bulls.

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.

Fig. 2.

Biplot of Principal Components (PC) 1 and 2. Two ellipses were drawn clustering the top and bottom 300 proven bulls ranked according to their Pro$ index values (left and right ellipses, respectively). Arrow size is directly proportional to the variable's contribution to each PC. Variables with arrows in opposite directions have different directions of variation in each PC. Proportion of variance (%) explained by each PC is shown inside the parentheses on their respective axis label. [Colour online.]

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.