- Research Article
- Open Access
- Published:

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

*Genetics Selection Evolution*
**volume 48**, Article number: 3 (2016)

## 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, 4–6]. 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 [10–13]; 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.

Because TC were not normally distributed (Fig. 1), a log-transformation was used such that LTTC = log_{10} (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.

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 M_{1} to M_{12} and are summarized in Table 2.

### Traditional animal model (AM)

Consider the following simple linear traditional animal model (M_{1}):

Here, *y*
_{
ijk
} is the *k*th 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; **x**′_{
j
} is the known incidence row vector of covariates connecting **β** to *y*
_{
ijk
}; *w*
_{
i
} is the random effect of CG *i* (\(i = 1, 2, \ldots ,\) 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:

and

where σ
^{2}_{w}
, σ
^{2}_{a}
, σ
^{2}_{c}
and σ
^{2}_{e}
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 (M_{1}) 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 M_{1} were used as covariate values in the following reaction norm model (M_{2}).

Here, \(\phi\) 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
}, **x**′_{
j
}
**β** and *e*
_{
ijk
} are defined as before.

Another two-step modeling strategy (M_{3}) that is very similar to M_{2} is given by Eq. (3):

In M_{3}, contemporary group effects are refitted as random effects rather than being treated as known covariates, such that M_{3} may be more flexible than M_{2} 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 (M_{4}):

Model M_{4} can be rewritten in matrix notation as below [17]:

where **y** = {*y*
_{
ijk
}} is the *nx1* vector of observations; **β** is the vector of fixed effects of order *p*; \({\mathbf{w = }} \, \left\{ {w_{i} } \right\}_{{i{ = }1}}^{{n_{w} }}\) is the vector of environmental effects; **a** = {*a*
_{
j
}}
^{q}_{
j=1}
is the vector of random genetic intercepts; **b** = {*b*
_{
j
}}
^{q}_{
j=1}
is the vector of random genetic slopes; **c** = {*c*
_{
j
}}
^{q}_{
j=1}
is the vector of random permanent environment intercepts; **d** = {*d*
_{
j
}}
^{q}_{
j=1}
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]:

For a homoscedastic residual specification such as for M_{1}, M_{2}, M_{3} and M_{4}, **R** = **I**σ
^{2}_{e}
, where σ
^{2}_{e}
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 (S_{1}) is defined by **R** = \(diag\left( {{\textbf{I}}_{{n_{i} }}\upsigma_{{{\text{e}}_{\text{i}} }}^{ 2} } \right)\), a diagonal matrix with elements equal to \(\upsigma_{{{\text{e}}_{\text{i}} }}^{ 2} \,= \,\upsigma_{\text{e}}^{ 2} \times\upeta^{{\hat{w}_{i} }}\) and \({\textbf{I}}_{{n_{i} }}\) denoting an identity matrix of order *n*
_{
i
}, where *n*
_{
i
} is the number of records in the *i*th 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 *i*th CG [26].

Based on S_{1}, we tested two two-step approaches (M_{6} and M_{7}) that used inferred values of *ŵ*
_{
i
} from M_{1} as if they were known and a one-step reaction norm model (M_{8}), where *w*
_{i} is an unknown covariate that is jointly inferred with the reaction norm and *η* parameters. Model M_{6} was a heteroscedastic residual extension of M_{2}, whereas M_{7} and M_{8} were heteroscedastic residual extensions of M_{3} and M_{4}, respectively.

Another heteroscedastic residual specification (S_{2}) 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 M_{1}, such that **R** = \(diag\left( {{\textbf{I}}_{{n_{k} }}\upsigma_{\text{e}}^{ 2} \upgamma_{k} } \right)\), 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, \ldots ,\) and 10. This specification was used to extend the two-step models M_{2}–M_{10} and M_{3}–M_{11} and the one-step model M_{4}–M_{12} with this particular heteroscedastic residual specification.

The last two models considered (M_{5} and M_{9}) were heteroscedastic residual animal models based on extending M_{1} with S_{1} and S_{2} 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:

where p(**β**) ∝ 1, σ
^{2}_{w}
is the environmental effect variance; σ
^{2}_{a}
and σ
^{2}_{b}
are the additive genetic variances due to the reaction norm intercept and slope, respectively; σ
^{2}_{c}
and σ
^{2}_{d}
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, \(\text{r}_{{\text{ab }}} = \upsigma_{\text{ab}} / \sqrt {\upsigma_{\text{a}}^{ 2} \times }\upsigma_{\text{b}}^{ 2}\) and \({\text{r}}_{\text{cd }} = \upsigma_{\text{cd}} / \sqrt {\upsigma_{\text{c}}^{ 2} \times }\upsigma_{\text{d}}^{ 2}\) 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., \(\upsigma_{\text{w}}^{ 2} | \upalpha_{\text{w}} ,\upbeta_{\text{w}} \sim {\text{IG\,(}}\upalpha_{\text{w}} = 1, \,\upbeta_{\text{w}} = 0 . 0 9 7 )\), where the mean of this distribution is:

Similarly, we specify \(\upsigma_{\text{e}}^{ 2} | {{\alpha }}_{\text{e}} ,\upbeta_{\text{e}} \sim {\text{IG\,(}}\upalpha_{\text{e}} = 1 , \upbeta_{\text{e}} = 0 . 0 7 2 8 )\).

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

Here, *v* = 4 represents a presumed known number of degrees of freedom, and \({\mathbf{T}}_{{\mathbf{0}}} = \left[ {\begin{array}{*{20}c} { 0. 0 1 8 1} & { 0. 0 1 0 2} \\ { 0. 0 1 0 2} & { 0. 0 1 6 1} \\ \end{array} } \right]^{ - 1}\) and \({\mathbf{T}}_{{\mathbf{1}}} = \left[ {\begin{array}{*{20}c} { 0.0156} & { 0. 0 0 8 7} \\ { 0. 0 0 8 7} & { 0 . 0 1 3 6} \\ \end{array} } \right]^{{ - 1}}\) are presumed scale matrices for additive genetic and permanent environmental effects, respectively, and \(E\left( {{\mathbf{G}}_{{\mathbf{0}}} } \right) = \frac{{{\mathbf{T}}_{{\mathbf{0}}}^{{{\mathbf{ - 1}}}} }}{v - p - 1}\) and \(E\left( {{\mathbf{U}}_{{\mathbf{0}}} } \right) = \frac{{{\mathbf{T}}_{{\mathbf{1}}}^{{ - {\mathbf{1}}}} }}{v - 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 (S_{1} or S_{2}) chosen, i.e.: \(\upeta | \upalpha_{\upeta} ,\upbeta_{\upeta} \sim {\text{p(}}\upeta |\upalpha_{\upeta} ,\upbeta_{\upeta} ) = {\text{ IG(}}\upalpha_{\upeta} ,\upbeta_{\upeta} )\), for S_{1}or \(\upgamma_{\text{k}} | \upalpha_{\upgamma} \sim {\text{p}}\left( {\upgamma_{\text{k}} |\upalpha_{\upgamma} } \right) = {\text{IG}}\left( {\upalpha_{\upgamma} ,\upalpha_{\upgamma} - 1} \right)\), \(k = 1, 2, \ldots , 10\) for S_{2} [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 (\({\text{var(}}\upalpha_{\gamma } ) = {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 M_{1} and M_{2} (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 M_{5}, M_{6}, M_{7} and M_{8} and for α_{γ} (S_{2}) in M_{9}, M_{10}, M_{11} and M_{12}. 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]:

where \(\rm {\overline{D}} (\uptheta ){\text{ }} = {\text{ }}E_{{\uptheta |y}} [D(\uptheta )]\) is the posterior expectation of Bayesian deviance; \({\text{p}}_{\text{D}} \,= \overline{\text{D}} \left(\uptheta \right) - {\text{D(}}\overline{\uptheta} )\) corresponds to the penalty for increasing model complexity, with θ being the vector of model parameters and \({\text{D(}}{\bar{\theta }} )\) 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:

Thus, the heritability (h
^{2}_{a}
) and repeatability (r) of TC for a specific environment was determined as:

and

where σ
^{2}_{c}
|w_{i} and σ
^{2}_{e}
|w_{i} are permanent environment and residual variances in environment *i*, respectively. For homoscedastic residual models (from M_{1} to M_{4}), σ
^{2}_{e}
|w_{i} is constant, i.e., σ
^{2}_{e}
|w_{i} = σ
^{2}_{e}
∀i. For heteroscedastic residual models, \(\upsigma_{{{\text{e}}_{\text{i}} }}^{ 2} | {\text{w}}_{\text{i}} \,= \,\upsigma_{\text{e}}^{ 2} \times\upeta^{{\hat{w}_{i} }}\) for M_{5}–M_{8}, and \(\upsigma_{{{\text{e}}_{\text{i}} }}^{ 2} | {\text{w}}_{\text{i}} \,= \,\upsigma_{\text{e}}^{ 2} \times\upgamma_{k:i}\), where *k*:*i* denotes the decile-based classification *k* for CG *i*, in models M_{9}, M_{10}, M_{11} and M_{12}.

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

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

### Estimated breeding values

An estimate of the breeding value of sire *j* for TC, specific to a given environment *i* was obtained by \({\hat{\text{g}}}_{\text{j}} | {\hat{\text{w}}}_{\text{i}} \,= \,\widehat{\text{a}}_{\text{j}} \,+ \,\widehat{\text{b}}_{\text{j}} \widehat{{\text{w}_{i} }}\) [10]. On the one hand, estimates of \(\hat{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 \(\hat{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 M_{1}, M_{5} and M_{9}, which were the only models that did not include G*E interactions with a linear reaction norm model, along with M_{7}, and M_{1}, had the highest or lowest DIC values. Comparison of DIC between models M_{1}, M_{5} and/or M_{9} 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 (M_{4} and M_{12}) had lower DIC values than the corresponding two-step reaction norm models, except for M_{10}. 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 M_{1}. 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 M_{12} had the lowest DIC value (Table 2). Recall that M_{12} 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 M_{1} estimated that CG posterior means (*ŵ*
_{
i
}) ranged from −0.849 to 0.805, which were considered fixed covariates for models M_{2}, M_{6} and M_{10} (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 M_{9}, M_{10}, M_{11} and M_{12}.

Posterior means *ŵ*
_{
i
} of *w*
_{
i
} were similar for all models, regardless of whether G*E interactions were considered, as in M_{3}, M_{4}, M_{7}, M_{8}, M_{11} and M_{12}, or not, as in M_{1} and M_{9} (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 (S_{2}) had substantially lower correlations with estimates from the heteroscedastic exponential function models (S_{1}), 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 M_{12} resulted in the highest (0.022 ± 0.04) estimate of the reference or intercept genetic variance (σ
^{2}_{a}
) compared with all other models, except for M_{8} (0.025 ± 0.03; Table 3). In addition, M_{12} showed the highest estimate of the genetic variance for slope (σ
^{2}_{b}
) compared to the two-step models, except for M_{3} and M_{8}, which had the same estimate (0.046 ± 0.022). Estimates of the variance components for reference permanent environment (PE) (σ
^{2}_{c}
) were similar among all models (ranging from 0.006 to 0.010). In agreement with σ
^{2}_{b}
, PE slopes (σ
^{2}_{d}
) 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 σ
^{2}_{a}
and σ
^{2}_{c}
, respectively.

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 (σ
^{2}_{e}
) were similar among models, ranging from 0.062 ± 0.001 to 0.074 ± 0.010, but they were slightly higher in traditional animal models M_{1} (0.072 ± 0.001), M_{5} (0.070 ± 0.001) and M_{9} (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 σ
^{2}_{e}
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 M_{5}, M_{6}, M_{7} and M_{8}, which modeled heteroscedastic residual variance as an exponential function (i.e., S_{1}). This function forced a gradual monotonic change in the residual variances over the CG classes, while M_{9}, M_{10}, M_{11} and M_{12} showed a more flexible pattern, perhaps reflecting the true residual variance behavior of the actual data.

Heritability estimates (\(\hat{h}^{2}\)) were generally higher for HBRNM and for M_{5} and M_{9}, compared to M_{1} (\(\hat{h}^{2}\) = 0.19 ± 0.04; Fig. 4a). Similar heritability estimates have been reported in the literature, using models such as M_{1} and logarithmic transformations of the observed data [1, 5]. With M_{12}, 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].

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 M_{12} 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.

### 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 M_{12}. 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.

## 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 (M_{12}), 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.

## References

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–24.

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–83.

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–82.

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–33.

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–7.

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–30.

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.

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.

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.

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–41.

Falconer DS, Mackay TF. Introduction to quantitative genetics. 4th ed. Harlow: Pearson Education Limited; 1996.

Lynch M, Walsh B. Genetics and analysis of quantitative traits. Sunderland: Sinauer Associates; 1998.

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–77.

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–9.

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–8.

Kirkpatrick M, Lofsvold D, Bulmer M. Analysis of the inheritance, selection and evolution of growth trajectories. Genetics. 1990;124:979–93.

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–7.

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–82.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.

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.

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.

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.

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–55.

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–23.

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.

Knap PW, Su G. Genotype by environment interaction for litter size in pigs as quantified by reaction norms analysis. Animal. 2008;2:1742–7.

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–79.

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.

Cardoso FF. Application of Bayesian inference in animal breeding using the Intergen program: manual of version 1.2. Brasilia: Embrapa Pecuária Sul; 2010.

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.

Brooks SP, Roberts GO. Convergence assessment techniques for Markov chain Monte Carlo. Stat Comput. 1998;8:319–35.

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.

Falconer DS. Selection in different environments: effects on environmental sensitivity (reaction norm) and on mean performance. Genet Res. 1990;56:57–70.

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–66.

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–7.

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.

## 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.

## Author information

### Authors and Affiliations

### Corresponding author

## Rights and permissions

**Open Access** This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

## About this article

### Cite this article

Mota, R.R., Tempelman, R.J., Lopes, P.S. *et al.* Genotype by environment interaction for tick resistance of Hereford and Braford beef cattle using reaction norm models.
*Genet Sel Evol* **48, **3 (2016). https://doi.org/10.1186/s12711-015-0178-5

Received:

Accepted:

Published:

DOI: https://doi.org/10.1186/s12711-015-0178-5

### Keywords

- Residual Variance
- Deviance Information Criterion
- Genetic Evaluation
- Genetic Merit
- Contemporary Group