This study considered CE scores from the first parity of calf birth events between 2002 and 2024 for Korean Holstein cattle. Phenotype data were provided by the Dairy Cattle Improvement Center (DCIC) of Korea. Under the Korean National Evaluation System, CE consists of four ordered categories (1 to 4) based on the observed level of difficulty at calf birth. The difficulty level was measured by the degree of assistance provided during calving. A CE of 1 indicated a normal calving event without help at birth (non-assisted). In contrast, a score of 4 indicated veterinary assistance (extreme difficulty at birth requiring surgical help). Intermediate scores (i.e., 2 and 3) indicated the total number of personnel assisted during calving: slightly assisted calving (assisted by one person) and difficult calving (assisted by two or more persons), respectively. We defined CE as a trait of calf (progeny), with an assumption that calf’s physical dimensions and birth weight play a major role in most calving problems. This trait definition would enable the evaluation of the direct contribution of sires in subsequent progeny generations.

We applied several data filters to obtain the final datasets. Animals (calves) were required to have valid parental information. Information of multiple births (twins and triplets) was discarded. Datasets only included observations associated with respective gestation lengths of 260 to 305 days of dams. More restrictions on the parity 1 dataset included dams’ calving ages residing within 20 to 42 months. We also removed farm records containing only normal calving to avoid any possible data recording bias by farmers. Using pedigree and birth information of animals and calving information of dams, we traced calf IDs for respective calving records. However, many male and some female calves’ identification numbers were missing, even though their valid sire and dam information and other essential information were available. Therefore, we assigned dummy IDs to those calves to avoid their exclusion from analysis. Afterwards, we put a restriction on the number of observations per HY level and prepared two different datasets such as DS5 (≥5 obs. per HY) and DS10 (≥10 obs. per HY) to study the effects of CG sizes on genetic analyses of animals. A total of 117,921 and 80,389 calf records were remained in DS5 and DS10 datasets, respectively, for further analyses. These final datasets (see more in Table 1) consisted of information on the sex of the calf (SEX), dam’s calving herd, dam’s calving year, dam’s calving season, and dam’s calving age (DCA). Four seasons of calving (i.e., spring, March to May; summer, June to August; autumn, September to November; and winter, December to February) were considered. We prepared the related animal pedigree files for each dataset from a pedigree database provided by another association, the Korea Animal Improvement Association (KAIA), which is responsible for managing the dairy cattle pedigree databases in the country. The pedigree for DS5 and DS10 datasets comprised 317,126 and 233,569 animals, in which animals were traced back up to 23 generations.

We estimated (co)variance components and genetic parameters for CE trait through animal models using the first parity dataset. We presented four univariate animal models for genetic parameter estimation. All four models included a random dam (maternal) effect, which was also assumed to be correlated with the random animal (direct) component. However, maternal permanent environment effect was ignored because of the single parity data structure. The SEX was a fixed effect across models, whereas the DCA (in days) was a fixed covariate effect in models to control for the effect of different ages of dams on CE. The CG effects e.g., dam’s calving herd, year, and season were considered as fixed effects but fitted differently across four models. To illustrate, Model 1 (M1) considered a combined fixed effect of the former three effects as herd-year-season (HYS). In model 2 (M2), the combined herd-year (HY) and year-season (YS) effects were fitted. Model 3 (M3) included a HY effect and a season (S) term as an independent effect. In model 4 (M4), herd (H), and YS effects were fitted. The reason for CG variations was to examine the influence of CG structures on variance component estimation. The BLUPF90+ software package [16] was used to estimate variance components through the average-information REML algorithm, genetic parameters, and their standard errors (SEs). The linear mixed model of the animal model with maternal effect in matrix notation was as follows:

where y was the vector related to CE; b was the vector of fixed effects, i.e., SEX, DAC, HYS or others (H, S, YS, HY etc.); d was the vector of random animal (direct) effect; m was the vector random maternal effect; and e was the vector of random residual effect. X, Z_d, and Z_m were design matrices relating effects to the phenotype. A covariance structure for random effects was assumed as follows:

where σd2 was the direct genetic variance, σm2 was the maternal genetic variance, σe2 was the residual variance, and σ_dm was the covariance between direct and maternal genetic effects. Therefore, the genetic covariance matrix (G₀) between d and m was:

Total phenotypic variance ( σp2) and different heritability estimates (direct, hd2; maternal, hm2; and total, T²), direct-maternal genetic correlation (r_dm) [17] were calculated using the above (co)variance components as follows:

and rdm=σdmσd2×σm2,

Here, hd2, hm2, and T² expressed the ratio of direct, maternal and total heritable variances available for response to selection over phenotypic variance, respectively.

Approximated SEs of genetic parameters were obtained from (co)variance components using the BLUPF90+ software package, in which a Monte Carlo method was implemented for computation of SE following the study of Houle and Meyer [18]. We plotted the distribution and genetic trends of predicted direct and maternal estimated breeding values (EBVs) using the ‘ggplot2’ package in R [19]. We also calculated the mean squared error (MSE) statistic of these models using predicted phenotype estimates for each model and dataset. In the following sections, we will mainly emphasize on direct and maternal heritability estimation and the genetic correlation between direct and maternal components. Total heritability estimates are provided for completeness purpose only.

In this study, we investigated two datasets of CE (DS5 and DS10) on calves, which were born from the first parity of Holstein cows. Table 2 presents (co)variance and genetic parameters estimates for direct and maternal genetic components of CE from four different animal models with a maternal effect. Across models, the estimated hd2 parameter ranged from 0.005±0.002 to 0.168±0.012 within DS5 animals and from 0.014±0.005 to 0.234±0.018 within DS10 animals. Estimates of maternal heritability ( hm2) were lower than hd2 from all models using both datasets. We found that hm2 ranged from 0.001±0.002 to 0.090±0.007 within studied animals. In both datasets, M4 model derived much higher hd2 estimates than first three models, which is also similar to hm2 estimates. This study also showed somewhat similar or slightly lower total heritability (T²) estimates than other heritability estimates from DS5 and DS10 analyses, where obtained T² ranges were from 0.007 to 0.018 across datasets. Like other h² estimates, first three models (M1 to M3) also derived slightly lower T² estimates than M4 model. Although most random (co)variance estimates were generally lower in this study, those from M1, M2, or M3 model appeared shrunk by some magnitude compared to the M4 model. Genetic correlation estimates between direct and maternal effects (r_dm) were mainly negative except for the slightly positive value (r_dm, 0.078±1.72) from M1 model using DS5 dataset. For DS10, r_dm estimates were moderate to strongly negative (−0.814±0.169 to −0.491±0.047). Overall, most results indicated negative correlation between direct and maternal genetic components of CE.

The knowledge of various genetic effects, such as direct and maternal effects and their relationship, is the key to CE improvement for any dairy cattle. In this study, we demonstrated that heritability (h²) estimates for CE in Korean Holsteins ranged from very small to low. In similar Korean Holsteins, a previous study using a sire-maternal grandsire (S-MGS) model for parity-1 progeny [15] revealed direct h² (0.11) and maternal h² (0.05) estimates, somewhat close to those of ours with M4 model estimates. Our results also partially agree with the study by Eaglen and Bijma [17] in Dutch Holstein-Friesian, showing their h² estimates (direct, ~0.08; maternal, ~0.04) in our estimation range. Mujibi and Crews [10] have reported h² estimates (direct, 0.14; maternal, 0.06) for Charolais cattle, also in agreement with ours. For those of very low h² estimates in this study, multiple studies from Iranian Holstein cattle [20,21] provide support with their very small direct (0.02 to 0.041) and maternal h² (0.002 to 0.012) estimates. Similar low h² for direct (0.03) and maternal CE (0.02) were also available from Eaglen et al [13]. However, Roughsedge et al [22] using a linear mixed model indicated a wide range of direct h² among beef breeds (0.13 to 0.35), aside from their agreeable maternal h² range (0.07 to 0.11). Differences among reports are mainly due to differences in fitted factors, model types, trait definitions, breeds, etc. According to some studies [8], linear models tend to yield lower estimates than threshold models. Salimi et al [20] have argued that their large phenotypic variances (or residual variances) compared to genetic variances are possibly arising from their recording methods and herd management practices, inevitably underestimating their population’s direct and maternal components. Eaglen and Bijma [17] have also argued that whether the model is an animal or S-MGS model, some h² parameters (e.g., maternal h²) are prone to inaccurate estimations despite having a sufficiently large dataset.

In the present study, the genetic correlation estimates between direct and maternal effects varied across models and datasets. Most models indicated a negative association between the two genetic components. The distribution of these effects among animals also supports this negative association (Figure 1), despite some dataset related variations. Overall, our r_dm estimates ranged from lowly negative to strongly negative. Previous reports provide good overall support for our observed negative correlation estimates. Alam et al [15] have used an S-MGS model in Korean Holstein and reported a genetic correlation (−0.68), which supported the present study. Salimi et al [20] and Ghiasi et al [21] have reported somewhat comparable r_dm (−0.41 to −0.43) estimates in Iranian Holstein. The r_dm estimates from the study of Eaglen and Bijma [17] also showed agreeable estimates using an animal model (−0.04 to −0.44). Some agreements in Holstein cattle by a near absence of correlation (i.e., weakly negative to weakly positive) are also available [23,24]. The r_dm estimate varied widely in the literature, and multiple factors could be involved in such variabilities. First, it could be due to differences in breed or population in studies such as beef cattle, in which correlation estimates are often highly negative [25]. Second, a possible estimation bias is also likely to genetic covariance between direct and maternal effects [26].

To further discuss the estimated genetic parameters, we argue that future genetic improvements for CE in Korean Holsteins using an animal model (with maternal effects) could be slower due to lower h² values, similar to any other reproductive trait. Given the increased calving difficulties in Korean Holstein cattle, the higher hd2 than hm2 indicates a greater contribution of sires for CE (through calf weight and calf dimension). In this study, we defined CE as a trait of the calf, which further indicates that present evaluation models will essentially allow relatively slender or less broad calves (any sex) to be born easily but finally cause more difficulties for female calves when they give birth as dams (due to reduction in pelvic dimensions) [13]. In this regard, a negative r_dm indicates the importance of the maternal component in CE evaluation. Our r_dm values suggest that selection improvement in the first parity animals could exert challenges due to such negative associations. Therefore, considering direct effects only for CE improvement (with an ignored negative correlation between direct and maternal genetic components) can eventually reduce the selection progress [27].

This study also revealed some additional challenges regarding data structure in Korean Holsteins. With animal models, genetic parameters were relatively sensitive to data structure changes in studied animals, indicating data connectedness problems with higher-order HY or HYS effects in the model. Our further investigation into datasets suggests that herds and dams are likely confounded with each other to some extent as dams and their first calving daughters hardly changed their herds. Therefore, accurately estimating genetic parameters for CE in Korean Holstein might be challenging. Such evaluations could inflate genetic parameters due to inaccurate estimation of herd and genetic effects. Hence, careful consideration is required for CE parameter estimations. An alternative model to counter such possible herd confounding errors could be S-MGS models, where one sire’s progenies are more likely spread over many herds, which can be a topic for future studies.

Figure 2 and 3 illustrate changes in average direct and maternal genetic merit (EBV) of CE within Holstein cattle based on their birth year. All four models showed no significant trends for direct and maternal EBV estimates. However, the first three models (M1 to M3; herd as a combined effect) and the last model (M4; herd as an independent effect) had noticeable differences in EBV estimates (also see Figure 1). Average animal genetic effects trends from the former three models are less sensitive to yearly differences. In contrast, the last model showed noticeable yearly variations irrespective of datasets.

Most animal models in this study demonstrated negligible genetic improvement of CE. Contrary to our reports, an earlier report in Korean Holsteins using the S-MGS model [15] presented some improvements in direct and maternal EBVs. Such prediction disagreements could appear due to differences in our model from the previous report. For the M4 model, animals’ direct EBVs are deemed negative for an extended period (Figure 2), which is desirable. Using the direct-maternal relationship in Figure 1, we could suggest that their maternal EBV estimates were mostly positive or higher (Figure 3), which is also desirable. Although such desirable scenarios from Figure 2 and 3 could appear to be some improvements in CE over time, these trends using M4 need careful consideration. The M4 model could also have a potential confounding concern between some non-genetic effects and direct genetic effects. The confounded herds and animals in the datasets and the M4 model’s inability to correct for such data structures could potentially over- or under-estimate animal predictions.

Table 2 also presents Akaike information criterion (AIC) estimates and MSE of predictions for each model fit obtained from BLUPF90+ analyses. We showed four linear mixed models following an order of more detailed to less detailed models and then compared their goodness-of-fit statistics. To illustrate, M1 was the most detailed model, which included the HYS effect (a 3-way interaction) and required higher computation requirements due to many levels fit, whereas M4 required the least computation requirements due to a relatively small number of levels fit. The M2 model was confirmed to have the best fit using DS10 and the second best fit model using the DS5 dataset. However, the M4 was a poorly fit model across datasets, despite the lowest model MSE estimates. Although M1 (with HYS) was the best-fit model using the larger dataset according to the AIC estimate, this model needs a careful interpretation due to the overall shrinkage of all variance components, especially with the maternal genetic component. Some possible reasons could be the lack of sufficient data per HYS level using M1 and the presence of disconnected data subsets within the primary dataset for specific factor combinations [28,29]. A possible fitting of data noises within disconnected subsets might not be unlikely due to the higher number of HYS levels with fewer records. The M4 model, which exhibits relatively higher (co)variances than other models, could potentially suffer from an ineffective separation of genetic effects from confounded herd effects. The genetic parameters estimated by M1 or M4 also require critical assessment. In this regard, M2 or M3 animal models or use of S-MGS model could be alternatives to account some of the challenges discussed above.

In contrast, for a routine national evaluation, designing a statistical model for maternally influenced traits should represent a careful balance of the model’s prediction ability and computational feasibility [12]. Note that the above factors might depend on several other factors, including the trait in question, dataset size, data recording biases, computing facility, and time availability for the task. From a simplicity standpoint, an animal model is preferable to other models (e.g., sire models) for routine evaluation of animals because more complicated models often have parameter convergence problems. Besides, an animal model is more appealing and practical, considering its ability to incorporate a large number of samples directly during animal evaluations. Therefore, we suggest considering models such as M2 or M3 in routine evaluation of Korean Holsteins, which could account for some existing data-related issues yet be practical.