COMPARATIVE STUDY OF CATTLE TICK RESISTANCE USING GENERALIZED LINEAR MIXED MODELS

 ABSTRACT: Comparison of tick resistance in Bos taurus indicus (Nelore) and Bos taurus taurus (Simmental and Caracu) subspecies was investigated utilizing generalized linear mixed models (GLMMs) with Poisson and Negative binomial distributions. Nelore animals (NE) are known to present greater resistance than B. t. taurus. Difference between tick resistance in Simmental (SI) and Caracu (CA) breeds has never been reported previously. Three artificial tick infestations were conducted to evaluate tick resistance in these breeds. The statistic point of the present study was to show alternative models for the evaluation of tick count data, the GLMMs. Analysis for tick resistance by GLMM with Negative binomial distribution has never been assessed previously. The analyses were performed by the use of the PROC GLIMMIX procedure of the SAS program. The results showed that GLMM with Negative binomial distribution is appropriated to evaluate tick count data with excess of zero observations avoiding overdispersion problems. Finally, considering multiple comparisons with the Bonferroni test, different pattern of tick infestation was observed for the studied breeds, suggesting that NE is the most resistant breed followed by CA.


Introduction
The cattle tick, Rhipicephalus microplus, is a blood-sucking parasite originally from Asia, however, it may be found largely distributed between the parallel 32ºN and 32ºS (JONSSON et al., 2000), affecting bovine herds from America, Africa, Asia, and Oceania.Tick infestations are a challenge for cattle livestock production since it causes significant economic losses such as a decrease in body weight, milk yield and diseases.
The tick R. microplus can be found across the Brazilian territory (GONZALES, 1993), and it is one of the major causes of detriment for cattle producers (PARIZI et al., 2009).Brazil has the biggest commercial cattle herd in size with nearly 212 million animals (IBGE, 2015).
The bovine genetic group may have an important impact on the level of tick infestation (WAMBURA et al., 1998).Zebu cattle (Bos taurus indicus) are resistant, while less productive, and have been used as an alternative for Brazilian producers to minimize parasitism problems.High proportion of B. t. taurus used in crossings generates an increase in tick susceptibility of the progenies (BONSMA, 1944;UTECH et al., 1978;UTECH and WHARTON, 1982;LEMOS et al., 1985;OLIVEIRA and ALENCAR, 1987;OLIVEIRA et al., 1989;OLIVEIRA and ALENCAR, 1990;SILVA et al., 2007;BIANCHIN et al., 2007;OLIVEIRA et al., 2013).Adaptation to tropical environmental conditions and its contact with parasites since early domestication (O'NEILL et al., 2010) may have contributed to development of natural resistance against tick in zebu cattle.
Identification of the factors influencing bovine tick resistance could help selection of the most resistant animals.For this purpose, the use of appropriated statistic methodologies are needed.Methods often used to evaluate tick count data are linear models using data log transformation (WHARTON et al., 1970;SEIFERT, 1971;MARUFU et al., 2011;IBELLI et al., 2012;BIEGELMEYER et al, 2015;MAPHOLI et al., 2016) or Poisson models (VAZQUEZ et al. 2009;PEÑAGARICANO et al. 2010;AYRES et al., 2013).
Linear models require normal distribution; consequently, it uses data transformation in analysis.Models with Poisson and negative binomial distributions plus random effects fit into generalized linear mixed models (GLMMs) (FITZMAURICE et al., 2004).The use of mixed models turns possible identification of covariance pattern for an individual due to the inclusion of animal random effect in the model.Usually studies for investigation of tick resistance use repeated measures on the same animal, that is, more than one observation are performed for an animal.This property promotes correlation structure between measures of the same animal in the data.GLMMs allow fitting this kind of data.Furthermore, GLMMs models allow explaining data overdispersion and using adequate distribution for the variable response, in this case, the tick count observations.The objective of this study was to compare tick resistance of Bos taurus indicus (Nelore) and Bos taurus taurus (Simmental and Caracu) subspecies raised in the southeast of Brazil utilizing generalized linear mixed models (GLMMs) with Poisson and negative binomial distributions.

Animals
Sixty heifers, with average weight of 340.2 ± 63.0 kg and average age of 450 days, were submitted to three artificial infestations with larvae of R. microplus.Animals were raised in the southeastern region of Brazil, an endemic area for the bovine tick.The animals were Nelore (NE, N=20), Simmental (SI, N=20) and Caracu (CA, N=20) breeds, belonging to B. t. indicus, B. t. taurus e B. t. taurus adapted, respectively.The criteria to select animals to study in this experiment was based on the measurement of skin thickness (ST), as described in the Skin thickness item (section 2.2).
The NE and CA animals remained at the same farm during the experiment, at the Centro Avançado de Pesquisa Tecnológica dos Agronegócios (APTA) -Bovinos de Corte, Instituto de Zootecnia de São Paulo located in Sertãozinho (21º10' south latitude and 48º5' west longitude), São Paulo, Brazil.The SI animals remained at other farm located in Avaré (23º00' south latitude and 48º5' west longitude), São Paulo, Brazil.This region is characterized by a tropical humid climate with conditions that increase the occurrence of the tick R. microplus throughout the year.Both farm treated the heifers with a commercial acaricide containing cypermethrin (Cypermil 5% Pour on, Ouro Fino, 100 ml / 100 kg).
The NE cattle belongs to the B. t. indicus group and is the most used breed in Brazil.The first NE animal arrived in Brazil in the decades of the nineteenth century with the first registered animal from India.Animals of this breed are characterized mainly by their tolerance and adaptation to the Brazilian tropical conditions.
The SI cattle belong to B. t. taurus group and they are originate from Simme Valley of Switzerland, country of continental climate.Animals from this breed can be found on all continents.The breed was introduced in Brazil in the first decades of the twentieth century and, nowadays, can be found across the Brazilian territory due to its good capability for beef and milk production, being frequently used in crossbreeding with zebu cattle.
The CA breed is originated from the crossing of Portuguese and Spanish breeds and belongs to the B. t. taurus group.Some animals are considerably adapted to the tropical climate, since the first ancestors of the CA breed were introduced in Brazil in the sixteenth century, before zebu cattle.Natural selection for CA cattle and selection of the strongest and more tolerant animals by producers were factors that had influence on cattle adaptation to the tropics (LIMA et al., 1992).

Skin thickness
Ten heifers with the smallest ST and 10 with the greatest ST of each breed were selected.The purpose of measuring the skin thickness was to evaluate the influence of either thick skin or thin skin on cattle's tick resistance.The means and standard deviation of ST according to the breeds are shown in Table 1.The highest measure of ST was observed for an animal from the CA breed (17.36).This breed also presented the highest ST mean in comparison to the other breeds.The minimum value of ST was observed for the SI breed (6.43), which presented the lowest mean value as well when compared to the breeds.Values of ST for different cattle breeds are not well described in the literature.
ST was measured with a millimetre calliper in the region posterior to the animal's scapula.Skinfold measurement is the method used and accepted for the evaluation of ST (DOWLING, 1955).The posterior scapular region has been indicated as the best site to measure this trait because the thickness of the skin is relatively uniform in this area (WESONGA et al., 2006).

Artificial infestations and tick counts
Three artificial infestations were realized successively at 14-days intervals, using larvae between 15 and 21 days old.The artificial infestations with R. microplus larvae started 30 days after the animals were treated with a commercial acaricide.The animals did not receive any further acaricides until the end of the experiment.
To obtain the infective larvae, only engorged R. microplus females with weights ranging from 160 to 300 mg were selected, since they are considered optimal for oviposition (BENNETT, 1974).The female ticks were washed and placed in sterile Petri dishes and incubated in a BOD at 27 ± 1•C and relative humidity above 80% for oviposition (15 days).The eggs were then weighed in aliquots of 1g, each containing about 20,000 eggs, and placed in 20-mL syringes (GONZALES, 1993).The syringes containing eggs were incubated again in the BOD under the same temperature and humidity conditions described above, until hatching of the larvae (14 days).Only syringes that showed hatching rates above 90% by visual inspection were used.
The animals were infested by securing them in a squeeze chute and emptying the entire contents of the syringes along the animal's lumbar region.The number of engorged females with at least 4.5 mm in diameter was counted on the left side of each animal from the 19th through the 23rd day after each infestation.The method of counting engorged females in one side of the animal's body was recommended by WHARTON and UTECH (1970) and is the method used to evaluate tick resistance in cattle.To avoid double counting of tick females, they were removed.Days of artificial infestation and days of tick count were displayed on a timeline to facilitate visualization of the experimental schedule (Figure 1).

Models
Tick count data was analysed by GLMMs with Poisson distribution and negative binomial distribution using the PROC GLIMMIX of the SAS software, version 9.3.The Gauss-Hermite Quadrature method was used for the estimation of the Maximum Likelihood, with number of quadrature points equal to 100.The total tick count (TTC) was considered as the response variable, comprising the sum of daily counts for each animal in each infestation.GLMMs models were used because the TTC variable is not normally distributed and due to the study being designed as repeated measures, assuming that the animal observations are correlated.
It was assumed that the tick infestation might be influenced by The linear predictor for the models is given by: where   is the response variable total tick count of the i-th animal at moment j, with  = 1, ⋯ , 60,   = 1,2,3;   is the expected number of tick on the animal i at the moment j,  is the scale parameter for the negative binomial model,   is the random intercept, and   is estimative of infestation variability of animal i.
The covariates of the models were defined as: Initially, six different models were tested.The models were characterized by incorporating the random intercept effect with different combination of fixed effects (B, ST, M) and interaction effects (B×ST, B×M) as well as different adjust for distribution (Poisson and negative binomial).The combination of fixed effects and interaction effect for each model are shown in table 2. The Akaike Information Criterion (AIC) was used to choose the most appropriate model.The smaller the value for AIC, the better is the adjustment of the model.Significant effects in the models were identified using the F test.The Bonferroni test was applied for multiple comparison of the estimated means (HSU, 1996).

Results
Differences in individual profiles of susceptibility were observed over time among the breeds (Figure 2).Visually, the NE heifers presented less number of ticks than B. t. taurus heifers, with low infestation remained on the three moments.Animals from NE breed exhibited similar patterns among themselves, the same happened with animals of SI breed.However, the amount of tick in animals SI was characterized as low in the first moment with an increase in the second and third moments.For CA breed, high tick resistance was observed for some animals, in addition, the pattern of tick infestation that can be observed for animals of this breed was low, intermediary and high.In general, the CA animals presented their infestation rates over time.Distribution of number of ticks by B and ST (Figure 3) showed that SI and CA breeds have the greater TTC and greater variability for ST was observed to the animals with thin skin.For the NE breed, lower variability was observed between the two groups of ST, however the animals with thick skin have the greater tick count.
The Model I had the lower AIC value among the models with Poisson distribution with value of 2187.40 (Table 2).For this model, the Pearson statistics was 5.61, indicating an overdispersion problem.From this point, we adjusted models with negative binomial distribution.The model VII had the lowest AIC value (1478.77).For the model VII, the Pearson statistics was 0.67, which can be considerate a value closest to the ideal (1).The model VII was classified as the best model based on the AIC criteria.In Figure 4, Q-Q plots of the selected models (model I and model VII) are shown.The Pearson conditional residuals for model VII (Figure 4B) fitted better to the line, indicating that the negative binomial distribution fitted better to the data.In addition, the Pearson statistic was closer to one for this model.The effects of the bread (B), Moment (M) and interactions were significant (P<0.05) in the chosen model (Table 3).Multiple comparisons for the Bonferroni test ( = 0.05) showed differences (P< 0.05) for the tick count among NE, SI and CA breeds (Table 4).The NE heifers presented the lowest tick infestation, followed by CA heifers with intermediary average, and SI animals with the highest average.The highest average for the SI breed is related with greater susceptibility to tick infestation when compared with the other breeds.Multiple comparisons by Bonferroni test showed differences (P< 0.05) in the interaction of B and M. The obtained results showed that SI animals presented different levels of infestation (P< 0.05) comparing the three moments.The NE animals presented low levels of infestation in the three moments ranging from 3.65 to 8.20.The same pattern can be observed in the individual profiles (Figure 2).Statistically, TTC increased over time for SI breed.The infestation levels were similar over time for the CA breed. 1 Means in the same row followed by the same lowercase letters do not differ from one another (P< 0.05).
2 Means in the same column followed by the same uppercase letter do not differ from one another (P< 0.05).

Discussion
The NE animals had, on average, TTC significantly lower (P< 0.05) than B. t. taurus animals.The greater resistance attributed to the B. t. indicus animals is well widespread and established in the literature (VILLARES, 1941;BONSMA, 1944;UTECH et al., 1978;UTECH and WHARTON, 1982;SILVA et al. 2007;BIANCHIN et al., 2007) especially for NE breed.Moreover, the mechanism involved in this resistance is not well understood.Speculation about the resistance of NE breed occurs around specific morphologic traits that could have an influence on the parasite responses.According to Ibelli et al. (2012), adaptative factors such as animal´s hair (length, weight, and density) contribute to tick resistance of NE.Despite little knowledge about the mechanism of NE resistance to tick, elaborated hypotheses in adaptive factors such as the ones reported by Ibelli et al. (2012) could justify low quantity of tick observed in the NE animals of this study.Some authors also observed that B. t. indicus animals present less R. microplus infestation when compared to the B. t. taurus and crossed animals from this two genetics group, as we can see in our results.Silva et al. (2013) reported differences comparing the breeds Gir and Hostein; Piper et al. (2009) comparing the breeds Brahman and Holstein-Friesian; Utech et al. (1978) comparing in the breeds Brahman and British cattle.
The superiority of NE breed has been shown in studies in Brazil in comparison with crossed cattle.Silva et al. (2007) compared tick resistance in NE, Simmental × Nelore, Canchim x Nelore and Angus × Nelore, and Oliveira et al. (2013) evaluated NE and bovines such as Three-Cross (1/2 Angus + 1/4 Canchim + 1/4 Nelore) with high percentage of B. t. taurus.Silva et al. (2013) analyzed factors of risk for tick infestation, using multiple logistic regression.They studied the influence of physiological state related with season of the year, breed, lactation number and milk production in dairy cows with tick resistance.The authors identified the influence of breed in the tick count more precisely when compared to other factors.In this case, they showed that animals are more sensitive in some physiological stages and B. t. indicus animals are more resistant than B. t. taurus animals, which was also supported by our findings.
Results presented here demonstrate clear differences between the breeds in the levels of the host's resistance to tick infestations.Therefore, the results obtained are in agreement with the literature and evidence that the B. t. indicus animals are the most resistant breed and can be indicated to be raised in tropical regions, especially, with a greater risk of tick R. microplus infestation, keeping the quantity of ticks lower and avoiding further problems.Otherwise comparative parameters between Caracu animals and other B. taurus had never been provided in the literature, and in the sense of the presence of some Caracu herds in Brazil focused on beef production, more information regarding on tick infestation in this breed are needed.
Our study also provided comparative parameters between the level of tick resistance in B. t. taurus, SI and CA.Higher means of TTC can be observed for the groups such as SI and CA compared to the NE breed.High tick count in SI animals was expected, mainly because this is a breed from Continental Europe, and animals B. t. taurus do not present favourable evidence to tick resistance.Rechav et al. (1990) and Rechav and Kostrzewski (1991) reported fragility for this breed in respect with tick infestations of the genres Boophilus decoloratus and Amblyomma hebraeu.Similar to our study, these authors verified that SI breed was less resistant comparing to B. t. indicus (Brahman).
The CA breed presented lower quantity of ticks in comparison to SI breed and greater quantity of ticks comparing to NE breed.In this way, the results of this experiment showed that even though CA belonging to B. t. taurus subspecies, it has some tick resistance in relation to the climate condition in Brazil.However, some level of susceptibility to the tick must be expected, which makes the CA breed far away from the B. t. indicus pattern of resistance.
Divergence in individual profiles of tick susceptibility between the studied groups was observed over time.The phenotypic data for CA animals showed extreme profiles of tick susceptibility, some animals showed resistance in the first infestation (low TTC) and other animals presented the TTC constant during all infestation.However, in general, SI animals showed tendency to increase TTC (P< 0.05) over time (Figure 2, Table 4).The animals in our experiment had previous tick contact, we could say that CA animals were able to develop a most effective response against the tick fixation comparing to SI animals.
Comparison between the results obtained in this study and the results obtained by Fraga et al. (2003) with CA breed showed differences in the observed TTC average.The difference can be explained mainly because of the use of the different methodologies of infestation.The authors used tick count of natural infestations and we used artificial infestations that could be a strong challenge to the animals compared to the natural conditions.
The ST effect was not statistically significant (P> 0.05; Table 3) in the analyses for model VII, however interaction between B and ST was detected as significant (P< 0.05; Table 3).Differences were detected in the TTC means between NE animals with thin skin and NE animals with thick skin.For NE, heifers with thick skin had more quantity of ticks compared to the ones with thin skin (P< 0.05).Wilkinson (1962) andRiek (1962) did not find any difference (P> 0.05) between ST and tick resistance in the Australian lllawarra Shorthorn, Shorthorns and Shorthorns × Zebu, respectively.The susceptibility to the tick presented low variability within the group of Nelore animals (Figure 2), and it is evidence that a broader study could better clarify the real association between tick susceptibility and ST in the animals of this breed.However, we did not find any difference for ST among animals B. t. taurus (P> 0.05), SI and CA.
Models with negative binomial distribution are not used very often in the methodology of studies similar to our study, when TTC is the response variable.The use of models with negative binomial distribution in the present study was due to the excess of zero during the tick counts, this excess may cause overdispersion problems.The Pearson statistics indicated the chosen model with negative binomial distribution (Model VII) could be an alternative to the Poisson model because it presented lower values of AIC and presented values for the Pearson statistics close to 1 (Table 2).In this way, GLMMs with negative binomial distribution could be an alternative for TTC data with an excess of number zero as shown in this study allowing the correct identification of factors influencing cattle's tick resistance. According

Figure 1 -
Figure 1 -Experimental schedule showing steps used in artificial infestations.

Figure 2 -
Figure 2 -Individual profiles of total tick count in Nelore, Simmental and Caracu heifers in the three analyzed moments of tick count.

Figure 3 -
Figure 3 -Boxplot representing the distribution of number of ticks according to breed and skin thickness.

B
Figure 4 -(A) Q-Q plot Conditional Pearson residuals for model I with Poisson distribution; (B) Q-Q plot Conditional Pearson residuals for model VII with Negative binomial distribution.

Table 2 -
Generalized linear mixed models fitted by Poisson and negative binomial distributions B = breed; ST = skin thickness; M = moment of the tick count;  ̂ = Pearson statistics;  ̂ = estimative of infestation variability for animal and  ̂= scale parameter of the model Negative binomial.

Table 3 -
Fixed effects for the chosen model (Model VII, with Negative binomial distribution)