Genetic analysis of twinning rate in Israeli Holstein cattle

Summary - Second and third parity twinning rate of Israeli Holsteins was analyzed by linear (LM) and threshold models (TM). Data were 124 553 calving records of daughters of 179 sires. Twinning rates were 4.8 and 6.9% for second and third parity, respectively. Heritabilities, as estimated by REML for the LM analysis, and the counterpart of REML for the TM analysis, were 2.2 and 10.1%. The correlation betwen LM and TM evaluations was 0.98. Both distributions of sire evaluations were positively skewed, but only the LM distribution differed significantly from normality. Heritability as estimated from the maternal grandsire effect on twinning rate, correlations between sire and maternal grandsire evaluations, intraclass correlations among half-brothers, and son-sire regressions were all consistent with the hypothesis of polygenic additive inheritance. Sire evaluations for twinning rate were economically favorably correlated with evaluations for dystocia. Breeding for twinning rate may be feasible by selection of sires with high repeatability evaluations, and will not result in significant undesirable correlated responses. cattle / twinning / threshold model / selection / additive inheritance

Twinning is a discrete observation with a binomial distribution. The heritability of twinning has been estimated on observed frequencies (Gaillard, 1969;Bar-Anan and Bowman, 1974) and on probit and arcsine values (Dempster and Lerner, 1950;Van Vleck, 1972) leading to an overall value of : 3%. The threshold model, which assumes an underlying normal distribution, has been applied to other dichotomous cattle traits such as dystocia and calf mortality (Gianola and Foulley, 1983;Weller et al, 1988;Weller and Gianola, 1989). It would seem that this model should also be appropriate for analyzing twinning rate.
Observations of cows with an exceptional rate of twinning and of sires with an exceptional prepotency for twinning have been reported (Bar-Anan and Bowman, 1974;Maijala and Syvajarvi, 1977;Morris, 1984). These were followed by attempts to select for twinning on the basis of high initial twinning rate of foundation dams and bulls (Morris, 1984). In a recent experiment, the first generation of daughters had a twinning frequency of 12.8% (Foulley et al, unpublished results, 1987).
The mode of inheritance of twinning in cattle has not been established. Two reports hypothesized that a single gene may be responsible for twinning in the New Zealand milking shorthorn (Morris and Day, 1986) and in the Israeli Friesian (Bar-Anan and Bowman, 1974). However, in an analysis of one million calvings in dairy cattle, Syrstad (1984) found no evidence to support this hypothesis. In sheep, a single dominant gene was shown to affect fecundity (Davis et al, 1982;Piper and Bindon, 1982).
The objectives of this study were to analyze twinning rate in the Israeli Holstein dairy population by linear and threshold models, and to clarify the mode of inheritance of twinning, and the genetic associations of twinning rate with other traits of economic importance.

MATERIALS AND METHODS
Records of second and third parity calvings of Israeli Holsteins from 1980 through 1987 were analyzed. First parity records were deleted because of the low incidence of twinning (Bar-Anan and Bowman, 1974). Records were deleted from the analysis if: 1) sire of cow or calf was unknown; 2) maternal grandsire (MGS) of the cow was unknown; 3) cow's freshening date or birthdate was unknown; 4) sire of cow had < 100 valid daughter records; and 5) MGS of cow had < 100 valid granddaughter records. The edited data set included 124 553 records, of which 77 631 (62%) were second parity calvings. Mean calving intervals between first and second, and second and third parity were 385 and 366 d, respectively. Basic statistics of the data set are given in table I.
Variance components and sire evaluations were computed by both linear model (LM) and threshold model (TM) analyses (Gianola and Foulley, 1983), using the algorithm of Misztal et al (1989). Variance components were estimated by REML for LM, and by the counterpart of REML for the TM analysis. The analysis model was as follows: where 5 i ;ki is 0 (single) or 1 (twin) calving of the lth cow in parity i, in herd-yearseason j, daughter of the kth sire; Pi is the effect of parity i(i = 2,3); HYS J is the effect of jth herd-year-season; S k is the effect of sire k; and ei jkl is the random residual associated with each record. HYS, sire and residual effects were random, while the effect of parity was fixed. Two calving seasons were defined within each year; calvings from March through September, and calvings from October through February. Since there were 7684 HYS, this effect was absorbed, and its variance component could not be estimated. Therefore the HYS variance component was arbitrarily set as 0.1 of the residual variance. Twinning rate was analyzed by a second model, in which the sire effect was replaced with the MGS effect. Sires and MGS were assumed to be unrelated in the analyses. Since no group of sire or MGS effects were included in the analysis models, the expectation of the evaluations was 0 for all analyses. For TM, 2 rounds of Fisher-scoring iteration were performed prior to variance component estimation by a modification of the expectation-maximization (EM) algorithm, and between each step of variance component estimation. Two rounds of EM iteration were performed between each step of Fisher-scoring iteration.
Estimation of LM variance components was also by EM. Iteration was continued for both methods until the change in the ratio of residual to sire variances was <1% of the previous value. At least 10 rounds of Fisher-scoring iteration were performed for the TM analyses. The heritability estimates for twinning rate were computed as 4 times the sire variance component and 16 times the grandsire variance component divided by the respective phenotypic variances. Phenotypic variances were computed as the sum of the sire (or MGS), HYS, and residual variance components. Variance components and solutions from the TM analyses are in arbitrary units. To facilitate comparisons between the 2 analyses, units in the TM analysis were set so that the residual variance component would be equal to the residual variance component of the corresponding LM analysis. Heritability on the underlying scale (h') was also estimated by the following equation (Dempster and Lerner, 1950): where hr is the estimate of heritability from the linear model analysis, p is the proportion of twins in the data set analyzed, and z is the normal ordinate for p.
Repeatabilities of sire evaluations from the LM analysis were computed as (Var Sire -PEV)/(Var Sire), where Var Sire is the REML estimate of the sire component of variance, and PEV is the prediction error variance, computed as the diagonal element of the inverse of the coefficient matrix. Repeatabilities of MGS evaluations were computed in ' a similar manner, with Var Sire replaced by the MGS variance component.
The Kolomogorov D statistic was used to test the deviation of the distributions of evaluations from normality (Snedecor and Cochran, 1967) and the g i statistic to test for skewness (Sokal and Rohlf, 1969). Intraclass correlations (t) of evaluations among half brothers and the regressions of evaluations of sons on sires (b) were computed.
Product moment correlations were computed between sire evaluations for twinning rate and for annualized milk, fat, and protein production [365 (total lactation yield/calving interval)], fat and protein percent, conception rate, dystocia, and calf mortality. Evaluations for dystocia and calf mortality were computed both for the sire of the calf and the cow calving. Only first parity single calvings were included in the analyses, with negative values indicating low incidence of dystocia. The sire evaluations for these traits were computed as described by Bar-Anan et al (1987), and Weller et ad (1988), and the data sets analyzed are given in Weller and Ron (1989).  (Weller et al, 1988;Weller and Gianola, 1989). The estimates of heritability on the underlying scale, derived from the equation of Dempster and Lerner (1950) are listed in parentheses after the TM estimates. These values are slightly lower than the direct TM estimates, and correspond to previous results for analysis of calving traits (Weller et al, 1988). For both the TM and LM analyses, the heritability estimates derived by the sire and MGS analyses were very similar.

RESULTS
The distributions of sire evaluations are given in figure 1 for the LM analysis, and in figure 2 for the TM analysis. The statistics for these distributions are given in table III. The ranges and standard deviations (SD) of evaluations for TM were twice the ranges and SD for the corresponding LM analyses, even though the TM evaluations were scaled to equal residual variances. As expected from the polygenic additive model of inheritance, the ranges for the sire evaluations were more than twice as large as the corresponding MGS ranges. Assuming additive genetic inheritance and equal repeatabilities, the SD of the sire evaluations should be twice that of the corresponding MGS evaluations. In fact, the SD of the sire evaluations were more than twice the corresponding MGS evaluations. This is not surprising, since the mean repeatability of the sire evaluations was nearly twice the mean repeatability of the MGS evaluations. The LM distributions were highly skewed, as evident both from figure 1 and the 91 statistics. Although g i values were lower for the TM analyses, both were positive, and skewness was significant for the sire distribution at p < 0.01. Skewness for the distribution of sires with repeatability > 65% was only marginally less. Both of the LM distributions deviated significantly from normality, as estimated by the D statistic, but the TM distributions did not.
The sire Shafan, No 781, had the highest sire evaluation for twinning rate; 5.9% and 9.1%, by the LM and TM analyses, respectively. The twinning rate of twinning rate of his granddaughters was 6.7%. His MGS evaluations were 0.5% and 0.8%, for LM and TM, respectively, but these evaluations are based only on 388 granddaughters (repeatability = 30%). Thus it is evident that this particular sire could be used to breed for increased twinning rate.
The correlations among TM and LM sire and MGS evaluations are presented in  (Weller et al., 1988;Weller and Gianola, 1989), correlations between TM and LM evaluations were greater than 0.95. The correlations between sires and MGS evaluations were 0.41-0.42.
The expectation of these correlations for evaluations with complete repeatability is 0.5. Therefore the results obtained for evaluations with moderate repeatability is in accord with the hypothesis of additive inheritance.
Coefficients of intraclass correlations and son-sire regressions of evaluations are presented in table V. Both the regressions and correlations were close to the expected values with complete repeatability. These results also support the hypothesis of polygenic additive inheritance.
Correlations between evaluations for twinning rate and dystocia are given in table VI. All other correlations between twinning and other traits of economic importance were within the range of -i to 1, and not significant. Correlations between sire evaluations for twinning rate and the sire of cow evaluation on dystocia were significant at P < 0.05. As expected, the TM and LM evaluations had very similar correlations with both dystocia traits. Although the correlation between the sire of calf evaluations for dystocia and twinning rate was slightly higher than the correlation between the sire of cow evaluations for dystocia and twinning rate, only the latter was significant, because of the greater number of sires with evaluations. These correlations are in fact positive on the economic scale, since negative values are favorable for dystocia.

DISCUSSION
The mean rate of twinning in the dairy cattle population of Israel is 5.6%, compared with 4.5% in this population 20 y ago (Bar-Anan and Bowman, 1974). This may be due to a difference in data recording or editing. However, within this period the milk yield per cow in Israel has increased by > 50% (Kalay, 1986).
The increase in twinning frequency may be due to improved feeding, therapy of anoestrous cows, and selection for milk. A positive genetic correlation between twinning rate and milk production was found in previous reports (Morris, 1984), but not in this study.
Linear model heritability estimates for twinning based on daughter and granddaughter groups are in accordance with many other studies (Morris, 1984). The 5-fold increase in heritability found in the threshold model analyses agrees with results for other dichotomous traits (Weller et al, 1988, Weller and.
The correlations between sires and MGS evaluations, the intraclass correlations, and the son-sire regressions all support the hypothesis of polygenic additive inheritance. This genetic variation may be used for selection.
Previous studies documented exceptional sires for twinning; daughters of the sires Zemed (Bar-Anan and Bowman, 1974), and Tahto (Maijala and Syvajarvi, 1977) had 11% and 12.2% twinning, respectively. A segregating major gene was proposed to explain these results. In the current study the twinning rate of 6 112 daughters of the sire Shafan was 10.6%. His evaluation for twinning rate was 4 SD units above the mean in the LM analysis, but only 3 SD units greater than the mean in the TM analysis. Under the .assumption of normality, the probability of observations > 3 SD units is 0.0013. Thus in sample of 179 sires, the expected frequency of sires > 3 SD units is 0.0013 (179) = 0.27. Thus a single observation in this range is consistent with the hypothesis of polygenic additive inheritance.
Although skewness was lower for the TM than for the LM distribution of evaluations, both the TM distributions for all sires, and for high repeatability sires were still significantly skewed. Deletion of Shafan reduced the g l value to 0.3, a value on the border of significance. Four other high repeatability sires, out of a total of 51, had evaluations > 2 SD units above the mean. These results are at variance with the threshold model assumption of an underlying normal distribution. Furthermore, in the previous analysis of calving traits in this population, none of the TM distributions were significantly skewed, even though the LM distributions were (Weller et al, 1988). Misztal et al (1988) demonstrated that if the distribution of the underlying variable is not normal, solutions for fixed and random effects can still be calculated as if normality held. Thus even in the present case TM is superior to LM.
The distribution of TM sire evaluations could be explained by a segregating gene with a substitution effect of ! 1 genetic SD unit and unequal allelic frequency. Since all analyses were based on sire or MGS evaluations, the possibility of a major recessive gene with a low frequency allele also cannot be excluded.
Sire evaluations for twinning rate were favorably correlated with evaluations for dystocia, in accordance with previous results on this population (Bar-Anan, 1971). Correlations between twinning rate and production traits were not significant. Because of low heritability, selection for twinning will not be practical on the sire-todam path. However, since sires of sires are often selected among high-repeatability proven sires, selection for twinning may be practised in the sire-to-sire path. Weller (1989) found that selection based on an index including milk production and fertility could increase conception rate by 5% in 10 y, but would result in an 8% reduction in the genetic gain for economically fat-corrected milk. Using the same method, and assuming the same reduction in genetic gain for production traits, selection for twinning rate in the male pathways would result in a genetic increase of 1.2% over this time period. In conclusion, the preponderance of evidence supports the hypothesis of polygenic additive inheritance for twinning rate, although the possibility of segregating major genes cannot be excluded. Breeding for twinning may be feasible by selection of sires with high repeatability evaluations, and will not result in significant undesirable correlated responses.