Conformation traits of Holstein cows and their association with a Canadian economic selection index

Abstract Pro$ is a Canadian economic selection index aimed to maximize profit by increasing production yields, while maintaining conformation and functional traits. Currently, there is an interest in understanding the individual contributions of conformation traits recorded in Canada to the overall economic value of a cow and whether they are equally important. We used multiple polynomial regression and principal component analysis to assess the association of 26 conformation traits with Pro$ using relative breeding values (RBVs) from 9351 proven bulls. The best reduced regression model explained 72.5% of the variance in Pro$, with heel depth and body depth having the highest and lowest effect on Pro$ values, respectively. Four traits classified as intermediate optimum traits, e.g., teat length, showed significant linear association with Pro$ instead of quadratic, whereas bone quality was not significantly associated with the index. Principal component analysis indicated that highly profitable bulls share similar RBV, with a subclustering of bulls of daughters with better mammary system versus better dairy strength and feet and legs. These results provide understanding of the individual contributions of conformation traits to Pro$ and give information to guide the Canadian dairy industry on how to best consider these traits in recording and genetic evaluation programs. Résumé Pro$ est un indice de sélection économique canadien qui vise à maximiser les profits en augmentant les rendements de production, tout en maintenant les caractéristiques de conformation et de fonction. Actuellement, il y a un intérêt pour la compréhension des contributions individuelles des caractéristiques de conformation enregistrées au Canada à la valeur économique générale d’un bovin et pour savoir s’ils sont d’importance égale. Nous avons utilisé la régression polynomiale multiple et l’analyse des composantes principales afin d’évaluer l’association de 26 caractéristiques de conformation à l’aide de Pro$ en utilisant les valeurs d’élevage relative (RBV — « Relative Breeding Values ») provenant de 9351 taureaux éprouvés. La meilleure régression réduite expliquait 72,5 % de la variance dans Pro$, avec la profondeur du talon et la profondeur corporelle ayant les plus grands et plus petits effets sur les valeurs Pro$, respectivement. Quatre caractéristiques classées comme caractéristiques optimales intermédiaires, p. ex. longueur de trayon, ont montré une association linéaire significative avec Pro$ au lieu de quadratique, tandis que la qualité d’os n’était pas associée de façon significative à l’indice. L’analyse de composantes principales a indiqué que les taureaux fortement profitables partagent des RBV semblables, avec un sous-regroupement de taureaux provenant de filles ayant un meilleur système mammaire contre une meilleure force laitière et meilleurs pieds et jambes. Ces résultats offrent une compréhension des contributions individuelles des caractéristiques de conformation à Pro$ et offre l’information pour guider l’industrie laitière canadienne sur comment mieux considérer ces caractéristiques pour enregistrement et programmes d’évaluation génétique. [Traduit par la Rédaction]


Introduction
Breeding programs across the globe are shifting their selection goals to more balanced indices that include traits of direct and indirect economic significance, such as conformation, fertility, and health (Cole et al. 2021). Selective breeding in Holstein cattle in Canada has utilized the Lifetime Performance Index (LPI), formerly known as Lifetime Profit Index, since 1991. LPI has been a key index for producers who want a herd with strong conformation and high solid yields, while maintaining functional traits. In the last two decades, Lactanet (Guelph, Canada) made three changes to LPI's main components (i.e., production, durability, and health & fertility), shifting their respective emphasis from 57: 38:5 (2001) to 54:36:10 (2005), 51:34:15 (2008), and 40:40:20 (2016), which is the current emphasis (Lactanet 2022).
Nonetheless, increasing requests from the industry led Lactanet to develop a second national index to help dairy producers who essentially have all their farm revenue based on milk sales (Van Doormaal 2015). Therefore, Lactanet introduced Pro$ in 2015, a new profit-based national economic index. The goal of Pro$ is to maximize profit by increasing production yields, while maintaining the level of functional and conformation traits. According to Lactanet (Van Doormaal 2015), sires with a Pro$ of $1000 are expected to produce daughters that have an average accumulated profit over their first six years of age that is $1000 higher than daughters of the average bull in Canada. Conformation traits have been of great interest in the dairy cattle industry for decades, with over two million classifications accumulated in Canada since the introduction of descriptive classification in the late 1960s and the linear scoring system in the early 1980s . Research studies in Canada made genetic evaluation of thurl placement possible in 2011 (Muir 2011), and the inclusion of a dedicated component for rump traits into the LPI selection index in 2019 (Beavers and van Doormaal 2019). Additionally, the needs of Canadian breeders led the Canadian Holstein breed association (Holstein Canada) to recently start recording locomotion, udder floor, and front legs view as research traits, with first genetic evaluations published in December 2020 (Fleming and van Doormaal 2020b).
Classification of dairy cows is a long-dated practice to help producers make mating and culling decisions. Classification focuses on the selection of healthier, more productive, and fertile animals for a longer herd life based on their conformation traits. According to Holstein Canada (Holstein Canada 2021), the Canadian Holstein classification system comprises 26 different linear and nonlinear descriptive traits that are measured or scored on first lactation animals. These traits are distributed in four scorecard sections: mammary system (MS), feet & legs (F&L), dairy strength (DS), and rump (R). After penalizing for defective characteristics (e.g., blind quarter, abnormal claw, weak back, and advanced tailhead), a final score is calculated considering the relevance of each of the scorecard sections for the dairy industry.
Given the extensive number of conformation traits being currently recorded, the dairy industry might consider stopping the genetic evaluation or even the recording of some of these traits in the future. Understanding the monetary contribution of conformation traits to an economic index that is expressed in dollar values, such as Pro$, may provide the dairy industry with more tangible and practical information to guide such decisions. Therefore, the present study used multiple polynomial regression and principal component analysis to assess the association of 26 conformation traits with the Pro$ selection index.

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

Multiple polynomial regression analysis
Pro$ was regressed on all conformation traits through a multiple polynomial regression. Intermediate-optimum traits (i.e., FTP, RTP, TL, UF, UD, FAN, HD, RLSV, FLV, RA, THP, CW, HFE, and ST) were considered to have a quadratic association with Pro$, therefore, their quadratic terms were also included in the regression models. In total, 26 linear and 14 quadratic terms were fitted. Regression analyses were carried out using the lm() function from the R Stats package (R Core Team 2020).

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

Principal component analysis
Principal component analysis (PCA) was carried out for all 26 conformation traits to explore clustering patterns in the data set using the prcomp() function implemented in the R Stats package (R Core Team 2020). The extraction and visualization of results were done through a biplot (variables and samples) using functions from the R factoextra package (Kassambara and Mundt 2020), to which two ellipses were drawn clustering the top and bottom 300 animals ranked ac-cording to their Pro$ values. Our main goals with this analysis were to investigate the strength and direction of correlations between conformation traits, and to identify a potential profile of conformation traits for top and bottom performing bulls based on their Pro$.

Results and discussion
An R Shiny WebApp was developed to display results from this research and it is freely available at https://alcantara.shin yapps.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).

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).
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. 2014Pryce et al. , 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.

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 2 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 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.] 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 intermediatesized udder is most desirable.

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.

Conclusion
With the growing number of traits being recorded and evaluated, this study provided a better understanding of the individual contributions of conformation traits to Pro$ and may guide the Canadian dairy industry on how best consider these traits in the recording and genetic evaluation programs in the future. Current genetically evaluated conformation traits in Canada jointly explain considerable variation in the Pro$ index (72%). The results presented herein showed that HD, BD, and DC had the highest impact on Pro$. On the other hand, BQ showed no association with Pro$ when all other conformation traits are considered.

Article information
History dates