Estimates of covariance components and genetic parameters for growth traits in Tswana cattle selected for increased weaning or eighteen-month weights

The objective of this study is to estimate variance-covariance components and genetic parameters for growth traits in Tswana cattle. Genetic analyses for average daily gain (ADGs) and growth traits were conducted using 7223 records of animals which were born between 1996 and 2013 from 1662 dams and 188 sires in 54 contemporary groups. Both univariate and bivariate animal models were used. Heritability estimates for growth traits ranged from 0.12±0.03 for BWT to 0.45±0.03 for EWT while those obtained for ADGs were 0.24±0.03 and 0.31±0.04 for ADG1 and ADG2, respectively. Permanent maternal environmental effects were significant for WWT and ADG1. Substantial maternal genetic effects were observed in BWT, WWT and ADG1. Genetic correlations among growth traits and ADGs ranged from 0.19±0.07 between BWT and ADG1 to 0.99±0.02 between WWT and ADG1. Phenotypic correlations among growth traits and ADGs ranged from 0.19±0.01 between BWT and ADG1 to 0.94±0.01 between WWT and ADG1. The two selected lines significantly improved annual genetic gain for all the growth traits and ADG1. Genetic gain for EWT was optimally enhanced in S2 than in S1. The control line exhibited significant annual genetic gain in WWT which was not anticipated and perhaps attributed to asymmetry of selection response. Substantial genetic variations were observed in all growth traits and ADGs suggesting that growth improvement can be attained through selection for growth rate. High genetic correlations between growth traits and ADGs indicated that selection for one of these traits may result in indirect correlated response on the other traits.


INTRODUCTION
Genetic parameter estimates are necessary for the implementation of genetic evaluations, subsequent assessment and quantification of genetic progress that can be achieved through on-going long term selection in economically based multi-trait breeding objectives (Enns and Nicoll, 2008). Therefore, the importance of genetic parameters such as heritability, correlations and repeatability in traits of economic importance in beef cattle cannot be underrated. On the other hand, these parameters are population specific and hence can by no means be generalised (Irgang et al., 1985). As a consequence, many studies have been conducted in recent years across the world with the aim of estimating both phenotypic and genetic parameters and to evaluate their importance in designing breeding programmes for different beef cattle breeds in different production environments (Chevraux and Bailey, 1977;Davis, 1987;Mercadante et al., 2003;Koch et al., 2004;Yilmaz et al., 2004;Enns and Nicoll, 2008;Boligon et al., 2010;Cervantes et al., 2010;Boligon et al., 2013;Shumbusho et al., 2013). This has led to validation of practically achieved progress with theoretically predicted direct and correlated responses for setting up of more efficient breeding programs. Vleck and Cundiff (1998) and Crews (2006) defined selection of beef animals for breeding to be based on breeding values for growth performance. The genetic potential for early growth traits is known to be influenced by both direct and maternal additive genetic merit (Willham, 1972;Robinson, 1996a). In addition to these genetic effects, Koots et al. (1994) and Robinson (1996b) emphasised that permanent environmental effects also need to be accounted for especially in tropical environments where accidents to the udder and exposure to diseases are likely to occur. Meyer et al. (1993) and Garrick (2010) also asserted that optimization of genetic progress should take the two into consideration. For traits expected to be influenced by maternal effects, genetic evaluation through best linear unbiased prediction (BLUP) involves partitioning of the phenotypic variances and covariances into their direct and maternal genetic, maternal permanent environmental, genetic correlations, and residual components. Due to the differences in the importance of direct and maternal genetic effects between the breeds, several studies have been conducted to estimate genetic and phenotypic parameters among weight traits of different breeds (Knights et al., 1984;Smith et al., 1989;Lee et al., 2000;Forni and Albuquerque, 2005).
In an established breeding programme, drawing a genetic trend is a useful tool to evaluate the impact of genetic improvement programme. The trend has to be monitored to check the validity of the predictions made and the direction of genetic change (Intaratham et al., 2008;Bosso et al., 2009). Currently, nothing has been reported on genetic change due to selection in indigenous Tswana cattle in Botswana. The objectives of this study were therefore 1) to estimate variancecovariance components and genetic parameters for growth traits, 2) to evaluate the importance of maternal effects on growth traits of Tswana cattle, and 3) to estimate annual genetic change in growth traits of Tswana cattle selected for increased weaning or 18-months weight relative to unselected control.

MATERIALS AND METHODS
The experimental ranch is located in Southern District about 180 kilometers in the South west of Gaborone, the capital city of Botswana. It covers an area of about 7628 ha with the recommended stoking rate of 10.3ha/Lsu and receives annual rainfall of 212 to 710 mm. The soil types are non-calcareous fine sandy soils; the vegetation is dominated by woody plants species such as Dichrostachys cinerea, Terminalia sericea, Acacia hebclada, Diospyros lycioides and Grewia flava which are mostly known as encroacher plants. However, the ranch condition has been largely characterized as fair in terms of pasture production due to high presence of palatable grass species such as Digitaria, Cynodon dactylon, Brachiaria nigropedata, Schimidtia pappophoroides and Urochloa trichopus. These are either increaser or decreaser species, and their abundance is said to suggest fair utilization. Animals are allowed to graze on natural pasture in a ranch for the better part of the year and only supplemented during dry seasons such as winter or drought periods following years of low rainfall. The supplementation is done using Urea-molasses block and dicalcium phosphate.
Data of selected Tswana cattle collected from 1995 to 2013, acquired from the Department of Agricultural Research (D.A.R.) in Botswana were used in this study. The data consisted of 2940 records for 7 months selection line (S1), 3034 for 18-months selection line (S2) and 1252 records for the unselected control line (S3). In both S1 and S2, mass selection was practiced based on the animals' weight index calculated at weaning (205 days) and at 18-month (540 days) of age, respectively. Within S1 and S2, indices were calculated based on pre-adjusted 205 and 540 weights, respectively, prior to the implementation of mass selection.
The data used in this study comprised the following: pedigree information, that is, calf identity number, sire identity number, dam identity number and associated important information such as birth date, sex of the calf, and selection line. The growth traits that were considered were birth weight (BWT); weaning weight (WWT) recorded at 7 months; yearling weight (YWT) recorded at 12 months; 18-month weight (EWT); pre-weaning average daily gain (ADG1); and post-weaning average daily gain (ADG2). ADG1 was calculated by dividing the change in body mass from birth to weaning by the number of days in the interval while ADG2 was calculated as the change in body mass from weaning to eighteen months divided by number of days in the interval. Dam age was grouped into three classes or categories namely young (≤5 yr), mature (>5 and ≤9 yr) and old (>9 yr). Weaning, yearling and eighteen-month weights were pre-adjusted to 205, 365 and 540 days respectively, therefore individual age differences at recording were not fitted as an effect in the analysis models. Since seasonal mating was practiced, all animals were born from late September to early January. However, there were fewer birth records in September and January than in other months. As a result, the birth seasons were regrouped as follows: those born in September and October were grouped as season one, those born in November as season two and those born in December and January as birth season three. After regrouping the birth month, contemporary group was then formed by concatenating birth year and birth season, and it was fitted as a sparse fixed effect in a mixed model. The summary of the final data used for genetic analyses in this study is presented in Table 1.
Mixed statistical models were used for genetic analysis of calf weight and growth rate traits. Fixed effects fitted were sex, dam age, line, contemporary group, sex by contemporary group and line by contemporary group interactions. An animal model integrating all pedigree information with performance data was used in the analysis of all the traits. In addition to the direct additive genetic effects, maternal genetic and permanent environmental effects on the trait studied were tested by fitting them as additional random effects. Because of the data size and depth of pedigree structure the correlation between direct and maternal effect was not estimable for all the traits. The full model considered for the analysis of BWT, WWT, YWT, EWT, ADG1 and ADG2 was: Where; Y = the observed growth trait of the calf (BWT, WWT, YWT, EWT, ADG1 or ADG2), X = the incidence matrix relating fixed effects to the observations of a trait, β = vector of fixed effects which were sex, dam age, selection line, contemporary group and their interaction, Z a = incidence matrix relating direct additive genetic effects to the observations of a trait, Z m = incidence matrix relating maternal additive genetic effects to the observations of a trait, Z c = incidence matrix relating permanent maternal environmental effects to the observations of a trait, u a = a vector of direct additive genetic effects, u m = a vector of maternal additive genetic effects, u c = a vector of permanent maternal environmental effects and ԑ = a vector of random residual effects.
The random effects in the mixed models were assumed to have the following variance covariance distributions: [µ d ', µ m ', µ pe ', e']' ~ N [(0', 0', 0', 0')', ∑] and the variance-covariance structure fitted in the analysis of a single trait genetic model was as follows: The variance-covariances estimates obtained from the fitted random effects were expected to appear as follows: Var(y) = σ 2 a + σ 2 m + σ 2 c + σ 2 e , Var(u a ) = Aσ 2 a , Var(u m ) = Aσ 2 m , Var(u c ) = Iσ 2 c and Var(e) = Iσ 2 e . A is the numerator relationship matrix and I is the identity matrix equal to the number of dams for Var(c) and number of animals for Var(e). These variance-covariances were later used to compute phenotypic variance (σ 2 p ), direct heritability (h 2 a ), maternal heritability (h 2 m ), the ratios for permanent maternal environmental effects (c 2 ) and random residual effects (e).
Analysis of data was carried out by fitting univariate and bivariate animal models that ignored one or two of the maternal effects to the data using the ASREML programme (Gilmour et al., 2015). Models were numbered as follows: Model 1: h 2 a ; Model 2: h 2 a and h 2 m ; Model 3: h 2 a , h 2 m , c 2 . The models were then compared using the log likelihood ratio (LR) test, which was calculated as minus twice the difference between the Log L of models under comparison. To declare the difference as significant, the likelihood ratio was compared with the tabulated χ 2 statistics with the number of degrees of freedom taken as the difference in the number of parameters fitted between the two models. Using the best model for each trait, phenotypic and genetic covariances, and correlations among traits were estimated through bivariate analyses.

RESULTS AND DISCUSSION
Based on the LR test, Model 3 was identified to be the best for BWT, WWT and ADG1. Model 2 was suitable for YWT and EWT while Model 1 was suitable for ADG2. Estimates of variance components and genetic parameters for each trait obtained from the three statistical models are presented in Tables 2 and 3, respectively. When one of the maternal effects was included in a model (Model 2 and Model 3) for BWT, WWT and ADG1, both direct genetic variance and the corresponding direct heritabilities were substantially reduced. The estimates obtained for both maternal genetic and environmental variances were significantly different from zero and showed substantial contribution to the overall phenotypic variances. The LR was significantly improved compared to the model fitting only direct genetic effect (Model 1). The maternal genetic variance estimate for YWT obtained from Model 2 was significantly different from zero. However, the estimates obtained for maternal effects when Model 2 and Model 3 were fitted for EWT and ADG2 were not significantly different from zero and did not show any important contribution to the phenotypic variance, despite the improved LR as compared to the model containing only direct genetic effect.
Although the inclusion of one of the maternal effects in Model 2 and Model 3 substantially reduced the direct genetic variances and the corresponding heritability estimates for BWT, WWT, YWT and ADG1 compared to Model 1. The two models did not show any considerable differences between each other in terms of both direct genetic variances and the corresponding heritability estimates. Besides, the magnitude of the estimates for both maternal genetic and maternal permanent environmental effects from Model 2 and Model 3 were more or less the same (Tables 2 and 3). There was a significant improvement on the LR, and a minor change in the magnitude of direct genetic variances and the corresponding heritabilities for EWT when one of the maternal effects was included in the analysis model (Model 2 and Model 3). The LR did not improve when model 3 was fitted for ADG2, while the magnitude of both maternal genetic and maternal permanent environmental variances and the respective heritability estimate from S.E. = standard error, h2a = heritability for direct genetic effects; h2m = heritability for maternal genetic effects; c2 = permanent maternal environmental proportion; ADG1 = pre-weaning average daily gain and ADG2 = post-weaning average daily gain. * = P<0.05 and *** = P<0.001.
both Model 2 and Model 3 were negligible and not different from zero. The range of direct heritability values obtained in the current study were consistent with those reported by Meyer (1997) for BWT and WWT in Australian beef cattle and Demeke et al. (2003) for BWT in Ethiopian indigenous breeds and their crosses. Direct heritability value obtained for ADG2 in the current study was similar to the value reported by Raphaka and Dzama (2010) for the same breed. On the other hand, they reported high direct heritability value for BWT and lower values for WWT, EWT and ADG1 than those obtained in the current study. In general, the current heritability estimates for various growth traits are within the range of values reported for other breeds (Groeneveld et al., 1998;Costa et al., 2011;Boligon et al., 2010). The magnitude of direct genetic variances and associated heritabilities currently obtained were generally high for WWT, YWT, EWT and ADGs which implied that generally there is sufficient genetic variation to implement genetic improvement through well planned selection.
The range of values obtained in the present study for maternal genetic heritabilities are within the range of estimates reported in literature for BWT (Groeneveld et al., 1998;Demeke et al., 2003;Iwaisaki et al., 2005), for WWT (Meyer, 1997;Costa et al., 2011;Boligon et al., 2012), for YWT and EWT (Groeneveld et al., 1998;Boligon et al., 2010;Costa et al., 2011), and for ADGs (Demeke et al., 2003;Raphaka and Dzama, 2010;Boligon et al., 2012). Permanent maternal environmental proportions from Model 3 were not significantly different from zero for BWT, YWT and EWT. The results showed that both maternal genetic and maternal environmental variances decreased by a large amount from WWT to EWT. However, a notable change observed on the magnitude of parameters and an improvement on LR between Model 1 and Model 3 indicated that the exclusion of maternal genetic effect may result, to some extent, in biased estimates of genetic predictions and genetic parameters for BWT, WWT, YWT and ADG1 in Tswana cattle. The current results are consistent with the findings by Thompson (1976) and Costa et al. (2011), who reported that maternal genetic effect remained present after weaning to yearling weight and mature weight. However, the inclusion of maternal effects did not substantially affect the estimates of direct genetic variances and the corresponding heritabilities for EWT and ADG2 in the current study. According to some reported findings, the opinions as to whether to consider the effects when computing estimated breeding values (EBV) for EWT and ADG2 seem to be contradictory (Meyer, 1999;Rumph et al., 2002;Costa et al., 2011). Despite that, Thompson (1976) emphasised that if a trait is affected by maternal effects, even to a smaller degree, these effects had to be considered in the model to avoid biased and inflated direct heritability estimates. Although Meyer (1992) argued that modern genetic evaluation of post-weaning weight traits tend to exclude maternal effects in the estimation of genetic parameters, hence estimation of direct genetic effects and direct heritabilities might be inflated. The current study suggests that the inclusion of maternal effects in EWT and ADG2 analyses may have little or no effect on the estimated direct genetic variances and the corresponding direct heritabilities in Tswana cattle.
The maternal permanent environmental effect explained roughly 3 to 8% of the variability in growth traits. The estimates for maternal permanent environmental effects followed the same pattern to maternal genetic effects. Estimates were large at younger age (WWT) and declined towards EWT. When compared with genetic parameters the maternal permanent environmental effects were generally smaller. However, the observed influence of this effect beyond weaning was not significant, suggesting that maternal permanent effects could even be excluded in the genetic analysis of YWT. As a result, the current findings conformed to the argument by Costa et al. (2011) that failure to account for maternal effects may result, to a certain extent but not greatly, in biased estimates of genetic predictions and genetic parameters in post-weaning growth traits. The same authors further stated that the resultant bias of EBV may be very small or even unnoticeable but the effects should be included in the analysis. Similar to the value obtained in Model 2 for maternal genetic heritability estimate, the zero value obtained for maternal permanent environmental effect in Model 3 for ADG2 may be attributed either to these effects being none influential on this trait or a statistical artefact due to insufficient records hence incomplete censoring.
Direct and maternal genetic covariances and correlations were not estimable for all growth traits and ADGs. According to Robinson (1996a) and Meyer (1997) the estimation of this effect is not simple and shows strong dependence on the data structure hence it should be set to zero. In addition to that, Boligon et al. (2012) also suggested that since direct and maternal genetic covariance is dependent on data structure, it should be set to zero when estimating genetic parameters for early growth traits from the data set with insufficient structure, that is, small proportion of cows with phenotypic information, small number of progenies and poorly linked pedigree relationships. Estimates of direct and maternal heritabilities for BWT, WWT, YWT, EWT, ADG1 and ADG2 obtained through bivariate analyses did not show any notable deviation from those obtained through univariate analyses.
Estimates of phenotypic and genetic correlations among BWT, WWT, YWT, EWT, ADG1 and ADG2 are presented in Table 4. Genetic correlation estimates were all positive and ranged from 0.19 to 0.98 (Table 4), indicating low to high genetic association between the traits. The positive correlations between traits suggest that selection for one trait may result in a correlated response in another trait. The genetic correlations between BWT and any of WWT, YWT, EWT and ADG1 were low to medium and ranged from 0.19 to 0.36. This therefore signalled low chances of increased BWT hence minimal risk associated with calving difficulty when selecting for increased WWT, YWT, EWT or ADG1. The current values obtained for correlations are consistent with those reported by Raphaka and Dzama (2010) for the same breed and Maiwashe et al. (2002) for Bonsmara breed. High correlation values obtained between WWT and each of ADG1, YWT, and EWT indicate that an increase in any of these weight traits can be attained through selection for WWT, hence the current results are in agreement with the findings of Burrow (2001) and Raphaka and Dzama (2010). Genetic correlation between YWT and EWT was also very high  signifying that selection for one of these traits can actually yield correlated response in the other. The inestimable genetic correlation between pre-weaning traits with ADG2 suggests that the estimates may be low due to compensatory growth effects (Raphaka and Dzama, 2010;Corbet et al., 2006).
Phenotypic correlations obtained in this study were generally lower than the corresponding genetic correlation estimates. This has also been observed in some reviews (Mohiuddin, 1993;Koots et al., 1994) and it was suggested that estimates of lower phenotypic correlations compared to the corresponding genetic correlations will arise where the genes governing the two traits are similar, but the environments affecting the expression of these traits have low correlation (Searle, 1961). As was the situation with the genetic correlations, phenotypic correlations were lowest between BWT and all the other traits. This may indicate that BWT will not be enhanced by selecting for increased WWT and postweaning weight traits, hence no threat of increasing chances of difficult calving can be encountered.
Genetic trends estimated by averaging predicted breeding values for growth traits and average daily gains of animals in selected and control lines are shown in Figures 1 to 6 and Table 5. Although the trends fluctuated over the entire period of selection, the distinct feature in genetic trends is exhibited in the two selected lines. The two lines showed positive trends in progression of average breeding values from 1995 to 2013 in various traits.     Regression coefficient not significantly different from zero, ** and *** regression coefficient significantly differ from zero at (P<0.01) and (P<0.001) respectively, BWT = birth weight, WWT = weaning weight, YWT = yearling weight, EWT = 18-month weight, ADG1 = pre-weaning average daily gain, and ADG2 = post-weaning average daily gain.
However, genetic trends for EWT of S1 did not show any significant progression and were not distinctive from those of the control line. Genetic trends for various traits in the control line were low and fluctuated below and above zero without any clear direction of change for the better part of the selection period as anticipated, although a significant positive progression was observed in WWT. However, trends for BWT of S2 and the control line surprisingly exhibited genetic progression toward a negative direction (Figure 1). Genetic trend estimates for both lines were not significantly different from zero ( Table  5).
The genetic trends for EBVs obtained currently for S2 and S3 were lower while the one for S1 was within the range of values; 0.061, 0.015, 0.024 and 0.090 kg/year as reported for BWT by Plasse et al. (2002) in selected Brahman cattle, Ferraz et al. (2000) in Santa Gertrudis, Assan (2013) in Tuli cattle and Tawah et al. (1993) in Gudali cattle respectively. The respective annual genetic trends obtained for WWT in S1, S2 and S3 were 0.99±0.12, 0.43±0.05 and 0.40±0.11 kg/year. The values obtained for WWT in the current study were within the range of values reported for selected Hereford and Angus cattle (0.63 to 1.03 kg/year) by Baker et al. (1991) Assan (2010). The annual fluctuations of estimated genetic trend values for BWT and WWT may be due to the circumstance that selection of reproducing animals was based on phenotypic indices which may not always be a true reflection of the allied breeding values.
Annual changes in breeding value reported for EWT in the current study for S2 was consistent with the previously reported values ranging from 0.19 to 0.72 kg/year by Enns and Nicoll (2008) in selected New Zealand Angus nucleus herd while the respective values obtained for S1 and S3 were higher and lower than the reported ranges. The respective annual genetic changes obtained for EWT were 0.23±0.11, 2.03±0.17 and 0.28±0.26 kg/year in S1, S2 and S3 ( Figure 5). The values obtained in the current study for EWT in S1 and S3 were within the previously reported range of 0.13 to 0.49 kg/year (Plasse et al., 2002;Wasike et al., 2006;Enns and Nicoll, 2008) while the one obtained for EWT in S2 was higher.
Genetic trends currently obtained for ADG1 and ADG2 were generally high and low respectively (Figures 3 and  6). The values currently obtained for ADG2 in S1, S2 and S3 were all not significantly different from zero. However, the values obtained in the current study for ADG1 in S2 and S3 were fairly closer to the reported values of 1.92 g/year by Beffa (2005) in Afrikaner cattle and 1.57 g/year by Tawah et al. (1993) in selected Wakwa cattle, respectively. The value currently obtained for S1 was way too high. The values currently obtained for ADG2 in S1 and S3 were far lower than the range of 0.13 g/year to 0.29 g/year previously reported for a selected New Zealand Angus nucleus herd (Enns and Nicoll, 2008) while the one currently obtained for S2 was almost similar to the value of 0.66 g/year reported by Tawah et al. (1993) for selected Gudali cattle. The deviations of annual mean breeding values for BWT, WWT, YWT and ADG1 between the two selection lines demonstrated effective improvement in S1 than in S2. The high annual breeding values obtained for both EWT and ADG2 in S2 were contrary to the findings of Rakwadi (2010) who reported no clear genetic improvement in Brahman, Tuli and Bonsmara breeds selected for EWT. The unselected control line exhibited generally lower annual breeding values for all traits than the two selected lines which may be due to sudden change in genotype of the selected lines over the years of selection (Yaeghoobi et al., 2011).
The current estimates for annual genetic gain indicated that in all traits, annual genetic gains were significantly enhanced among selected lines compared to the unselected control line. The directions of genetic change were desirable for all the traits except for BWT in S2 and S3. Both significant genetic progression in WWT of the control line and genetic change towards negative direction in BWT of the S2 and S3 were not anticipated and the cause could not be well identified in the current study. However, it may be attributed to asymmetry of selection response due to maternal effects or random drift (Falconer, 1989). High genetic gain for EWT in S2 demonstrates that selection for increased eighteen-month weight enhanced breeding value for EWT than selection for increased weaning weight. This clearly indicates the potential for improvement in the rate of change in annual genetic gain if an effective selection program was to be implemented.

Conclusion
Direct genetic variances and corresponding heritability estimates obtained in the current study for WWT, YWT, EWT, ADG1 and ADG2 were high, demonstrating the breed's potential for genetic improvement if planned selection has to be carried out. Maternal genetic effects had significant influence on BWT, WWT, YWT and ADG1, while maternal permanent environmental effects had significant influence on BWT, WWT and ADG1, hence the effects should always be accounted for when estimating breeding values and genetic parameters for these growth traits and pre-weaning ADG in Tswana cattle. High genetic correlations among growth traits suggested that selection for any of the growth traits from weaning weight to 18-month weight and average daily gains can result in correlated responses in other traits.
Low genetic correlation between birth weight and other growth traits and average daily gains indicated the low chances of increased birth weight due to selection for increased growth traits and average daily gains hence minimal risks of dystocia. The two selected lines significantly enhanced annual genetic gain for all growth traits and ADG1, but not for ADG2. However, genetic gain for EWT was optimally improved in S2 than in S1. The control line exhibited significant annual genetic gain in WWT which was not anticipated and perhaps attributed to asymmetry of selection response. In general the study revealed that selection has yielded positive response over the years hence presents an opportunity for planned selection program to improve growth traits of the Tswana cattle breed.