Variance components and correlations of female fertility traits in Chinese Holstein population

Background The objective of the present study was to estimate (co)variance components of female fertility traits in Chinese Holsteins, considering fertility traits in different parities as different traits. Data on 88,647 females with 215,632 records (parities) were collected during 2000 to 2014 from 32 herds in the Sanyuan Lvhe Dairy Cattle Center, Beijing, China. The analyzed female fertility traits included interval from calving to first insemination, interval from first to last insemination, days open, conception rate at first insemination, number of inseminations per conception and non-return rates within 56 days after first insemination. Results The descriptive statistics showed that the average fertility of heifers was superior to that of cows. Moreover, the genetic correlations between the performances of a trait in heifers and in cows were all moderate to high but far from one, which suggested that the performances of a trait in heifers and cows should be considered as different but genetically correlated traits in genetic evaluations. On the other hand, genetic correlations between performances of a trait in different parities of cows were greater than 0.87, with only a few exceptions, but variances were not homogeneous across parities for some traits. The estimated heritabilities of female fertility traits were low; all were below 0.049 (except for interval from calving to first insemination). Additionally, the heritabilities of the heifer interval traits were lower than those of the corresponding cow interval traits. Moreover, the heritabilities of the interval traits were higher than those of the threshold traits when measuring similar fertility functions. In general, estimated genetic correlations between traits were highly consistent with the biological categories of the female fertility traits. Conclusions Interval from calving to first insemination, interval from first to last insemination and non-return rates within 56 days after first insemination are recommended to be included in the selection index of the Chinese Holstein population. The parameters estimated in the present study will facilitate the development of a genetic evaluation system for female fertility traits to improve the reproduction efficiency of Chinese Holsteins.


Background
Female fertility traits are the most economically functional traits affecting the efficiency of the dairy industry. Cows with excellent fertility resume cycles soon after calving and become pregnant when properly inseminated [1]. According to previous studies, poor reproduction has become a worldwide problem and the most probable cause for culling [2,3]. Due to the unfavorable correlations between yield and female fertility traits, intensive selection for milk production over the past five decades has led to a strong decline in cow reproduction, particularly in highproducing populations [4,5]. Moreover, the heritabilities of most fertility traits are typically below 5% [6][7][8], thus, selection for female fertility is not as efficient as selection for production traits. Therefore, the reproductive decline in modern dairy cattle cannot be ignored, as the costs of extra inseminations, fertility treatments and involuntary culling are increasing.
Although the heritabilities of reproduction traits are relatively low, their additive genetic variances have remained sufficiently high to perform effective selection. Nordic countries (including Denmark, Finland and Sweden) are pioneers in applying fertility traits in dairy cattle breeding and have implemented fertility traits in their selection index since 1994. The Nordic Cattle Genetic Evaluation has shown that the genetic trend for the female fertility index of Holstein bulls has increased by approximately 1% per year since 2000 (http://www.sweebv.info/ba52nycknav.aspx). Furthermore, the results of routine genetic evaluations in Canada (https://www.cdn.ca/) and the United States (http:// www.holsteinusa.com/) have shown that the decline in the genetic trend of fertility traits has flattened and even slightly reversed after increasing the weight of fertility traits in the selection index [9]. In a recent test run of the multiple across-country evaluation (MACE) by Interbull in August 2016, 21 countries participated in the genetic evaluation of female fertility traits [10]. In China, however, an insemination and calving record collection system has not yet been implemented on the national scale, therefore, female fertility traits are still not included in the current selection index for Chinese Holsteins. Preliminary results from a genetic evaluation of five fertility traits (age at first insemination, number of inseminations per conception, days from calving to first insemination, days open and calving interval) in Chinese Holsteins have been recently published [11]. However, that dataset was relatively small and did not include several internationally recommended female fertility traits, such as conception rate and non-return rate at first insemination. Moreover, the main limitation of the previous analysis is that the performances of a trait in heifers and cows were simply treated as the same trait, which differs from various previous genetic evaluation studies [6,12]. The physiological status of cows changes considerably after first calving [13], thus, the performances of heifers and cows have typically been considered biologically distinct [14,15]. Nevertheless, the genetic relationships among female fertility traits across multiple parities in dairy cattle remain less well investigated.
With respect to the models and algorithms typically used in the genetic evaluation of female fertility traits, many studies have shown that linear models perform similarly to threshold models but with less computational demand, consequently, linear models have been widely used in routine genetic evaluations [12,16]. In addition, animal and sire models are currently the most frequently implemented models in routine genetic evaluations. Compared with sire models, animal models are theoretically superior for estimating variance components and breeding values. Several previous studies have reported that animal models provide more accurate predictions than sire models [17][18][19].
The aim of the present study was to assess the genetic relationships across multiple parities, estimate the variance and covariance components, and calculate the heritabilities and genetic correlations of female fertility traits in Chinese Holsteins. The results of this study will facilitate the development of a routine genetic evaluation system for female fertility traits, resulting in a more balanced dairy cattle breeding program in China.

Data
Field records of insemination and calving from 2000 to 2014 were collected from 32 herds in the Sanyuan Lvhe Dairy Cattle Center, Beijing, China. The following traits were recorded: interval from calving to first insemination (ICF), days open (DO), interval from first to last insemination (IFL), conception rate at first insemination (CR), number of inseminations per conception (AIS) and nonreturn rate within 56 days after first insemination (NRR56). For binary traits, CR was coded as 0 (failure) or 1 (success), and NRR56 was coded as 1 when there was no subsequent insemination within 56 days after the first insemination and 0 otherwise. All traits recorded before first calving were considered as heifer traits and were coded as parity 0. Traits measured in first-lactation cows were coded as parity 1, and the same strategy was used for the remaining lactations. Among these traits, ICF and DO were available only for cows. All the remaining traits were measured for both heifers and cows. In this context, h and c were used as suffixes in trait abbreviations to distinguish heifers and cows.
Most fertility traits included censored records, as mistakes in recording and some culling events arose during data collection or fertility events were still in progress. The last insemination with a positive pregnancy test result or with the following calving record was considered a verified conception. In the present data, 5.92% of heifers and 13.92% of cow-lactations showed no verified conception from the last insemination. In such cases, the IFL and AIS were calculated according to the last insemination, and the expectation of the trait was added as a penalty. Based on uncensored data, the expectations of IFLh, IFLc, AISh and AISc were 21 days, 43 days, 1.5, and 2.0, respectively. Similarly, for cows without ascertained conception, the penalty for DO was the same as that for IFL. To further reduce the proportion of censored records, we only used data for females with birth years from 1999 to 2013, insemination years from 2000 to 2014, and calving years from 2000 to 2014. To reduce possible selection bias, cows without heifer records were also removed from the data. Moreover, the data in the analysis contained only females inseminated by AI bulls. Other criteria for data editing included age at first insemination between 270 and 900 days, age at first calving between 500 and 1100 days, gestation length between 240 and 300 days, ICF between 20 and 230 days, IFL less than 365 days and AIS less than eight times per parity. Records with values smaller than the lower limit were deleted as outliers, while those with values larger than the upper limit were set as the upper limit.
Based on the distribution of records across parities for each trait, only data from heifers and the first three lactations were used in the analysis. After editing, the total number of cows in the final dataset was 88,647, and the total number of records (parities) was 215,632. The numbers of observations available for heifers and the first three parities of cows were 88,647, 60,022, 41,858 and 25,105, respectively. Heifer records comprised 41.11% of the entire dataset.
Full pedigree data were provided by the Dairy Association of China (Beijing, China), and each animal was traced back as many generations as possible. Unknown parents were assigned into 16 phantom groups, based on sex and birth year. In total, the pedigree included 179,939 females and 4310 males, with birth years from 1930 to 2013. Considering the number of heifers/cows in the present study, the size of pedigree was small. This was determined by the data structure. In the current study, the phenotypic data covered a long period, and the females came from 32 herds during birth year 1999 to 2013. Thus, the phenotypic data already covered a large number of female ancestors (such like mother, grandmother and so on). Therefore, the number of animals in the pedigree was only about two times as large as the number of animals with fertility records, even though the heifers/cows were traced back as many generations as possible.

Models
Five linear Gaussian animal models were used to estimate the (co)variance components of female fertility traits. First, data across multiple parities were analyzed using a multi-trait model in which the performances of a trait in different parities were treated as different traits. For each of the traits defined both in heifers and cows (IFL, CR, AIS and NRR56), a four-trait model was used. For traits only defined in cows (ICF and DO), a threetrait model was used.
In addition, cow traits (ICF, DO, IFLc, CRc, AISc and NRR56c) were further analyzed by pooling data over parities and using a single-trait model treating records in multiple parities as repeated measurements.
Correlations between a pair of heifer traits and a pair of cow traits were estimated using a two-trait model. For a pair of heifer traits (IFLh, CRh, AISh and NRR56h), the following two-trait model was used: For a pair of cow traits (ICF, DO, IFLc, CRc, AISc and NRR56c), the following two-trait model was used: In the above models, y i was the vector of observations (i = 0, 1, 2, 3 representing heifers and parities 1, 2 and 3 of cows in models (1) and (2); i = 1, 2 representing fertility traits 1 and 2 in models (4) and (5)); β i was the vector of fixed effects for the ith trait, including age group (in terms of age at first insemination, divided into 270-439 days, 440-469 days, 470-499 days, 500-529 days and 530-900 days), herd-year of calving (for ICF and DO) or herd-year of first insemination within parity (for IFL, CR, AIS and NRR56), year-month of calving (for ICF and DO) or year-month of first insemination within parity (for IFL, CR, AIS and NRR56), AI technician, sexed semen (coded as 1 for the use of sexed semen and 0 for the use of non-sexed semen), and parity (available only for models (3) and (5); a i was the vector of random additive genetic effects for the ith trait; pe i was the vector of permanent environmental effects for the ith trait (available only for models (3) and (5)); and e i was the vector of random residual effects for the ith trait. X i , Z i and W i were the incidence matrices connecting β i , a i and pe i to y i . It was assumed that In the above distributions, A is the matrix of additive genetic relationships between individuals in the pedigree; I is the identity matrix; σ 2 a i , σ 2 pe i and σ 2 e i were the additive genetic variance, permanent environmental variance and residual variance of the ith trait, respectively; σ a i a j , σ pe i pe j and σ e i e j (i ≠ j) were the additive genetic covariance, permanent environmental covariance and residual covariance between the ith trait and the jth trait, respectively.
The (co)variance components in the above models were estimated using average information-restricted maximum likelihood (AI-REML) implemented in the DMU package [20]. The asymptotic standard errors of the variances and covariances were obtained from the average information matrix. Standard errors of heritabilities and correlations were calculated according to an expansion of the Taylor series, as previously described [21].

Results
The descriptive statistics for female fertility traits regarding heifers and cows are presented in Table 1. In general, all fertility capacities of the heifers were superior to those of the cow: IFL were 47.4 days shorter, CR was 27.0% higher, NRR56 was 14.8% higher, and AIS was 1.035 less. For all traits of both heifers and cows, large phenotypic variations were observed.
Estimates of the variance components and heritabilities of fertility traits using a multi-trait model and a singletrait model are shown in Table 2. Estimates of heritability of female fertility traits were generally low, ranging from 0.005 (first parity of cows for NRR56) to 0.102 (first parity of cows for ICF). The approximate standard errors for the heritability estimates ranged from 0.001 to 0.011. The estimated heritabilities all differed significantly from zero. As shown in Table 2, binary traits (CR and NRR56) did not exhibit large differences in terms of variance components between heifers and cows. However, for IFL and AIS, the additive variances of heifers (56.6 for IFLh and 0.018 for AISh) were only approximately one = fifth those of cows (257.8 for IFLc and 0.095 for AISc in an analysis of data pooled over parities). For cow traits, additive variances for ICF, DO and IFLc in the first parity were higher than those in the second and third parities, but CRc, AISc and NRR56c exhibited no clear trends. In addition, the three interval traits of cows (ICF, IFL and DO) had higher heritabilities (ranging from 0.027 to 0.083) than the other three traits (ranging from 0.006 to 0.026). On the other hand, the heritabilities of the four heifer traits did not differ dramatically (ranging from 0.008 to 0.013).
The genetic and residual correlations between parities of the performances of a trait using the multi-trait model are presented in Table 3. Generally, the genetic correlations between the performances in heifers and in cows were moderate but far from one, ranging from 0.25 for AIS between heifers and cows at the third parity to 0.74 for NRR56 between heifers and cows at the second parity. However, the genetic correlations between parities of cows were generally high, except for those between the third parity and the other two parities which were 0.63 and 0.66 and had large standard errors, the other genetic correlation ranged from 0.87 to 1. On the other hand, the residual correlations between heifers and cows and between parities of cows were very low, ranging from −0.002 to 0.14, indicating small permanent effect on fertility traits. Since all fertility traits had low heritabilities, the phenotypic correlations (not presented) were just very slightly higher than the residual correlations. Table 4 presents the genetic and phenotypic correlations between heifer traits. The genetic correlations between IFLh, CRh, AISh and NRR56h were strong. The absolute values of these genetic correlations ranged from 0.63 (negative between IFLh and NRR56h) to 0.94 (negative between CRh and AISh). The phenotypic correlations between these traits were moderate to high. The absolute ð Þ = additive genetic variance of ith parity using multi-trait models (model 1); σ 2 ei i ¼ 0; 1; 2; 3 ð Þ = residual variance of ith parity using multi-trait models (model 1); ĥ2 i = estimated heritability of ith parity using multi-trait models (model 1); σ 2 a123 = additive genetic variance of a cow trait over parities using single-trait models (model 4); σ 2 pe123 = permanent environmental variance of a cow trait over parities using single-trait models; σ 2 e123 = residual variance of a cow trait over parities using single-trait models (model 4); ĥ2 123 = estimated heritability of acow trait over parities using single-trait models b IFL interval from first to last insemination, CR conception rate at first insemination, AIS number of inseminations per conception, NRR56 non-return rate at 56 days after first insemination, ICF interval from calving to first insemination, DO days open  The genetic and phenotypic correlations between cow traits, estimated using data pooled over three parities, are shown in Table 5. There were strong correlations between DO, IFLc, CRc and AISc. The absolute values of genetic correlations between the four traits ranged from 0.82 (positive between DO and AISc) to 0.98 (negative between DOc and IFLc), and the phenotypic correlations ranged from 0.55 (negative between DO and CRc) to 0.97 (positive between DOc and IFLc). ICF had high genetic correlation and moderate phenotypic correlation with DO, while NRR56c had moderate genetic and phenotypic correlations with CRc and AISc.

Discussion
This study investigated the relationships among parities and variance components in each parity for female fertility traits in Chinese Holsteins. In addition, the present study estimated genetic and phenotypic correlations between heifer traits and between cow traits. The results obtained from the present study will be useful for genetic evaluations in dairy cattle breeding programs, particularly for developing a routine genetic evaluation system for female fertility traits in Chinese Holstein population.
The average performances of fertility traits in the present study were similar to those of other populations in general [6,12,22]. It was also observed that all the fertility performances of heifers were better than those of cows, consistent with previous studies [23]. The differences in the reproductive physiology between maiden heifers and cows could be caused by first calving and the negative energy balance resulting from high milk yield [13]. In addition, the larger amount of censored data for cows (13.92%) than for heifers (5.92%) could imply that cows had poorer fertility conditions than heifers.
The heritabilities obtained in the present study were consistent with the previous estimates for fertility traits in various populations [6,12,24]. For example, Veerkamp and Beerda [25] reviewed the heritabilities estimated in 17 studies and reported that the average heritability estimates for IFL, CR, AIS, NRR56, ICF and DO were less than 5%. Similarly, in the present study, the estimates of heritabilities of IFL, CR, AIS and NRR56 in various parities ranged from 0.005 to 0.034. Furthermore, the estimates of heritabilities of ICF (0.069-0.102) and DO (0.037-0.049) were higher than those of the other traits, which was consistent with previous studies [12,24]. IFLh interval from first to last insemination in heifers, CRh conception rate at first insemination in heifers, AISh number of inseminations per conception in heifers, NRR56h non-return rate within 56 days after first insemination in heifers Table 5 Genetic (above the diagonal) and phenotypic (below the diagonal) correlations between cow traits using a two-trait model For traits defined in both heifers and cows, variance components for IFL and AIS differed considerably between heifers and cows, and heritability estimates of the two heifer traits were much lower than those of the corresponding cow traits. These results were consistent with previous observations in Canadian Holsteins [6], Brown Swiss [23] and US Holsteins [26]. Moreover, genetic correlations between the performances of heifers and cows were moderate but far from one, consistent with previous studies [6,14,23,26]. These results further indicate that traits in heifers and cows should be treated as different but genetically correlated traits in genetic evaluations. In a breeding program, both heifer and cow traits should be considered. In addition, information for heifer traits is useful for early selection, as heifer traits can be measured during a relatively early period of life.
Genetic correlations between the performances of a trait across parities of cows were high (over 0.87) with the exception for CRc, for which the genetic correlation between the first and second parities was 0.66 and between the second and third parities was 0.63 but with a relatively large standard error. The results were consistent with previous literature [15,23,27]. Despite the high genetic correlation, the additive genetic variances of ICFc and DOc in the first parity were substantially larger than those in later parities. The additive genetic variance of IFLc in the first parity was also clearly larger than those in later parities. The differences in additive genetic variances led to differences in heritabilities between the first parity and the other two parities, but the differences in heritabilities were smaller than the differences in additive genetic variances. Differences in additive genetic variances and heritabilities between parities of cows for fertility traits have also been reported in previous studies [15,23]. The results from the present study suggest that even if there is a high genetic correlation between parities of cows, the performance of a trait in the first parity should be treated different from those in later parities in genetic evaluations, especially for ICF, DO and IFL. To increase reliability and reduce bias in routine genetic evaluations for female fertility traits, Nordic countries have treated parities from 1 to 3 as separate traits instead of repeated observations since May 2015 (http:// www.nordicebv.info).
Although the genetic correlations among traits varied substantially, moderate or strong correlations were obtained between traits in groups with similar definitions of fertility functions. Therefore, four groups of traits defined based on reproductive biology highly consistent with the patterns of genetic correlations. The first group of traits measured the abilities of the maiden heifers to conceive and maintain pregnancy, including IFLh, CRh, AISh and NRR56h, which had genetic correlations ranging from 0.63 to 0.94 in absolute values. The second trait group measured the abilities of lactating cows to recycle after calving, indicated by ICF in the present study. The third group of traits described the capacities of cows to conceive and maintain pregnancy, including IFLc, CRc, AISc and NRR56c. Most genetic correlations between these traits were moderate to high, ranging from 0.51 to 0.95 in absolute values, except for the correlations of NRR56 with IFLc (−0.15) and CRc (0.31). The fourth group of traits was DO, which measured the combined abilities to recycle, conceive and maintain pregnancy. This grouping of traits is consistent with the current Interbull concepts of grouping fertility traits.
The traits IFL, CR, AIS and NRR56 all reflect the ability to conceive and maintain pregnancy. According to the estimated genetic parameters, the expected responses to selection given the same selection intensity can be calculated. Among these traits, selection for IFL would achieve the largest genetic progress in the ability to conceive and maintain pregnancy in both heifers and cows. Since IFL had a very high genetic correlation with CR and AIS, it is not necessary to include CR and AIS in a fertility index that contains IFL. NRR56, which measures the capacity for conception and early embryo loss, had a low genetic correlation with IFL for cows. It would be reasonable to include NRR56 in the fertility index. There is no need to include DO in the fertility index given its component traits, ICF and IFL, are already in the index. Thus, IFL, NRR56 and ICF are recommended to be included in the selection index of the Chinese Holstein population.
The present study also estimated genetic and phenotypic correlations between a pair of traits in each parity, using a two-trait model (results not shown). Compared with the results from the analysis using pooled data over parities and treating them as repeated records, the patterns of estimated correlations in each parity were somewhat different from those estimated using the repeatability model (model 5). Moreover, the estimated correlations were not stable across different parities. One possible reason could be bias caused by selection on previous parities, especially the selection for heifer fertility. To avoid this bias, a multi-trait model, which used data from all parities and treated fertility traits from different parities as different traits, was implemented to estimate correlations between a pair of traits for each parity. However, this multi-trait model had difficulty reaching convergence. Therefore, only correlations estimated using data pooled over parities and treating them as repeated records were presented in this study.
For traits with binary phenotypes, a linear model was used in the present study due to low computational requirements. The results showed that the heritability estimates of binary traits were generally lower than those of interval traits measuring conception abilities, consistent with previous studies [6,8]. Many sophisticated models, such as the threshold model and the generalized linear mixed model, have been proposed for the genetic evaluation of categorical fertility traits. Differences between linear and non-linear models were negligible when phenotypes were polychotomous or the incidences of the binary response were not extreme, e.g. between 25 and 75% [28]. Sun and Su [7] reported that probit and logit models performed slightly better than linear models in terms of the prediction of breeding values for non-return rate and success in first insemination, but the correlations between estimated breeding values from different models were high (up to 0.99). Furthermore, a study assessing NRR56 and ICF using a bivariate threshold-linear and a bivariate-linear model showed that no difference was observed between two models in predictive ability [29]. Nevertheless, due to its intensive computational costs and possible numerical difficulties, the threshold model was not recommended for practical use [29]. Moreover, a simulation study comparing computing aspects of linear and non-linear models illustrated that the computing cost of a threshold model analysis was three to five times larger than that of a linear model [7]. Hence, linear models are widely used in practical routine genetic evaluations of categorical fertility traits in many countries, which is also the reason that a linear model was used in the present study. In fact, four types of generalized linear mixed models, including a probit animal model, a probit sire model, a logit animal model and a logit sire model, were also applied for CR and NRR56 (results not shown). The estimates of heritability from these models were similar with those of the linear animal models. Therefore, only the results of the linear animal models were presented here.
One large challenge in genetic evaluations of fertility interval traits is the presence of censoring, which reflects data truncated through culling or fertility events still in progress. Regarding censored data as uncensored or simply excluding them from the dataset would favor bulls with many daughters having censored records [30], and would fail to account for selection bias [14]. A simple strategy to reduce bias is to add a penalty to censored records. In previous studies of Danish Holsteins, a penalty of 21 days was added to censored ICF, IFL and DO, and a penalty of 1 was added to censored AIS, assuming that cows failing to become pregnant would conceive when provided an extra estrus cycle [19,24,31]. In the present study, a penalty of the expected value for a fertility trait was added to a censored record, assuming that the success of an insemination was independent of the result of previous inseminations. The methods for adding a penalty to censored records are simple and offer a potentially sub-optimal strategy to manage censored data. Other methods, such as proportional hazard models [24,32], may be more efficient for handling data that include censored records. However, such models also have some limitations, e.g., high computational demand or extra work to obtain animal model solutions. Moreover, choosing the optimum editing criteria using data truncation with different threshold limits could also be an alternative strategy for dealing with censored data [30].

Conclusions
The parameters estimated in the present study will facilitate the development of a genetic evaluation system for fertility traits to improve the female fertility abilities of Chinese Holsteins. The performances of a trait in heifers and cows were significantly different in genetic and should be considered as distinct but correlated traits in genetic evaluations. In addition, even though high genetic correlations between the performances of a trait in different parities of cows, additive genetic variances and heritabilities differed between the first parity and the later parities, which suggests that the performances of a trait in the first and the later parities should be treated as different traits in genetic evaluations. Generally, the estimated heritabilities of all female fertility traits were low, but selection to improve fertility is promising since the genetic variations are large. IFL, NRR56 and ICF are recommended to be included in the selection index of the Chinese Holstein population.