PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of gseBioMed CentralBiomed Central Web Sitesearchsubmit a manuscriptregisterthis articleGenetics, Selection, Evolution : GSEJournal Front Page
 
Genet Sel Evol. 2016; 48: 3.
Published online 2016 January 14. doi:  10.1186/s12711-015-0178-5
PMCID: PMC5518165

Genotype by environment interaction for tick resistance of Hereford and Braford beef cattle using reaction norm models

Abstract

Background

The cattle tick is a parasite that adversely affects livestock performance in tropical areas. Although countries such as Australia and Brazil have developed genetic evaluations for tick resistance, these evaluations have not considered genotype by environment (G*E) interactions. Genetic gains could be adversely affected, since breedstock comparisons are environmentally dependent on the presence of G*E interactions, particularly if residual variability is also heterogeneous across environments. The objective of this study was to infer upon the existence of G*E interactions for tick resistance of cattle based on various models with different assumptions of genetic and residual variability.

Methods

Data were collected by the Delta G Connection Improvement program and included 10,673 records of tick counts on 4363 animals. Twelve models, including three traditional animal models (AM) and nine different hierarchical Bayesian reaction norm models (HBRNM), were investigated. One-step models that jointly estimate environmental covariates and reaction norms and two-step models based on previously estimated environmental covariates were used to infer upon G*E interactions. Model choice was based on the deviance criterion information.

Results

The best-fitting model specified heterogeneous residual variances across 10 subclasses that were bounded by every decile of the contemporary group (CG) estimates of tick count effects. One-step models generally had the highest estimated genetic variances. Heritability estimates were normally higher for HBRNM than for AM. One-step models based on heterogeneous residual variances also usually led to higher heritability estimates. Estimates of repeatability varied along the environmental gradient (ranging from 0.18 to 0.45), which implies that the relative importance of additive and permanent environmental effects for tick resistance is influenced by the environment. Estimated genetic correlations decreased as the tick infestation level increased, with negative correlations between extreme environmental levels, i.e., between more favorable (low infestation) and harsh environments (high infestation).

Conclusions

HBRNM can be used to describe the presence of G*E interactions for tick resistance in Hereford and Braford beef cattle. The preferred model for the genetic evaluation of this population for tick counts in Brazilian climates was a one-step model that considered heteroscedastic residual variance. Reaction norm models are a powerful tool to identify and quantify G*E interactions and represent a promising alternative for genetic evaluation of tick resistance, since they are expected to lead to greater selection efficiency and genetic progress.

Background

The cattle tick is a parasite that adversely affects beef cattle production in tropical areas such as Brazil. Retail beef markets are imposing restrictions on meat, ensuring that it is free of chemical residues that are perceived as having negative impacts on environment, public health and human welfare. Therefore, to remain competitive in foreign beef markets, Brazil must aim at complying with these higher standards.

To ensure market competitiveness, one strategy might be to increase the contribution of the Bos taurus breeds to Brazilian herds because they are more advantageous in terms of productive traits [1], such as carcass yield, gain weight, meat quality and sexual precocity compared to Bos indicus breeds. However, Bos taurus breeds tend to have greater susceptibility to tick infestation than Bos indicus breeds [2, 3]. Hence, selection of animals for tick resistance would be useful to reduce the need for chemical control while also increasing productivity.

Evidence for additive genetic variability of tick counts in cattle includes reported heritability estimates, which range from 0.05 to 0.42 [2, 46]. Genetic evaluations for tick counts are routinely performed in countries such as Australia and South Africa, which have a similar climate as Brazil and where the cattle tick is also present. Examples of breeds with such evaluations include breeds such as Bonsmara and Belmont Red [7], and Brahman and Hereford-Shorthorn) [8]. In Brazil, the Conexão Delta G (Delta G Connection) company has used a genetic improvement program based on selection for tick resistance in Hereford and Braford cattle since 2003 [9].

These and other research studies and genetic evaluations [7, 9] have not considered genotype by environment interactions (G*E). Failing to consider G*E interactions in genetic evaluations can adversely affect breeding programs if relative genetic merit is affected by the environment [1013]; specifically, animals that are identified as top breeders in one environment may not be ideal in other environments. This issue is further exacerbated if progeny are raised in environments that differ from that of their parents [13]. In addition, most current genetic evaluation systems assume homogeneous residual variances across environments, although evidence of residual heteroscedacity has been reported, which is defined as heterogeneity of residual variances across contemporary groups, for traits such as milk yield [14] and post-weaning gain [10, 15].

Linear reaction norm models capture a simple form of G*E interactions. They are based on the use of covariance functions [16] that allow for the prediction of the relative genetic merit of animals as a function of gradual linear changes in an environmental covariate. Sometimes this environmental covariate is not known with certainty and must be estimated from the data; Su et al. [17] demonstrated how this inference uncertainty could be formally accounted for by using Bayesian methods. If G*E interactions are important for tick resistance, reaction norm models could be used to fine-tune genetic improvement for tick resistance in Brazilian beef cattle. Because G*E interactions contribute to heterogeneous genetic variability across environments, if heteroscedastic residual variability across environments is ignored, inferences on G*E interactions based on reaction norm models could be biased.

The objective of this study was to infer upon G*E interactions based on models with different assumptions regarding the nature of genetic and residual variation and with different approaches to account for uncertainty on environmental gradient.

Methods

Tick count data

Data used in this current study were obtained from a breeding program conducted by Conexão Delta G (Delta G Connection). Data included records of tick counts (TC) on Hereford and Braford beef cattle from eight herds from the Rio Grande do Sul state, Brazil. TC were obtained on each animal from 326 to 729 days of age using the method described by Wharton and Utech [18], for which all engorged female ticks larger than 4.5 mm were counted on the entire left side of the animal when average management group infestation, i.e., animals under the same feeding and sanitary management, exceeded 20 ticks per animal. Up to three such counts were obtained for each animal, with each count separated by a minimum of 30 days, as described in other studies [5, 19, 20]. A total of 241, 1934 and 2188 animals for which, respectively, one, two and three TC were recorded. The average age during the evaluation period was 524 ± 65 days, and the overall mean TC was equal to 35.0 with a standard deviation of 42.2 (ranging from 0 to 532).

The 4363 animals with records were born between 2008 and 2011 and originated from 604 sires and 3966 dams, with 10 generations of pedigree depth. A total of 11,967 animals remained after pedigree pruning (i.e., removing any terminal ancestors that occur only once in the pedigree file). Pedigree information was incomplete due to the use of multiple-sire matings; 36 % of the animals only had their dam known. For animals with TC, this increased to 65 %. Similar pedigree structures from this same population have been used in other studies [20], and they have not affected the results of genetic evaluation. A detailed breakdown of the pedigree structure is in Table 1.

Table 1
Pedigree structure as defined by parentage certainty and pedigree completeness

Because TC were not normally distributed (Fig. 1), a log-transformation was used such that LTTC = log10 (TC + 1.001), which was used as the response variable [1, 20]. The constant 1.001 was included because some TC were equal to 0 [1, 20]. Skewness and kurtosis tests were performed and ensured the normality of the residuals from the fitted models.

Fig. 1
Distribution of tick counts

Contemporary groups (CG) were defined as groups of animals within the same herd, year of birth, season of birth (April–July; August–November and December–March), sex and from the same management group. From 11,316 observations, we selected 10,673 records pertaining to 146 CG with at least five animals and with each LTTC record being within 3.5 standard deviation (SD) from its specific CG. Connectedness among the CG was determined by each having more than 10 genetic links in the dataset, using the AMC software [21]. Estimates of CG effects on LTTC were assumed to be the environmental covariates for a linear reaction norm model because they are the most appropriate entities used to describe the environmental conditions for beef cattle production [10, 22, 23].

Statistical models

Twelve analyses based on different models and/or inferential methodologies and specifications on residual variability were conducted on the data. These analyses are described below as M1 to M12 and are summarized in Table 2.

Table 2
Statistical model implemented for analysis of tick counts, including approach, contemporary group effect, heteroscedasticity specification, deviance criterion information (DIC) value with respective model ranking

Traditional animal model (AM)

Consider the following simple linear traditional animal model (M1):

yijk=xjβ+wi+aj+cj+eijk.
1

Here, y ijk is the kth LTTC record of animal j from CG i, β is the vector of fixed effects that includes an overall intercept, linear regression coefficients for Nellore breed proportion, heterozygosity and recombination loss (predetermined by Cardoso et al. [9]), as well as linear and quadratic regression coefficients on age of the animal; xj is the known incidence row vector of covariates connecting β to y ijk; w i is the random effect of CG i (i = 1, 2, …,  146 levels); a j is the random additive genetic effect of animal j; c j is the random permanent environmental effect of animal j; and e ijk is the residual error.

The following distributional assumptions were specified:

w=wiN(0,Iσw2),

a=ajN(0,Aσa2),

c=cjN(0,Iσc2)

and

e=eijkN(0,Iσe2),

where σw2, σa2, σc2 and σe2 represent variance components due to CG, additive genetics, permanent environment and residual terms, respectively. Here, A represents the numerator of the relationship matrix between the animals in the pedigree, and I is the identity matrix.

Hierarchical bayesian reaction norm models (HBRNM)

Two somewhat different approaches were used to estimate environmental sensitivities of animals. One approach was based on a commonly used two-step model [24, 25], in which in the first step, the regular animal model (M1) from Eq. (1) is used to estimate CG effects ŵ i. The second step consists of using these ŵ i estimates as if they were “known” environmental covariates in a linear reaction norm model. More specifically, posterior means of ŵ i obtained from M1 were used as covariate values in the following reaction norm model (M2).

yijk=xjβ+ϕw^i+aj+bjw^i+cj+djw^i+eijk.
2

Here, ϕ is an overall linear regression coefficient of y ijk on ŵ i; a j is the additive genetic intercept of animal j pertaining to genetic merit for an average environment (ŵ i = 0); b j is the random additive genetic effect of the reaction norm slope of animal j on ŵ i; c j is the non-genetic (e.g., permanent environmental effect) intercept of animal j, as defined for an average environment (ŵ i = 0); and d j is the random permanent environmental effect of the reaction norm slope of animal j on ŵ i. Note that y ijk, xj β and e ijk are defined as before.

Another two-step modeling strategy (M3) that is very similar to M2 is given by Eq. (3):

yijk=xjβ+wi+aj+bjw^i+cj+djw^i+eijk.
3

In M3, contemporary group effects are refitted as random effects rather than being treated as known covariates, such that M3 may be more flexible than M2 for modeling CG effects. Nevertheless, ŵ i was again used as a “known” covariate in the random regression portion of the model.

Including ŵ i as if it is a “known” covariate in the second step of this approach is clearly a limitation that may understate statistical uncertainty and lead to biased predictions on animal genetic merit. These biases may be due to genetic trend, differences in environmental covariate values across CG, or both [10, 17]. An appealing one-step approach that avoids these limitations of the two-step approach was proposed by Su et al. [17]. This approach is purely Bayesian in that the covariate associated with the reaction norm is treated as unknown, which allows inferences for all unknowns together within a one-step linear reaction norm model (M4):

yijk=xjβ+wi+aj+bjwi+cj+djwi+eijk.
4

Model M4 can be rewritten in matrix notation as below [17]:

y=Xβ+Pw+Zaa+Zbb+Zcc+Zdde
5

where y = {y ijk} is the nx1 vector of observations; β is the vector of fixed effects of order p; w=wii=1nw is the vector of environmental effects; a = {a j}j=1q is the vector of random genetic intercepts; b = {b j}j=1q is the vector of random genetic slopes; c = {c j}j=1q is the vector of random permanent environment intercepts; d = {d j}j=1q is the vector of random permanent environment slopes; and e is the nx1 vector of residuals. Matrices X, P, Z a and Z c are known incidence matrices, where the column address of matrices Z b and Z d has exactly one element equal to the environmental covariate (w i or an estimate of w i) for that CG in the row address of the observation, with all other elements in that row equal to 0.

Prior distributional specifications

To infer environmental sensitivities using a hierarchical Bayesian model, three stages are required: the first stage defines the distribution of the observed data conditional on all other parameters [17]:

y|β,w,a,b,c,d,RN(Xβ+Pw+Zaa+Zbb+Zcc+Zdd,R).
6

For a homoscedastic residual specification such as for M1, M2, M3 and M4, R = Iσe2, where σe2 is the residual variance and I is the identity matrix. However, as previously noted, it might be important to model residual heteroscedasticity. We propose two alternative strategies for this. The first heteroscedastic residual specification (S1) is defined by R = diagIniσei2, a diagonal matrix with elements equal to σei2=σe2×ηw^i and Ini denoting an identity matrix of order n i, where n i is the number of records in the ith CG. Here, η is an unknown scaling parameter that characterizes the degree of heterogeneity of residual variance across environments, and ŵ i is the solution for the ith CG [26].

Based on S1, we tested two two-step approaches (M6 and M7) that used inferred values of ŵ i from M1 as if they were known and a one-step reaction norm model (M8), where w i is an unknown covariate that is jointly inferred with the reaction norm and η parameters. Model M6 was a heteroscedastic residual extension of M2, whereas M7 and M8 were heteroscedastic residual extensions of M3 and M4, respectively.

Another heteroscedastic residual specification (S2) was based on residual variance subclasses determined by a decile-based classification of ŵ i, following Cardoso and Tempelman [10]. That is, CG were ordered into one of 10 categories based on decile delimiters of ŵ i obtained from M1, such that R = diagInkσe2γk, where the order n k denotes the number of records delimited by deciles k − 1 and k, and was 1157, 1174, 1047, 765, 1188, 1192, 1208, 918, 1150 and 874, respectively, for k = 1, 2, …,  and 10. This specification was used to extend the two-step models M2–M10 and M3–M11 and the one-step model M4–M12 with this particular heteroscedastic residual specification.

The last two models considered (M5 and M9) were heteroscedastic residual animal models based on extending M1 with S1 and S2 heteroscedastic residual specifications, and were used as control models to determine the consequences of failing to model G*E interactions versus failing to model residual heteroscedasticity [10].

The second stage of HBRNM is represented by the prior distributions of the location parameters, as follows:

β ∼ p(β), 
7

w|σw2N(0,Iσw2),
8

abN00,σa2σabσabσb2A,
9

cdN00,σc2σcdσcdσd2I,
10

where p(β[proportional, variant] 1, σw2 is the environmental effect variance; σa2 and σb2 are the additive genetic variances due to the reaction norm intercept and slope, respectively; σc2 and σd2 are permanent environment variances due to reaction norm intercept and slope, respectively; σab is the genetic covariance between reaction norm intercept and slope; and σcd is the permanent environment covariance between reaction norm intercept and slope. Then, rab=σab/σa2×σb2 and rcd=σcd/σc2×σd2 are the corresponding genetic and permanent environment correlations between intercept and slope, respectively.

Finally, the third stage of HBRNM was based on specifying an inverted gamma (IG) distribution for the variance of the contemporary group effects, i.e., σw2|αw,βwIG\,(αw=1,βw=0.097), where the mean of this distribution is:

Ewi|αw,βw=αwβw.
11

Similarly, we specify σe2|αe,βeIG\,(αe=1,βe=0.0728).

Likewise, an inverted Wishart distribution (IW) prior distribution was specified for the permanent environment and additive genetic covariance matrices, as follows:

G0=σa2σabσabσb2IW(T0,v),
12

U0=σc2σcdσcdσd2IW(T1,v).
13

Here, v = 4 represents a presumed known number of degrees of freedom, and T0=0.01810.01020.01020.0161-1 and T1=0.01560.00870.00870.0136-1 are presumed scale matrices for additive genetic and permanent environmental effects, respectively, and EG0=T0-1v-p-1 and EU0=T1-1v-p-1 are the prior means for v > p+1, where p is the number of parameters In the models with heterogeneous residual variances, additional hierarchical specifications were required, depending on the nature of the function (S1 or S2) chosen, i.e.: η|αη, βη ∼ p(η|αη, βη) =  IG(αη, βη), for S1or γkγ ∼ p(γkγ) = IG(αγ, αγ - 1), k = 1, 2, …, 10 for S2 [10, 27, 28]. We specified αη = −1 and βη = 0, where the prior p(αγ) on αγ was a gamma with shape and scale hyperparameter values of 0.03 and 0.01, respectively [10]. This assumption leads to a prior mean of αγ equal to 3 [E(αγ) = 3] and a large prior variance (var(αγ) = 300) [27].

Due to the absence of relevant previous knowledge, flat or highly dispersed prior densities were assumed for all parameters of all models, and hyperparameters for variance components priors were specified on the basis of REML estimates obtained by M1 and M2 (not shown).

Bayesian inference

Bayesian analyses were conducted to sample all parameters from their fully conditional posterior distributions. Gibbs sampling was generally used except for the w i’s and η in M5, M6, M7 and M8 and for αγ (S2) in M9, M10, M11 and M12. MCMC sampling of these parameters required a random walk Metropolis–Hastings step because their full conditional posterior distributions were unrecognizable (see Cardoso and Tempelman [10] for further details).

Monte Carlo Markov chain (MCMC) based inferences were implemented using the INTERGEN software [29] by saving every 10th cycle from a total of 1,000,000 cycles, after 100,000 cycles of burn-in. Global convergence was checked using Geweke’s Z criterion [30] applied to the conditional distribution of the data, as proposed by Brooks and Roberts [31]. In addition, visual inspection of trace plots was conducted, and a minimum effective sample size of 100 for all unknown parameters was obtained.

Model comparison

The deviance information criterion (DIC) was used to compare model fit and model complexity [32]:

DIC=D¯θ+pD=2D¯(θ)-D(θ¯),
14

where D¯(θ)=Eθ|y[D(θ)] is the posterior expectation of Bayesian deviance; pD=D¯θ-D(θ¯) corresponds to the penalty for increasing model complexity, with θ being the vector of model parameters and D(θ¯) being the Bayesian deviance as a function of the posterior mean of the parameters. Smaller values of DIC thereby indicate better-fitting models, while taking a penalty for model complexity into consideration.

Variance components and genetic parameters

The additive genetic variance of TC for a specific environment i with effect w i was obtained as follows:

σa2|wi= varaj+ bjwi=σa2+ wi2σb2+ 2wiσab.
15

Thus, the heritability (ha2) and repeatability (r) of TC for a specific environment was determined as:

ha2|wi=σa2|wiσa2wi+σc2wi+σe2|wi,
16

and

r|wi=σa2|wi+σc2|wiσa2wi+σc2wi+σe2|wi,
17

where σc2|wi and σe2|wi are permanent environment and residual variances in environment i, respectively. For homoscedastic residual models (from M1 to M4), σe2|wi is constant, i.e., σe2|wi = σe2∀i. For heteroscedastic residual models, σei2|wi=σe2×ηw^i for M5–M8, and σei2|wi=σe2×γk:i, where k:i denotes the decile-based classification k for CG i, in models M9, M10, M11 and M12.

The genetic covariance of TC between two environmental gradients based on covariate values w i and w i was calculated as:

covaaj+ bjwi,aj+ bjwi=σa2+wi+ wiσab+ wiwiσb2,
18

so that the corresponding correlation between TC in two specific environments was calculated as described below:

raaj+ bjwi,aj+ bjwi=covaaj+ bjwi,aj+ bjwi(σa2+ wi2σb2+ 2wiσab)(σa2+ wi2σb2+ 2wiσab).
19

Estimated breeding values

An estimate of the breeding value of sire j for TC, specific to a given environment i was obtained by g^j|w^i=a^j+b^jwi^ [10]. On the one hand, estimates of b^j close to 0 indicate that ĝ j is relatively constant across various environments (ŵ i) such that sire j has an environmentally robust genetic merit. On the other hand, an environmentally sensitive genetic merit has a large estimate b^j, meaning its relative performance should substantially change on the environmental gradient [33].

The sire breeding value estimates were compared based on the ranking of the animals obtained by AM and HBRNM for low, medium and high environmental levels. These values were defined by the 10, 50 and 90th percentiles for ŵ i. Potential differences in re-rankings of sires for selection based on these models were also determined by the Spearman correlation between estimated breeding values. Spearman correlations were obtained for all animals and for the top 10 % (60) of sires with 12 or more progeny between low, medium and high environmental levels under different fitted models.

Results and discussion

Model comparison

Models M1, M5 and M9, which were the only models that did not include G*E interactions with a linear reaction norm model, along with M7, and M1, had the highest or lowest DIC values. Comparison of DIC between models M1, M5 and/or M9 implies that considering heterogeneity of residual variance across environments is important for modeling LTTC. However, these DIC improvements from homoscedastic to heteroscedastic residual models were small compared to the improvements in DIC when going from regular animal to linear reaction norm models. This suggests that modeling G*E interactions is more important than modeling heterogeneous residual variances (Table 2).

The two one-step reaction norm models (M4 and M12) had lower DIC values than the corresponding two-step reaction norm models, except for M10. Thus, treating all CG effects as uncertain when modeling G*E interactions based on reaction norms seems to be important. This observation is in agreement with the findings of Su et al. [17], who demonstrated by simulation that jointly estimating all unknown parameters is more reliable than using previously estimated environmental effects from a simple animal model as known covariates. DIC can be used to compare any type of model (not necessarily nested models) [10, 13, 22, 23]. However, when fitting two-step models, the reported DIC values come from the second step because we could not account for the uncertainty about ŵ i estimates from the first step model M1. This limitation might have yielded downwards-biased p D and DIC values for two-step models, but even so, their fit was much poorer compared to their counterpart one-step models (Table 2).

Model M12 had the lowest DIC value (Table 2). Recall that M12 allows for residual variance groupings into decile-based subclasses, which agrees with the findings of Cardoso and Tempelman [10], who reported this same model as the best-fitting in the characterization of post-weaning gain in Angus cattle.

Inferences on contemporary group effects

Model M1 estimated that CG posterior means (ŵ i) ranged from −0.849 to 0.805, which were considered fixed covariates for models M2, M6 and M10 (Fig. 2). Going from the 0–1st to the 9–10th deciles, corresponding values of ŵ i were equal to −0.424, −0.224, −0.121, −0.032, 0.032, 0.107, 0.182, 0.240 and 0.316, respectively. Following Cardoso and Tempelman [10], these values were used as the cutoff points for the decile-based heteroscedastic residual subclasses defined in M9, M10, M11 and M12.

Fig. 2
Distribution of the frequencies of environmental gradient estimates (posterior means for contemporary group effects) based on different models

Posterior means ŵ i of w i were similar for all models, regardless of whether G*E interactions were considered, as in M3, M4, M7, M8, M11 and M12, or not, as in M1 and M9 (Fig. 2); Pearson correlations among these estimates between methods always exceeded 0.99, which means they were also not influenced by homoscedastic versus heteroscedastic residual modeling. These results do not agree with Cardoso and Tempelman [10] for post-weaning gain in Angus cattle, for which estimates ŵ i from the model with the decile-based heteroscedastic classification function (S2) had substantially lower correlations with estimates from the heteroscedastic exponential function models (S1), or even the conventional animal models. Furthermore, every model resulted in negative skewness on the ŵ i, ranging from −0.521 to −0.415.

Inferences on variance components and genetic parameters

The one-step model M12 resulted in the highest (0.022 ± 0.04) estimate of the reference or intercept genetic variance (σa2) compared with all other models, except for M8 (0.025 ± 0.03; Table 3). In addition, M12 showed the highest estimate of the genetic variance for slope (σb2) compared to the two-step models, except for M3 and M8, which had the same estimate (0.046 ± 0.022). Estimates of the variance components for reference permanent environment (PE) (σc2) were similar among all models (ranging from 0.006 to 0.010). In agreement with σb2, PE slopes (σd2) were also significant (ranging from 0.015 to 0.084). These results show that the one-step approach confirmed the presence of G*E interactions. Biegelmeyer [20], in a study on tick resistance in Hereford and Braford beef cattle reported similar estimates, i.e., 0.012 and 0.022 for σa2 and σc2, respectively.

Table 3
Posterior means and 95 % posterior probability intervals reported as (2.5, 97.5th) posterior percentiles of dispersion parameters estimated for tick counts of Hereford and Braford cattle by different models

Estimates of the correlation between intercept and slope for the additive genetic and permanent environment effects were characterized by a great deal of uncertainty, as shown by the widths of their respective 95 % posterior probability interval (PPI; Table 3). This large uncertainty differs from those of previous studies [10, 25, 34], which estimated large and positive correlations. These differences may in part be caused by the fact that the correlation estimates depend upon the scale used for ŵ i or because the biological nature of tick counts is different from that of production traits.

Residual variance estimates (σe2) were similar among models, ranging from 0.062 ± 0.001 to 0.074 ± 0.010, but they were slightly higher in traditional animal models M1 (0.072 ± 0.001), M5 (0.070 ± 0.001) and M9 (0.074 ± 0.010), which confirms the importance of considering G*E interactions in genetic evaluations for Hereford and Braford beef cattle (Table 3). Cardoso and Tempelman [10] also reported that HBRNM resulted in lower estimates of σe2 than AM. However, despite the similarity of the residual variances across the various reaction norm models, Fig. 3 illustrates the need to consider residual heteroscedasticity. The first decile class was particularly deviant from the other classes. This unexpected, very large residual variance at the lowest extreme of the CG effects boundary may be due to data artifacts or a non-obvious biological condition associated with low tick infestation levels. Similar results were demonstrated by Cardoso and Tempelman [10], with residual variances being remarkably decreased at the extremes of the CG average performance. Figure 3 also explains the poor fit of models M5, M6, M7 and M8, which modeled heteroscedastic residual variance as an exponential function (i.e., S1). This function forced a gradual monotonic change in the residual variances over the CG classes, while M9, M10, M11 and M12 showed a more flexible pattern, perhaps reflecting the true residual variance behavior of the actual data.

Fig. 3
Posterior mean residual variances along the environmental gradient for each 10th percentile of tick counts

Heritability estimates (h^2) were generally higher for HBRNM and for M5 and M9, compared to M1 (h^2 = 0.19 ± 0.04; Fig. 4a). Similar heritability estimates have been reported in the literature, using models such as M1 and logarithmic transformations of the observed data [1, 5]. With M12, average heritability estimates were higher, which also indirectly indicates the better fit of one-step versus two-step models that consider residual heteroscedasticity. Other studies in beef cattle also found higher average heritability estimates for weaning weight and 450-day weight, respectively, using HBRNM compared to AM [10, 35]. Therefore, greater response to selection is expected when using reaction norm models that model heterogeneity of residual variances across CG. Considering that using data from animals with unknown sires could lead to lower heritability estimates, we found that our heritability estimates were similar to those previously reported in the literature [1, 5, 20, 35].

Fig. 4
Heritabilities and repeatabilities of tick counts for the 10 and 90th percentiles of the environmental gradient

Estimates of repeatability varied along the environmental gradient (ranging from 0.18 to 0.45) and were, in general, higher under high levels of tick infestation (Fig. 4b). These results demonstrate the particular importance of modeling permanent environmental effects in harsh environments, where more resistant animals are more likely to maintain a consistent performance.

Posterior means of the genetic correlations [see Eq. (19)] between breeding values along the environmental gradients for Hereford and Braford LTTC that were obtained by the best-fitting model M12 demonstrated a large plateau above 0.80 (Fig. 5). Furthermore, estimated genetic correlations decreased as the tick infestation level increased, with negative correlations between extreme environmental levels, i.e., between more favorable (low infestation) and harsh environments (high infestation). Similar results that demonstrate differences in genetic correlations between breeding values along environmental levels, mainly between high challenge conditions and favorable environments, have been reported in the literature [10, 13, 25]. However, Ambrosini [36] estimated small differences for Nellore yearling weight, with genetic correlations between breeding values along the environmental gradient ranging from 0.78 to 1.00.

Fig. 5
Estimates of genetic correlations between tick counts in different environmental conditions obtained by model M12

Inferences on genetic merit

A low genetic correlation between breeding values in extreme infestation environments (Fig. 5) could indicate that different animals would be selected when using the reaction norm model M12. However, Spearman rank correlations among genetic values obtained by different models were always higher than 0.85 (Table 4), which indicates that rankings of animals would be similar and, thus, substantial losses on selection precision might not be observed when using a traditional animal model.

Table 4
Spearman rank correlations among posterior means of genetic values for tick counts of Hereford and Braford cattle at different tick infestation levels obtained by different models

Conclusions

Hierarchical Bayesian reaction norm models can be used to describe the presence of genotype by environment interactions for tick resistance in Hereford and Braford beef cattle. The model that best fitted tick counts in Brazilian climates was a one-step model that considered heteroscedastic residual variance based on ten discrete classes of deciles of average CG performance (M12), and hence, this model should be considered as the preferred model for genetic evaluation of this population. However, other functions on residual variance and other classes of models can be evaluated as viable approaches. Reaction norm models are a powerful tool to identify and quantify genotype by environment interactions and present a promising alternative for genetic evaluation of tick resistance, since they are expected to lead to greater selection efficiency and genetic progress.

Authors’ contributions

RRM performed the statistical analyses, conceived and designed the study, assisted in the interpretation of the results and contributed to the preparation of the final manuscript. RJT conceived and designed the study, assisted in the interpretation of the results and contributed to the preparation of the final manuscript. IA performed computational assistance, assisted in the interpretation of the results and contributed to the preparation of the final manuscript. PSL assisted in the interpretation of the results and contributed to the preparation of the final manuscript. FFS assisted in statistical analyses and in the interpretation of the results and contributed to the preparation of the final manuscript. FFC provided the database and the software for the analysis, conceived and designed the study, assisted in the interpretation of the results, contributed to the preparation of the final manuscript and directed the overall project. All authors read and approved the final manuscript.

Acknowledgements

The authors thank Delta G Connection for providing the data used in this research and CAPES and CNPq for granting the scholarships. This project was supported by CNPq—National Council for Scientific and Technological Development grant 478992/2012-2, Embrapa—Brazilian Agricultural Research Corporation grants 02.09.07.004 and 01.11.07.002.07, and Agriculture and Food Research Initiative Competitive Grant no. 2011-67015-30338 from the USDA National Institute of Food and Agriculture.

Competing interests

The authors declare that they have no competing interests.

Contributor Information

Rodrigo R. Mota, rb.vfu@atom.ogirdor.

Robert J. Tempelman, ude.usm@amlepmet.

Paulo S. Lopes, rb.vfu@sepolp.

Ignacio Aguilar, yu.gro.aini@raliugai.

Fabyano F. Silva, rb.vfu@acesnofonaybaf.

Fernando F. Cardoso, rb.aparbme@osodrac.odnanref.

References

1. Ayres DR, Pereira RJ, Boligon AA, Silva FF, Schenkel FS, Roso VM, et al. Linear and Poisson models for genetic evaluation of tick resistance in cross-bred Hereford × Nellore cattle. J Anim Breed Genet. 2013;130:417–424. doi: 10.1111/jbg.12036. [PubMed] [Cross Ref]
2. Prayaga KC, Henshall JM. Adaptability in tropical beef cattle: genetic parameters of growth, adaptive and temperament traits in a crossbred population. Anim Prod Sci. 2005;45:971–983. doi: 10.1071/EA05045. [Cross Ref]
3. Silva AM, Alencar M, Regitano LCA, Oliveira MCS. Infestação natural de fêmeas bovinas de corte por ectoparasitas na região Sudeste do Brasil. R Bras Zootec. 2010;39:1477–1482. doi: 10.1590/S1516-35982010000700012. [Cross Ref]
4. Burrow HM. Variances and covariances between productive and adaptive traits and temperament in a composite breed of tropical beef cattle. Livest Prod Sci. 2001;70:213–233. doi: 10.1016/S0301-6226(01)00178-6. [Cross Ref]
5. Budeli MA, Nephawe KA, Norris D, Selapa NW, Bergh L, Maiwashe A. Genetic parameter estimates for tick resistance in Bonsmara cattle. S Afr J Anim Sci. 2009;39:321–327.
6. Ibelli AMG, Ribeiro ARB, Giglioti R, Regitano LCA, Alencar MM, Chagas ACS, et al. Resistance of cattle of various genetic groups to the tick Rhipicephalus microplus and the relationship with coat traits. Vet Parasitol. 2012;186:425–430. doi: 10.1016/j.vetpar.2011.11.019. [PubMed] [Cross Ref]
7. Corbert NJ, Sheperd RK, Burrow HM, van der Westhuizen J, Strydom DJ, Bosman DJ. Evaluation of Bonsmara and Belmont Red cattle breeds in South Africa. 1. Productive performance. Aust J Exp Agr. 2006;46:199–212. doi: 10.1071/EA05223. [Cross Ref]
8. Shyma KP, Gupta JP, Singh V. Breeding strategies for tick resistance in tropical cattle: a sustainable approach for tick control. J Parasit Dis. 2015;39:1–6. doi: 10.1007/s12639-013-0294-5. [PMC free article] [PubMed] [Cross Ref]
9. Cardoso FF, Gomes CCG, Sollero BP, Oliveira MM, Roso VM, Piccoli ML, et al. Genomic prediction for tick resistance in Braford and Hereford cattle. J Anim Sci. 2015;93:1–13. doi: 10.2527/jas.2014-8832. [PubMed] [Cross Ref]
10. Cardoso FF, Tempelman RJ. Linear reaction norm models for genetic merit prediction of Angus cattle under genotype by environment interaction. J Anim Sci. 2012;90:2130–2141. doi: 10.2527/jas.2011-4333. [PubMed] [Cross Ref]
11. Falconer DS, Mackay TF. Introduction to quantitative genetics. 4. Harlow: Pearson Education Limited; 1996.
12. Lynch M, Walsh B. Genetics and analysis of quantitative traits. Sunderland: Sinauer Associates; 1998.
13. Corrêa MBB, Dionello NJL, Cardoso FF. Caracterização da interação genótipo-ambiente e comparação entre modelos para ajuste do ganho pós-desmama de bovinos Devon via normas de reação. Rev Bras Zootec. 2009;38:1468–1477. doi: 10.1590/S1516-35982009000800010. [Cross Ref]
14. Torres RA, Bergmann JAG, Costa CN, Pereira CS, Valente J, Penna VM, et al. Heterogeneidade de variância e avaliação genética de bovinos da raça Holandesa no Brasil. Rev Bras Zootec. 2000;29:1050–1059. doi: 10.1590/S1516-35982000000400015. [Cross Ref]
15. Carvalheiro R, Fries LA, Schenkel FS, Albuquerque LG. Efeitos da heterogeneidade de variância residual entre grupos de contemporâneos na avaliação genética de bovinos de corte. Rev Bras Zootec. 2002;31:1680–1688. doi: 10.1590/S1516-35982002000700010. [Cross Ref]
16. Kirkpatrick M, Lofsvold D, Bulmer M. Analysis of the inheritance, selection and evolution of growth trajectories. Genetics. 1990;124:979–993. [PubMed]
17. Su G, Madsen P, Lund MS, Sorensen D, Korsgaard IR, Jensen J. Bayesian analysis of the linear reaction norm model with unknown covariates. J Anim Sci. 2006;84:1651–1657. doi: 10.2527/jas.2005-517. [PubMed] [Cross Ref]
18. Wharton RH, Utech KBW. The relation between engorgement and dropping of Boophilus microplus (Canestrini) (Ixodidae) to the assessment of tick numbers on cattle. Aust J Entomol. 1970;9:171–182. doi: 10.1111/j.1440-6055.1970.tb00788.x. [Cross Ref]
19. Frisch JE, O’Neill CJ. Comparative evaluation of beef cattle breeds of African, European and Indian origins. 2. Resistance to cattle ticks and gastrointestinal nematodes. Anim Sci. 1998;67:39–48. doi: 10.1017/S1357729800009772. [Cross Ref]
20. Biegelmeyer P. Resistência genética à infestação natural e artificial por Rhipicephalus (Boophilus) microplus em bovinos das raças Hereford e Braford. Pelotas: Universidade Federal de Pelotas; 2012. p. 96p.
21. Roso VM, Schenkel FS. AMC-a computer programme to assess the degree of connectedness among contemporary groups. In Proceedings of the 8th World Congress on Genetics Applied to Livestock Production: 13–18 August 2006; Belo Horizonte. 2006.
22. Cardoso LL, Braccini Neto J, Cardoso FF, Cobuci JA, Biassus IO, Barcellos JOJ. Hierarchical Bayesian models for genotypex environment estimates in post-weaning gain of Hereford bovine via reaction norms. Rev Bras Zootec. 2011;40:294–300. doi: 10.1590/S1516-35982011000200009. [Cross Ref]
23. Mattar M, Silva LO, Alencar MM, Cardoso FF. Genotype x environment interaction for long-yearling weight in Canchim cattle quantified by reaction norm analysis. J Anim Sci. 2011;89:2349–2355. doi: 10.2527/jas.2010-3770. [PubMed] [Cross Ref]
24. Calus MPL, Groen AF, de Jong G. Genotypex environment interaction for protein yield in Dutch dairy cattle as quantified by different models. J Dairy Sci. 2002;85:3115–3123. doi: 10.3168/jds.S0022-0302(02)74399-3. [PubMed] [Cross Ref]
25. Kolmodin R, Strandberg E, Madsen P, Jensen J, Jorjani H. Genotype by environment interaction in Nordic dairy Cattle studied using reaction norms. Acta Agr Scand A-AN. 2002;52:11–24.
26. Knap PW, Su G. Genotype by environment interaction for litter size in pigs as quantified by reaction norms analysis. Animal. 2008;2:1742–1747. doi: 10.1017/S1751731108003145. [PubMed] [Cross Ref]
27. Cardoso FF, Rosa GJ, Tempelman RJ. Multiple-breed genetic inference using heavy-tailed structural models for heterogeneous residual variances. J Anim Sci. 2005;83:1766–1779. doi: 10.2527/2005.8381766x. [PubMed] [Cross Ref]
28. Kizilkaya K, Tempelman RJ. A general approach to mixed effects modeling of residual variances in generalized linear mixed models. Genet Sel Evol. 2005;37:31–56. doi: 10.1186/1297-9686-37-1-31. [PMC free article] [PubMed] [Cross Ref]
29. Cardoso FF. Application of Bayesian inference in animal breeding using the Intergen program: manual of version 1.2. Brasilia: Embrapa Pecuária Sul; 2010.
30. Geweke J. Evaluating the accuracy of sampling-based approaches to the calculation of posterior moments. Federal Reserve Bank of Minneapolis, Research Department Staff Report 148; 1991.
31. Brooks SP, Roberts GO. Convergence assessment techniques for Markov chain Monte Carlo. Stat Comput. 1998;8:319–335. doi: 10.1023/A:1008820505350. [Cross Ref]
32. Spiegelhalter DJ, Best NG, Carlin BP, Van Der Linde A. Bayesian measures of model complexity and fit. J R Stat Soc Series B Stat Methodol. 2002;64:583–639. doi: 10.1111/1467-9868.00353. [Cross Ref]
33. Falconer DS. Selection in different environments: effects on environmental sensitivity (reaction norm) and on mean performance. Genet Res. 1990;56:57–70. doi: 10.1017/S0016672300028883. [Cross Ref]
34. Shariati MM, Su G, Madsen P, Sorensen D. Analysis of milk production traits in early lactation using a reaction norm model with unknown covariates. J Dairy Sci. 2007;90:5759–5766. doi: 10.3168/jds.2007-0048. [PubMed] [Cross Ref]
35. Pégolo NT, Oliveira HN, Albuquerque LG, Bezerra LAF, Lôbo RB. Genotype by environment interaction for 450-day weight of Nelore cattle analyzed by reaction norm models. Genet Mol Biol. 2009;32:281–287. doi: 10.1590/S1415-47572009005000027. [PMC free article] [PubMed] [Cross Ref]
36. Ambrosini DP, Carneiro PLS, Neto JB, Malhado CHM, Martins Filho R, Cardoso FF. Interação genótipo × ambiente para peso ao ano em bovinos Nelore Mocho no Nordeste do Brasil. Brasília: Pesqui Agropecu Bras; 2012;47:1489-95.

Articles from Genetics, Selection, Evolution : GSE are provided here courtesy of BioMed Central