# Additional file 1: Description of the used within herd model of Map Spread

## A Equations for the within-herd dynamic

The description of the model we used corresponds to the within-herd model described in [2] and based on [11]. The set of parameters used are described in tables 1, 2 and 3.

### Notations

In the following equations, $X_{(t,a)}$ is the number of animals in health state $X$ at time $t$ and age $a$. Age is given in weeks until the entry in the adult stage (from 1 to cal, with intermediary stages $w$ for weaning age, $y$ for young heifers age, $h$ for heifers age and $u$ for maximal age in the susceptible compartments), then by age group ($P_1$ to $P_{5+}$). Some variables can have a prefix: "b" for births in health states $X$ ($bX$), "N" for animals transiting between two health states (NX) and "s" for exits (mortality and culling) ($sX$). After entering the adult stage, flows corresponding to aging are noted using a superscript ng, whereas those remaining in the same age group are noted using a superscript sg. $N_{(t,a)}$ is the number of animals of age $a$ at time $t$. Average duration in health states are noted by $v_X$. The remaining terms used are defined when introduced.

### Equations for the updating of variables describing health states

In this section we introduce the equations for the updating of variables corresponding to the health states, for a given herd $i$.

#### Susceptible (S) and No more Susceptible (R)

$$S_{(t+1,a=1)} = bS_{(t)} - NT_{(t+1,a)}$$

$$S_{(t+1,a \in [2;52])} = [S_{(t,a-1)} - sS_{(t,a-1)}] - NT_{(t+1,a)}$$

$$R_{(t+1,53)} = S_{(t,52)} - sS_{(t,52)}$$

$$R_{(t+1,a \in [54;cal])} = R_{(t,a-1)} - sR_{(t,a-1)}$$

$$R_{(t+1,P_1)} = R_{(t,P_1)}^{sg} - sR_{(t,P_1)}^{sg} + R_{(t,cal)} - sR_{(t,cal)}$$

$$R_{(t+1,P_i \in [P_2;P_4])} = R_{(t,P_i)}^{sg} - sR_{(t,P_i)}^{sg} + R_{(t,P_{i-1})}^{ng} - sR_{(t,P_{i-1})}^{ng}$$

$$R_{(t+1,P_{5+})} = R_{(t,P_{5+})} - sR_{(t,P_{5+})} + R_{(t,P_4)}^{ng} - sR_{(t,P_4)}^{ng}$$
# Which phenotypic traits of resistance should be improved in cattle to control paratuberculosis dynamics in a dairy herd: a modelling approach

Racem Ben Romdhane, Gaël Beaunée, Guillaume Camanes, Raphaël Guatteo, Christine Fourichon, Pauline Ezanno

► To cite this version:

Racem Ben Romdhane, Gaël Beaunée, Guillaume Camanes, Raphaël Guatteo, Christine Fourichon, et al.. Which phenotypic traits of resistance should be improved in cattle to control paratuberculosis dynamics in a dairy herd: a modelling approach. Veterinary Research, 2017, 48 (1), pp.62. 10.1186/s13567-017-0468-8. [hal-01621477](https://hal.science/hal-01621477)

HAL Id: hal-01621477

https://hal.science/hal-01621477

Submitted on 23 Oct 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
Transiently infected (T)

$$T_{(t+1,a=1)} = bT_{(t)} + nT_{(t+1,a)}$$

$$T_{(t+1,a \in [2;52])} = [T_{(t,a-1)} - sT_{(t,a-1)}] - nL_{(t+1,a)} + nT_{(t+1,a)}$$

$$T_{(t+1,a \in [53;cal])} = [T_{(t,a-1)} - sT_{(t,a-1)}] - nL_{(t+1,a)}$$

Latently infected (L)

$$L_{(t+1,a=2)} = nL_{(t+1,a)}$$

$$L_{(t+1,a \in [3;h])} = [L_{(t,a-1)} - sL_{(t,a-1)}] + nL_{(t+1,a)}$$

$$L_{(t+1,a \in [h+1;cal])} = [L_{(t,a-1)} - sL_{(t,a-1)}] - nIs_{(t+1,a)} + nL_{(t+1,a)}$$

$$L_{(t+1,P_1)} = [L_{(t,P_1)}^{sg} - sL_{(t,P_1)}^{sg} + L_{(t,cal)} - sL_{(t,cal)} + T_{(t,cal)} - sT_{(t,cal)}] - nIs_{(t+1,P_1)}$$

$$L_{(t+1,P_i \in [P_2;P_4])} = [L_{(t,P_i)}^{sg} - sL_{(t,P_i)}^{sg} + L_{(t,P_{i-1})}^{ng} - sL_{(t,P_{i-1})}^{ng}] - nIs_{(t+1,P_i)}$$

$$L_{(t+1,P_{5+})} = [L_{(t,P_{5+})} - sL_{(t,P_{5+})} + L_{(t,P_4)}^{ng} - sL_{(t,P_4)}^{ng}] - nIs_{(t+1,P_{5+})}$$

Moderate shedding (Is)

$$Is_{(t+1,a=h+1)} = nIs_{(t+1,a)}$$

$$Is_{(t+1,a \in [h+2;cal])} = [Is_{(t,a-1)} - sIs_{(t,a-1)}] - nIc_{(t+1,a)} + nIs_{(t+1,a)}$$

$$Is_{(t+1,P_1)} = [Is_{(t,P_1)}^{sg} - sIs_{(t,P_1)}^{sg} + Is_{(t,cal)} - sIs_{(t,cal)}] - nIc_{(t+1,P_1)} + nIs_{(t+1,P_1)}$$

$$Is_{(t+1,P_i \in [P_2;P_4])} = [Is_{(t,P_i)}^{sg} - sIs_{(t,P_i)}^{sg} + Is_{(t,P_{i-1})}^{ng} - sIs_{(t,P_{i-1})}^{ng}] - nIc_{(t+1,P_i)} + nIs_{(t+1,P_i)}$$

$$Is_{(t+1,P_{5+})} = [Is_{(t,P_{5+})} - sIs_{(t,P_{5+})} + Is_{(t,P_4)}^{ng} - sIs_{(t,P_4)}^{ng}] - nIc_{(t+1,P_{5+})} + nIs_{(t+1,P_{5+})}$$

High shedding and clinically affected (Ic)

$$Ic_{(t+1,a=h+2)} = nIc_{(t+1,a)}$$

$$Ic_{(t+1,a \in [h+3;cal])} = [Ic_{(t,a-1)} - sIc_{(t,a-1)}] + nIc_{(t+1,a)}$$

$$Ic_{(t+1,P_1)} = [Ic_{(t,P_1)}^{sg} - sIc_{(t,P_1)}^{sg} + Ic_{(t,cal)} - sIc_{(t,cal)}] + nIc_{(t+1,P_1)}$$

$$Ic_{(t+1,P_i \in [P_2;P_4])} = [Ic_{(t,P_i)}^{sg} - sIc_{(t,P_i)}^{sg} + Ic_{(t,P_{i-1})}^{ng} - sIc_{(t,P_{i-1})}^{ng}] + nIc_{(t+1,P_i)}$$

$$Ic_{(t+1,P_{5+})} = [Ic_{(t,P_{5+})} - sIc_{(t,P_{5+})} + Ic_{(t,P_4)}^{ng} - sIc_{(t,P_4)}^{ng}] + nIc_{(t+1,P_{5+})}$$
<figure><img src="image_2.png" /></figure>

# Which phenotypic traits of resistance should be improved in cattle to control paratuberculosis dynamics in a dairy herd: a modelling approach

Racem Ben Romdhane<fnref n="1" />*, Gaël Beaunée², Guillaume Camanes<fnref n="1" />, Raphaël Guatteo<fnref n="1" />, Christine Fourichon<fnref n="1" /> and Pauline Ezanno<fnref n="1" />

## Abstract

Paratuberculosis is a worldwide disease causing production losses in dairy cattle herds. Variability of cattle response to exposure to *Mycobacterium avium* subsp. *paratuberculosis* (Map) has been highlighted. Such individual variability could influence Map spread at larger scale. Cattle resistance to paratuberculosis has been shown to be heritable, suggesting genetic selection could enhance disease control. Our objective was to identify which phenotypic traits characterising the individual course of infection influence Map spread in a dairy cattle herd. We used a stochastic mechanistic model. Resistance consisted in the ability to prevent infection and the ability to cope with infection. We assessed the effect of varying (alone and combined) fourteen phenotypic traits characterising the infection course. We calculated four model outputs 25 years after Map introduction in a naïve herd: cumulative incidence, infection persistence, and prevalence of infected and affected animals. A cluster analysis identified influential phenotypes of cattle resistance. An ANOVA quantified the contribution of traits to model output variance. Four phenotypic traits strongly influenced Map spread: the decay in susceptibility with age (the most effective), the quantity of Map shed in faeces by high shedders, the incubation period duration, and the required infectious dose. Interactions contributed up to 12% of output variance, highlighting the expected added-value of improving several traits simultaneously. Combinations of the four most influential traits decreased incidence to less than one newly infected animal per year in most scenarios. Future genetic selection should aim at improving simultaneously the most influential traits to reduce Map spread in cattle populations.

## Introduction

Bovine paratuberculosis or Johne's disease (JD) is a bacterial infection caused by *Mycobacterium avium* subsp. *paratuberculosis* (Map). It mainly affects domestic ruminants. Paratuberculosis has a worldwide distribution with a high prevalence, herd prevalence being around 50% in Europe [1]. The progressive evolution of the infection leads to a chronic diarrhoea, an emaciation and death. This infection is responsible for significant weight losses, a decrease in milk production, an increase in mortality,

and the early culling of infected animals, inducing economic losses [2]. Infectious animals shed bacteria in their faeces, milk, and colostrum. Susceptible animals are infected by ingesting Map or in utero. Calves are known to be the age group most susceptible to infection [3].

Individual response to a given exposure to Map differs among animals. Within-herd prevalence is usually low, with 2.8–27% of infected animals [4, 5]. Field observations have reported substantial variation in individual response to Map exposure: among birth cohorts which are assumed to have been similarly exposed to Map, some are later shown to be infected/infectious, while others remain not infected/infectious. In addition, the following observations have been made following experimental

<footnote>*Correspondence: racem.ben-romdhane@oniris-nantes.fr</footnote>

<footnote n="1">¹ BIOEPAR, INRA, ONIRIS, 44307 Nantes, France</footnote>

<footnote>Full list of author information is available at the end of the article</footnote>
## Equations describing flows

New incoming (and outgoing) flows in each health states are mainly drawn using binomial laws.

### Births (bX)

At each time step t, births are calculated with regards to the health state of the dam. These births are then distributed into S and T states

$$bS_{(t)} = bS_{(t)}^{R} + bS_{(t)}^{L} + bS_{(t)}^{Is} + bS_{(t)}^{Ic} \quad \text{et} \quad bT_{(t)} = bT_{(t)}^{L} + bT_{(t)}^{Is} + bT_{(t)}^{Ic},$$

where $bS_{(t)}^X$ et $bT_{(t)}^X$ ($X \in R, L, Is, Ic$) represent the number of births at time t from cows in health state X:

$$bS_{(t)}^{R} \sim \text{Bin} \left( R_{(t,cal)} + \sum_{i=1}^{i=5} [R_{(t,P_i)}^{ng}] ; 1 - \sigma_B \right)$$

$$bS_{(t)}^{L} = nbV_{(t)}^{L} - bT_{(t)}^{L}, \quad bT_{(t)}^{L} \sim \text{Bin} \left( nbV_{(t)}^{L} ; p_L * \varphi_{p_{LI}s} \right)$$

$$bS_{(t)}^{Is} = nbV_{(t)}^{Is} - bT_{(t)}^{Is}, \quad bT_{(t)}^{Is} \sim \text{Bin} \left( nbV_{(t)}^{Is} ; p_{Is} * \varphi_{p_{LI}s} \right)$$

$$bS_{(t)}^{Ic} = nbV_{(t)}^{Ic} - bT_{(t)}^{Ic}, \quad bT_{(t)}^{Ic} \sim \text{Bin} \left( nbV_{(t)}^{Ic} ; p_{Ic} * \varphi_{p_{Ic}} \right)$$

In equations above, $nbC_{(t)}^X$ is the number of female calves alive born at time t. It is obtained from $nbV_{(t)}^X$, the number of female calves born at time t, from cows in the health state X:

$$nbV_{(t)}^{L} \sim \text{Bin} \left( L_{(t,cal)} + T_{(t,cal)} + \sum_{i=1}^{i=5} [L_{(t,P_i)}^{ng}] ; 1 - \sigma_B \right)$$

$$nbV_{(t)}^{Is} \sim \text{Bin} \left( Is_{(t,cal)} + \sum_{i=1}^{i=5} [Is_{(t,P_i)}^{ng}] ; 1 - \sigma_B \right)$$

$$nbV_{(t)}^{Ic} \sim \text{Bin} \left( Ic_{(t,cal)} + \sum_{i=1}^{i=5} [Ic_{(t,P_i)}^{ng}] ; 1 - \sigma_B \right)$$

We note that at the age of moving in the adult group (cal), all the animals in the health state T enter the health state L ($L_{(t,cal)} = L_{(t,cal)} + T_{(t,cal)}$). $\varphi_{px}$ is a factor, varying from 100% to 0%, of decrease in probability of in utero transmission where x correspond to the dam giving birth in health states L and Is ($\varphi_{p_{LI}s}$) or Ic ($\varphi_{p_{Ic}}$).

### Change in age group ($X^x$)

$$X_{(t,P_i)} = X_{(t,P_i)}^{sg} + X_{(t,P_i)}^{ng}$$

$$X_{(t,P_i)}^{ng} \sim \text{Bin} \left( X_{(t,P_i)} ; \frac{1}{\tau_{aa}} \right) ; \text{ where } X = \{R, L, Is, Ic\},$$

and $\tau_{aa}$ is the average time spent in each of age group $P_1$ to $P_4$.
infection of similar aged calves with similar infectious dose of Map: (1) a wide range of paratuberculosis lesion severity have subsequently been observed [3], (2) different quantities of Map are shed in their faeces [6], and (3) different antibody responses have been detected, suggesting a variable duration of the latency period (being the period between infection and later detection by direct or indirect tests) [7]. The duration of the incubation period (which is defined as the period between infection and clinical signs) varied greatly between animals, ranging from 4 months to 15 years [8–12]. The amount of bacteria shed by infectious cattle is also highly variable, some being high shedders, while others are low shedders. Both intermittent and continuous shedding has been observed.

Individual resistance to paratuberculosis is assumed to be highly variable among, and expresses as different courses of infection. The phenotype of cattle resistance to paratuberculosis can be divided into (1) the ability to prevent infection and (2) the ability to cope with infection. This resistance in response to Map exposure involves different mechanisms and individual characters. Each of these characters will be denoted thereafter as phenotypic traits, a phenotype being defined by combined phenotypic traits. At the population scale, the distribution of phenotypic traits among individuals will influence the level of herd immunity, and therefore impair Map spread.

Strategies to control Map spread within dairy cattle herds usually consist in two main actions: hygiene improvement to reduce environmental and food contamination by Map, and a test-and-cull strategy to identify and remove infected animals. These control measures are not sufficient to control Map spread at herd and regional scales [13–15]. Vaccines against paratuberculosis have also been developed. Available vaccines decrease shedding of Map by infectious animals and decrease clinical signs of the disease [14, 16]. However, they do not prevent the infection of susceptible animals. In addition, most licensed vaccines show a cross reaction with tuberculosis diagnostic tests [17]. Therefore, the use of vaccination is restricted in many countries.

The observed variability of the individual response to Map exposure could support the development of innovative control measures applied at population scale if the most resistant animals can be selected. Several studies demonstrated a heritability of resistance to paratuberculosis in cattle ranging from 0.01 to 0.23 [17–21]. Recent studies highlighted an association between genetic markers and the course of Map infection [22–26]. Other genome markers were associated with Map shedding in faeces, presence of Map in several tissues, and seropositivity, in animals from comparable herds regarding paratuberculosis infection and of the same age group. Therefore, these animals were assumed to have been

exposed in a similar way [25, 27, 28]. This highlights the potential to select for cattle more resistant to paratuberculosis. However, there are still gaps of knowledge concerning the phenotypic traits of resistance that would be the most relevant to improvements in the control of Map spread at population scale.

Modelling is the most appropriate approach to investigate the dynamics of complex systems such as within-herd Map transmission. Observational and experimental studies are both difficult to implement and expensive regarding the long evolution of paratuberculosis. In addition, a modelling approach allows us to overpass the lack of knowledge on genetic resistance of cattle to paratuberculosis by assuming improved phenotypic traits as if they were already selected for. Simulations then provide information on how such modifications of phenotypic traits would influence Map spread. Only one recent study investigated the potential effectiveness of hypothetical genetic selection as a strategy to control paratuberculosis at herd scale [29]. The authors assessed the effect of varying three phenotypic traits of resistance: (1) length of the susceptibility period, (2) level of the susceptibility to infection (expressed as the dose of Map required resulting in infection), and (3) duration of the latency period. Each tested phenotypic trait has been tested one-at-a-time and ranked by the time required to reach eradication. Modelling predictions showed that, when only genetic selection is implemented, eradication takes hundreds of years. However, this study did not investigate the potential progress in disease control when combining variations in several traits. In addition, other traits also could influence Map spread including intensity of shedding by infectious animals, in utero transmission, and progress of the infection course through different infection stages.

Our objective was to identify which phenotypic traits of resistance to paratuberculosis have the strongest influence on Map spread within a dairy cattle herd. The purpose was to identify ranges of phenotypic trait variations and trait combinations that limit Map spread in the herd. We assessed three categories of phenotypic traits both one-at-a-time and in combination, including: infection susceptibility, delays in the infection course, and shedding levels.

## Materials and methods

### Overall study design and model choice

A modelling approach was used to predict the effect of varying phenotypic traits of resistance to paratuberculosis on Map spread in a dairy cattle herd. We compared a situation where phenotypic traits were set at current observed levels with situations reached after a successful hypothetical genetic selection of more resistant animals
## Exits (sX)

The mortality of calves during the first week of life is applied at birth and defined above in the section concerning births.

From age 1 to cal, mortality and culling rates $\sigma_x$ are defined as:

* $a \in \{1; 2\} \rightarrow \sigma_x = \sigma_{c1}$

* $a \in [3; w] \rightarrow \sigma_x = \sigma_{c2}$

* $a \in [w + 1; cal] \rightarrow \sigma_x = \sigma_{c3}$

Then, exits following death and culling write as:

for $a = 1$:

$$sX_{(t,a)} \sim \text{Bin}(X_{(t,a)}; \sigma_x), \quad \text{where } X = \{S, T\},$$

for $a \in [2; 4]$:

$$sX_{(t,a)} = \text{sale} + \text{death}, \quad \text{where } X = \{S, T, L\}$$

$$\text{sale} \sim \text{Bin}(X_{(t,a)}; \sigma_m)$$

$$\text{death} \sim \text{Bin}(X_{(t,a)} - \text{sale}; \sigma_x)$$

for $a \in [5; cal - 11]$:

$$sX_{(t,a)} \sim \text{Bin}(X_{(t,a)}; \sigma_x) \quad ; \quad X = \{S, R, T, L\}$$

for $a = cal - 10$, we consider management by heifers (safe management, keep all female calves):

$$Cows_{(t)} = \sum_{i=1}^{5} R_{(t,P_i)} + \sum_{i=1}^{5} L_{(t,P_i)} + \sum_{i=1}^{5} Is_{(t,P_i)} + \sum_{i=1}^{5} Ic_{(t,P_i)}$$

$$Heifers_{(t)} = \sum_{a=1}^{u} S_{(t,a)} + \sum_{a=u+1}^{cal} R_{(t,a)} + \sum_{a=1}^{cal} T_{(t,a)} + \sum_{a=1}^{cal} L_{(t,a)} + \sum_{a=1}^{cal} Is_{(t,a)} + \sum_{a=1}^{cal} Ic_{(t,a)}$$

If the number of heifers, Heifers(t), is greater than $K_g$ or that the number of cows, Cows(t), is greater than $K_v$, we consider the sale of heifers:

$$sX_{(t,a)} = \text{sale} + \text{death} \quad ; \quad X = \{R, T, L, Is, Ic\}$$

$$\text{sale} \sim \text{Bin}(X_{(t,a)}; \exp(-\sigma_h \cdot (Cows_{(t)}/K_v)^6) \cdot ((Heifers_{(t)}/K_g)^6))$$

$$\text{death} \sim \text{Bin}(X_{(t,a)} - \text{sale}; \sigma_x)$$

where $K_v$ is the capacity of the holding in number of cows (see Table 1) and $K_g = \sigma_P * K_v * (cal - h)$ is the capacity of the holding in number of heifers.

Otherwise, we do not consider the sale of heifers:

$$sX_{(t,a)} \sim \text{Bin}(X_{(t,a)}; \sigma_x)$$
in response to Map exposure. For each change of a trait, the resistance level was simulated as constant over time assuming that this average level had been reached in the population after a (not modelled) selection period. Several scenarios were simulated where one or several phenotypic traits were varied. The scenarios were compared regarding Map spread in the herd.

Several models have been published that represent Map spread within a dairy cattle herd (reviewed in [30], and more recently [29, 31–38]). We selected a stochastic compartmental model that offers an up-to-date description of Map spread within a dairy cattle herd. This model takes into account all of the major processes involved (according to the most recent literature) and allowed us to represent phenotypic traits of resistance corresponding to all of our hypotheses of interest. This model adequately combines demographical and infection dynamics, and accounts for herd structure, all these processes having been shown to highly influence Map spread [32]. The chosen model is mechanistic: each step and mechanism of the infection course is represented by a model parameter. This allowed us to simulate changes in phenotypic traits of resistance by minimal changes in the model.

## Main features of the model

The within-herd transmission model and the corresponding equations are fully described in Marcé et al. [32] and Beaunée et al. [15]. And a detailed description of the model is presented in Additional file 1.

The main modelling assumptions are the following: the herd population dynamics reflects the one of a typical western Europe Holstein herd with 5 age groups (unweaned calves, weaned calves, young heifers, bred heifers, and cows), a high renewal rate of cows (one-third per year), and no males kept in the herd. The within-herd contact structure varies seasonally between housing and pasturing periods. The infection dynamics is represented by successive health states (Figure 1): animals initially susceptible (state S) are assumed to be no longer susceptible (state R) after a susceptibility period of duration *u*. Susceptibility decreases with age, assuming an exponential decay coefficient *h*. The possible infection of adults is neglected as it rarely occurs (it has only been demonstrated in adults following sudden exposure to a highly contaminated environment [39]). Infection can occur when a susceptible animal is in contact with a sufficient infectious dose per animal *α* (explicit indirect transmission), and then becomes transiently infectious (state T) for an average duration *v*<sub>T</sub>. Then, infected animals enter a latent state (state L), during which shedding is neglected. After this latent period of average duration *v*<sub>L</sub>, they become moderate shedders (state Is). For some animals, the evolution of the infection leads to a persistently

high shedding and most of the animals are likely to have reduced milk yield or clinical signs called here high shedding and clinically affected state (Ic) after an average duration *v*<sub>Is</sub> in the moderate shedding state. Animals are assumed to be culled on average 6 months after entering Ic state.

Susceptible animals can be infected through five transmission routes [32]: (1) contact with bacteria present in the general environment of the farm contaminated by all shedders, (2) contact with bacteria present in the local environment of calves contaminated by shedding calves, (3) in utero transmission from infected cows to their foetus, and (4) ingestion of contaminated milk or (5) colostrum from infectious cows.

## Phenotypic traits of resistance to paratuberculosis

In this study we assessed the effect of varying 14 phenotypic traits of resistances to paratuberculosis and combination of them on Map transmission in the herd. Each of the tested scenarios corresponds to a variation, or combinations in variation, in one or more phenotypic trait of resistance to paratuberculosis.

The phenotype of cattle resistance to paratuberculosis is classically divided into (1) resistance to infection defined as the ability to prevent infection when exposed to a given dose of Map, and (2) tolerance to infection defined as the ability to cope with infection when infected [40, 41]. On the one hand, animals are considered to be resistant (ability to prevent infection) if they show a decrease in susceptibility to infection if they are no longer susceptible at a younger age, if they need to be exposed to a higher dose of Map to be infected, or if

<figure><img src="image_3.png" /><figcaption>**Figure 1** Paratuberculosis infection course and phenotypic traits of interest to reduce Map spread at herd scale: in bracket: the corresponding parameters coding for them in the model. Boxes, disease states; S, susceptible; R, no longer susceptible; T, transiently infectious; L, latent; Is, moderate shedder; Ic, high shedder or clinically affected animals. Green compartments: non infected states, orange compartments: infected states, dashed arrows: shedding, solid arrows: transitions between states, blue (large) arrows: changes in individual phenotypic traits that could limit Map spread at herd scale.</figcaption></figure>
for $a \in [cal - 9; cal]$:

$$sX_{(t,a)} = Bin(X_{(t,a)}; \sigma_x)$$

for $a \in [P_1; P_5]$:

$$sX_{(t,P_i)}^x \sim Bin(X_{(t,P_i)}; \sigma_{P_i}), \quad \text{where } X = \{R, L, Is\},$$

$$sIc_{(t,P_i)}^x \sim Bin\left(Ic_{(t,P_i)}^x; \frac{1}{v_{Ic}}\right)$$

where $v_{Ic}$ corresponds to the average time spent in the health state $Ic$.

## New infections ($S \to T$, except for in-utero transmission)

We have:

$$NT_{(t+1,a)} = inf_{(t+1,a)}^c + inf_{(t+1,a)}^m + inf_{(t+1,a)}^l + inf_{(t+1,a)}^g$$

Superscripts correspond to different possible routes of transmission, respectively colostrum (c), milk (m), local environment (l) and global environment (g). New infections by in-utero transmission are accounted for through births.

By age, the possible routes of infection are:

<figure><img src="image_1.png" /></figure>

## Transmission through colostrum

It is consideblackcol that calves drink colostrum from their mothers during the first three days before drinking milk:

$$inf_{(t+1,a=1)}^c = \sum_{1}^{bS_{(t)}^{Is}} \left[ Bern \left( 1 - exp \left( -\frac{\beta_l q_c^{Is}}{\alpha} \right) \right) \right] + \sum_{1}^{bS_{(t)}^{Ic}} \left[ Bern \left( 1 - exp \left( -\frac{\beta_l q_c^{Ic}}{\alpha} \right) \right) \right],$$

with

$$q_c^{Is} \sim Bern(sh_{Is}) \times [3 \times b (10^5 \cdot Beta(8; 8) + 1 + 10^3 \cdot Beta(1; 25))] * \varphi Milk_{Is},$$

$$q_c^{Ic} \sim Bern(sh_{Ic}) \times [3 \times b (10^5 \cdot Beta(8; 8) + 10^{(3+10 \cdot Beta(50;200))})] * \varphi Milk_{Ic}.$$

Where $\varphi Milk_x$ is a factor, varying from 100% to 0%, of decrease in quantity of bacteria shed in colostrum and milk by an animal in health state $x$ (here $x$ could be $Is$ or $Ic$)
they show a faster decrease in susceptibility with age than less resistant animals. On the other hand, animals are considered tolerant (ability to cope with infection once infected) if they show longer latency and incubation periods, and a lower shedding level when infectious than less tolerant animals. In addition, foetuses of the latter may have a lower chance to be infected in utero.

We accounted for all of the potential mechanisms involved in an increased resistance to paratuberculosis because we assumed they can all contribute to Map spread at herd scale (Figure 1). The ability to remain non-infected was composed of four components: (1) a shorter susceptibility period for calves, (2) a faster decrease in age-related susceptibility, and (3) a higher infectious dose of Map needed to be infected after birth. The ability to cope with infection was represented by a longer latency period before the onset of moderate shedding, a longer incubation period before high shedding and clinical signs, a decrease in the amount of Map shed through the different transmission routes and a decrease in the probability of in utero transmission. Overall, we studied 14 parameters coding for the identified phenotypic traits of resistance to paratuberculosis (Table 1).

Based on the literature, we defined a realistic variation of resistance levels to simulate within observed values for the investigated traits. The reference value was the worst one. Changes were simulated from reference to the most favourable value observed value, Indeed, calves susceptibility can sharply reduce, and animals are no longer susceptible, as soon as their first week of life [53–56]. Some susceptible animals have been shown to need a dose of bacteria as high as $10^{12}$ to become infected [57]. After a transient shedding period, infected animals can have a barely detectable level of shedding for about 4 years (208 weeks) [10, 11]. Infected animals can show clinical signs of the disease up to more than 9 years after infection (468 weeks). Concerning the probability of transmission of Map in utero from infected dam to its foetus and the quantities of bacteria shed through different routes, only partial information was available. Hence, we chose to test for extreme values by assuming that animals can stop shedding completely with no further in utero transmission of the infection. Nevertheless, it has been shown that high shedders and clinically affected animals can shed as few as $10^8$ bacteria/kg of faeces, which corresponds to $1/100^{th}$ of the reference value that we have assumed in our model [58, 59].

## Initial conditions and model outputs

Map spread was initiated by the introduction of a moderate shedding cow into a fully naïve herd of 260 animals. We assumed that herd renewal is mainly driven by internal demographic processes (no further introduction),

which is typical of western Europe farming systems. Map spread was predicted over 25 years. To obtain accurate outputs from the stochastic model, we ran 500 repetitions for each of the tested scenario. A scenario represented one phenotype of interest. Each phenotype was defined by a set of values of 14 parameters.

Four model outputs described Map spread within a herd (Table 2). All outputs were calculated at the end of the simulation, $t = 25$ years after Map introduction. The first output was the cumulative incidence calculated as the mean cumulative number of newly infected animals over the 25 years of simulation. The second output was the infection persistence defined as the proportion of runs where the infection persisted until 25 years after Map introduction, i.e. where there was at least one infected animal of state *T*, *L*, *Is*, or *Ic*, or bacteria in the environment. The third output was the prevalence of infected animals calculated as the median prevalence of infected animals in the population 25 years after Map introduction for runs where the infection persisted. Finally, the fourth output was the prevalence of affected animals calculated as the median prevalence of high shedding and clinically affected animals in the population 25 years after Map introduction for runs where the infection persisted. Outputs related to prevalence were calculated only if Map persistence was higher than 6% (30 runs out of 500) in order to provide a sufficient number of runs to estimate medians.

## Simulation protocol and output analysis

First, we performed a univariate simulation study: each of the traits of interest was varied one-at-a-time, assuming they varied independently (Table 1). Second, we performed a multivariate simulation study: combinations of phenotypic traits were studied to test for a potential enhanced effect of simultaneously improving several phenotypic traits simultaneously. The R programming language [42] was used for data analyses. Results obtained in the univariate simulation study revealed that some parameters—when analysed one-at-a-time—did not influence model outputs. Instead of keeping numerous parameters or removing some of them expected not to be influential, we grouped in this second step non-influential parameters when they are untangled in the same trait or when involved in a given transmission route. This decrease in the number of considered parameters without losing information eased the interpretation of the multivariate simulation study results. Ranges of variation of phenotypic traits were represented by five possible values per trait (including the reference value) combined in the multivariate simulation study using a complete factorial design, leading to 390 625 scenarios (Table 1). Five levels of variation per trait appeared to be a good
## Transmission through milk

Regarding the age, we have:

$$a = 1: \quad inf_{(t+1,a)}^m \sim Bin \left( bS; 1 - \exp \left( -\frac{\beta_l q_l \frac{4}{7}}{\alpha} \right) \right),$$

$$a \in \{2; w\}: \quad inf_{(t+1,a)}^m \sim Bin \left( [S_{(t,a-1)} - sS_{(t,a-1)}]; \left[1 - \exp\left(-e^{(-\gamma(a-1)) \cdot \frac{\beta_l q_l}{\alpha}}\right)\right] \right),$$

with

$$q_l = \frac{7 \times d \times (Q_{Is(t)}^{milk} + Q_{Ic(t)}^{milk})}{MilkTot_{(t)}},$$

where

$$Q_{Is(t)}^{milk} = 7 \times \varepsilon \times g_{Is} \times \left( \sum_{1}^{nbExcr_{(t)}^{Is}} (10^5 \cdot \text{Beta}(8;8)) + \sum_{1}^{nbLac_{(t)}^{Is}} (1 + 10^3 \cdot \text{Beta}(1;25)) \right) * \varphi \text{Milk}_{Is},$$

$$Q_{Ic(t)}^{milk} = 7 \times \varepsilon \times g_{Ic} \times \left( \sum_{1}^{nbExcr_{(t)}^{Ic}} (10^5 \cdot \text{Beta}(8;8)) + \sum_{1}^{nbLac_{(t)}^{Ic}} 10^{(3+10 \cdot \text{Beta}(50;200))} \right) * \varphi \text{Milk}_{Ic},$$

$$MilkTot_{(t)} = 7 \times \varepsilon \times (nbLac_{(t)}^{R} + g_L \cdot nbLac_{(t)}^{L} + g_{Is} \cdot nbLac_{(t)}^{Is} + g_{Ic} \cdot nbLac_{(t)}^{Ic}).$$

with

$$nbExcr_{(t)}^{Is} \sim Bin \left( nbLac_{(t)}^{Is}; sh_{Is} \right),$$

$$nbExcr_{(t)}^{Ic} \sim Bin \left( nbLac_{(t)}^{Ic}; sh_{Ic} \right),$$

$$nbLac_{(t)}^{R} \sim Bin \left( \sum_{i=1}^{i=5} R_{(t,P_i)}, prop \right),$$

$$nbLac_{(t)}^{L} \sim Bin \left( \sum_{i=1}^{i=5} L_{(t,P_i)}, prop \right),$$

$$nbLac_{(t)}^{Is} \sim Bin \left( \sum_{i=1}^{i=5} Is_{(t,P_i)}, prop \right),$$

$$nbLac_{(t)}^{Ic} \sim Bin \left( \sum_{i=1}^{i=5} Ic_{(t,P_i)}, prop \right).$$

## Local transmission (in collective pens, $a \in [1; u]$)

$$inf_{(t+1,a)}^l \sim Bin \left( S_{(t,a-1)} - sS_{(t,a-1)}; p_{inf(t+1,a)}^i \right),$$
<table><caption>Table 1 Parameters coding for the phenotypic traits of resistance: definition and values</caption><thead><tr><th rowspan="2">Parameters</th><th rowspan="2">Definition</th><th rowspan="2">Reference value</th><th rowspan="2">Univariate simulations: [min-max]</th><th colspan="4">Multivariate simulations: tested values</th><th rowspan="2">Source</th></tr><tr><th>#1</th><th>#2</th><th>#3</th><th>#4</th></tr></thead><tbody><tr><td>u</td><td>Susceptibility period duration</td><td>52 weeks</td><td>[1–52]</td><td>–</td><td></td><td></td><td></td><td>[51–53]</td></tr><tr><td>h</td><td>Decay in susceptibility with age (coefficient)</td><td>0.1</td><td>[0.1–1]</td><td>0.2</td><td>0.3</td><td>0.4</td><td>0.5</td><td>[54]</td></tr><tr><td>a</td><td>Required infectious dose to be infected</td><td>10<sup>6</sup> bacteria</td><td>[10<sup>6</sup>–10<sup>12</sup>]</td><td>1.5 × 10<sup>6</sup></td><td>2 × 10<sup>6</sup></td><td>2.5 × 10<sup>6</sup></td><td>3 × 10<sup>6</sup></td><td>[55]</td></tr><tr><td>v<sub>L</sub></td><td>Duration of latent state</td><td>52 weeks</td><td>[52–208]</td><td>–</td><td></td><td></td><td></td><td>[10, 11]</td></tr><tr><td>v<sub>L</sub> + v<sub>Is</sub></td><td>Duration before high shedding and clinically affected state</td><td>156 weeks</td><td>[156–468]</td><td>234</td><td>312</td><td>390</td><td>468</td><td>[10–12]</td></tr><tr><td>v<sub>T</sub> with v<sub>T</sub> + v<sub>L</sub> = constant</td><td>Duration of transiently infectious state with constant duration before moderate shedding state</td><td>vT = 25 weeks (vT + vL = 77)</td><td>[1–25]</td><td>–</td><td></td><td></td><td></td><td>[8–11]</td></tr><tr><td>v<sub>Is</sub> with v<sub>L</sub> + v<sub>Is</sub> = constant</td><td>Duration of moderate shedding state with constant duration before high shedding or clinically affected state</td><td>vIs = 104 weeks (vL + vIs = 156 weeks)</td><td>[60–104]</td><td>95</td><td>86</td><td>77</td><td>68</td><td>[8–12]</td></tr><tr><td>φMilk<sub>X</sub></td><td>Factor of decrease of Map shed in milk by animals in health state X</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>φMilk<sub>Is</sub></td><td>Moderate shedding state (Is)</td><td>100%</td><td>[0–100]</td><td>50%</td><td>10%</td><td>5%</td><td>0%</td><td>[56]</td></tr><tr><td>φMilk<sub>Is</sub></td><td>High shedding and clinically affected state (Is)</td><td>100%</td><td>[0–100]</td><td></td><td></td><td></td><td></td><td>[57]</td></tr><tr><td>φFeces<sub>X</sub></td><td>Factor of decrease of Map shed in faeces by animals in health state X</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>φFeces<sub>T</sub></td><td>Transient state (T)</td><td>100%</td><td>[0–100]</td><td>50%</td><td>10%</td><td>5%</td><td>0%</td><td>[9]</td></tr><tr><td>φFeces<sub>Is</sub></td><td>Moderate shedding state (Is)</td><td>100%</td><td>[0–100]</td><td></td><td></td><td></td><td></td><td>[58]</td></tr><tr><td>φFeces<sub>Is</sub></td><td>High shedding or clinically affected state (Is)</td><td>100%</td><td>[0–100]</td><td>66%</td><td>50%</td><td>40%</td><td>33%</td><td>[59, 60]</td></tr><tr><td>φP<sub>X</sub></td><td>Factor of decrease of probability of in utero transmission for cows in health state X</td><td></td><td></td><td></td><td></td><td></td><td></td><td>[2, 61]</td></tr><tr><td>φP<sub>LIs</sub></td><td>Latent and moderate shedding states (LIs)</td><td>100%</td><td>[0–100]</td><td>50%</td><td>10%</td><td>5%</td><td>0%</td><td></td></tr><tr><td>φP<sub>Is</sub></td><td>High shedding or clinically affected state (Is)</td><td>100%</td><td>[0–100]</td><td></td><td></td><td></td><td></td><td></td></tr></tbody></table>

compromise between parameter space exploration and number of scenarios to investigate interactions. A complete factorial design was required to assess all interaction orders.

We performed a cluster analysis of the multivariate scenarios based on two of our model outputs. Scenarios were grouped to minimise outputs variability within a cluster and maximise this variability among clusters. The aim was to identify and characterize groups of scenarios. We build clusters using the two model outputs

available on all model repetitions of each scenario: cumulative incidence and infection persistence 25 years after Map introduction, after they were standardized into variables of comparable scales. To define the appropriate number of clusters, we studied the sum of squared distances between each scenario and the centroid of its corresponding cluster (called the sum of squared error or the within-group sum-of-squares) for different number of clusters [43–45]. Clusters were built using k-means clustering method (kmeans function from R package
where

during housing period:

$$p_{inf(t+1,a)}^i = 1 - \exp \left( -e^{(-\gamma(a-1))} \cdot \frac{\beta_c E_{(t+1)}^i}{\alpha N_{(t+1)}^i} \right), \text{ and}$$

during grazing period, regarding the age,

$$a \in [1; 26] : \quad p_{inf(t+1,a)}^i = 1 - \exp \left( -e^{(-\gamma(a-1))} \cdot \frac{\beta_c E_{(t+1)}^i}{\alpha N_{(t+1)}^i} \right),$$

$$a \in [27; u] : \quad p_{inf(t+1,a)}^i = 1 - \exp \left( -e^{(-\gamma(a-1))} \cdot \frac{\beta_o E_{(t+1)}^i}{\alpha N_{(t+1)}^i} \right).$$

$N_{(t+1)}^i$ is the total number of animals in environment $i$, across all health states, and $E^i$ represents the quantity of bacteria in the environment, with $i$ corresponding to the specific area.

## Global transmission

This occurs up to the age allowing to go to the pasture (26 weeks - 6 months) during the grazing period, and up to the age limit for sensitivity ($u$) during the housing period:

$$inf_{(t+1,a)}^g \sim Bin \left( S_{(t,a)} - sS_{(t,a)}; p_{inf(t+1,a)}^g \right),$$

where

$$p_{inf(t+1,a)}^g = 1 - \exp \left( -\exp[-h(a-1)] \cdot \frac{\beta_g E_{(t+1)}^g}{\alpha N_{(t+1)}^g} \right), \quad \text{with} \quad E_{(t+1)}^g = \sum_{i=1}^{i=5} E_{(t+1)}^{INTi}$$

$E^{INTi}$ represents the quantity of bacteria in the environment where $INTi$ corresponds to a specific area (see the section below about the dynamics of bacteria in the environments).

## New latently infected (T → L)

For $a < cal$:

$$\mathrm{NL}_{(t+1,a \in [2;cal-1])} \sim \mathrm{Bin} \left( T_{(t,a)} - sT_{(t,a)}; \frac{1}{v_T} \right).$$

After age $cal$, there are no more animals in $T$ state:

$$\mathrm{NL}_{(t+1,P_1)} = T_{(t,cal)} - sT_{(t,cal)}.$$
<table><caption>Table 2 Contribution of the four most influential phenotypic traits to the model output variance</caption><thead><tr><th></th><th>Parameters</th><th>Cumulative incidence</th><th>Persistence</th><th>Prevalence of infected animals</th><th>Prevalence of affected animals</th></tr></thead><tbody><tr><td rowspan="4">Principal effect</td><td>h</td><td>0.25</td><td>0.34</td><td>0.30</td><td>0.22</td></tr><tr><td>vL + vls</td><td>0.16</td><td>0.16</td><td>0.18</td><td>0.22</td></tr><tr><td>φFeceslc</td><td>0.14</td><td>0.17</td><td>0.14</td><td>0.09</td></tr><tr><td>α</td><td>0.13</td><td>0.17</td><td>0.12</td><td>0.08</td></tr><tr><td rowspan="6">First order</td><td>h:vL + vls</td><td>0.07</td><td>0.05</td><td>0.09</td><td>0.12</td></tr><tr><td>h:φFeceslc</td><td>0.06</td><td>0.04</td><td>0.06</td><td>0.03</td></tr><tr><td>h:α</td><td>0.06</td><td>0.04</td><td>0.05</td><td>0.03</td></tr><tr><td>vL + vls:φFeceslc</td><td>0.03</td><td>0.02</td><td>0.03</td><td>0.04</td></tr><tr><td>α:vL + vls</td><td>0.03</td><td>0.02</td><td>0.03</td><td>0.03</td></tr><tr><td>h:vL + vls</td><td>0.03</td><td>2 × 10<sup>-3</sup></td><td>0.02</td><td>0.03</td></tr><tr><td rowspan="5">Second order</td><td>h:vL + vls</td><td>0.03</td><td>10<sup>-3</sup></td><td>0.02</td><td>0.03</td></tr><tr><td>h:vL + vls:φFeceslc</td><td>0.03</td><td>0.02</td><td>0.01</td><td>6 × 10<sup>-3</sup></td></tr><tr><td>α:φFeceslc</td><td>0.02</td><td>2 × 10<sup>-4</sup></td><td>5 × 10<sup>-3</sup></td><td>10<sup>-3</sup></td></tr><tr><td>h:α:φFeceslc</td><td>0.02</td><td>3 × 10<sup>-4</sup></td><td>2 × 10<sup>-3</sup></td><td>4 × 10<sup>-3</sup></td></tr><tr><td>α:vL + vls:φFeceslc</td><td>8 × 10<sup>-3</sup></td><td>3 × 10<sup>-3</sup></td><td>4 × 10<sup>-4</sup></td><td>7 × 10<sup>-5</sup></td></tr><tr><td>Third order</td><td>h:vL + vls:φFeceslc</td><td>3 × 10<sup>-3</sup></td><td>3 × 10<sup>-3</sup></td><td></td><td></td></tr><tr><td></td><td>Residuals</td><td>0.31</td><td>0.15</td><td>0.27</td><td>0.38</td></tr></tbody></table>

Contribution was estimated from the ANOVA. In italic, values above 5%.

"FactoMineR" [46]). A descriptive analysis was performed to characterize clusters for phenotypic traits using *catdes* function ("FactoMineR" package [46]). This step aimed to identify if tested variations in phenotypic traits are uniformly distributed in regards of cluster or if some values are over represented in a given cluster. Besides, we performed an ANOVA to quantify the contribution of each trait to the variance of each of the four model outputs. Each trait contribution to the model output variance (κ) was calculated as:

$$ \kappa = \text{principal effect of the trait} $$

$$ + \sum_{i=1}^{m} \frac{\text{ith order interaction effect involving the trait}}{i+1}, $$

with *i* the interaction order and *m* the highest interaction order in which the trait was involved. A second ANOVA was performed on the influential phenotypic traits to quantify contribution of each trait (principal effect) and each combination of them (interaction) to the variance of each of the four model outputs. Factors (individual traits or combinations of traits) were influential if they contributed to more than 5% of the variance of at least one of the four model outputs.

In order to identify the most effective combinations of variation in phenotypic traits to decrease Map spread, we used the cumulative incidence output as an indicator of a successful Map control at herd scale. For each combination, the cumulative incidence was plotted and visually described (Figure 6). In addition, we chose two thresholds to evidence the most effective combinations with emphasis on the ones with the lowest variations in

parameters (Figures 6, 7): (1) 25 newly infected animals over the 25 years of simulation, interpreting such a level of one newly infected animal per year as an infection under control, and (2) half this threshold, i.e. 12 newly infected animals over the 25 years of simulation.

## Results

In the univariate simulation study the variation of six phenotypic traits influenced at least one model output resulting in decrease in Map spread (Figure 2): a shorter susceptibility period (u), an increase in the decay in susceptibility with age (h), an increase in the required infectious dose (α), a longer latent state (vL), a delayed occurrence of the high shedding or clinically affected state (vL + vIs), and a decrease in the quantity of Map shed in faeces by high shedders or clinically affected animals (φFecesIc).

Eight of the traits investigated in the univariate simulation study did not influence Map spread dynamics. These traits were: a shorter transiently infectious state when assuming a constant duration before the moderate shedding state (vT with vT + vL = constant), a shorter latent state when assuming a constant duration before the high shedding and clinically affected state (vIs with vL + vIs = constant), a decrease in the quantity of Map shed in milk by moderate shedders (φMilkIs) and high shedders and clinically affected animals (φMilkIc), a decrease in the quantity of Map shed by transiently infectious animals (φFecesT), and by moderate shedders (φFecesIs), and a decrease in the probability of in utero transmission by latent infectious animals and moderate
## New subclinically infected (L → Is)

For heifers:

$$NI s_{(t+1,a \in [h+1;cal])} \sim Bin \left( L_{(t,a)} - sL_{(t,a)}; \frac{1}{v_L} \right).$$

For cows:

$$NI s_{(t+1,P_i)} \sim Bin \left( n; \frac{1}{v_L} \right),$$

with, regarding the age,

$$P_1 \rightarrow n = [L_{(t,P_1)}^{sg} - sL_{(t,P_1)}^{sg} + L_{(t,cal)} - sL_{(t,cal)} + T_{(t,cal)} - sT_{(t,cal)}],$$

$$\{P_2; P_4\} \rightarrow n = [L_{(t,P_i)}^{sg} - sL_{(t,P_i)}^{sg} + L_{(t,P_{i-1})}^{ng} - sL_{(t,P_{i-1})}^{ng}],$$

$$P_{5+} \rightarrow n = [L_{(t,P_{5+})} - sL_{(t,P_{5+})} + L_{(t,P_4)}^{ng} - sL_{(t,P_4)}^{ng}].$$

## New clinically infected (Is → Ic)

For heifers:

$$NI c_{(t+1,a \in [h+1;cal])} \sim Bin \left( Is_{(t,a)} - sIs_{(t,a)}; \frac{1}{v_{Is}} \right).$$

For cows:

$$NI c_{(t+1,P_i)} \sim Bin \left( n; \frac{1}{v_{Is}} \right),$$

with

$$P_1 \rightarrow n = [Is_{(t,P_1)}^{sg} - sIs_{(t,P_1)}^{sg} + Is_{(t,cal)} - sIs_{(t,cal)}],$$

$$\{P_2; P_4\} \rightarrow n = [Is_{(t,P_i)}^{sg} - sIs_{(t,P_i)}^{sg} + Is_{(t,P_{i-1})}^{ng} - sIs_{(t,P_{i-1})}^{ng}],$$

$$P_{5+} \rightarrow n = [Is_{(t,P_{5+})} - sIs_{(t,P_{5+})} + Is_{(t,P_4)}^{ng} - sIs_{(t,P_4)}^{ng}].$$

## New resistant (S → R)

At age u, the transition from compartment S to compartment R is done in a deterministic way.

## Dynamics of bacteria in the environments (E)

The composition of the environments according to the season is the following:

$$\text{Grazing : } a = \frac{\text{INT1}}{1 \cdots \cdots w} \frac{\text{INT2}}{26} \frac{\text{EXT1} + \text{EXT2}}{\cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots h} \frac{\text{EXT3}}{\cdots \cdots \cdots cal}$$
<figure><img src="image_4.png" /><figcaption>**Figure 2** Changes in model outputs resulting from univariate variations of phenotypic traits of resistance to paratuberculosis. A cumulative incidence, B persistence, C prevalence of infected animals, D prevalence of high shedders and affected animals at the end of simulations (25 years). "Ref" corresponding to the reference value of the phenotypic trait. See Table 1 for parameter definitions and tested values, Table 2 for output definitions. Vertical and horizontal solid lines give reference values, dotted lines give associated 5th and 95th percentiles.</figcaption></figure>

shedders ($\phi P_{LIs}$), and by high shedders and clinically affected animals ($\phi P_{Ic}$).

We chose traits to be included in the multivariate simulation study in light of these results, noting that it was not possible to evaluate interactions among traits could have been evaluated with such a univariate analysis. Among traits highlighted as influential, we kept all except *u* that was redundant with *h*. Among other traits, we kept $v_{Is}$ (assuming $v_L + v_{Is} = \text{constant}$) and we grouped traits related to Map shedding in milk and colostrum ($\phi \text{Milk}$), to in utero transmission ($\phi P$), and to Map shedding in faeces by transiently infectious animals and moderate shedders ($\phi \text{Feces}_{TIs}$).

The cluster analysis of multivariate scenarios identified seven groups of scenarios from current (A; assumed as the worst) to the best control of Map spread (G; Figures 3A and B). This analysis highlighted three distinct dynamics (Figure 3): clusters A and B represented low control with a decrease in cumulative incidence, a slight

decrease in infection persistence, and an almost as high prevalence of infected animals. Clusters C, D, and E represented a good control with a low cumulative incidence, persistence and prevalence of infected animals, but with the occurrence of high shedders and clinically affected animals. Clusters F and G represented complete control with a very low cumulative incidence and persistence. Up to 80% of the scenarios were in these most favourable clusters F and G (Figure 3A).

The descriptive analysis (Figure 4) of clusters showed that the dynamics toward the most favourable clusters was mainly driven by four out of the eight traits: increasing the decay in susceptibility with age (*h*), lengthening the incubation period ($v_L + v_{Is}$), decreasing the quantity of Map shed in faeces by high shedders or clinically affected animals ($\phi \text{Feces}_{Ic}$), and increasing the required infectious dose ($\alpha$). The ANOVA (Figure 5) evidenced that these four phenotypic traits contributed most to the variance of model outputs, and allowed us to rank
<figure><img src="image_2.png" /></figure>

Dynamics of bacteria in the environments (E) are defined below:

$E_{(t+1)}^{\text{INT1}} = E_{(t)}^{\text{INT1}}(1 - \mu_g^{\text{int}}) + Q_{(t)}^{\text{TNS}}$. If pens are empty, it becomes $E_{(t+1)}^{\text{INT1}} = E_{(t+1)}^{\text{INT1}}(1 - \mu_{cp})$, where Q represents the quantity of bacteria shed.

During grazing period, we have:

$$ E_{(t+1)}^{\text{INT2}} = E_{(t)}^{\text{INT2}}(1 - \mu_g^{\text{int}}) + Q_{(t)}^{\text{Ts1}} \\ \text{If } \sum_{a=w+1}^{a=26} SR_{t,a} + \sum_{a=w+1}^{a=26} T_{t,a} + \sum_{a=w+1}^{a=26} L_{t,a} == 0 \text{ then } E_{(t+1)}^{\text{INT2}} = E_{(t+1)}^{\text{INT2}}(1 - \mu_{cp}) $$

$$ E_{(t+1)}^{\text{INT3}} = E_{(t)}^{\text{INT3}}(1 - \mu_g^{\text{int}}) $$

$$ E_{(t+1)}^{\text{INT4}} = E_{(t)}^{\text{INT4}}(1 - \mu_g^{\text{int}}) $$

$$ E_{(t+1)}^{\text{INT5}} = E_{(t)}^{\text{INT5}}(1 - \mu_g^{\text{int}}) $$

$$ E_{(t+1)}^{\text{EXT1}} = E_{(t)}^{\text{EXT1}}(1 - \mu_g^{\text{ext}}) + Q_{(t)}^{\text{Ts2}} $$

$$ E_{(t+1)}^{\text{EXT2}} = E_{(t)}^{\text{EXT2}}(1 - \mu_g^{\text{ext}}) + Q_{(t)}^{\text{Ty}} $$

$$ E_{(t+1)}^{\text{EXT3}} = E_{(t)}^{\text{EXT3}}(1 - \mu_g^{\text{ext}}) + Q_{(t)}^{\text{TH}} + Q_{(t)}^{\text{ISH}} + Q_{(t)}^{\text{ICH}} $$

During housing period, we have:

$$ E_{(t+1)}^{\text{INT2}} = E_{(t)}^{\text{INT2}}(1 - \mu_g^{\text{int}}) + Q_{(t)}^{\text{Ts1}} + Q_{(t)}^{\text{Ts2}} \\ \text{If } \sum_{a=w+1}^{a=52} SR_{t,a} + \sum_{a=w+1}^{a=52} T_{t,a} + \sum_{a=w+1}^{a=52} L_{t,a} == 0 \text{ then } E_{(t+1)}^{\text{INT2}} = E_{(t+1)}^{\text{INT2}}(1 - \mu_{cp}) $$

$$ E_{(t+1)}^{\text{INT3}} = E_{(t)}^{\text{INT3}}(1 - \mu_g^{\text{int}}) + Q_{(t)}^{\text{Ty}} $$

$$ E_{(t+1)}^{\text{INT4}} = E_{(t)}^{\text{INT4}}(1 - \mu_g^{\text{int}}) + Q_{(t)}^{\text{TH}} + Q_{(t)}^{\text{ISH}} + Q_{(t)}^{\text{ICH}} $$

$$ E_{(t+1)}^{\text{INT5}} = E_{(t)}^{\text{INT5}}(1 - \mu_g^{\text{int}}) + Q_{(t)}^{\text{IS}} + Q_{(t)}^{\text{IC}} $$

$$ E_{(t+1)}^{\text{EXT1}} = 0, \quad E_{(t+1)}^{\text{EXT2}} = 0, \quad E_{(t+1)}^{\text{EXT3}} = 0. $$

Shed quantities of bacteria are defined, regarding the health states and the age, by:

unweaned calves T:

$$ Q_{(t)}^{\text{TNS}} = \sum_{a=w}^{a=w} \left[ 7 \cdot f_1 \cdot 10^6 \cdot \sum_{(t,a)}^{T_{(t,a)}} \text{Beta}(8.8; 19) \right] * \varphi \text{faeces}_T, $$

weaned calves T, without access to grazing:

$$ Q_{(t)}^{\text{Ts1}} = \sum_{a=w+1}^{a=26} \left[ 7 \cdot f_2 \cdot 10^6 \cdot \sum_{(t,a)}^{T_{(t,a)}} \text{Beta}(8.8; 19) \right] * \varphi \text{faeces}_T, $$
<figure><img src="image_5.png" /><figcaption>**Figure 3** Model outputs 25 years after initial Map introduction for all of the multivariate scenarios: scenarios were clustered by cumulative incidence and persistence. A cumulative incidence and persistence per scenario and proportion of scenarios per cluster; B prevalence of infected animals and of affected animals per scenario where persistence was higher than 6%; C boxplots of the cumulative incidence for each cluster; D evolution over time of the prevalence of infected animals for the centroids of the seven clusters (A–G). Solid lines show output reference values, the dashed line represents the threshold of 30 runs where infection persists, asterisks indicate centroids of clusters. Total number of scenarios is 380 625.</figcaption></figure>

phenotypic traits from the most to the less influential. The increase in the decay in susceptibility with age (h) contributed the most to the variance of the four model outputs, while a decrease in the quantity of Map shed in faeces by high shedders and clinically affected animals ($\phi Feces_{Ic}$), a longer incubation period ($v_L + v_{Is}$), and an increase in the required infectious dose (α) led to almost equivalent contributions to model output variances.

Combined variations of phenotypic traits of resistance contributed to decrease Map spread dynamics. Interactions among traits showed contributions to model output variance ranging from 0.007% to up to 12% (Table 2). The interaction between increased decay in susceptibility with age (h) and a lengthened incubation period ($v_L + v_{Is}$) contributed to 12% of the variance of the prevalence of affected animals, and was also the most contributing interaction for other model outputs. In addition, h was involved in all of the contributing interactions, thus having both the highest principal and interaction effects.

The combined variation in the four most influential phenotypic traits of resistance to paratuberculosis

decreased the collative incidence to <1 newly infected animal over 25 years of simulation when set at their highest tested level. Over the 625 scenarios combining variations of the four most influential phenotypic traits, 537 scenarios resulted in decrease the cumulative incidence from 617 newly infected animals when phenotypic traits were set to their current values to 25 newly infected animals over the 25 years of simulation, and 473 scenarios allowed a cumulative incidence to be reached of 12 newly infected animals over the 25 years of simulation (Figure 6). Fourteen of the tested scenarios allowed to achieve a good control of the disease dynamics in the herd (<25 newly infected animals over the 25 years of simulation) with one of the four most influential traits at its reference value and the other trait at value 1 or value 2 (Figure 7). Some of the tested scenarios allowed to decrease the cumulative incidence to 25 newly infected animals over the 25 years of simulation or less were based on moderate variations of phenotypic traits (Figure 7). For example, a combined variation of the four most influential traits at their first tested level (#1) led to
weaned calves T, with access to grazing:

$$Q_{(t)}^{\text{Ts2}} = \sum_{a=27}^{a=y} \left[ 7 \cdot f_2 \cdot 10^6 \cdot \sum_{t=a}^{T_{(t,a)}} \text{Beta}(8.8; 19) \right] * \varphi \text{faeces}_{T},$$

young heifers T:

$$Q_{(t)}^{\text{TY}} = \sum_{a=y+1}^{a=h} \left[ 7 \cdot f_Y \cdot 10^6 \cdot \sum_{t=a}^{T_{(t,a)}} \text{Beta}(8.8; 19) \right] * \varphi \text{faeces}_{T},$$

heifers T:

$$Q_{(t)}^{\text{TH}} = \sum_{a=h+1}^{a=cal} \left[ 7 \cdot f_A \cdot 10^6 \cdot \sum_{t=a}^{T_{(t,a)}} \text{Beta}(8.8; 19) \right] * \varphi \text{faeces}_{T},$$

heifers Is:

$$\text{if } \sum_{a=h+1}^{a=cal} (Is_{(t,a)}) > 0 \quad : \quad Q_{(t)}^{\text{ISH}} = \sum_{a=h+1}^{a=cal} \left[ 7 \cdot f_A \cdot 10^{(4+10 \times \sum^{Is(t,a)} \text{Beta}(2.65;17))} \right] * \varphi \text{faeces}_{Is}$$

else $Q_{(t)}^{\text{ISH}} = 0$,

cows Is:

$$\text{if } \sum_{i=1}^{i=5} (Is_{(t,P_i)}) > 0 \quad : \quad Q_{(t)}^{\text{Is}} = \sum_{i=1}^{i=5} \left[ 7 \cdot f_A \cdot 10^{(4+10 \times \sum^{Is(t,P_i)} \text{Beta}(2.65;17))} \right] * \varphi \text{faeces}_{Is}$$

else $Q_{(t)}^{\text{Is}} = 0$,

heifers Ic:

$$\text{if } \sum_{a=h+1}^{a=cal} (Ic_{(t,a)}) > 0 \quad : \quad Q_{(t)}^{\text{ICH}} = \sum_{a=h+1}^{a=cal} \left[ 7 \cdot f_A \cdot 10^{(8+10 \times \sum^{Ic(t,a)} \text{Beta}(2;17))} \right] * \varphi \text{faeces}_{Ic}$$

else $Q_{(t)}^{\text{ICH}} = 0$,

cows Ic:

$$\text{if } \sum_{i=1}^{i=5} (Ic_{(t,P_i)}) > 0 \quad : \quad Q_{(t)}^{\text{IC}} = \sum_{i=1}^{i=5} \left[ 7 \cdot f_A \cdot 10^{(8+10 \times \sum^{Ic(t,P_i)} \text{Beta}(2;17))} \right] * \varphi \text{faeces}_{Ic}$$

else $Q_{(t)}^{\text{IC}} = 0$.

$\varphi \text{faeces}_x$, in the equation above is a factor of decrease in quantity of bacteria shed by an animal in health state $x$. $x$ could be $T$, $Is$ or $Ic$.
<figure><img src="image_6.png" /><figcaption>**Figure 4** Distribution of scenarios among tested values for each phenotypic trait per cluster (A–G). See Table 1 for parameter definitions and values, and Figure 3 for cluster definition.</figcaption></figure>

a cumulative incidence of 19 newly infected animals. On the other hand, improving a single trait, even a fivefold increase in the decay in susceptibility with age (#4), the most influential trait, was not sufficient to reach accumulative incidence of 25 newly infected animals or less over the 25 years of simulations. Interesting examples of combined moderate variations of traits allowing decreasing the cumulative incidence were highlighted. First, halving the decay in susceptibility with age (h) (#2) together with a 50% increase in duration before entering the high shedding and clinically affected state ($V_L + V_{Is}$) (#1), and a 34% decrease in the quantity of bacteria shed by high shedders or clinically affected animals ($\phi FecesIc$) (#1) results in threshold being reached of 23 newly infected animals over the 25 years of simulation. Second, tripled threefold increase in the decay in susceptibility with age (h) (#2) combined with a doubling of the required infectious dose ($\alpha$) (#2), and a 34% decrease in the quantity of bacteria shed by high shedders or clinically affected animals ($\phi FecesIc$) (#1) resulted in a cumulative incidence of 13 newly infected animals.

## Discussion

Variations of 4 of 14 phenotypic traits strongly reduced Map spread within a dairy cattle herd: the decay in susceptibility with age, this being the most influential trait, the quantity of Map shed in faeces by high shedders and clinically affected animals, the duration of the incubation period, and the required infectious dose. Combining these phenotypic traits was the sole way to effectively control Map spread at the herd scale. Most tested combinations of these influential phenotypic traits allowed the cumulative incidence to be reduced to <25 newly infected animals over the 25 years of simulation, which was interpreted here as an infection under control. Interestingly, such a low level of cumulative incidence could not be reached when varying a single phenotypic trait.

The increase in the decay in susceptibility with age is largely related to a shorter susceptibility period. We also highlighted the required infectious dose as an influential phenotypic trait. Our results concerning these traits are in agreement with van Hulzen et al. [29], who in a theoretical study also identified that an earlier resistance
B Parameters related to population dynamics

<table><caption>Table 1: Parameters for management and population dynamics used in the *Mycobacterium avium subsp. paratuberculosis* (Map) infection dynamics model within a structublackcol dairy herd (reproduced from [2], Table 1).</caption><thead><tr><th>Notation</th><th>Value</th><th>Definition</th><th>Source</th></tr></thead><tbody><tr><td>σ<sub>B</sub></td><td>0.07</td><td>Mortality rate of calves at birth</td><td>*, [17]</td></tr><tr><td>σ<sub>m</sub></td><td>0.206</td><td>Exit rate of male calves, weeks 2 to 4 (per week)</td><td></td></tr><tr><td>σ<sub>c1</sub></td><td>0.015</td><td>Death rate of female calves, weeks 1 and 2 (per week)</td><td>[17]</td></tr><tr><td>σ<sub>c2</sub></td><td>0.0035</td><td>Death rate of female calves, weeks 3 to weaning (per week)</td><td>[9]</td></tr><tr><td>σ<sub>c3</sub></td><td>0.00019</td><td>Death rate of heifers from weaning to entry in adult group (per week)</td><td>†</td></tr><tr><td>σ<sub>h</sub></td><td>0.011</td><td>Sale rate of bblackcol heifers 10 weeks before 1st calving</td><td>†</td></tr><tr><td>σ<sub>Ai</sub></td><td>0.27,0.25,0.31, 0.31,0.62</td><td>Yearly culling rate of cows in adult group i: 1, 2, 3, 4 and above 5 respectively (%)</td><td>*, [1]</td></tr><tr><td>w</td><td>10</td><td>Weaning age (weeks)</td><td>[12]</td></tr><tr><td>y</td><td>52</td><td>Age when entering the young heifer group (weeks)</td><td></td></tr><tr><td>h</td><td>91</td><td>Age when entering the heifer group (weeks)</td><td>*</td></tr><tr><td>cal</td><td>130</td><td>Age when entering the adult group (weeks)</td><td>*, †</td></tr><tr><td>τ<sub>aa</sub></td><td>56.3</td><td>Mean time spent in adult age groups 1 to 4 (weeks)</td><td>*, †</td></tr><tr><td>b</td><td>5</td><td>Quantity of colostrum fed to calves (L/day for 3 days)</td><td>†</td></tr><tr><td>d</td><td>7</td><td>Quantity of milk fed to calves after 3 days (L/day/calf)</td><td>†</td></tr><tr><td>prop</td><td>0.85</td><td>Proportion of lactating cows</td><td>*</td></tr><tr><td>ε</td><td>25</td><td>Quantity of milk or colostrum produced (L/day/cow)</td><td>*</td></tr><tr><td>f<sub>1</sub></td><td>0.5</td><td>Quantity of feaces produced by a non-weaned calf (kg/day)</td><td>†</td></tr><tr><td>f<sub>2</sub></td><td>5.5</td><td>Quantity of feaces produced by a weaned calf (kg/day)</td><td>†</td></tr><tr><td>f<sub>Y</sub></td><td>10</td><td>Quantity of feaces produced by a heifer (kg/day)</td><td>†</td></tr><tr><td>f<sub>A</sub></td><td>30</td><td>Quantity of feaces produced by a cow (kg/day)</td><td>†</td></tr><tr><td>Graz</td><td>[14 – 26]</td><td>Grazing period (1 being the first week of the year)</td><td>†</td></tr><tr><td>K<sub>c</sub></td><td>110</td><td>Number of cows above which the heifer selling rate increases</td><td>†</td></tr></tbody></table>

<footnote>* Agricultural statistics.</footnote>

<footnote>† Based on expert opinion.</footnote>
<figure><img src="image_7.png" /></figure>

Model outputs

Increase in decay in susceptibility with age (h)

Decrease in duration before high shedding and clinically affected state ($v_L+v_{Is}$)

Increase in required infectious dose ($\alpha$)

Decrease in quantity of Map shed in faeces by high shedders and clinically affected animals ($\phi Feces_{Is}$)

Other phenotypic traits

<caption>**Figure 5** Total contribution of phenotypic traits to model output variance. Contribution includes the principal effect of a given factor and interaction effects in which this factor was involved divided by the number of factors involved.</caption>

acquisition would be crucial when it comes to control paratuberculosis using genetic selection. However, there is nowadays no available knowledge to implement a genetic selection on these traits. These traits are not easily measurable in field conditions.

A decrease in the quantity of Map shed in faeces by high shedders and clinically affected animals, which was also identified as an influential phenotypic trait, might be achieved thanks to genetic selection. Currently, it has been shown that genetic markers could be associated with the occurrence of shedding versus no shedding at all by animals in infected herds [23, 27, 47]. More precise knowledge is needed concerning our ability to select cattle that will shed less Map in faeces while in their last stage of infection.

While van Hulzen et al. [28] identified the increase in duration of the latency period as an effective phenotypic trait in controlling paratuberculosis through genetic selection, we highlighted that an increase in this latency period (this being the period between infection and the occurrence of a moderate detectable shedding) without delaying the start of the high shedding or clinically affected state did not influence Map spread dynamics in

the herd. We have shown that it will be more interesting to lengthen the incubation period, as this delays the occurrence of the high shedding or clinically affected state.

Phenotypic traits identified as influencing Map spread dynamics at the herd scale also are related to control measures currently implemented in infected herds in the field [14]. Therefore, a valuable interaction can be expected between routine control plans and innovative control through genetic selection.

The variation of several other traits did not influence Map spread dynamics: decrease in duration of transiently infectious state with a constant duration before moderate shedding state, decrease in quantity of Map shed in milk and colostrum irrespective of the animal infection state, decrease in quantity of Map shed in faeces by transiently infectious animals and moderate shedders, and decrease in probability of in utero transmission irrespective of animal infection state. A decrease in duration of moderate shedding state (from 104 to 60 weeks) did not influence Map spread dynamics. The range of variation modelled for this trait was lower than for other traits due to limitations inherent to the compartmental model.
C Parameters related to infection dynamics

<table><caption>Table 2: Parameters for infection and transmission used in the *Mycobacterium avium subsp. paratuberculosis* (*Map*) infection dynamics model within a structublackcol dairy herd (reproduced from [2] - Table 2).</caption><thead><tr><th>Notation</th><th>Value</th><th>Definition</th><th>Source</th></tr></thead><tbody><tr><td>p<sub>X</sub></td><td>p<sub>L</sub> = 0.149<br>p<sub>Is</sub> = 0.149<br>p<sub>Ic</sub> = 0.65</td><td>Probability of in utero transmission for cow in health state X<br>X = latently infected (L)<br>X = subclinically infected (Is)<br>X = clinically affected (Ic)</td><td>[4, 26]</td></tr><tr><td>u</td><td>52</td><td>Maximal age in the susceptible compartment (weeks)</td><td>[6, 23]</td></tr><tr><td>γ</td><td>0.1</td><td>Susceptibility follows an exponential decrease: exp(-γ(age - 1))</td><td>[27]</td></tr><tr><td>v<sub>X</sub></td><td>v<sub>T</sub> = 25<br>v<sub>L</sub> = 52<br>v<sub>Is</sub> = 104<br>v<sub>Ic</sub> = 26</td><td>Mean time spent in health state X (weeks)<br>X = transiently infectious (T)<br>X = latently infected (L)<br>X = subclinically infected (Is)<br>X = clinically affected (Ic)</td><td>[21]<br>[16, 14]<br>[13]<br>†</td></tr><tr><td>sh<sub>X</sub></td><td>sh<sub>L</sub> = 0<br>sh<sub>Is</sub> = 0.4<br>sh<sub>Ic</sub> = 0.9</td><td>Probability of shedding in colostrum or milk for a cow in health state X<br>X = latently infected (L)<br>X = subclinically infected (Is)<br>X = clinically affected (Ic)</td><td>[20, 19]</td></tr><tr><td>α</td><td>10<sup>6</sup></td><td>Map infectious dose</td><td>[3]</td></tr><tr><td>β<sub>l</sub></td><td>5 × 10<sup>-4</sup> × 7</td><td>Transmission rate if ingestion of an infectious dose (per week)</td><td>‡</td></tr><tr><td>β<sub>c</sub></td><td>5 × 10<sup>-5</sup> × 7</td><td>Transmission rate if one infectious dose is present in the local environment (per week)</td><td>[21]</td></tr><tr><td>β<sub>g</sub></td><td>9.5 × 10<sup>-7</sup> × 7</td><td>Transmission rate if one infectious dose is present in the global environment (per week)</td><td>[21]</td></tr><tr><td>β<sub>o</sub></td><td>5 × 10<sup>-6</sup> × 7</td><td>Transmission rate if one infectious dose is present in the pasture (per week)</td><td>‡</td></tr><tr><td>g<sub>X</sub></td><td>g<sub>Is</sub> = 1 - 0.11<br>g<sub>Ic</sub> = 1 - 0.25</td><td>Decrease in milk production for cattle in health state X (per week)<br>X = subclinically infected (Is)<br>X = clinically affected (Ic)</td><td>[15]</td></tr><tr><td>μ<sub>k</sub></td><td>μ<sub>g</sub><sup>int</sup> = 0.4<br>μ<sub>g</sub><sup>ext</sup> = 1/14<br>μ<sub>cp</sub> = 0.17</td><td>Removal rate of Map from environment k<br>all the environments (per week)<br>all the environments (per week)<br>collective pens (when empty)</td><td>[7, 24]</td></tr></tbody></table>

† Expert opinions.

‡ Parameters' values are assumed.
<figure><img src="image_8.png" /><figcaption>**Figure 6** Effect of combined variations of the four most influential phenotypic traits on cumulative incidence. Cumulative incidence was calculated 25 years after initial Map introduction. The external dashed circle corresponds to the lowest cumulative incidence (0.35 newly infected animals) obtained among tested scenarios (log₁₀(0.35) = -0.46), the internal dashed circles correspond to thresholds of 25 (log₁₀(25) = 1.39) and 12 (log₁₀(12) = 1.09) newly infected animals. Asterisk corresponds to the cumulative incidence (log₁₀(617) = 2.79) for the scenario with current values of phenotypic traits. *h*, decay in susceptibility with age; *α*, increased required infectious dose; *v*<sub>L</sub> + *v*<sub>IS</sub>, increased duration before high shedding or clinically affected state; *φFeces*<sub>IC</sub>, decreased quantity of bacteria shed by high shedders or clinically affected animals. Each leaf is one scenario with branches representing variations of the four traits. Scenarios are presented by increasing level of variation in each trait successively (*h*, *α*, *v*<sub>L</sub> + *v*<sub>IS</sub>, and, *φFeces*<sub>IC</sub>). Tested values are given in Table 1.</figcaption></figure>

Nevertheless, as no effect was evidenced with a reduction of one-third of that duration, we assumed that this trait was not highly influential over the simulated range.

The four traits identified as influential are well described in the literature therefore we can assume that their tested ranges of variation were realistic. We assumed extreme ranges of variation for traits for which information was missing. The other traits assessed were not influential even with such extreme, non-realistic, variations. Using a different set of variation in the investigated traits is not expected to change our conclusions concerning which traits influence Map spread within dairy cattle herds.

As our objective was to assess Map spread in herds in which phenotypic traits would have been improved, we did not account for the long time needed [29] to reach such targeted levels of phenotypic traits by a potential genetic selection. On the one hand, recent studies identified several genetic markers associated with resistance (reviewed in [48, 49]), but genes and mechanisms

responsible for the tested phenotypic traits are still unknown. Further genetic studies of resistance of cattle to paratuberculosis are required, especially to identify genes and mechanisms involved in these relevant phenotypic traits to allow potential future selection of more resistant cattle. In addition, diagnostic tests currently available in the field do not allow identifying animals having the phenotypic traits identified here. Concerning future genetic selection, tests more sensitive during the early stage of the infection would be needed to distinguish infected animals from others and to better quantify the individual duration of incubation periods. On the other hand, there is a risk of a negative association between phenotypic traits of resistance to paratuberculosis and other traits of economic importance. For example, it has been shown that genetic markers associated with susceptibility to paratuberculosis could be associated to lactation persistence [50].

Our model represents a typical western Europe farming system for dairy cattle herds. Demographic processes have been shown to highly influence the disease dynamics [32], and therefore, different farming systems could change the influence of the studied phenotypic traits on Map spread dynamics in the herd. We assumed a closed herd without introduction of animals from other herds. Animal exchanges between herds could reintroduce

<figure><img src="image_9.png" /><figcaption>**Figure 7** Threshold of variation in influential parameters and combinations necessary to reach a low cumulative incidence.</figcaption></figure>

Squares, triangles and dots represent the threshold value for the parameter "duration before high shedding and clinically affected state" needed to reach the cumulative incidence <25 over 25 years of simulation at the given value of the three other influential parameters (required infectious dose, Map shedding in faeces by high shedders and clinically affected animals, and decay in susceptibility). The ten empty positions corresponds to combinations in variations where "duration before high shedding and clinically affected state" have to be more than doubled (>V2) to have a low cumulative incidence (<25). Combinations represented here account only for threshold of variations below V2. Tested values (Vref, V1, and V2) are given in Table 1.
<table><caption>Table 3: Summary of published data and assumed distributions of the quantities of *Mycobacterium avium* subsp. *paratuberculosis* (*Map*) shed, depending on the health state (X) and the route of transmission (r) in the *Map* infection dynamics model within a structublackcol dairy herd (reproduced from [2] - Table 3).</caption><thead><tr><th rowspan="2">Route of transmission (r)</th><th rowspan="2">Health state (X)</th><th colspan="4">Literature</th><th rowspan="2">Assumed distributions</th></tr><tr><th>Minimal value</th><th>Maximal value</th><th>Mean value</th><th>Source</th></tr></thead><tbody><tr><td rowspan="2">Map direct shedding in milk and colostrum (Map/L)</td><td>Subclinically infected</td><td>2.2 × 10<sup>4</sup></td><td>8.8 × 10<sup>4</sup></td><td>5 × 10<sup>4</sup></td><td>[20]</td><td>10<sup>5</sup> × Beta(8; 8)</td></tr><tr><td>Clinically affected</td><td>-</td><td>-</td><td>5 × 10<sup>4</sup></td><td>[5]</td><td>10<sup>5</sup> × Beta(8; 8)</td></tr><tr><td rowspan="2">Map indirect shedding in milk and colostrum (faecal contamination) (Map/L)</td><td>Subclinically infected</td><td>0</td><td>2 × 10<sup>10</sup></td><td>40</td><td>[10, 22]</td><td>1 + 10<sup>3</sup> × Beta(1; 25)</td></tr><tr><td>Clinically affected</td><td>700</td><td>2 × 10<sup>10</sup></td><td>14 × 10<sup>4</sup></td><td>[10, 22]</td><td>10<sup>(3+10×Beta(50;200))</sup></td></tr><tr><td rowspan="3">Map shedding in faeces (Map/kg)</td><td>Transiently infectious</td><td>6 × 10<sup>4</sup></td><td>6.3 × 10<sup>5</sup></td><td>3 × 10<sup>5</sup></td><td>[21]</td><td>10<sup>6</sup> × Beta(8; 19)</td></tr><tr><td>Subclinically infected</td><td>10<sup>4</sup></td><td>10<sup>15</sup></td><td>2.6 × 10<sup>6</sup></td><td>[18]</td><td>10<sup>(4+10×Beta(2.65;17))</sup></td></tr><tr><td>Clinically affected</td><td>10<sup>8</sup></td><td>10<sup>15</sup></td><td>10<sup>10</sup></td><td>[25, 8]</td><td>10<sup>(8+10×Beta(2;17))</sup></td></tr></tbody></table>
Map in free herds and thus influence disease dynamics. However, it is not expected to modify our conclusions as regards the identification of crucial phenotypic traits to better manage infected herds. Indeed, a single Map introduction can lead to infection persistence in 40% of the cases under current situation as regards phenotypic traits [32], with a huge cumulative incidence reached after 25 years if no control is applied. Animal movements are not expected to modify significantly this finding. However, the occurrence of animal movements might increase the cumulative incidence under controlled situations with improved traits.

This study highlighted four phenotypic traits of resistance of cattle to paratuberculosis influencing Map spread within a dairy herd: decay in susceptibility with age, quantity of Map shed in faeces by high shedders and clinically affected animals, duration of the incubation period, and required infectious dose. A combination of these traits strongly contributes to limit Map spread. Further genetic study should aim at better identifying cattle genes involved in these traits in order to allow their potential selection.

## Additional file 1

**Additional file 1. Description of the used within herd model of Map spread.**

### Competing interests

The authors declare that they have no competing interests.

### Authors' contributions

Conceived and designed the study: RBR, CF, PE, and RG. Adapted the model: GC, RBR, GB, and PE. Performed the simulations: RBR. Performed statistical analyses of simulated data: RBR and PE. Interpreted the results and wrote the paper: RBR, CF, and PE. All authors read and approved the final manuscript.

### Acknowledgements

The authors acknowledge the funding from the INRA metaprogram "Integrated Management of Animal Health" for the "GISA-PICSAR" project, from the Regional Council of Pays de la Loire, from APIS-GENE and GDS France for the "PARADIGM" project, with additional support (computational resources and expertise) from The French Research Agency, program investments for the future, Grant #ANR-10-BINF-07 (MIHMES), and the European fund for the regional development of Pays de la Loire. The funders had no role in the production, analysis, interpretation and writing of the data. Alain Joly (Animal Health Services of Brittany), and Elisabeta Vergu and Sandie Arnoux (INRA) are thanked for their helpful comments on model analyses and output interpretation. Pr. Simon MORE is thanked for his help in improving the language editing of the manuscript.

### Author details

### Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Received: 24 February 2017 Accepted: 22 August 2017

Published online: 10 October 2017

## References

1. Nielsen SS, Toft N (2009) A review of prevalences of paratuberculosis in farmed animals in Europe. Prev Vet Med 88:1–14

2. Whittington RJ, Windsor PA (2009) In utero infection of cattle with *Mycobacterium avium* subsp. paratuberculosis: a critical review and meta-analysis. Vet J 179:60–69

3. Mortier R, Orsel K, Barkema HW, Atkins G, De Buck J (2011) Age and dose dependent susceptibility to *Mycobacterium avium* subsp. paratuberculosis in dairy cattle. WCDS Adv Dairy Technol 23:375

4. Good M, Clegg T, Sheridan H, Yearsly D, O'Brien T, Egan J, Mullowey P (2009) Prevalence and distribution of paratuberculosis (Johnne's disease) in cattle herds in Ireland. Ir Vet J 62:597–606

5. Raizman EA, Wells SJ, Muñoz-Zanzi CA, Tavornpanich S (2011) Estimated within-herd prevalence (WHP) of *Mycobacterium avium* subsp. paratuberculosis in a sample of Minnesota dairy herds using bacterial culture of pooled fecal samples. Can J Vet Res 75:112–116

6. Mortier RA, Barkema HW, Orsel K, Wolf R, De Buck J (2014) Shedding patterns of dairy calves experimentally infected with *Mycobacterium avium* subspecies paratuberculosis. Vet Res 45:71

7. Mortier RA, Barkema HW, De Buck J (2015) Susceptibility to and diagnosis of *Mycobacterium avium* subspecies paratuberculosis infection in dairy calves: a review. Prev Vet Med 121:189–198

8. Stewart DJ, Vaughan JA, Stiles PL, Noske PJ, Tizard ML, Prowse SJ, Michalski WP, Butler KL, Jones SL (2007) A long-term bacteriological and immunological study in Holstein-Friesian cattle experimentally infected with *Mycobacterium avium* subsp. paratuberculosis and necropsy culture results for Holstein-Friesian cattle, Merino sheep and Angora goats. Vet Microbiol 122:83–96

9. van Roermund HJW, Bakker D, Willemsen PTJ, de Jong MC (2007) Horizontal transmission of *Mycobacterium avium* subsp. paratuberculosis in cattle in an experimental setting: calves can transmit the infection to other calves. Vet Microbiol 122:270–279

10. Mitchell RM, Medley GF, Collins MT, Schukken YH (2012) A meta-analysis of the effect of dose and age at exposure on shedding of *Mycobacterium avium* subspecies paratuberculosis (MAP) in experimentally infected calves and cows. Epidemiol Infect 140:231–246

11. Nielsen SS (2008) Transitions in diagnostic tests used for detection of *Mycobacterium avium* subsp. paratuberculosis infections in cattle. Vet Microbiol 132:274–282

12. Matthews HT (1947) On Johnne's disease. Vet Rec 59:397–401

13. Ezanno P, van Schaik G, Weber MF, Heesterbeek JA (2005) A modeling study on the sustainability of a certification-and-monitoring program for paratuberculosis in cattle. Vet Res 36:811–826

14. Bastida F, Juste RA (2011) Paratuberculosis control: a review with a focus on vaccination. J Immune Based Ther Vaccines 9:8

15. Beaunée G, Vergu E, Ezanno P (2015) Modelling of paratuberculosis spread between dairy cattle farms at a regional scale. Vet Res 46:111

16. Kalis CH, Hesselink JW, Barkema HW, Collins MT (2001) Use of long-term vaccination with a killed vaccine to prevent fecal shedding of *Mycobacterium avium* subsp. paratuberculosis in dairy herds. Am J Vet Res 62:270–274

17. Behr MA, Collins DM (2010) Paratuberculosis: organism, disease. Control, vol 1. CABI, Wallingford, UK, Cambridge, USA, p 1–388

18. Kirkpatrick BW, Shook GE (2011) Genetic susceptibility to paratuberculosis. Vet Clin N Am Food Anim Pract 27:559–571

19. van Hulzen KJ, Nielen M, Koets AP, de Jong G, van Arendonk JA, Heuven HC (2011) Effect of herd prevalence on heritability estimates of antibody response to *Mycobacterium avium* subspecies paratuberculosis. J Dairy Sci 94:992–997

20. Küpper J, Brandt H, Donat K, Erhardt G (2012) Heritability estimates for *Mycobacterium avium* subspecies paratuberculosis status of German Holstein cows tested by fecal culture. J Dairy Sci 95:2734–2739

21. Zare Y, Shook GE, Collins MT, Kirkpatrick BW (2014) Short communication: heritability estimates for susceptibility to *Mycobacterium avium* subspecies paratuberculosis infection defined by ELISA and fecal culture test results in Jersey cattle. J Dairy Sci 97:4562–4567

22. Purdie AC, Plain KM, Begg DJ, de Silva K, Whittington RJ (2011) Candidate gene and genome-wide association studies of *Mycobacterium avium* subsp. paratuberculosis infection in cattle and sheep: a review. Comp Immunol Microbiol Infect Dis 34:197–208

<footnote n="1">¹ BIOEPAR, INRA, ONIRIS, 44307 Nantes, France. ² MAIAGE, INRA, 78352 Jouy-en-Josas, France.</footnote>
Bibliography

[1] Beaudeau, F. (1999). *Review of culling and replacement practices in dairy cattle herds, Master Thesis. Institut Supérieur des Productions Animales*. PhD thesis, Rennes, France.

[2] Beaunée, G., Vergu, E., and Ezanno, P. (2015). Modelling of paratuberculosis spread between dairy cattle farms at a regional scale. *Veterinary Research*, 46:1-13.

[3] Begg, D. J. and Whittington, R. J. (2008). Experimental animal infection models for Johne's disease, an infectious enteropathy caused by Mycobacterium avium subsp. paratuberculosis. *The Veterinary Journal*, 176(2):129-145.

[4] Benedictus, A., Mitchell, R. M., Linde-Widmann, M., Sweeney, R., Fyock, T., Schukken, Y. H., and Whitlock, R. H. (2008). Transmission parameters of Mycobacterium avium subspecies paratuber-culosis infections in a dairy herd going through a control program. *Preventive Veterinary Medicine*, 83(3-4):215-227.

[5] Giese, S. B. and Ahrens, P. (2000). Detection of Mycobacterium avium subsp paratuberculosis in milk from clinically affected cows by PCR and culture. *Veterinary Microbiology*, 77(3-4):291-297.

[6] Hagan, W. A. (1938). Age as a Factor in Susceptibility to Johne's Disease. *Cornell Veterinarian*, 28:34-40.

[7] Jørgensen, J. B. (1977). Survival of Mycobacterium paratuberculosis in slurry. *Nordisk Veterinaer Medicin*, 29(6):267-270.

[8] Jørgensen, J. B. (1982). An improved medium for culture of Mycobacterium paratuberculosis from bovine faeces. *Acta veterinaria Scandinavica*, 23(3):325-335.

[9] Jégou, V., Porhiel, J., and Brunschwig, P. (2006). Risk management factors affecting mortality among dairy calves herds in 80 herds in Brittany. In *Proc. Journées Bovines Nantaises*, Nantes.

[10] Magnusson, M., Christiansson, A., Svensson, B., and Kolstrup, C. (2006). Effect of Different Premilking Manual Teat-Cleaning Methods on Bacterial Spores in Milk. *Journal of Dairy Science*, 89(10):3866-3875.

[11] Marcé, C., Ezanno, P., Seegers, H., Pfeiffer, D. U., and Fourichon, C. (2011). Predicting fadeout versus persistence of paratuberculosis in a dairy cattle herd for management and control purposes: a modelling study. *Veterinary Research*, 42(36).

[12] Marcé, C., Guatteo, R., Bareille, N., and Fourichon, C. (2010). Dairy calf housing systems across Europe and risk for calf infectious diseases. *animal*, 4(09):1588-1596.

[13] Matthews, H. T. (1947). On Johne's disease. *The Veterinary Record*, 59(31):397-401.
23. Kirkpatrick BW, Shi X, Shook GE, Collins MT (2011) Whole-Genome association analysis of susceptibility to paratuberculosis in Holstein cattle. *Anim Genet* 42:149–160

24. van Hulzen KJ, Koets AP, Nielen M, Hoeboer J, van Arendonk JA, Heuven HC (2012) Genetic variation for infection status as determined by a specific antibody response against *Mycobacterium avium* subspecies paratuberculosis in milk of Dutch dairy goats. *J Dairy Sci* 95:6145–6151

25. Alpay F, Zare Y, Kamalludin MH, Huang X, Shi X, Shook GE, Collins MT, Kirkpatrick BW (2014) Genome-wide association study of susceptibility to infection by *Mycobacterium avium* subspecies paratuberculosis in Holstein cattle. *PLoS One* 9:e111704

26. Zanella R, Settles ML, McKay SD, Schnabel R, Taylor JF, Whitlock RH, Schukken Y, Van Kessel JS, Smith JM, Neibergs HL (2011) Identification of loci associated with tolerance to Johne's disease in Holstein cattle. *Anim Genet* 42:28–38

27. Settles M, Zanella R, McKay SD, Schnabel RD, Taylor JF, Whitlock R, Schukken Y, Van Kessel JS, Smith JM, Neibergs H (2009) A whole genome association analysis identifies loci associated with *Mycobacterium avium* subsp. paratuberculosis infection status in US holstein cattle. *Anim Genet* 40:655–662

28. Neibergs HL, Settles ML, Whitlock RH, Taylor JF (2010) GSEA-SNP identifies genes associated with Johne's disease in cattle. *Mamm Genome* 21:419–425

29. van Hulzen KJ, Koets AP, Nielen M, Heuven HC, van Arendonk JA, Klinkenberg D (2014) The effect of genetic selection for Johne's disease resistance in dairy cattle: results of a genetic-epidemiological model. *J Dairy Sci* 97:1762–1773

30. Marcé C, Ezanno P, Weber MF, Seegers H, Pfeiffer DU, Fourichon C (2010) Invited review: modeling within-herd transmission of *Mycobacterium avium* subspecies paratuberculosis in dairy cattle: a review. *J Dairy Sci* 93:4455–4470

31. Lu Z, Schukken YH, Smith RL, Grohn YT (2010) Stochastic simulations of a multi-group compartmental model for Johne's disease on US dairy herds with test-based culling intervention. *J Theor Biol* 264:1190–1201

32. Marcé C, Ezanno P, Seegers H, Pfeiffer DU, Fourichon C (2011) Predicting fadeout versus persistence of paratuberculosis in a dairy cattle herd for management and control purposes: a modelling study. *Vet Res* 42:36

33. Marcé C, Ezanno P, Seegers H, Pfeiffer DU, Fourichon C (2011) Within-herd contact structure and transmission of Mycobacterium avium subspecies paratuberculosis in a persistently infected dairy cattle herd. *Prev Vet Med* 100:116–125

34. Robins J, Bogen S, Francis A, Westhoek A, Kanarek A, Lenhart S, Eda S (2015) Agent-based model for Johne's disease dynamics in a dairy herd. *Vet Res* 46:68

35. Martcheva M, Lenhart S, Eda S, Klinkenberg D, Momotani E, Stabel J (2015) An immuno-epidemiological model for Johne's disease in cattle. *Vet Res* 46:69

36. Smith RL, Schukken YH, Gröhn YT (2015) A new compartmental model of *Mycobacterium avium* subsp. paratuberculosis infection dynamics in cattle. *Prev Vet Med* 122:298–305

37. Koets AP, Gröhn YT (2015) Within- and between-host mathematical modeling of *Mycobacterium avium* subspecies paratuberculosis (MAP) infections as a tool to study the dynamics of host-pathogen interactions in bovine paratuberculosis. *Vet Res* 46:60

38. Al-Mamun MA, Smith RL, Schukken YH, Gröhn YT (2016) Modeling of *Mycobacterium avium* subsp. paratuberculosis dynamics in a dairy herd: an individual based approach. *J Theor Biol* 408:105–117

39. Pradhan AK, Mitchell RM, Kramer AJ, Zurakowski MJ, Fyock TL, Whitlock RH, Smith JM, Hovingh E, Van Kessel JA, Kams JS, Schukken YH (2011) Molecular epidemiology of *Mycobacterium avium* subsp. paratuberculosis in a longitudinal study of three dairy herds. *J Clin Microbiol* 49:893–901

40. Råberg L, Sim D, Read AF (2007) Disentangling genetic variation for resistance and tolerance to infectious diseases in animals. *Science* 318:812–814

41. Schneider DS, Ayres JS (2008) Two ways to survive infection: what resistance and tolerance can teach us about treating infectious diseases. *Nat Rev Immunol* 8:889–895

42. R Core Team (2016) R: A language and environment for statistical computing. Vienna, Austria. R Foundation for Statistical. http://www.R-project.org/

43. Everitt BS, Hothorn T (2010) Cluster analysis: classifying Romano-British pottery and exoplanets. In: A handbook of statistical analyses using R, 2nd ed. CRC Press, Taylor and Francis Group edition, Boca Raton, p 315–348. http://www.bagualu.net/wordpress/wp-content/uploads/2015/10/A_Handbook_of_Statistical_Analyses_Using_R_Second_Edition.pdf

44. Thinsungnoen T, Kaoungku N, Durongdumronchai P, Kerdprasop K, Kerdprasop N (2015) The clustering validity with silhouette and sum of squared errors. In: Proceedings of the 3rd international conference on industrial application engineering, p 44–51. https://www2.ia-engineers.org/iciae/index.php/iciae2015/paper/viewFile/576/380

45. MacQueen J (1967) Some methods for classification and analysis of multivariate observations. In: Proceedings of the fifth berkeley symposium on mathematical statistics and probability, vol 1. University of California Press, Berkeley, California, p 281–297. https://projecteuclid.org/download/pdf_1/euclid.bmsp/1200512992

46. Husson F, Josse J, Le S, Mazet J (2016) FactoMineR: multivariate exploratory data analysis and data mining. R package version 1.32. https://CRAN.R-project.org/package=FactoMineR

47. Zare Y, Shook GE, Collins MT, Kirkpatrick BW (2014) Genome-wide association analysis and genomic prediction of *Mycobacterium avium* subspecies paratuberculosis infection in US Jersey cattle. *PLoS One* 9:e88380

48. McSpadden K, Caires K, Zanella R (2013) The effect of *Mycobacterium avium* subspecies paratuberculosis exposure on animal health. *Acta Sci Vet* 41:1095

49. Pauciullo A, Küpper J, Brandt H, Donat K, Iannuzzi L, Erhardt G (2015) Wingless-type MMTV integration site family member 2 (WNT2) gene is associated with resistance to MAP in faecal culture and antibody response in Holstein cattle. *Anim Genet* 46:122–132

50. Carvajal AM, Huircan P, Lepori A (2013) Single nucleotide polymorphisms in immunity-related genes and their association with mastitis in Chilean dairy cattle. *Genet Mol Res* 12:2702–2711

51. Rankin JD (1962) The experimental infection of cattle with *Mycobacterium bohnei*. IV. Adult cattle maintained in an infectious environment. J Comp Pathol 72:113–117

52. Hagan WA (1938) Age as a factor in susceptibility to Johne's disease. Cornell Vet 28:34–40

53. Whitlock RH, Buergelt C (1996) Preclinical and clinical manifestations of paratuberculosis (including pathology). *Vet Clin N Am Food Anim Pract* 12:345–356

54. Windsor PA, Whittington RJ (2010) Evidence for age susceptibility of cattle to Johne's disease. *Vet J* 184:37–44

55. Begg DJ, Whittington RJ (2008) Experimental animal infection models for Johne's disease, an infectious enteropathy caused by *Mycobacterium avium* subsp. paratuberculosis. *Vet J* 176:129–145

56. Sweeney RW, Whitlock RH, Rosenberger AE (1992) *Mycobacterium paratuberculosis* isolated from fetuses of infected cows not manifesting signs of the disease. Am J Vet Res 53:477–480

57. Giese SB, Ahrens P (2000) Detection of *Mycobacterium avium* subsp. paratuberculosis in milk from clinically affected cows by PCR and culture. *Vet Microbiol* 77:291–297

58. Rossiter CA, Burhans WS (1996) Farm-specific approach to paratuberculosis (Johne's disease) control. *Vet Clin N Am Food Anim Pract* 12:383–415

59. Whittington RJ, Reddacliff LA, Marsh I, McAllister S, Saunders V (2000) Temporal patterns and quantification of excretion of *Mycobacterium avium* subsp. paratuberculosis in sheep with Johne's disease. Aust Vet J 78:34–37

60. Jørgensen JB (1982) An improved medium for culture of *Mycobacterium paratuberculosis* from bovine faeces. Acta Vet Scand 23:325–335

61. Benedictus A, Mitchell RM, Linde-Widmann M, Sweeney R, Fyock T, Schukken YH, Whitlock RH (2008) Transmission parameters of *Mycobacterium avium* subspecies paratuberculosis infections in a dairy herd going through a control program. *Prev Vet Med* 83:215–227
[14] Nielsen, S. S. (2008). Transitions in diagnostic tests used for detection of Mycobacterium avium subsp. paratuberculosis infections in cattle. *Veterinary Microbiology*, 132(3-4):274-282.

[15] Nielsen, S. S., Enevoldsen, C., and Toft, N. (2006). Milk production losses associated with bovine paratuberculosis diagnosed from repeated testing. In *11th international Symposia of Veterinary Epidemiology and Economics*, Cairns.

[16] Nielsen, S. S. and Ersbøll, A. K. (2006). Age at Occurrence of Mycobacterium avium Subspecies paratuberculosis in Naturally Infected Dairy Cows. *Journal of Dairy Science*, 89(12):4557-4566.

[17] Rio, O. (1999). *Frequency and risks of mortality and health disorders of calves in dairy cattle herds*. PhD thesis, Nantes.

[18] Rossiter, C. A. and Burhans, W. S. (1996). Farm-specific approach to paratuberculosis (Johne's disease) control. *The Veterinary clinics of North America. Food animal practice*, 12(2):383-415.

[19] Streeter, R. N., Hoffsis, G. F., Bechnielsen, S., Shulaw, W. P., and Rings, M. (1995). Isolation of Mycobacterium-Paratuberculosis From Colostrum and Milk of Subclinically Infected Cows. *American Journal of Veterinary Research*, 56(10):1322-1324.

[20] Sweeney, R. W., Whitlock, R. H., and Rosenberger, A. E. (1992). Mycobacterium paratuberculosis isolated from fetuses of infected cows not manifesting signs of the disease. *American Journal of Veterinary Research*, 53(4):477-480.

[21] van Roermund, H. J. W., Bakker, D., Willemsen, P. T. J., and de Jong, M. C. M. (2007). Horizontal transmission of Mycobacterium avium subsp. paratuberculosis in cattle in an experimental setting: Calves can transmit the infection to other calves. *Veterinary Microbiology*, 122(3-4):270-279.

[22] Vissers, M. M. M., Driehuis, F., Te Giffel, M. C., De Jong, P., and Lankveld, J. M. G. (2006). Improving Farm Management by Modeling the Contamination of Farm Tank Milk with Butyric Acid Bacteria. *Journal of Dairy Science*, 89(3):850-858.

[23] Whitlock, R. H. and Buergelt, C. (1996). Preclinical and clinical manifestations of paratuberculosis (including pathology). *Veterinary Clinics of North America-Food Animal Practice*, 12(2):345-356.

[24] Whittington, R. J., Marshall, D. J., Nicholls, P. J., Marsh, I. B., and Reddacliff, L. A. (2004). Survival and dormancy of Mycobacterium avium subsp. paratuberculosis in the environment. *Applied and Environmental Microbiology*, 70(5):2989-3004.

[25] Whittington, R. J., Reddacliff, L. A., Marsh, I., McAllister, S., and Saunders, V. (2000). Temporal patterns and quantification of excretion of Mycobacterium avium subsp paratuberculosis in sheep with Johne's disease. *Australian Veterinary Journal*, 78(1):34-37.

[26] Whittington, R. J. and Windsor, P. A. (2009). In utero infection of cattle with Mycobacterium avium subsp. paratuberculosis: A critical review and meta-analysis. *The Veterinary Journal*, 179(1):60-69.
[27] Windsor, P. A. and Whittington, R. J. (2010). Evidence for age susceptibility of cattle to Johne's disease. *The Veterinary Journal*, 184(1):37-44.