|Home | About | Journals | Submit | Contact Us | Français|
Tick-borne fever (TBF) is stated as one of the main disease challenges in Norwegian sheep farming during the grazing season. TBF is caused by the bacterium Anaplasma phagocytophilum that is transmitted by the tick Ixodes ricinus. A sustainable strategy to control tick-infestation is to breed for genetically robust animals. In order to use selection to genetically improve traits we need reliable estimates of genetic parameters. The standard procedures for estimating variance components assume a Gaussian distribution of the data. However, tick-count data is a discrete variable and, thus, standard procedures using linear models may not be appropriate. Thus, the objectives of this study were twofold: 1) to compare four alternative non-linear models: Poisson, negative binomial, zero-inflated Poisson and zero-inflated negative binomial based on their goodness of fit for quantifying genetic variation, as well as heritability for tick-count and 2) to investigate potential response to selection against tick-count based on truncation selection given the estimated genetic parameters from the best fit model. Our results showed that zero-inflated Poisson was the most parsimonious model for the analysis of tick count data. The resulting estimates of variance components and high heritability (0.32) led us to conclude that genetic determinism is relevant on tick count. A reduction of the breeding values for tick-count by one sire-dam genetic standard deviation on the liability scale will reduce the number of tick counts below an average of 1. An appropriate breeding scheme could control tick-count and, as a consequence, probably reduce TBF in sheep.
Sheep farming in Norway is based on grazing unfenced rangeland and mountain pastures in summer. Grazing such pastures do however provide welfare and production challenges with losses to predators, diseases and accidents. There is a yearly national economic loss of approximately 100 million NOK (125,000 sheep and lambs lost every year) while it is also a severe animal welfare issue [1,2].
Tick-borne fever (TBF) is stated as one of the main disease challenges in Norwegian sheep farming during the grazing season, particularly along the coast of south-western Norway . It is proposed to be an explanatory factor of the observed increase in lamb loss during the last decades . TBF is caused by the bacterium Anaplasma phagocytophilum, an intragranulocytic alpha-proteobacterium belonging to the family Rickettsiaceae, that is transmitted by Ixodes ticks [5, 6, 7]. A.phagocytophilum infects neutrophils and survives for several months by avoiding bacteriacidal defence mechanisms in immune-competent sheep [8, 9]. Lamb losses as high as 30% in a flock due to A. phagocytophilum infection are reported  and the Norwegian Food Safety Authority considers restrictions of grazing on pastures with high losses due to the severe welfare problems. Ticks and tick-borne diseases have, for a long time, also been a major concern to livestock producers in tropical and sub-tropical regions of the world [11, 12, 13].
Tick-infestation in sheep is commonly controlled by chemical treatment or insecticide, known as acaricides. The frequent use of such treatment is however costly and associated with development of resistance against such treatment, particularly in one-host ticks [14, 15]. An alternative strategy to control tick-infestation is to breed for genetically resistant or more robust animals. Theoretical arguments are raised as potential problems in selection for resistance, but a number of studies support that it can be sustainable, feasible and desirable . Genetic variation in resistance is shown in many farmed species, where numerous diseases are involved  i.e., gastrointestinal nematode infections , mastitis , foot rot , ectoparasites, i.e., flies and lice [21, 22] and scrapie  in sheep. Various levels of host resistance to tick-infestation are found to occur in different breeds of cattle and have been implemented in breeding schemes [24–26]. Individual variation in response to A. phagocytophilum infection in sheep is evident and shown by Granquist et al.  and Stuen et al. . Furthermore, genetic variation in lamb survival on tick exposed pastures has also been reported . This variation in response to infection might include genetic variation in resistance or susceptibility to tick-borne infections. The risk of being infected is likely to increase with number of ticks infested as approximately 8.8% of the ticks are infected with the bacterium . Hence, tick-count on lambs may, to some degree, reflect the susceptibility of A. phagocytophilum infection in lambs. However, to our knowledge, genetic variation of tick-count of the Ixodes ricinus on sheep has not been studied. For cattle, host resistance to ticks has been shown to be under genetic control  and tick count is suggested to be an appropriate method to measure tick resistance or susceptibility [32, 33]. Resistance of cattle to ticks is heritable and responsive to selection  and heritability estimates for resistance to ticks on cattle range from 0.05 to 0.42 [35–37].
In order to use selection to genetically improve traits we need reliable estimates of genetic parameter. The standard procedures using linear models for estimating variance component assume a Gaussian distribution of the data . However, tick-count data is a discrete variable and, thus, such standard procedures may not be appropriate. One potential alternative could be to use nonlinear or generalized linear models with discrete link functions, such as the Poisson or negative binomial distributions. Data with a Poisson distribution have equal mean and variance whereas negative binomial allows the data to be overdispersed. Moreover, the nature of tick-count data expresses a higher-than-expected incidence of zero scores, but this can be accommodated with zero-inflated versions [39–41] of the above described distributions.
In addition, response to selection for tick-count in sheep has not been studied, and our effort is to study the potential of selection for tick-count to consider selecting for it in a breeding program. So far, only one theoretical study on prediction of response to selection for Poisson distributed traits is published . Thus, such prediction may deviate from the response to the truncation selection assuming the estimated genetic effects here or when the distribution of tick-count data deviates from the conventional Poisson distribution.
Therefore, the objectives of this study were two folds: 1) to compare four different statistical models, i.e., Poisson, negative binomial, zero-inflated Poisson, and zero-inflated negative binomial based on the their goodness of fit for quantifying within-breed genetic variation, as well as, heritability for tick-count in lambs on tick-infested pastures, and 2) to investigate potential response to selection against tick-count based on truncation selection given the estimated genetic parameters from the best fit model.
Furthermore, following Foulley and Im ’s methodology, we derived a heritability (observed scale) while accounting for additional random effects in the non-linear model which it has not been accounted for when deriving the heritability in the previous study .
All procedures were in accordance with ethical standards as regulated by Regulation for use of animals in trials (FOR-2015-06-18-761) which is in accordance with EU directive 2010/63/EU on protection of animals that are used for scientific purposes. The regulation is administrated by the Norwegian Food Safety Authorities (www.mattilsynet.no). The procedure of counting ticks on lambs did not need official approval as of current regulations at start of study. All sheep included in the study were managed, handled and cared for in the same manner as in regular sheep production systems. Catching and holding sheep for examination is a common and normal procedure in sheep production, e.g. for shearing, treatment, health check. Counting ticks on sheep did not cause any additional pain or discomfort apart from the routine care of being caught and held for visual examination. Handling sheep for tick count was conducted by experienced sheep handlers to keep handling stress at a minimum. The sheep handlers were not certified veterinarians, but experienced sheep handlers.
The study was conducted in 2011, 2012 and 2013 on 6 sheep farms in tick endemic areas on the west coast of Norway in Vestnes (62°37'16.5"N 7°05'23.0"E), Rauma (62°35'15.2"N 7°41'10.0"E) and Fræna (62°51'13.3"N 7°09'14.5"E) municipalities where ticks and losses to TBF had been observed earlier. Presence of ticks was confirmed by examining the pastures for questing ticks with the cloth lure method  at the same day as counting ticks on the lambs. The lambs and their mothers were not treated with acaricides until after the last tick count was conducted. Tick counts were recorded on 555 lambs (Norwegian White Sheep breed) grazing on 12 different fenced pastures during spring. Lambs were sired by 87 rams that were mated with 283 dams. Tick-counts were repeatedly recorded once on the same lambs on approximately 1 and 2 weeks after turn out to pasture from the winter indoor feeding period. Lambs were turned out on spring pasture between 2–18 May on the various farms and pastures in the study. First tick count was conducted 7 days after turn out on spring pasture; i.e. between 9–25 May. Second tick count was conducted 6–7 days after first tick count; i.e. 16 May–1 June. Ewes and lambs were turned out on summer range pasture in early June for approximately 12–15 weeks. Tick counts were observed on the head, armpit, and groin of the lamb. The average number (and standard deviation) of tick counts were 1.35 (2.06). The distribution of ticks on the different sites and body locations is presented in Table 1. The raw distribution of tick counts is presented in Fig 1, ranging between 0 and 21, with 49% of 0 tick, 46% of 1–5 ticks, and 5% of observations with greater than 5 ticks.
Data were analysed using Poisson and negative binomial sire-dam models and their zero-inflated versions. Zero-inflated models assumed the data were generated from a mixture of two distributions although the population membership was not observed. The first level of the hierarchy, the probability mass function of the response Yi (tick-count) was:
where Yi has a probability mass function corresponding to a Poisson or negative binomial distribution defined by a set of parameters θi, which is assigned with a probability (1 − η) and a degenerate distribution supported at zero with probability η. Further, the standard or non zero-inflated version of the models is generated by setting η = 0.
The mean and the variance of the zero-inflated are:
For the Poisson model θi = (λi), the probability mass function is:
with mean and variance E(Yi|λi) = Var(Yi|λi) = λi.
For the negative binomial model θ = (ω, κ) and the probability mass function is:
where ω is the mean and κ is the shape parameter which approaches to the Poisson distribution as κ → ∞, and with E(Y|ωi, κ) = ωi and .
At the second level of hierarchy, we define and for the Poisson and negative binomial model respectively. Further n is the number of tick count observations. The underlying linear model for ln λ and ln ω are:
where bλ(bω) is the vectors of fixed systematic effects (record-farm and pasture-farm with 12 levels each). The uλ(uω), pλ(pω), and mλ(mω) are the vectors of sire-dam genetic random effects, permanent environmental random effects and maternal environmental random effects, respectively. The X, Zs, Zd, W and T are the corresponding incidence matrices. The link function depends on the distribution of the parameters: Poisson, negative binomial, or the zero-inflated version of those. Prior distributions for u, p and m were the following multivariate Gaussian (MVN) distributions:
In addition prior distributions of the variance components and systematic effects was bounded uniform. Moreover, prior distributions of κ in the negative binomial model was also assumed bounded uniform. Finally, the prior distribution for η in the zero-inflated models was uniform between 0 and 1.
Models were implemented using a Gibbs Sampler McMC algorithm  that involved an updating sampling scheme from the conditional posterior distributions of all the unknowns in the model. Here, conditional posterior distributions of location parameters (b, u, p and m) where unknown and single-site Metropolis-Hasting updates with random walk proposal were implemented . Further, conditional posterior distributions for variance components were scaled inverted chi-square distributions. The Gibbs sampler was implemented with a single long chain of 1,250,000 after discarding the first 250,000. Convergence of the McMC chains was explored by visual inspection of the chain.
Models were compared by means of the pseudo log-marginal probability of data (logCPO) and the Deviance Information Criterion (DIC).
Log-marginal probability (logCPO): If we consider the data vector y = (yi,y−i), where yi is the ith datum and y-i is the vector of data with ith datum deleted, the conditional predictive distribution has a probability density equal to:
where θ is the vector of parameters. Therefore, p(yi|y−i) can be interpreted as the probability of each datum given the rest of the data, and it is known as the conditional predictive ordinate (CPO) for the ith datum. The pseudo log-marginal probability of the data is then:
A Monte Carlo approximation of the logCPO suggested by Gelfand  is:
and N is the number of the McMC draws, and θj is the jth draw from the posterior distribution of the corresponding parameter. The higher the value of the LogCPO, the better the model fit to the data.
Deviance Information Criterion: The Deviance Information Criterion (DIC) was defined by Spiegelhalter et al. . It compares the global quality of 2 or more models accounting for model complexity. For a particular model M, the DIC is defined as:
where is the posterior expectation of the deviance D(θM), and is the deviance evaluated at the posterior mean estimate of the parameter vector (θM). The computation of DIC is composed by two terms, i.e., is a measure of model fit and is related with the effective number of parameters. Models with smaller DIC exhibit a better global fit after accounting for model complexity. Differences in DIC of more than 7 units are considered important by Spiegelhalter et al. .
For the sire-dam model, the estimated variance component for sires was equal to the variance component for dams and equal to one quarter of the additive genetic variance or . Therefore, the additive genetic variance for tick-count on the liability scale was estimated as .
Further, the phenotypic variance () was calculated, following Foulley and Im  as:
and the residual variance () was calculated as:
Under the Poisson distribution:
And under the negative binomial distribution,
where the nr is the non-residual components and the .
Finally, and following Foulley and Im , the additive genetic variance on the observed scale () can be achieved by:
Therefore, the heritability on the observed scale is computed as:
After the model comparison, the response to selection for tick-count was calculated based on the best-fitted model. In order to quantify the potential application of the procedure for breeding for low tick number, we computed the expected number of ticks for the progeny of the individuals with the bottom 5, 10 and 20% breeding values. For comparison and theoretical interest, top 5, 10 and 20% individuals were also selected. Expected response was calculated for the total population. Given the non-linear nature of the model, calculations were performed for several reference points of average tick number (0.5, 1.0, 1.5, 2.0, 2.5 and 3.0). For each group, we calculated the average posterior mean estimate of the breeding value () and, later on, the expectation of the progeny of the selected group computed as:
where η is the posterior mean estimate and , where r is the reference point (0.5, 1.0, 1.5, 2.0, 2.5 and 3.0).
The most frequent model for analysis of animal breeding data is the standard linear model , or the threshold model for categorical data . In this study, we discarded its implementation given the clear discrepancy with the Gaussian distribution and the wide range of categories of data (Fig 1). Thus, we restricted our comparison to Poisson  and negative binomial  models and its zero-inflated versions [41,46] using a Bayesian approach. The interpretation of zero-inflated models in the tick count data implies that the records are generated from a mixture of two distributions. First, a Bernoulli distribution that determines whether the individual has been in contact with ticks or not and, second, a Poisson (P) or negative binomial (NB) distribution that affects tick count only if the individual has been in contact with ticks.
The results of the variance component estimation and the model comparison are presented in Table 2. Both LogCPO and DIC criteria indicated that zero-inflated versions were better in goodness of fit than nonzero-inflated ones. Moreover, the zero-inflated Poisson (ZIP) model (logCPO = -1435.33: DIC = 2805.54) fitted the data slightly better than zero-inflated negative binomial (ZINB) model (logCPO = -1435.66: DIC = 2807.58), as the model is slightly more parsimonious avoiding the estimation of the parameter κ. In fact, considering the DIC, a reduction of more than 7 units usually indicates a significant improvement of goodness of fit . In this study, we found that ZINB model results in considerably lower DIC by 46.36 units than NB model. Similarly, ZIP model also provides 14.82 lower units of DIC than conventional P model. Hence, the most fitted model in this study was ZIP model.
These results can be confirmed based on the average posterior mean estimates of λ, ω and κ under different models. First, we will compare the P and NB models. The average location parameter of negative binomial distribution (ω = 1.359) was very similar to the average Poisson parameter (λ = 1.349), but the variance of the negative binomial distribution was greater than λ in Poisson distribution as the posterior mean estimate of the κ parameter (12.89) was far away from infinity. This indicated the presence of over-dispersion in the dataset. Consequently, based on logCPO, NB model (logCPO = -1441.58) fitted the data better than the P model did (logCPO = -1460.21), indicating that accounting for over-dispersion can improve goodness of fit (Table 2). Although we implemented Bayesian process, our results based on logCPO is in line with the simulation study by Tempelman and Gianola  who used marginal maximum likelihood methodology to compare P and NB models and found that estimates under NB model were less biased and yielded lower mean square error. However, DIC in our study indicated the opposite because the Poisson model (DIC = 2820.4) fitted the data better than the NB model (DIC = 2853.9). These contradicting results may be due to the effect of excessive presence of zeroes in tick counts, which it could not be appropriately accounted for with conventional Poisson and NB models. Nevertheless, when the link function takes into account the possibility of an excess of zeroes by the zero-inflated version the distributions, the average posterior estimate of ω moved towards λ (1.421), while the κ increased substantially, suggesting that the variance of the zero-inflated negative binomial distribution was moving closer to the ω. This yielded similar distribution to ZIP, and the ZIP model was selected as the best when using logCPO and DIC approaches. Both the ZIP and ZINB models yielded a posterior mean estimate of 0.042 for the probability of zero (η), proving that even a low probability of zero ticks improves considerably the goodness of fit of models.
The model comparisons for count data has been studied empirically and reported mainly in cattle and sheep [53–55]. The use of Poisson model was often the first choice for fitting the count data in comparison to a linear or threshold model. Nevertheless, most of the analysis was dedicated to litter size data, and they provide a better adjustment for threshold and linear models than for Poisson, such as in the study of Olesen et al. . However, the raw distribution of tick count data (Fig 1) diverges clearly from the distribution expected for litter size. In fact, the only precedent in analysing tick count data was a study of crossbreed Hereford x Nellore cattle  suggesting that the Poisson model fits better than linear model. However, there was also a clear excess of zeroes (48.6%) in that data set, indicating that zero inflated distributions could have performed even better as shown by the model comparison analysis described above.
We also performed some additional model comparisons with the same models (P, NB, ZIP and ZINB) while fixing in order to determine the genetic determinism of tick count. The results are presented in S1 Table. We found that DIC decreased by 7.24 units for P model, 12.38 units for NB model, 6.82 units for ZIP model and 7.42 units for ZINB and also the logCPO estimates were worst for the model with null sire-dam genetic variation. These results confirmed that was relevant and different from zero for tick counts in lambs.
Different models did neither noticeably influence the estimates of variance components on the observed nor the liability scales (Table 2). In Fig 2, the posterior distributions of the variance components () and heritability of the observed scale calculated using Eq 20 and for the selected model (ZIP) are presented. The average of over McMC samples was 0.114 with a standard deviation of 0.055, confirming the results of the model comparison with the reduced models () presented above. At the liability scale, the estimate of additive genetic variance () under model ZIP was 0.456, whereas the posterior estimate of additive genetic variance () on the observable scale was higher (0.844). The reason for lower variation at the liability scale is the logarithm transformation of the location parameters of both Poisson and negative Binomial distribution.
The mean posterior estimate of h2 was 0.320 with a posterior standard deviation of 0.162. For ticks in sheep there is, to date, only a published estimate of h2 for body louse (Bovicola ovis) in Romney sheep. The h2 for body louse score (loge[louse score +1]) varied from low to moderate, depending on the time of measurements, i.e., January 2002: h2 = 0.13, March 2002: h2 = 0.17 and November 2004: h2 = 0.37 . The genetic parameters for tick count in sheep are still lacking in the literature. In cattle, the h2 of tick count has been estimated and range widely; from low estimates of h2 (0.11 to 0.13) for the cattle tick Boophilus microplus in crossbred cattle [57,58] to moderate and high estimates of h2 (0.17 to 0.42) [59–61]. Our estimate of h2 for tick count is in the range of the previous studies in sheep and cattle, indicating that accuracy in predicting genetic effects for tick count is relatively high and the tick load can be reduced through selective breeding.
The estimate of h2 for tick count in both sheep and cattle varied due to different statistical models and trait definitions. For the trait definition, we used the sum of tick counts from three positions, which are very likely to be in contact to ticks during grazing season. With this approach, it may not represent the total number of ticks on each lamb as applied in the other studies. Yet, variation of ticks in the hot spots per se should still reflect the variation in susceptibility or a general robustness to tick. Counting ticks of the whole body is more labour-intensive and time consuming than just counting at only three positions. Hence, a future study should quantify the genetic correlation between tick counts from these three positions and the whole body to verify if they can be considered the same trait for which it may be selected using a selection index.
In addition, the results proved that the maternal care may be a very important factor that determines the variation in tick count. The mean of maternal environmental variance over McMC replicates (Fig 2) was high ( = 0.183) and even higher than . We did not quantify maternal genetic variance but in cattle, it has shown to be very low (0.004 or only 2% of phenotypic variance). The maternal environmental variance may also reflect that the ewes take their lambs to different parts of the pastures with different levels of tick infestation. Maternal care is a temporary environmental factor that can be improved and taken into account by farming management to reduce tick infection.
In order to quantify the relevance of the genetic variation in a practical breeding scheme, we calculated the expected number of tick count for the top and bottom individuals of the population. The results are presented in Table 3. Given the non-linear nature of the model, the expected selection response will depend on a reference point defined by the trait mean. Thus, a selection of 5% downward would imply a reduction of ticks from 0.064 to 0.383 depending on the reference point (0.5 to 3 ticks) and a selection of 5% upward will increase the number of ticks from 0.084 to 0.502. The asymmetry of selection response is also associated with the non-linear nature of the model, which is coherent with the theoretical predictions of Foulley . Thus, even when the breeding values were symmetric in the log scale the response of the observed scale is asymmetric. To illustrate this fact, we present the distribution of the simulated breeding values given the posterior mean estimate of in the ZIP model on the observed scale with reference means of 0.5 and 1.5 (Fig 3).
Given the mean of tick count of this studied population of 1.350, the reduction of the breeding values on the liability in one sire-dam genetic standard deviation (σsd = 0.337), will imply a reduction of the number of ticks of approximately 0.37 leading to an average tick count lower than 1. Such expected genetic gain will considerably reduce the need for treatment with insecticide and thereby likely reduce the incidence of tick-borne infections.
To apply a selective breeding program, a trait needs to be routinely measured. This could imply a challenge as counting tick is labour intensive. Furthermore, to get accurate tick counts, sheep farmers cannot use acaricides. This is still not practically feasible at most commercial sheep farms with tick-infested pastures. In addition, the maternal environmental effects may be important as the ewes’ care of their lambs may be a crucial factor determining the number of ticks on their lambs. Future studies may also investigate maternal genetic effects on tick-count to be included in a selection index. However, a more efficient system to record the trait has to be developed to select accurately. The use of genomic selection may unravel the genetic potential of this trait and enables the ability to predict breeding values without the need to record phenotypes of selection candidates on a routine basis . Thus, genomic selection may provide more accuracy to selection on maternal traits. However, knowledge about the genetic architecture of this trait is needed to design a breeding scheme using genomic selection against tick load in lambs.
Zero-inflated Poisson model is the most parsimonious for fitting tick count data. Genetic variation and heritability estimates for number of ticks on Norwegian lambs suggest potential for improving tick resistance through selective breeding. A reduction of the breeding values by one sire-dam genetic standard deviation on the liability scale may reduce the number of tick counts below an average of 1. The maternal care may be an important environmental factor affecting the number of ticks on lambs, and should be accounted for. Further, knowledge on the correlation between tick infestation and risk of developing tick-borne disease needs to be established.
This study was part of the TICKLESS project (project number 207737) which was funded by the Norwegian Foundation for Research Levy on Agricultural Products (FFL) and the Agricultural Agreement Research Funds (JA), the Møre og Romsdal Council and County Governor in Møre og Romsdal. We thank the farmers and technicians at NIBIO for collecting data of tick count on lambs.
This study was part of the TICKLESS project, which was funded by the Norwegian Foundation for Research Levy on Agricultural Products (FFL) and the Agricultural Agreement Research Funds (JA) (https://www.slf.dep.no/no/fou-midler/landbruks-og-matforskning), grant nr 207737, together with Møre and Romsdal Council (www.mrfylke.no), grant nr 127/2009, and the County Governor in Møre and Romsdal (www.fylkesmannen.no/More-og-Romsdal/), grant nr 2010/5782/OTLO/. PSL, LG and IO received the funding. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Data are available from the Norwegian Institute of Bioeconomy Research (NIBIO) for researchers who meet the criteria for access to confidential data. NIBIO initiated and owned the project “TICKLESS”. The sheep pedigree data in the TICKLESS project belongs to the Norwegian Sheep Recording System and was provided to the TICKLESS project only for this study after an application. The tick count data were collected from different commercial sheep farms in Norway and belongs to NIBIO. Hence, from a commercial point of view, the tick count data itself is confidential. On the research point of view, pedigree data is owned by the Norwegian National Sheep Recording System Norway and is used only within NIBIO for research purposes. The first author (Panya Sae-Lim) signed the confidentiality agreement to perform the analysis. Future interested researchers are able to obtain the dataset of tick counts in the exact same way as the first author (Nofima) did by contacting Norwegian Institute of Bioeconomy Research (NIBIO), PO. Box 115, NO-1431 Ås or firstname.lastname@example.orgI. The pedigree information is available from Norwegian Sheep Recording System (email@example.comM).