Non-additive genetic effects for fertility traits in Canadian Holstein cattle (Open Access publication )

The effects of additive, dominance, additive by dominance, additive by additive and dominance by dominance genetic effects on age at first service, non-return rates and interval from calving to first service were estimated. Practical considerations of computing additive and dominance relationships using the genomic relationship matrix are discussed. The final strategy utilized several groups of 1000 animals (heifers or cows) in which all animals had a non-zero dominance relationship with at least one other animal in the group. Direct inversion of relationship matrices was possible within the 1000 animal subsets. Estimates of variances were obtained using Bayesian methodology via Gibbs sampling. Estimated non-additive genetic variances were generally as large as or larger than the additive genetic variance in most cases, except for non-return rates and interval from calving to first service for cows. Non-additive genetic effects appear to be of sizeable magnitude for fertility traits and should be included in models intended for estimating additive genetic merit. However, computing additive and dominance relationships for all possible pairs of individuals is very time consuming in populations of more than 200 000 animals.


INTRODUCTION
Genetic evaluations of dairy cattle are generally based on additive genetic models as either sire or animal models [8,12]. Total genetic values of animals may also contain non-additive components. Non-additive genetic variation, namely dominance and epistasis, are the result of interactions of alleles and loci. Interactions of alleles at the same locus result in dominance effects and interactions of two alleles at different loci result in additive by additive variation. Interactions of more than two alleles in different loci result in different levels of epistatic variations such as additive by dominance and dominance by dominance, and many others [19]. Although non-additive genetic effects are not directly transmitted from parents to offspring, they are important for traits closely related to fitness or having low heritability [2]. If non-additive genetic effects are significant, then ignoring them in genetic evaluation models could bias estimates of additive genetic effects.
Selection to improve production and conformation traits [11] has led to a decline in some fertility traits due to negative genetic correlations with production [1,14,18]. Fertility traits generally have low heritability and could have large non-additive genetic effects. Canada has recently adopted new genetic evaluations for a number of fertility traits [9], and part of that project had the objective to study non-additive genetic effects.
In order to separate non-additive genetic variances from additive genetic variance, data should contain individuals having non-zero dominance relationships to each other. Van Tassell et al. [21] suggested a minimum of 20% full sibs in the population to be successful in estimating non-additive genetic variation. Extensive utilization of multiple ovulation and embryo transfer (MOET), embryo splitting and nuclear cloning have produced groups of close relatives which share additive and non-additive genetic effects [20]. Misztal et al. [13] pointed out that ignoring non-additive genetic effects for animals with dominance relationships tended to inflate additive genetic variances. Even though the number of full sibs may be small compared to the total population, such animals may be the elite animals of the population, and their influence could be disproportionally greater.
The objective of this study was to estimate several non-additive genetic variances including dominance (D), additive-by-additive (AA), additive-bydominance (AD) and dominance-by-dominance (DD), together with the additive genetic variance (A), for age at first service and non-return rate for heifers, and non-return rate and interval from calving to first service for cows.

Pedigree and fertility information
A pedigree file containing 1 811 940 Holsteins, with the eldest born in 1961 and the youngest born in 2002, was obtained from the Canadian Dairy Network. Animals with unknown parents were deleted (11%), and sires with only one progeny (15 821) were deleted (0.8%). Sire-by-dam combinations appearing more than once gave 47 104 full sib families ranging in size from 2 to 14 animals per family, with two animals per family as the most frequent. A total of 106 640 full sibs were found in 1 606 468 animals or 6.6%. A total of 22 611 animals in 106 640 had the same birth year within family, a group which likely included natural twins and full sibs from embryo transfers. Fertility records on 1 376 934 Holsteins were obtained from the Canadian Dairy Network. Records were divided into heifers (486 012) and cows (507 351). A total of 243 907 heifers appeared also as a cow. Traits recorded for heifers were age at first service (AFS) and non-return rate (NRRH). Traits for cows were non-return rate (NRRC) and interval from calving to first service (CTFS), both with repeated records for cows. Descriptive statistics are in Table I. The 106 640 full sibs from the pedigree file were matched with the fertility data giving a total of 49 629 full sibs with fertility records. All herdmates of those full sibs were retrieved (154 387 animals), giving a total of 203 400 animals (138 596 heifers and 64 804 cows) and 389 085 records.
The total number of animals in the pedigrees for the 203 400 animals was 598 230. Among the 203 400 animals with data, 18 782 (9%) appeared to be possible embryo transfers and among those 15 991 were raised in the same herd.

Construction of additive and dominance matrices
The intention of the study was to calculate and save all of the additive and dominance relationships (using the algorithm of Smith and Maki-Tanila, [16]) for the 598 230 animals in the pedigree file. The intention was to proceed as in Schaeffer [15] using shortcuts to avoid inversion of the non-additive covariance matrices. However, with 598 230 animals in the pedigree file there would be 3.56 × 10 11 possible additive and dominance relationships to calculate. This would have required many terabytes to store the non-zero coefficients as well as several years of computing time. Therefore, the data file was split into four groups of about 50 000 animals each plus their complete pedigree information. A genomic matrix was constructed for each group and additive and dominance matrices were obtained (approximately one month to calculate all relationships, about one terabyte of storage space, and coefficients were stored in binary format to save space). Additive and dominance matrices were checked for the average number of non-zero off-diagonal elements. The percentage of non-zero off-diagonals for dominance matrices for all groups, showed an overall average of only 14%. Due to the very small percentage of non-zero offdiagonals, some non-additive matrices such as dominance-by-dominance and additive-by-dominance would have been equally sparse. Processing the relationships many times, as would be necessary with Gibbs sampling or iteration on data, would have required many months of computing. Thus, another strategy had to be found.
Van Tassell et al. [21] recommended that about 20% of the population should be full sibs in order to separate non-additive variances from the additive variance. In order to increase the percentage of non-zero dominance relationships in the covariance matrix, animals were sampled as follows. Animals were divided into heifers and cows. Two subsets of 1000 heifers and two subsets of 1000 cows were extracted from the latest years from each of the four groups of 50 000 animals. To be in a subset an animal had to have fertility data and a non-zero dominance relationship with at least one other animal in the subset. The additive and dominance relationships among the 1000 animals were retrieved from the already computed relationships in the four groups of 50 000 (plus pedigree animals) from which they were sampled. This gave an additive relationship matrix (A) of order 1000 and a dominance relationship matrix (D) of order 1000. The AA, AD, and DD matrices were obtained as the Hadamard products of the A and D matrices for the 1000 animals. The matrices included inbred animals and appropriate inbreeding coefficients in A and D [16]. After the sampling, the percentage of non-zero off-diagonal elements in A and D matrices ranged from 31% to 50% for dominance matrices and over 69% to 96% for additive matrices. Additive and dominance matrices for all subsets appeared sensibly fuller with non-zero off-diagonals, implying that epistatic effect matrices would be fuller with non-zero off-diagonals. The size of these matrices allowed direct inversion and the procedures given by Henderson [6] could be followed directly.

Models and methodology
Each fertility trait had a different model. Besides A, D, AD, DD, and AA genetic effects in each model, other effects included region-year of birth-season of birth (RYS, fixed effect), month of insemination (M, fixed effect), AI technician (T, fixed effect), service sire by year of insemination (SY, random effect), parity-age of calving-month of previous calving (APC, fixed effect), parity-age of calving-month of first insemination (API, fixed effect) and permanent environmental effects (for cow traits) (PE, random effect). Models for each trait are in Table II.
Due to the method of sampling subsets of 1000 animals, herd effects were not modeled because complete contemporary groups were not selected and slightly more than 500 herds were represented in each 1000 animal datasets, giving approximately two animals per herd, none of them belonging to the same full sib family. The assumption was that variability due to herd effects would go into the residual variance and would not bias the other parameter estimates.
Additive and non-additive variances were estimated applying Bayesian estimation via Gibbs sampling. Additive, dominance, epistatic, permanent environmental and residual effects were assumed to follow normal distributions, and their variances were assumed to follow inverted Chi square distributions. In order to estimate each gene interaction variance, a preliminary run, per trait per group was made, where each component was estimated one at a time, starting with the additive effect. To avoid the influence of the starting values on inferences and due to the small number of animals considered, priors were entered with degrees of belief. In each preliminary run, each effect was given a degree of belief equal to zero, while previously estimated interaction variances were given 500 degrees of belief. In this way the prior value needed to start sampling, was given zero credibility letting the data provide all information needed to reach the posterior distribution. The high degrees attributed to previously estimated variances influenced their new estimation to remain close to their estimated values allowing the new parameters to be separated from the residual. After all five gene interaction variances were estimated, a final run with all five parameters together, using the preliminary estimates as priors each with one degree of belief, was made. All gene interaction variances were estimated together with equal degrees of belief to see how much individual estimation changed when all parameters were estimated simultaneously. The low degrees of belief were necessary to get as much information as possible from the data structure, avoiding the influence of the priors. For each run a total of 40 000 samples were made. Burn-in was achieved after 25 000 samples based upon visual assessment of sample value behavior. All samples after burn-in were used to calculate means and standard deviations of estimates for each component. The estimates from the eight subsets for heifers and eight subsets for cows were weighted by their posterior standard deviations and averaged. Phenotypic variance was the sum of all variances and the ratios of each component to the phenotypic variance were calculated.

RESULTS
Estimates of non-additive genetic variances, permanent environmental and residual variances are in Table III for all four fertility traits. Estimates of variances for the service sire-year of insemination effects for non return rates were not presented since these were very close to zero and represented less than 0.1% of the total phenotypic variance.
In general, non-additive variance estimates were as large as or greater than the additive variance estimates. Posterior standard deviations of estimates were almost as large as the estimates, which reflects the small subset sizes and low heritabilities for these traits.
Dominance genetic variances were greater than the additive genetic variance for age to first service, heifer non return rate, and interval from calving to first service. All non-additive components were greater than the additive variance 10 Only additive h 2 ** 0.14 * As found by Jamrozik et al. [9] using an additive animal model. ** Results from our analysis including only additive genetic effect.
for age to first service and heifer non-return rate. Permanent environmental variances for cow traits were very low. There were roughly 1.5 records per cow in the eight cow subsets, and only single records for heifers.

DISCUSSION
The results for fertility traits were similar to those of Miglior et al. [10] regarding production traits. Miglior et al. [10] found AA and D estimates that were as large as or larger than the A component estimates, indicating that nonadditive variances may be larger than commonly believed. Hoeschele [7] found Restricting days open to no more than 150 days reduced the D variance. Fuerst and Solkner [5] showed that D and AA components were greater than the A component for calving interval. The interval from calving to first service is a component of both days open and calving interval, and showed similar results to [5] and [7] with the addition that the AD variance estimate was larger than the other non-additive components.
Comparing genetic variance estimates between heifer and cow non-return rate, non-additive genetic variance estimates were similar in value, but the additive component was much greater for cows than for heifers. The heritabilities in the narrow sense for heifers and cows were 0.005 and 0.11, respectively. This may be due to the lower mean non-return rate for cows compared to heifers. Another possible reason may be that non-return rate in cows is influenced by other factors that regulate ovarian activity and may have a heritability value greater than that of non-return rate.
Permanent environmental effects for both cow traits were low. The low estimates of permanent environmental variances were a direct consequence of the limited number of repeated records per animal in analysis (1000 cows with an average of 1.5 records each). Nevertheless, similar estimates of permanent environmental variances were also found by Fatehi et al. [3]. Fatehi et al. [3] used a single trait animal model on age at first service, non-return rate, for heifers and cows, and calving to first service interval on 477 748 Holstein animals. Thus, the small subset sizes were not a cause for this result.
Heritability in the broad sense was calculated as the ratio between the sum of additive and non-additive variances over the phenotypic variance of the trait. Broad sense heritability was high for all traits highly influenced by nonadditive genetic effects, i.e. age at first service (0.74) and interval from calving to first service (0.34), even for non-return rate for heifers heritability in the broad sense, although low in value, was higher than the additive genetic effect (0.049). Heritability in the broad sense was only slightly higher than additive genetic effect for non-return rate for cows (0.13) where all non-additive effects were very low and almost negligible in value. In order to get an estimate of broad sense heritability, genetic evaluation must account for non-additive genetic effects in estimating breeding value. Broad sense heritabilities, as given in this research, are intended to be general, providing an idea of the effect of the inclusion of non-additive genetic effects and are not intended to be fully exhaustive, since non-additive genetic covariances were not taken into account.
Heritability in the narrow sense (i.e. additive genetic variance to phenotypic variance) was lower when accounting for non-additive genetic effects than using an additive animal model. This phenomenon has been already reported by several authors [3,5,7,10,13] and confirmed by our analysis, Table III. Whenever gene interactions are omitted from the model their variance gets split between the additive and the residual effect therefore determining the additive effect to be overestimated.
An analysis including a regression on inbreeding was made utilizing data, as well as additive and dominance relationships, from the second subset for both heifers and cows. Animals belonging to this group showed an average inbreeding level of 1.5%. The results from this analysis revealed that inbreeding depression was quite important for the two temporal traits, age at first service and calving to first service, accounting for 1 and 0.89 days per 1% increase in inbreeding, respectively. The effect of inbreeding depression for the heifer and cow non-return rate was, instead, negligible and very close to zero (−0.0002 and −0.006, respectively). Similar values of inbreeding depression were reported in different studies [5,7,10,11].
A comparison of variance estimates with and without regression for inbreeding is available in Tables IV and V. Accounting for inbreeding depression reduced posterior standard deviations in magnitude for all four traits and parameters. Estimates from the two runs were very close in value for age at first service and calving to first service while estimates for non-return rates for heifer and cow showed a slight difference. The same trend of relevance for non-additive genetic components, for all traits, was still found even including a regression on inbreeding. Due to the smaller value of posterior standard deviations obtained, inclusion of inbreeding depression is recommended.
There were no previous studies of the four fertility traits, and none for the levels of non-additive genetic effects considered in this study. The sampling of subsets of 1000 animals to increase the percentage of non-zero relationships in the D matrix may have generated some bias, although results obtained for dominance and additive by additive epistatic effects were similar to those available in the literature [5,7,10,11]. Not accounting for herd effects in the subsets may have had some effect on estimates. One possibility would be to estimate the herd effects in an additive genetic model for the four groups of 50 000 animals and then use deviations from the herd solutions as the observations in the groups of 1000 animals.
The consequences of this study on genetic evaluations for fertility traits, and maybe other traits, are that the ratio of the variance explained by nonadditive genetic effects to phenotypic variance appears larger than heritability in the narrow sense for age at first service, heifer non-return rate and calving to first service. To ignore non-additive genetic effects may cause additive genetic effects to be overestimated and possibly biased, as seen by comparison of the results in Table III with Jamrozik et al. [9] and with several studies on this matter [3-5, 10, 11, 13, 17]. However, the calculation and storage of non-zero additive and dominance relationships for reasonably large data sets will always be a bottleneck because there are no easy ways to calculate the inverses of the non-additive genetic covariance matrices as there is for the additive genetic covariance matrix. Applying the model of Henderson [6] assumes the population is randomly mating, and this is clearly not the case. Thus, with non random mating there are covariances created between additive and dominance genetic effects and most likely between other epistatic effects [7]. For these reasons non-additive genetic effects will likely continue to be ignored in genetic evaluation models. Studies should be conducted to determine the risks associated with that scenario in terms of reduced genetic change. This will likely depend on the percentage of non-zero dominance relationships that exist in the population. As time progresses the proportion of animals with dominance relationships are likely to increase, just as inbreeding increases, particularly in Holstein dairy cattle where the number of effective sires is estimated to be 30 on a worldwide basis. Perhaps the importance of non-additive genetic effects may increase in the next decade such that they cannot be ignored.