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

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. [ Table Omitted - see PDF ]

Table 1

Pedigree structure as defined by parentage certainty and pedigree completeness

With tick counts

Without tick counts

Total

Both parents known

1515

1917

3432

Both parents unknown

4

3807

3811

Only sire known

7

4

11

Only dam known

2837

1876

4713

Total

4363

7604

11,967

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. [ Table Omitted - see PDF ]

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

Mx

Two-step

CG

HET

DIC value

Ranking

M1

N/A

R

S0

4828.60

12

M2

Y

C

S0

3736.80

4

M3

Y

R

S0

4010.93

8

M4

N

R

S0

3590.84

3

M5

N/A

R

S1

4507.90

11

M6

Y

C

S1

3863.92

7

M7

Y

R

S1

4258.18

10

M8

N

R

S1

3819.11

5

M9

N/A

R

S2

4129.39

9

M10

Y

C

S2

3549.11

2

M11

Y

R

S2

3823.49

6

M12

N

R

S2

3114.77

1

M x model number x; two-step (Y yes, N no, N/A non-applicable because the model is not a reaction norm model); CG specification on contemporary group effects (C covariate, R random classification effect); HET heterogeneous residual variance; S 0 homoscedastic residual variance; S 1 exponential function on heteroscedastic residual variance; S 2 discrete subclasses based on classification function on heteroscedastic residual variance

Traditional animal model (AM)

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

[Math Processing Error]

(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; x? j is the known incidence row vector of covariates connecting ? to y ijk ; w i is the random effect of CG i ([ Equation Omitted - see PDF ] 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:

[Math Processing Error]

[Math Processing Error]

[Math Processing Error]

and

[Math Processing Error]

where ? w 2 , ? a 2 , ? c 2 and ? e 2 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).

[Math Processing Error]

(2)

Here, [ Equation Omitted - see PDF ] ? 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 (M3) that is very similar to M2 is given by Eq. (3):

[Math Processing Error]

(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):

[Math Processing Error]

(4)

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

[Math Processing Error]

(5)

where y = {y ijk } is the nx1 vector of observations; ? is the vector of fixed effects of order p; [ Equation Omitted - see PDF ][Math Processing Error] is the vector of environmental effects; a = {a j } j=1 q is the vector of random genetic intercepts; b = {b j } j=1 q is the vector of random genetic slopes; c = {c j } j=1 q is the vector of random permanent environment intercepts; d = {d j } j=1 q 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]:

[Math Processing Error]

(6)

For a homoscedastic residual specification such as for M1, M2, M3 and M4, R = I? e 2 , where ? e 2 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 = [ Equation Omitted - see PDF ][Math Processing Error] , a diagonal matrix with elements equal to [ Equation Omitted - see PDF ][Math Processing Error] and [ Equation Omitted - see PDF ] 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 = [ Equation Omitted - see PDF ][Math Processing Error] , 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 [ Equation Omitted - see PDF ] 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:

[Math Processing Error]

(7)

TeX

w|?2w?N(0,I?2w),

(8)

[Math Processing Error]

(9)

[Math Processing Error]

(10)

where p(?) ? 1, ? w 2 is the environmental effect variance; ? a 2 and ? b 2 are the additive genetic variances due to the reaction norm intercept and slope, respectively; ? c 2 and ? d 2 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, [ Equation Omitted - see PDF ] rab =?ab/?2a×?????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????2b and [ Equation Omitted - see PDF ] rcd =?cd/?2c×?????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????2d 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., [ Equation Omitted - see PDF ] , where the mean of this distribution is:

TeX

(11)

Similarly, we specify [ Equation Omitted - see PDF ] .

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

TeX

(12)

TeX

(13)

Here, v = 4 represents a presumed known number of degrees of freedom, and [ Equation Omitted - see PDF ] and [ Equation Omitted - see PDF ] are presumed scale matrices for additive genetic and permanent environmental effects, respectively, and [ Equation Omitted - see PDF ] and [ Equation Omitted - see PDF ] 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.: [ Equation Omitted - see PDF ] , for S1or [ Equation Omitted - see PDF ] , [ Equation Omitted - see PDF ] 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 ([ Equation Omitted - see PDF ] ) [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]:

TeX

(14)

where [ Equation Omitted - see PDF ] is the posterior expectation of Bayesian deviance; [ Equation Omitted - see PDF ] corresponds to the penalty for increasing model complexity, with ? being the vector of model parameters and [ Equation Omitted - see PDF ] 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:

TeX

(15)

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

TeX

(16)

and

TeX

(17)

where ? c 2 |wi and ? e 2 |wi are permanent environment and residual variances in environment i, respectively. For homoscedastic residual models (from M1 to M4), ? e 2 |wi is constant, i.e., ? e 2 |wi = ? e 2 ?i. For heteroscedastic residual models, [ Equation Omitted - see PDF ] for M5-M8, and [ Equation Omitted - see PDF ] , 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:

TeX

(18)

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

TeX

(19)

Estimated breeding values

An estimate of the breeding value of sire j for TC, specific to a given environment i was obtained by [ Equation Omitted - see PDF ] [10]. On the one hand, estimates of [ Equation Omitted - see PDF ] 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 [ Equation Omitted - see PDF ] , 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. [ Table Omitted - see PDF ]

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

Parameter

Models

M1

M2

M3

M4

M5

M6

? a 2

0.019

0.019

0.020

0.020

0.018

0.019

(0.011, 0.026)

(0.012, 0.025)

(0.013, 0.026)

(0.013, 0.026)

(0.011, 0.025)

(0.013, 0.025)

? b 2

N/A

0.032

0.046

0.036

N/A

0.036

(0.006, 0.074)

(0.011, 0.094)

(0.008, 0.087)

(0.009, 0.069)

? c 2

0.010

0.008

0.008

0.007

0.010

0.008

(0.003, 0.018)

(0.003, 0.014)

(0.003, 0.015)

(0.002, 0.014)

(0.003, 0.017)

(0.003, 0.015)

? d 2

N/A

0.063

0.053

0.084

N/A

0.04

(0.022, 0.098)

(0.009, 0.093)

(0.031, 0.123)

(0.009, 0.075)

rab

N/A

?0.23

?0.19

?0.14

N/A

?0.28

(?0.69, 0.42)

(?0.65, 0.43)

(?0.65, 0.55)

(?0.73, 0.40)

rcd

N/A

?0.09

?0.07

?0.11

N/A

0.33

(?0.65, 0.62)

(?0.71, 0.74)

(?0.69, 0.60)

? w 2

0.099

N/A

0.097

0.101

N/A

(-0.47, 0.88)

(0.079, 0.126)

(0.076, 0.123)

(0.080, 0.129)

N/A

? e 2

0.072

0.063

0.064

0.062

0.070

0.064

(0.069, 0.074)

(0.060, 0.065)

(0.061, 0.066)

(0.060, 0.065)

(0.067, 0.072)

(0.061, 0.067)

M7

M8

M9

M10

M11

M12

? a 2

0.018

0.025

0.021

0.020

0.021

0.022

(0.012, 0.025)

(0.019, 0.030)

(0.013, 0.028)

(0.014, 0.027)

(0.014, 0.028)

(0.014, 0.028)

? b 2

0.038

0.046

N/A

0.031

0.035

0.046

(0.013, 0.065)

(0.021, 0.072)

(0.009, 0.057)

(0.010, 0.063)

(0.009, 0.098)

? c 2

0.009

0.006

0.010

0.010

0.009

0.009

(0.004, 0.015)

(0.002, 0.010)

(0.004, 0.017)

(0.004, 0.016)

(0.003, 0.016)

(0.003, 0.016)

? d 2

0.023

0.015

N/A

0.021

0.020

0.055

(0.005, 0.051)

(0.003, 0.039)

(0.004, 0.049)

(0.004, 0.050)

(0.006, 0.106)

rab

?0.16

?0.45

N/A

?0.28

?0.17

?0.07

(?0.08, 0.94)

(?0.68, ?0.18)

(?0.67, 0.29)

(?0.61, 0.45)

(?0.57, 0.58)

rcd

0.63

0.35

N/A

0.53

0.39

0.30

(?0.08, 0.94)

(?0.51, 0.87)

(?0.36, 0.93)

(?0.58, 0.91)

(?0.50, 0.89)

? w 2

0.097

0.113

0.098

N/A

N/A

0.097

(0.076, 0.123)

(0.088, 0.144)

(0.077, 0.125)

(0.076, 0.124)

? e 2

0.066

0.071

0.074

0.068

0.069

0.065

(0.063, 0.068)

(0.068, 0.074)

(0.059, 0.099)

(0.055, 0.089)

(0.056, 0.089)

(0.053, 0.085)

? a 2 reaction norm intercept genetic variance; ? b 2 reaction norm slope genetic variance; ? c 2 reaction norm intercept permanent environment variance; ? d 2 reaction norm slope permanent environment variance; r ab genetic correlation between intercept and slope of the reaction norm; r cd permanent environment correlation between intercept and slope of the reaction norm; ? w 2 contemporary group effect (environmental) variance; ? e 2 residual variance

M 1 linear animal model; M 2 two-step linear reaction norm model; M 3 two-step linear reaction norm model with the random contemporary group (CG) effect being re-estimated; M 4 one-step linear reaction norm model with homoscedastic residual variance; M 5 linear animal model with exponential function on heteroscedastic residual variance; M 6 two-step linear reaction norm model with exponential function on heteroscedastic residual variance; M 7 two-step linear reaction norm model with exponential function on heteroscedastic residual variance and with the random CG effect being re-estimated; M 8 one-step linear reaction norm model with exponential function on heteroscedastic residual variance; M 9 linear animal model with classification function grouped in discrete subclasses on heteroscedastic residual variance; M 10 two-step linear reaction norm model with classification function grouped in discrete subclasses on heteroscedastic residual variance; M 11 two-step linear reaction norm model with classification function grouped in discrete subclasses on heteroscedastic residual variance and the random CG effect being re-estimated; M 12 one-step linear reaction norm model with classification function grouped in discrete subclasses on heteroscedastic residual variance; N/A not applicable

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 (? e 2 ) 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 ? e 2 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. [ Table Omitted - see PDF ]

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

Model

M1

M4

M4

M4

M10

M10

M10

M12

M12

M12

M9

Environmental level

(Ov.)

(LTI)

(MTI)

(HTI)

(LTI)

(MTI)

(HTI)

(LTI)

(MTI)

(HTI)

(Ov.l)

M1 (Ov.)

0.97

0.98

0.96

0.98

0.99

0.98

0.94

0.97

0.94

0.98

M4 (LTI)

0.96

0.97

0.91

0.96

0.97

0.97

0.98

0.96

0.89

0.94

M4 (MTI)

0.99

0.97

0.98

1.00

1.00

1.00

0.94

0.99

0.96

0.99

M4 (HTI)

0.96

0.91

0.98

0.99

0.98

0.98

0.88

0.97

0.99

0.98

M10 (LTI)

0.98

0.96

1.00

0.99

1.00

1.00

0.94

0.99

0.97

0.99

M10 (MTI)

0.99

0.97

1.00

0.98

1.00

1.00

0.95

0.99

0.96

0.98

M10 (HTI)

0.99

0.97

1.00

0.98

1.00

1.00

0.94

0.99

0.96

0.99

M12 (LTI)

0.92

0.98

0.93

0.85

0.92

0.93

0.92

0.95

0.88

0.93

M12 (MTI)

0.98

0.97

0.99

0.97

0.99

1.00

0.99

0.94

0.98

0.99

M12 (HTI)

0.96

0.91

0.98

1.00

0.98

0.98

0.98

0.86

0.97

0.98

M9 (Ov.)

0.99

0.96

0.99

0.98

1.00

0.99

0.99

0.92

0.99

0.98

Correlations between all animal above the diagonal and between the most used sires (60) below the diagonal

LTI = 10th (?0.424); MTI = 50th (0.032); HTI = 90th (0.316) percentiles of the contemporary group effects (environmental gradient)

M 1 linear animal model based on homoscedastic residual variance; M 4 one-step linear reaction norm model with homoscedastic residual variance; M 9 linear animal model with classification function grouped in discrete subclasses on heteroscedastic residual variance; M 10 two-step linear reaction norm model with classification function grouped in discrete subclasses on heteroscedastic residual variance; M 12 one-step linear reaction norm model with classification function grouped in discrete subclasses on heteroscedastic residual variance. LTI low tick infestation; MTI medium tick infestation; HTI high tick infestation; Ov. overall

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.

Declarations

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.

Open AccessThis 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.

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.