Home | About | Journals | Submit | Contact Us | Français |

**|**PLoS One**|**v.12(3); 2017**|**PMC5336382

Formats

Article sections

- Abstract
- Introduction
- Materials and methods
- Results and discussion
- Conclusions
- Supporting information
- References

Authors

Related links

PLoS One. 2017; 12(3): e0172711.

Published online 2017 March 3. doi: 10.1371/journal.pone.0172711

PMCID: PMC5336382

Ulrike Gertrud Munderloh, Editor^{}

University of Minnesota, UNITED STATES

**Conceptualization:**LG IO.**Data curation:**LG.**Formal analysis:**LV PSL.**Funding acquisition:**LG IO.**Investigation:**LG.**Methodology:**LV PSL.**Project administration:**LG.**Resources:**LG LV.**Software:**LV PSL.**Supervision:**LG IO.**Validation:**LG IO PSL LV.**Visualization:**LV PSL.**Writing – original draft:**PSL.**Writing – review & editing:**IO LG LV.

* E-mail: on.amifon@mil-eas.aynap

Received 2016 August 16; Accepted 2017 February 8.

Copyright © 2017 Sae-Lim et al

This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

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 [3]. It is proposed to be an explanatory factor of the observed increase in lamb loss during the last decades [4]. 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 [10] 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 [16]. Genetic variation in resistance is shown in many farmed species, where numerous diseases are involved [17] i.e., gastrointestinal nematode infections [18], mastitis [19], foot rot [20], ectoparasites, i.e., flies and lice [21, 22] and scrapie [23] 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. [27] and Stuen et al. [28]. Furthermore, genetic variation in lamb survival on tick exposed pastures has also been reported [29]. 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 [30]. 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 [31] 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 [34] 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 [38]. 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 [42]. 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 [43]’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 [43].

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 [44] 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 *Y*_{i} (tick-count) was:

$$\mathrm{Pr}\left({Y}_{i}={y}_{i}|\eta ,{\theta}_{i}\right)=\{\begin{array}{c}\eta +f\left({Y}_{i}=0|{\theta}_{i}\right)\left(1-\eta \right),\phantom{\rule{0.25em}{0ex}}{y}_{i}=0,\\ f\left({Y}_{i}={y}_{i}|{\theta}_{i}\right)\left(1-\eta \right),\phantom{\rule{0.25em}{0ex}}{y}_{i}=1,2,\mathrm{\dots ..}\phantom{\rule{0.25em}{0ex}}0\le \eta \le 1\end{array}$$

(1)

where *Y*_{i} 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:

(2)

and

(3)

For the Poisson model *θ*_{i} = (*λ*_{i}), the probability mass function is:

$$f\left({Y}_{i}={y}_{i}|{\lambda}_{i}\right)=\frac{{\lambda}_{i}^{{y}_{i}}\mathrm{exp}\left(-{\lambda}_{i}\right)}{{y}_{i}!},\phantom{\rule{0.25em}{0ex}}\lambda >0,\phantom{\rule{0.25em}{0ex}}{Y}_{i}=0,1,\mathrm{\dots .},$$

(4)

with mean and variance *E*(*Y*_{i}|*λ*_{i}) = *Var*(*Y*_{i}|*λ*_{i}) = *λ*_{i}.

For the negative binomial model *θ* = (*ω*, *κ*) and the probability mass function is:

$$f\left({Y}_{i}={y}_{i}|\omega ,\kappa \right)=\frac{\mathrm{\Gamma}\left(\kappa +{y}_{i}\right)}{\mathrm{\Gamma}\left(\kappa \right){y}_{i}!}{\left(\frac{\omega}{\kappa +\omega}\right)}^{{y}_{i}}{\left(\frac{\kappa}{\kappa +\omega}\right)}^{\kappa},$$

(5)

where *ω* is the mean and *κ* is the shape parameter which approaches to the Poisson distribution as *κ* → ∞, and with *E*(*Y*|*ω*_{i}, *κ*) = *ω*_{i} and $var\left(Y|{\omega}_{i},\kappa \right)={\omega}_{i}+\frac{{\omega}_{i}^{2}}{\kappa}$.

At the second level of hierarchy, we define $\mathrm{ln}\phantom{\rule{0.25em}{0ex}}\lambda ={\left\{\mathrm{l}\mathrm{n}\phantom{\rule{0.25em}{0ex}}{\lambda}_{i}\right\}}_{i=1}^{n}$ and $\mathrm{ln}\phantom{\rule{0.25em}{0ex}}\omega ={\left\{\mathrm{ln}\phantom{\rule{0.25em}{0ex}}{\omega}_{i}\right\}}_{i=1}^{n}$ 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:

ln*λ* = **X****b**_{λ} + (**Z**_{s} + **Z**_{d})**u**_{λ} + **W****p**_{λ} + **T****m**_{λ},

(6)

and

ln*ω* = **X****b**_{ω} + (**Z**_{s} + **Z**_{d})**u**_{ω} + **W****p**_{ω} + **T****m**_{ω},

(7)

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**, **Z**_{s}, **Z**_{d}, **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:

$$\mathbf{u}\sim MVN\left(\mathbf{0},\mathbf{A}{\sigma}_{u}^{2}\right),\phantom{\rule{0.25em}{0ex}}\mathbf{p}\sim MVN\left(\mathbf{0},\mathbf{I}{\sigma}_{p}^{2}\right)\phantom{\rule{0.25em}{0ex}}\mathrm{and}\phantom{\rule{0.25em}{0ex}}\mathbf{m}\sim MVN\left(\mathbf{0},\mathbf{I}{\sigma}_{m}^{2}\right).$$

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 [45] 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 [46]. 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** = (*y*_{i},**y**_{−i}), where *y*_{i} is the *i*^{th} datum and **y**_{-i} is the vector of data with *i*^{th} datum deleted, the conditional predictive distribution has a probability density equal to:

(8)

where *θ* is the vector of parameters. Therefore, *p*(*y*_{i}|**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 *i*^{th} datum. The pseudo log-marginal probability of the data is then:

$$\sum _{i}\mathrm{ln}}\phantom{\rule{0.15em}{0ex}}p\left({y}_{i}|{\mathbf{y}}_{-i}\right).$$

(9)

A Monte Carlo approximation of the logCPO suggested by Gelfand [47] is:

$$\sum _{i}\mathrm{ln}}\phantom{\rule{0.15em}{0ex}}\widehat{p}\left({y}_{i}|{\mathbf{y}}_{-i}\right),$$

(10)

where

$$\widehat{p}\left({y}_{i}|{\mathbf{y}}_{-i}\right)=N{\left[{\displaystyle \sum _{j=1}^{N}\frac{1}{p\left({y}_{i}|{\theta}^{j}\right)}}\right]}^{-1},$$

(11)

and *N* is the number of the **McMC** draws, and *θ*^{j} is the *j*^{th} 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. [48]. It compares the global quality of 2 or more models accounting for model complexity. For a particular model *M*, the DIC is defined as:

$$DI{C}_{M}=2{\overline{D}}_{M}-D\left({\overline{\theta}}_{M}\right),$$

(12)

where ${\overline{D}}_{M}$ is the posterior expectation of the deviance *D*(*θ*_{M}), and $D\left({\overline{\theta}}_{M}\right)=-2\mathrm{l}\mathrm{o}\mathrm{g}\left(p\left(\mathit{y}|{\theta}_{M}\right)\right)$ 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., ${\overline{D}}_{M}$ is a measure of model fit and ${\overline{D}}_{M}-D\left({\overline{\theta}}_{M}\right)$ 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. [48].

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 ${\sigma}_{u}^{2}={\sigma}_{sd}^{2}=\frac{1}{4}{\sigma}_{a}^{2}$. Therefore, the additive genetic variance for tick-count $\left({\sigma}_{a}^{2}\right)$ on the liability scale was estimated as $4{\sigma}_{u}^{2}$.

Further, the phenotypic variance (${\sigma}_{y}^{2}$) was calculated, following Foulley and Im [39] as:

$${\sigma}_{y}^{2}={E}_{\theta}\left[var\left(Y|\theta \right)\right]+{var}_{\theta}\left[E\left(Y|\theta \right)\right]$$

(13)

and the residual variance (${\sigma}_{r}^{2}$) was calculated as:

$$\begin{array}{l}{\sigma}_{r}^{2}={E}_{\theta}\left[var\left(Y|\theta \right)\right]\\ {\sigma}_{r}^{2}=\left(1-\eta \right)var\left(Y|\theta \right)+\eta \left(1-\eta \right){\left[E\left(Y|\theta \right)\right]}^{2}\\ {\sigma}_{r}^{2}=\left(1-\eta \right)var\left(Y\right)+\eta \left(1-\eta \right){\left[E\left(Y\right)\right]}^{2}.\end{array}$$

(14)

Under the Poisson distribution:

(15)

And under the negative binomial distribution,

$$E\left(Y\right)={\omega}_{i}\phantom{\rule{0.25em}{0ex}}\mathrm{and}\phantom{\rule{0.25em}{0ex}}var\left(Y\right)={\omega}_{i}+\frac{{\omega}_{i}^{2}}{\kappa},$$

(16)

$${\sigma}_{nr}^{2}=va{r}_{\theta}\left[E\left(Y|\theta \right)\right]={\left(1-\eta \right)}^{2}E{\left(Y\right)}^{2}\left[\text{exp}\left({\sigma}_{t}^{2}\right)-1\right],$$

(17)

where the *nr* is the non-residual components and the ${\sigma}_{t}^{2}={2\sigma}_{sd}^{2}+{\sigma}_{p}^{2}+{\sigma}_{m}^{2}$.

Finally, and following Foulley and Im [39], the additive genetic variance on the observed scale (${\sigma}_{\mathrm{g}}^{2}$) can be achieved by:

$${\sigma}_{\text{g}}^{2}=\frac{{\left[cov\left(\text{g},a\right)\right]}^{2}}{{\sigma}_{a}^{2}}=\frac{{\left[\left(1-\eta \right)E\left(Y\right){\sigma}_{a}^{2}\right]}^{2}}{{\sigma}_{a}^{2}}$$

(18)

$${\sigma}_{\text{g}}^{2}={\left(1-\eta \right)}^{2}E{\left(Y\right)}^{2}{\sigma}_{a}^{2}$$

(19)

Therefore, the heritability on the observed scale is computed as:

$${h}^{2}=\frac{{\left(1-\eta \right)}^{2}{\left[E\left(Y\right)\right]}^{2}{\sigma}_{a}^{2}}{\left(1-\eta \right)var\left(Y\right)+\eta \left(1-\eta \right){\left[E\left(Y\right)\right]}^{2}+{\left(1-\eta \right)}^{2}{\left[E\left(Y\right)\right]}^{2}\left[exp\left({\sigma}_{t}^{2}\right)-1\right]}$$

(20)

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 (${\overline{u}}_{SEL}$) and, later on, the expectation of the progeny of the selected group computed as:

(21)

where *η* is the posterior mean estimate and ${\theta}_{i}=\mathrm{e}\mathrm{x}\mathrm{p}(\mathrm{log}\left(r\right)+{\overline{u}}_{SEL})$, 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 [49], or the threshold model for categorical data [50]. 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 [51] and negative binomial [52] 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 [48]. 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 [52] 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. [55]. 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 [54] 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 ${\sigma}_{sd}^{2}=0$ 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 ${\sigma}_{sd}^{2}$ 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 (${\sigma}_{sd}^{2},{\sigma}_{m}^{2},{\sigma}_{pe}^{2}$) and heritability of the observed scale calculated using Eq 20 and for the selected model (ZIP) are presented. The average of ${\sigma}_{sd}^{2}$ over **McMC** samples was 0.114 with a standard deviation of 0.055, confirming the results of the model comparison with the reduced models (${\sigma}_{sd}^{2}=0$) presented above. At the liability scale, the estimate of additive genetic variance (${\sigma}_{a}^{2}=4{\sigma}_{sd}^{2}$) under model ZIP was 0.456, whereas the posterior estimate of additive genetic variance (${\sigma}_{g}^{2}$) 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 *h*^{2} was 0.320 with a posterior standard deviation of 0.162. For ticks in sheep there is, to date, only a published estimate of *h*^{2} for body louse (*Bovicola ovis*) in Romney sheep. The *h*^{2} for body louse score (log_{e}[louse score +1]) varied from low to moderate, depending on the time of measurements, i.e., January 2002: *h*^{2} = 0.13, March 2002: *h*^{2} = 0.17 and November 2004: *h*^{2} = 0.37 [56]. The genetic parameters for tick count in sheep are still lacking in the literature. In cattle, the *h*^{2} of tick count has been estimated and range widely; from low estimates of *h*^{2} (0.11 to 0.13) for the cattle tick *Boophilus microplus* in crossbred cattle [57,58] to moderate and high estimates of *h*^{2} (0.17 to 0.42) [59–61]. Our estimate of *h*^{2} 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 *h*^{2} 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 (${\sigma}_{m}^{2}$ = 0.183) and even higher than ${\sigma}_{sd}^{2}$. 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 [42]. 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 ${\sigma}_{sd}^{2}$ 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 [62]. 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.

(PDF)

Click here for additional data file.^{(199K, pdf)}

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 Availability

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 on.oibin@nesnah.regnI. The pedigree information is available from Norwegian Sheep Recording System (on.ailamina@datsyl.tiraM).

1. NIBIO, (2016) National statistics on sheep lost on range pasture in Norway. Available: http://www.skogoglandskap.no/filearchive/OBB_Tapsprosent_1970-2015.pdf. Accessed 14 July 2016.

2. Ulvund MJ (2012) Important sheep flock health issues in Scandinavia/northern Europe. Small Rum Res
106/1: 6–10.

3. Stuen S, (2003) Anaplasma phagocytophilum infection in sheep and wild ruminants in Norway. A study on clinical manifestation, distribution and persistence. Dr Philosophiae thesis Norwegian School of Veterinary Science, Oslo, Norway.

4. Stuen S (2007) *Anaplasma phagocytophilum*-the most widespread tick-borne infection in animals in Europe. Vet Res Commun
31: 79–84. doi: 10.1007/s11259-007-0071-y
[PubMed]

5. Foggie A (1951) Studies on the infectious agent of tick-borne fever in sheep. J.Pathol.Bacteriol. 63, 1–15. [PubMed]

6. Dugat T, Lagrée AC, Maillard R, Boulouis HJ, Haddad N (2015) Opening the black box of Anaplasma phagocytophilum diversity: current situation and future perspectives. Front Cell Infect Microbiol; 5:61
doi: 10.3389/fcimb.2015.00061
[PMC free article] [PubMed]

7. Stuen S, Granquist EG, Silaghi C (2013) Anaplasma phagocytophilum—a widespread multi-host pathogen with highly adaptive strategies. Front Cell Infect Microbiol.;3:31
doi: 10.3389/fcimb.2013.00031
[PMC free article] [PubMed]

8. Granquist EG, Stuen S, Lundgren AM, Braten M, Barbet AF (2008) Outermembrane protein sequence variation in lambs experimentally infected with Anaplasma phagocytophilum. Infect.Immun. 76, 120–126. doi: 10.1128/IAI.01206-07
[PMC free article] [PubMed]

9. Woldehiwet Z (2010). The natural history of Anaplasma phagocytophilum. Vet.Parasitol. 167, 108–122. doi: 10.1016/j.vetpar.2009.09.013
[PubMed]

10. Stuen S, Kjølleberg K (2000). An investigation of lamb deaths on tick pastures in Norway In: Kazimirova M, Labunda M, Nuttall PA (Eds.). Slovak Academy of Sciences, Bratislava: 111–115.

11. Mota R, Lopes P, Tempelman R, Silva F, Aguilar I, Gomes C, Cardoso F (2016) Genome-enabled prediction for tick resistance in Hereford and Braford beef cattle via reaction norm models. J. Anim Sci
94: 1834–1843. doi: 10.2527/jas.2015-0194
[PubMed]

12. Passafaro TL, Carrera JPB, dos Santos LL, Raidan FSS, dos Santos DCC, Cardoso EP, Leite RC, Toral FLB (2015) Genetic analysis of resistance to ticks, gastrointestinal nematodes and *Eimeria spp*. in Nellore cattle. Vet Parasitol
210: 224–234. doi: 10.1016/j.vetpar.2015.03.017
[PubMed]

13. Jongejan F, Uilenberg G (2004) The global importance of ticks. Parasitology, 129: S3–S14. [PubMed]

14. Muyobela J, Nkunika POY, Mwase ET (2015) Resistance status of ticks (Acari: Ixodidae) to amitraz and cypermethrin acaricides in Isoka District, Zambia. Trop Anim Health Pro
47: 1599–1605. [PubMed]

15. Rajput ZI, Hu SH, Chen WJ, Arijo AG, Xiao CW (2006) Importance of ticks and their chemical and immunological control in livestock. J Zhejiang Univ- Sc B
7: 912–921. [PMC free article] [PubMed]

16. Stear M, Bishop SC, Mallard BA and Raadsma H (2001) Review: The sustainability, feasibility and desirability of breeding livestock for disease resistance. Res Vet Science
7:1
1–7. [PubMed]

17. Bishop SC, Axford RF, Nicholas FW, Owen JB (2010) Breeding for disease resistance in farm animals. 3rd ed. CABI.

18. Stear M, Murray M (1994) Genetic resistance to parasitic disease: particularly of resistance in ruminants to gastrointestinal nematodes. Vet Parasitol
54: 161–176. [PubMed]

19. Rupp R, Bergonier D, Dion S, Hygonenq MC, Aurel MR, Robert-Granié C, Foucras G (2009) Response to somatic cell count-based selection for mastitis resistance in a divergent selection experiment in sheep. J Dairy Sci
92: 1203–1219. doi: 10.3168/jds.2008-1435
[PubMed]

20. Raadsma H, Dhungyel O (2013) A review of footrot in sheep: new approaches for control of virulent footrot. Livest Sci
156: 115–125.

21. Pfeffer A., Morris C.A., Green R.S., Wheeler M., Shu D., Bisset S.A., Vlassoff A., 2007.
Heritability of resistance to infestation with the body louse, Bovicola ovis, in Romney sheep bred for differences in resistance or resilience to gastro-intestinal nematode parasites. International Journal for Parasitology
37, 1589–1597. doi: 10.1016/j.ijpara.2007.05.008
[PubMed]

22. Raadsma H.W., 1991.
Fleece rot and body strike in Merino sheep. V. Heritability of liability to body strike in weaner sheep under flywave conditions. Aust. J. Agric. Res. 42, 279–293.

23. Dawson M, Hoinville L, Hosie B, Hunter N (1998) Guidance on the use of PrP genotyping as an aid to the control of clinical scrapie. Vet Rec
142: 623–625. [PubMed]

24. Lemos A, Teodoro R, Oliveira G, Madalena F (1985) Comparative performance of six Holstein-Friesian× Guzera grades in Brazil. 3. Burdens of *Boophilus microplus* under field conditions. Anim Prod
41: 187–191.

25. Prayaga K (2003) Evaluation of beef cattle genotypes and estimation of direct and maternal genetic effects in a tropical environment. 2. Adaptive and temperament traits. Crop Pasture Sci
54: 1027–1038.

26. Utech K, Wharton R, Kerr J (1978) Resistance to *Boophilus microplus* (Canestrini) in different breeds of cattle. Crop Pasture Sci
29: 885–895.

27. Granquist EG, Stuen S, Crosby L, Lundgren AM, Alleman AR, Barbet AF (2010) Variant-specific and diminishing immune responses towards the highly variable MSP2(P44) outer membrane protein of *Anaplasma phagocytophilum* during persistent infection in lambs. Vet Immunol Immunop
133: 117–124. [PMC free article] [PubMed]

28. Stuen S, Grøva L, Granquist EG, Sandstedt K, Olesen I, Steinshamn H (2011) A comparative study of clinical manifestations, haematological and serological responses after experimental infection with *Anaplasma phagocytophilum* in two Norwegian sheep breeds. Acta Vet Scand
53. [PMC free article] [PubMed]

29. Grøva L (2011) Tick-borne fever in sheep-production loss and preventive measures: Norwegian University of Life Sciences. Department of Animal and Aquacultural Sciences.

30. Mysterud A, Easterday WR, Qviller L, Viljugrein H, Ytrehus B (2013) Spatial and seasonal variation in the prevalence of *Anaplasma phagocytophilum* and *Borrelia burgdorferi* sensu lato in questing *Ixodes ricinus* ticks in Norway. Parasite Vector
6: 1. [PMC free article] [PubMed]

31. Hewetson R (1972) The inheritance of resistance by cattle to cattle tick. Aust Vet J
48: 299–303. [PubMed]

32. Neto LRP, Jonsson NN, Michael J, Barendse W (2011) Molecular genetic approaches for identifying the basis of variation in resistance to tick infestation in cattle. Vet Parasitol
180: 165–172. doi: 10.1016/j.vetpar.2011.05.048
[PubMed]

33. Regitano L, Martinez M, Machado M (2006) Molecular aspects of bovine tropical adaptation. Instituto Prociência. pp. 16–06.

34. Turner LB, Harrison BE, Bunch RJ, Neto LRP, Li Y, Barendse W (2010) A genome-wide association study of tick burden and milk composition in cattle. Anim Prod Sci
50: 235–245.

35. Wharton R, Utech K, Turner H (1970) Resistance to the cattle tick, *Boophilus microplus* in a herd of Australian Illawarra Shorthorn cattle: its assessment and heritability. Crop Pasture Sci
21: 163–181.

36. Burrow H (2001) Variances and covariances between productive and adaptive traits and temperament in a composite breed of tropical beef cattle. Livest Prod Sci
70: 213–233.

37. Porto-Neto LR, Reverter A, Prayaga KC, Chan EK, Johnston DJ, Hawken RJ, Fordyce G, Garcia JF, Sonstegard TS, Bolormaa S (2014) The genetic architecture of climatic adaptation of tropical cattle. PloS one
9: e113284
doi: 10.1371/journal.pone.0113284
[PMC free article] [PubMed]

38. Lynch M, Walsh B (1998) Genetics and analysis of quantitative traits: Sinauer Associates, Sunderland, MA.

39. Lambert D (1992) Zero-Inflated Poisson Regression, with an Application to Defects in Manufacturing. Technometrics
34: 1–14.

40. Naya H, Urioste JI, Chang Y-M, Rodrigues-Motta M, Kremer R, Gianola D (2008) A comparison between Poisson and zero-inflated Poisson regression models with an application to number of black spots in Corriedale sheep. Genet Sel Evol
40: 379–394. doi: 10.1051/gse:2008010
[PMC free article] [PubMed]

41. Rodrigues-Motta M, Gianola D, Heringstad B, Rosa G, Chang Y (2007) A zero-inflated Poisson model for genetic analysis of the number of mastitis cases in Norwegian Red cows. J Dairy Sci
90: 5306–5315. doi: 10.3168/jds.2006-898
[PubMed]

42. Foulley JL (1993) Prediction of selection response for Poisson distributed traits. Genet Sel Evol
25: 297–303.

43. Foulley J, Im S (1993) A marginal quasi-likelihold approach to the analysis of Poisson variables with generalized linear mixed models. Genet Sel Evol
25: 101.

44. Vassallo M, Pichon B, Cabaret J, Figureau C, Pérez-Eid C (2000) Methodology for Sampling Questing Nymphs of *Ixodes Ricinus* (Acari: Ixodidae), the Principal Vector of Lyme Disease in Europe. J Med Entomol
37: 335–339. [PubMed]

45. Gelfand AE, Smith AF (1990) Sampling-based approaches to calculating marginal densities. J Am Stat Assoc
85: 398–409.

46. Varona L, Sorensen D (2010) A genetic analysis of mortality in pigs. Genetics
184: 277–284. doi: 10.1534/genetics.109.110759
[PubMed]

47. Gelfand AE (1996) Model determination using sampling-based methods. Markov chain Monte Carlo in practice: pp145-161.

48. Spiegelhalter DJ, Best NG, Carlin BP, Van Der Linde A (2002) Bayesian measures of model complexity and fit. J R Stat Soc B
64: 583–639.

49. Henderson C (1984) Applications of Linear Models in Animal Breeding. In: Schaeffer LR, editor. 3rd ed.
Guelph University.

50. Sorensen D, Andersen S, Gianola D, Korsgaard I (1995) Bayesian inference in threshold models using Gibbs sampling. Genet Sel Evol
27: 1.

51. Foulley J, Gianola D, Im S (1987) Genetic evaluation of traits distributed as Poisson-binomial with reference to reproductive characters. Theor and Appl Genet
73: 870–877. [PubMed]

52. Tempelman RJ, Gianola D (1996) A mixed effects model for overdispersed count data in animal breeding. Biometrics: 265–279.

53. Matos C, Thomas D, Gianola D, Perez-Enciso M, Young L (1997) Genetic analysis of discrete reproductive traits in sheep using linear and nonlinear models: II. Goodness of fit and predictive ability. J Anim Sci
75: 88–94. [PubMed]

54. Ayres DR, Pereira RJ, Boligon AA, Silva FF, Schenkel FS, Roso VM, Albuquerque LG (2013) Linear and Poisson models for genetic evaluation of tick resistance in cross-bred Hereford x Nellore cattle. J Anim Breed Genet
130: 417–424. doi: 10.1111/jbg.12036
[PubMed]

55. Olesen I, Perez-Enciso M, Gianola D, Thomas D (1994) A comparison of normal and nonnormal mixed models for number of lambs born in Norwegian sheep. J Anim Sci
72: 1166–1173. [PubMed]

56. Pfeffer A, Morris CA, Green RS, Wheeler M, Shu D, Bisset SA, Vlassoff A (2007) Heritability of resistance to infestation with the body louse, *Bovicola ovis*, in Romney sheep bred for differences in resistance or resilience to gastro-intestinal nematode parasites. Int J Parasito
37: 1589–1597. [PubMed]

57. Prayaga KC, Henshall JM (2005) Adaptability in tropical beef cattle: genetic parameters of growth, adaptive and temperament traits in a crossbred population. Aust J Exp Agr
45: 971–983.

58. Ayres DR, Pereira RJ, Boligon AA, Baldi F, Roso VM, Albuquerque LG (2014) Genetic parameters and investigation of genotype × environment interactions in Nellore × Hereford crossbred for resistance to cattle ticks in different regions of Brazil. J Appl Genet
56: 107–113. doi: 10.1007/s13353-014-0238-5
[PubMed]

59. Henshall JM (2004) A genetic analysis of parasite resistance traits in a tropically adapted line of *Bos taurus*. Aust J Agr Res
55.

60. Budeli M, Nephawe K, Norris D, Selapa N, Bergh L, Maiwashe A (2009) Genetic parameter estimates for tick resistance in Bonsmara cattle. S Afr J Anim Sci
39: 321–327.

61. Ali AA, O’Neill CJ, Thomson PC, Kadarmideen HN (2012) Genetic parameters of infectious bovine keratoconjunctivitis and its relationship with weight and parasite infestations in Australian tropical *Bos taurus* cattle. Genet Sel Evol
44: 1–11. [PMC free article] [PubMed]

62. Meuwissen THE, Hayes BJ, Goddard ME (2001) Prediction of Total Genetic Value Using Genome-Wide Dense Marker Maps. Genetics
157: 1819–1829. [PubMed]

Articles from PLoS ONE are provided here courtesy of **Public Library of Science**

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |