Frede Jared O'connell

Aakmann Tøgersen

Nicolas C Friggens

Peter Løvendahl

Søren Højsgaard

Combining Cattle Activity and Progesterone Measurements Using Hidden Semi-Markov Models

Keywords: Dairy cow, EM-algorithm, Oestrus detection, Online data, Streaming data, Time series. 1. INTRODUCTION

Hourly pedometer counts and irregularly measured concentration of the hormone progesterone were available for a large number of dairy cattle. A hidden semi-Markov was applied to this bivariate time-series data for the purposes of monitoring the reproductive status of cattle. In particular, the ability to identify oestrus is investigated as this is of great importance to farm management. Progesterone concentration is a more accurate but more expensive method than pedometer counts, and we evaluate the added benefits of a model that includes this variable. The resulting model is biologically sensible, but validation is difficult. We utilize some auxiliary data to demonstrate the model's performance.

understanding of the reproductive cycle in cattle. This paper extends some initial ideas by [START_REF] Friggens | The Potential of On-Farm Fertility Profiles: In-Line Progesterone and Activity Measurements[END_REF]. Automated oestrus detection using pedometers has been previously explored using a variety of methodologies, see [START_REF] Firk | Automation of Oestrus Detection in Dairy Cows: A Review[END_REF] for a review.

In this paper we investigate hidden semi-Markov models as a possible methodology for monitoring the reproductive status of cattle based on proxy measurements. Progesterone level is a reliable indicator of oestrus (Friggens et al. 2008), however a financial cost is associated with each measurement taken. Activity measurements are essentially free after the initial set-up overhead, but the data is much noisier leading to missed oestruses and false positives. Univariate models for each measure and a bivariate model using both measures were estimated with the data to assess the utility of each variable and the improvements that can be gained by using both.

In Section 2, we give some biological background and describe the available data. In Section 3, we introduce hidden semi-Markov models and the methods of parameter estimation and prediction and show how they were applied to our data. Section 4 contains results and Section 5 summarizes our work and covers future areas of investigation.

BACKGROUND

THE OVARIAN CYCLE IN CATTLE

Cows are a polyoestrous animal, that is, once started, the oestrus cycle will continue until pregnancy or illness occurs (as opposed to the seasonal oestrus of some animals). The oestrus cycle takes approximately 21 days (but with a large variation) and can be divided into four stages, which can be further grouped into two longer stages (see Table 1). The use of days is an insufficiently granular timescale for the oestrus stage, as oestrus lasts between 6 and 30 hours [START_REF] Ball | Reproduction in Cattle[END_REF]. Ovulation occurs in the day following oestrus, so identifying oestrus allows a farm manager to know when to artificially inseminate a cow.

Following ovulation a structure called the corpus luteum forms in the ovary. The corpus luteum produces the hormone progesterone and remains in the ovary until a few days prior to the next ovulation (the luteal phase), at which time it degenerates rapidly. That is, progesterone is high during the luteal phase and low during the follicular phase. Progesterone directly reflects the biological processes occurring in the ovary and is a very useful indicator of reproductive status. Progesterone can be automatically measured from milk samples, with milking occurring every 6-20 hours (in robotic milking systems). A typical progesterone profile after calving can be viewed in Figure 1 (bottom). Figure 1. Top: The raw hourly activity counts from a pedometer against time since calving (grey), the black line is as 24-hour centred moving average. Bottom: The progesterone concentration over the same time and cow. The dashed vertical lines are times when an artificial insemination occurred (an indicator of oestrus). Note that these occur after a drop in progesterone and during a brief spike in activity levels.

In the period leading up to ovulation the cow will try to attract the attention of a bull by standing to be mounted, mounting other cows or being mounted by other cows. This is sometimes referred to as standing heat and is traditionally how a stockman would identify a cow that is about to ovulate, allowing them to proceed with artificial insemination or bring the cow to a bull. This behaviour leads to an increase in the number of counts on a pedometer the cow is wearing, and so can be exploited for automated detection of oestrus. Having a stockman manually detect oestrus in farms with hundreds of cows is expensive and inaccurate, so automated systems are of great interest. These spikes in activity (together with the corresponding drops in progesterone) are illustrated in Figure 1 (top).

DATA

The data was collected at the Danish Cattle Research Centre. Activity data consisted of counts per hour from ALPRO R pedometers by DeLaval. Progesterone is measured in nanograms per millilitre (ng/mL) at milking times (Friggens et al. 2008). Progesterone and activity data from 58 cows were analyzed. These cows were selected from a database of 424 cows using filtering rules described below.

Three criteria were applied to individual cows to determine whether they were suitable for this analysis. Firstly, a cow had to have become pregnant for us to consider it. The only sure case of an oestrus occurring is when a cow has become pregnant and we used this in our model validation. Secondly, only cows with both progesterone and activity data available were used in this work, since we wish to evaluate the benefits of measuring both variables. Thirdly, for reasons described in Section 3.4, series starting at the time of a calving event were desirable, so cows without a single measurement of either variable during the first 21 days after calving were also discarded. After filtering the data, we were left with 58 cows and data spanning 128,578 hours with 112,775 activity counts and 1,917 progesterone measurements. The median time between progesterone measurements was 49 hours.

We rounded the time of progesterone measurements to the nearest hour and merged these with activity measurements to create a time series of bivariate data, suitable to be analyzed as a Markov chain (in respect to the timing of observations). Most hours had only an activity count available, and in these cases progesterone was treated as a missing value. Useful auxiliary information was the date of artificial inseminations and the dates when a cow gave birth. Combining these allowed us to determine when a successful artificial insemination occurred and hence when we can be certain that oestrus occurred.

ACTIVITY INDEX

Cows are more active during the day and more passive during the night (diurnal cycle) and the level of activity does not remain the same over a cow's lifespan. Hence, we cannot treat activity levels as a stationary process. We choose to remove the effect of this diurnal cycle and any long-term trend. In short, we want to see how active a cow is at a given hour relative to its recent levels of activity. We applied a Holt-Winters model [START_REF] Holt | Forecasting Trends and Seasonals by Exponentially Weighted Moving Averages[END_REF][START_REF] Winter | Forecasting Sales by Exponentially Weighted Moving Averages[END_REF]) using an implementation by [START_REF] Meyer | Naive Time Series Forecasting Methods[END_REF] to achieve this. Rather than looking at raw activity we worked with x t = log(activity t + 1), to stabilize variance. A histogram of the transformed variable can be seen in Figure 2 (left).

The Holt-Winters seasonal method is a simple technique for seasonally adjusting data using exponential moving averages. Let x t be a univariate time series (in our case the transformed activity measure). Forecasting h steps ahead involves variables a t and s t with and the updating formulae

ât = α(x t -ŝt-p ) + (1 -α) ât-1 and ŝt = β(x t -ât ) + (1 -β)ŝ t-p . (2.2)
Here we have a local component, a t , and a seasonal component, s t , with period, p, set to 24 to remove the diurnal trend. The parameters α and β are estimated by minimizing the sum of squared errors of one-step-ahead predictions.

If we assume that the Holt-Winters model captures a typical cow behaviour, then the residuals of the one-step-ahead predictions should be large when the cow deviates from normal behaviour. That is, when a cow becomes atypically active (a possible standing heat event) we should observe a large positive residual. Henceforth, we call the one-step-ahead residuals the activity index. A histogram of the activity index can be viewed in Figure 3 (centre). We used the activity index as it does not introduce any lag and does not unnecessarily smooth the data (smoothing is already implicit in the hidden semi-Markov models). Visual inspection of the activity index along with the timing of artificial inseminations and drops in progesterone suggests they were correlated, confirming that this index was a reasonable variable to use (Figure 3,right).

Alternative methods for oestrus detection solely using pedometer readings are in active development, see [START_REF] Jonsson | Oestrus Detection in Dairy Cows Using Likelihood Ratio Tests[END_REF] and references therein. [START_REF] Roelofs | Various Behavioral Signs of Estrous and their Relationship with Time of Ovulation in Dairy Cattle[END_REF] provide a method using simple statistical tests. [START_REF] Firk | Automation of Oestrus Detection in Dairy Cows: A Review[END_REF] review a number of methods. The aim of this research was to investigate the relation between progesterone and activity, rather than developing an absolute best technique for oestrus detection with activity readings. Furthermore, different activity measures (or indeed any variable) can be easily "plugged in" to hidden Markov or semi-Markov models at a later date.

METHODOLOGY

We give a brief summary of Markov chains, hidden Markov and semi-Markov models, or HMMs and HSMMs respectively. For a complete introduction we refer the reader to [START_REF] Rabiner | A Tutorial on Hidden Markov Models and Selected Applications in Speech Recognition[END_REF]. Implementation of all models was done using a combination of R (R Development Core Team 2008) and the C programming language. We have released this implementation as an R package called mhsmm (O'Connell and Højsgaard 2009b) and a paper has been submitted on its use (O'Connell and Højsgaard 2009a). We only consider the discrete time case here, as activity measures are hourly counts and the times of progesterone measurements were rounded to the nearest hour.

DISCRETE MARKOV CHAINS

A discrete Markov chain is a random process (in discrete time) taking discrete values (states) from the state space S, that is, S t ∈ S = {1, . . . , J } for t = 1, 2, . . . , T . The process S t is a Markov chain if it has the Markov property,

P (S t+1 = s t+1 |S 0 = s 0 , S 1 = s 1 , . . . , S t = s t ) = P (S t+1 = s t+1 |S t = s t ), (3.1)
for any s 0 , s 1 , . . . , s t+1 ∈ S. That is, the state at any given time t + 1 depends only on the previous states through the state at time t. Let p ij = P (S t+1 = j |S t = i) with the properties J j =1 p ij = 1 and p ij ≥ 0 so that P = (p ij ) is the transition matrix of the Markov chain. To fully specify the model we require the distribution of the initial state π i = P (S 0 = i).

For later use we notice that under this model, the probability of spending u consecutive time steps in state i is

d i (u) = P (S t+u+1 = i, S t+u = i, S t+u-1 = i, . . . , S t+2 = i|S t+1 = i, S t = i) (3.2) = p u-1 ii (1 -p ii ). (3.3)
We call this the duration density and this is inherently geometrically distributed for any Markov chain. This may be an unreasonable assumption for many real-world processes; see Section 3.3 for a discussion.

HIDDEN MARKOV MODELS

Suppose we can only observe a variable X t related to the state S t but not the state itself. The conditional distribution of the observed variable X t conditioned on the unobserved (or hidden) state S t is referred to as the emission distribution. We refer to the parameters defining such a process as a hidden Markov model, henceforth referred to as an HMM.

These models have been used for a variety of different applications, such as speech recognition [START_REF] Rabiner | A Tutorial on Hidden Markov Models and Selected Applications in Speech Recognition[END_REF], weather modelling [START_REF] Hughes | A Non-Homogeneous Hidden Markov Model for Precipitation Occurrence[END_REF] and DNA sequence analysis [START_REF] Krogh | A Hidden Markov Model that Finds Genes in E. Coli DNA[END_REF].

In addition to the usual Markov chain parameters, π and P, an HMM also requires an emission distribution to be defined:

b i (x t ) = P (X t = x t |S t = i).
(3.4)

For example b i (x) may be a multivariate Gaussian distribution. We assume that the emitted variables are conditionally independent given the underlying state, that is,

P (X t = x t , X t-1 = x t-1 |S t = i) = P (X t = x t |S t = i). (3.5)
A full HMM is hence specified by θ = (π, P, b). The Baum-Welch algorithm is the original procedure for estimating the parameters of a HMM [START_REF] Baum | A Maximization Technique Occurring in the Statistical Analysis of Probabilistic Functions of Markov Chains[END_REF]. This technique was later grouped with a more general class of algorithms for incomplete data, named the Expectation-Maximization (EM) algorithm [START_REF] Dempster | Maximum Likelihood from Incomplete Data Via the EM Algorithm[END_REF]. [START_REF] Rabiner | A Tutorial on Hidden Markov Models and Selected Applications in Speech Recognition[END_REF] provides a very clear overview.

HIDDEN SEMI-MARKOV MODELS

A severe limitation of standard HMMs is the implicit geometrically distributed duration time of states in the model (Equation (3.3)), as many of the real-world problems we wish to model do not have geometric sojourn times. For example, the reproductive states of cows are not a memoryless process since a follicular stage is likely to occur after 18 days in a luteal stage. A potential solution to this issue is to explicitly estimate d(u), producing what is referred as a hidden semi-Markov model, henceforth called an HSMM. [START_REF] Ferguson | Hidden Markov Analysis: An Introduction, Hidden Markov Models for Speech[END_REF] was the first to propose such models along with an algorithm to fit them; [START_REF] Rabiner | A Tutorial on Hidden Markov Models and Selected Applications in Speech Recognition[END_REF] provides a good summary. [START_REF] Guédon | Estimating Hidden Semi-Markov Chains From Discrete Sequences[END_REF] developed a more efficient algorithm and a method to deal with right censoring which we have implemented.

In a standard discrete Markov chain the state duration density has the form

d i (u) = p u
ii (1p ii ), whereas in a semi-Markov model, we model d i (u) explicitly. That is, in addition to the parameters already defined for the HMM, we now have explicit probabilities for our state duration density, θ = (π, P, b, d). [START_REF] Guédon | Estimating Hidden Semi-Markov Chains From Discrete Sequences[END_REF] utilizes a non-parametric probability mass function of d(u) in his derivations, but then proposes an ad hoc method to introduce parametric distributions (for example, Poisson). The latter procedure has been followed here.

Parameter Estimation

The complete data likelihood of a HSMM is

P (X = x, S = s; θ) = π s * 1 d s * 1 (u 1 ) R r=2 p s * r-1 s * r d s * r (u r ) p s * R-1 s * R D s * R (u R ) T t=1 b s t (x t ) (3.6)
where s * r is the rth visited state and u r is the time spent in that state. Guédon proposed using the survivor function

D i (u) = v≥u d i (v) (3.7)
for the sojourn time in the last state, so we do not have to assume the process is leaving a state immediately after time T . This is particularly important for online applications where we wish to estimate the most recent state so we can monitor some process. As we have not observed the state sequence, maximizing this likelihood constitutes an incomplete data problem. A local maximum can be found using the Expectation-Maximization, or EM, algorithm. We briefly outline the procedure below.

The EM algorithm involves iterating over two steps until convergence. In the E-step, we calculate the expected complete data likelihood given the value of the parameters at iteration k and the observed data: (k) .

Q θ θ (k) = E log P (X = x, S = s; θ) X = x; θ
(3.8)

This term is typically broken down into a sum of terms involving subsets of the parameters. The M-step then involves choosing θ (k+1) as the values that maximize Q(θ |θ (k) ). These steps are repeated until convergence. We now provide the specific terms calculated for an HSMM.

The E-step For HSMMs, the E-step involves estimating three separate terms: the probability of being in state i at time t, γ t (i) = P (S t = i|X = x; θ), (3.9) the probability that the process left state i at time t and entered state j at t + 1, ξ t (i, j ) = P (S t = i, S t+1 = j |X = x; θ), (3.10) and finally the expected number of times a process spent u time steps in state i,

η i,u = P (S u = i, S u-v = i, v = 1, . . . , u|X = x; θ) + T t=1 P (S t+u+1 = i, S t+u-v = i, v = 0, . . . , u -1, S t = i|X = x; θ). (3.11)
These values can be calculated via a dynamic programming method known as the forwardbackward algorithm. It should be noted that the derivations and algorithms involved are complex; we refer the reader to [START_REF] Guédon | Estimating Hidden Semi-Markov Chains From Discrete Sequences[END_REF] for details.

The M-step The starting and transition probabilities are implicit from the definition of the Markov chain,

π i = γ 1 (i), p ij = T -1 t=1 ξ t (i, j ) T -1 t=1 i =j ξ t (i, j )
.

(3.12) Guédon provides derivations for d i (u) as a non-parametric probability mass function, but then proposes an ad hoc solution for using parametric distributions with η i,u as weights which we have followed. We assume that the sojourn times are Gamma distributed, that is, U r |S r = i ∼ (c i , λ i ). We estimated the parameters as follows: the likelihood for the Gamma distribution can be maximized with respect to its parameters by solving

log( ĉi ) -ψ( ĉi ) = log( ūi ) -log u i (3.13)
where ψ( ) is the digamma function. We use

ūi = u η i,u u u η i,u , log u i = u η i,u log(u) u η i,u (3.14)
and then solve the equation using Newton's method [START_REF] Choi | Maximum Likelihood Estimation of the Parameters of the Gamma Distribution and Their Bias[END_REF]. The scale parameter is estimated as λi = ĉi / ūi . We assume X t are normally distributed given S = i, that is,

X t |S t = i ∼ N(μ i , σ 2 i ).
We can re-estimate the Gaussian parameters via the equations:

μi = T t=1 γ t (i)x t T t=1 γ t (i) , σ 2 i = T t=1 γ t (i)(x t -μi ) 2 T t=1 γ t (i)
.

(3.15) [START_REF] Sansom | Fitting Hidden Semi-Markov Models to Breakpoint Rainfall Data[END_REF] provide derivations for these re-estimating equations. In the case where X t is multivariate, for our application we assumed that the individual components of X i were conditionally independent given the current state and so take the joint density as the product of the individual densities.

Prediction: The Viterbi Algorithm

Given a model, we are interested in the most likely sequence of states, given the sequence of observations. That is, we wish to find the sequence of states that maximizes P (S|X, θ ). Calculating P (S|X, θ ) for every possible state sequence is not computationally feasible. A dynamic programming technique known as the Viterbi algorithm [START_REF] Forney | The Viterbi Algorithm[END_REF], can be used to maximize P (S|X, θ ) in a computationally feasible manner. It should be noted that this is distinct from taking argmax i γ t (i) as the state at time t, which is the individually most likely state. We start by defining the variable

δ t (j ) = max s 1 ,...,s t-1 P (S t+1 = j, S t = j, S 1 = s 1 , . . . , S t-1 = s t-1 , X 1 = x 1 , . . . , X t-1 = x t-1 ; θ) (3.16)
with the special case at the last time step

δ T (j ) = max s 1 ,...,s T -1 P (S T = j, S 1 = s 1 , . . . , S T -1 = s T -1 , X 1 = x 1 , . . . , X T = x T ; θ),
(3.17) where we are only concerned with the probability of being in state j at time T , not with leaving the state at T + 1. The variable δ can be calculated via induction with the initialization step (3.18) the recursion step (3.19) and the terminating step

δ 1 (j ) = b j (x 1 )d j (1)π j ,
δ t (j ) = b j (x t ) max max 1≤u≤t u-1 v=1 b j (x t-v ) d j (u) max i =j p ij δ t-u (i) ; t v=1 b j (x t-v ) d j (t + 1)π i ,
δ T (j ) = b j (x T ) max max 1≤u≤T u-1 v=1 b j (x T -v ) D j (u) max i =j p ij δ T -u (i) ; T v=1 b j (x T -v ) D j (T )π i . (3.20)
The likelihood of the optimum state sequence is then just max i δ T (i). To obtain the optimal sequence of states we must retain the values of u and i that maximize δ t (j ) at each step. We record the most likely value of u, that is the most likely sojourn time,

υ t (i) = t -argmax u,t max 1≤u≤t u-1 v=1 b i (x t-v ) d j (u) max i =j p ij δ t-u (j ) ; t v=1 b i (x t-v ) d j (t + 1)π j , (3.21)
and the most likely preceeding state,

ψ t (j ) = argmax i =j p ij δ υ t (i) (j ) , 2 ≤ t ≤ T . (3.22)
Finally, we apply a backtracking algorithm to find the optimal state sequence:

1. Set t = T 2. s t = argmax i (δ t (i)) 3. s υ t (s t )+1 , . . . , s t-1 = s t 4. s υ t (s t ) = ψ t (s t ) 5. Set t = υ t (s t )
6. If t = 1 stop, else return to 3

Missing Data

Missing data occurs in both of the observed variables for our application. Pedometer readings for a given hour may fail to be recorded due to hardware problems or a cow being away from the stable. There are many missing progesterone measurements due to the fact it is measured at a much sparser resolution than activity (less than 2% of activity recordings have a matching progesterone measurement). As such, we require a method for handling missing observations, which is straightforward for HSMMs. If an observation x t is missing, we set b j (x t ) = 1 for all j , that is, there is no contribution to the likelihood at that time step.

APPLICATION TO CATTLE DATA

Three models were investigated: two univariate models using each variable and one model using both variables. The states identifiable by each combination of variables differ. For each univariate model we used three states, but with different definitions. The activityonly model contained states S Activity = {post-partum anoestrus, standing heat, not standing heat} (Figure 4,top). For progesterone measurements we had S Progesterone = {post-partum anoestrus, luteal, follicular} (Figure 4,bottom). Standing heat occurs near the beginning of the follicular state but the follicular state is much longer. The bivariate model used six states, the combination of activity and progesterone allowing us to break the follicular stage up into three different sub-stages in addition to the PPA and luteal states (Figure 5).

We define t = 1 as the hour immediately following a calving. This is a useful way to arrange the data as we know that the cow will be in post-partum anoestrus at this time. Gaussian emission distributions were used for each variable. We worked with log(progesterone + 1) rather than raw progesterone to stabilize the variance. Gamma distributions were used for the distribution of sojourn times for each state. Starting values for distribution parameters were found via manual inspection of the data and biological knowledge of the ovarian cycle as detailed in Section 2.1.

RESULTS

Validation for this data is difficult as the true timing of ovulation is never known. One proxy is when a successful artificial insemination occurs, that is, an insemination that results in a cow giving birth approximately 280 days later. We take two approaches here, first investigating the ability of the activity-only model to predict the states of a progesteroneonly model. Secondly, we validate the bivariate model using successful insemination times as proxies for ovulation.

PREDICTION OF FOLLICULAR STATE BY ACTIVITY-ONLY MODEL

We estimated parameters for two separate univariate models for each variable using the methods described in Section 3. The Viterbi algorithm was applied to find the optimal state sequences. We then looked at how well the state sequence from the activity model aligned with that from the progesterone model, treating the state sequence from progesterone as validation data. According to biological knowledge, standing heat should occur near the beginning of the follicular stage [START_REF] Ball | Reproduction in Cattle[END_REF]. We considered a predicted standing heat from the activity model that was bounded (in time) by a predicted follicular state from the progesterone model to be a correct prediction. We do not claim that the state sequence from progesterone is true, but a closer reconstruction of the biological states of the cow.

Two possible prediction errors can occur in this framework; false positive (prediction of standing heat outside a follicular state) and false negative (failure to predict standing heat within a follicular stage). An example of some correctly predicted follicular stages is shown in Figure 6. The results are summarized in Table 2.

OESTRUS DETECTION

We compared the bivariate model and the activity-only univariate model as online systems by using information about successful artificial inseminations and applying crossvalidation. Cross-validation was used as follows. We iterated through a set of 58 cows with successful pregnancies, excluding one cow's full data sequence and estimating the HSMM parameters using the remaining data. We then applied the Viterbi algorithm to the left-out cows' sequence up until the time of successful insemination. Finally, we compared the time of the most recent predicted (if any) oestrus event by the models to the artificial insemination event. Artificial insemination is thought to be best applied within 24 hours following oestrus so we expect a positive difference between AI and the last detected oestrus hour.

A kernel density estimate of the bivariate models time differences can be viewed in Figure 7 (left) and a comparison of the two models in Figure 7 (right). The results for each model are summarized in Table 3 where we can see that the bivariate model had slightly We have limited our analysis to insemination events that led to a pregnancy. The majority of inseminations do not lead to a pregnancy and we have excluded these failed inseminations from the analysis. Inseminations may fail due to incorrect timing, technical problems with the process or cow infertility. We can only be certain an oestrus has occurred if a cow becomes pregnant and hence have used these events to validate our model.

DISCUSSION

We have investigated and compared univariate and bivariate HSMMs for two different variables related to cattle reproduction. Approximately 70% of follicular states were predicted by the activity-only model. This is less accurate than results reported by [START_REF] Roelofs | Various Behavioral Signs of Estrous and their Relationship with Time of Ovulation in Dairy Cattle[END_REF]. However those detections rates were for behavioral oestrus events, whereas our data also contains silent oestruses. The false positive rate of 14% is in line with their and other results [START_REF] Jonsson | Oestrus Detection in Dairy Cows Using Likelihood Ratio Tests[END_REF]. The results from our bivariate model indicate that the introduction of progesterone can reduce the number of false positives. The bivariate model was slightly better in terms of predicting the timing of oestrus. In 9 out of 58 cases the activity-only model failed to detect the relevant oestrus.

Progesterone has a financial cost associated with each measurement taken due to the chemicals used, while activity recordings are essentially free after initial overhead. Whether the reduction in false positives and negatives and slightly improved timing justifies these costs is a matter for further investigation. It would also be useful to investigate how the frequency of progesterone measurements affects accuracy. Another possibility is using activity as an indicator as to when to begin progesterone monitoring. For example, progesterone measurements are not useful during post-partum anoestrus. We found an activity-only model could often identify when a cow has left post-partum anoestrus (although we have not formally reported the results here); this could be a trigger to begin progesterone monitoring.

Figure 2 .

 2 Figure 2. Left: Histogram of the raw activity data across the entire population. Right: Histogram of the raw progesterone data across the entire population.

Figure 3 .

 3 Figure 3. Left column: Histogram of log(activity + 1) for two individual cows. Centre: One-step-ahead forecast residuals from a Holt-Winters model for the same cows (the activity index). Right: The activity index in grey with a 24-hour rolling mean overlayed in black, the dashed vertical line is when artificial insemination occurred.

Figure 4 .

 4 Figure 4. Top: The underlying states for a univariate activity model. Bottom: The underlying states for a univariate progesterone model.

Figure 5 .

 5 Figure 5. The hidden states used for the bivariate HSMM. The use of both variables allows us to break the follicular stage up into three different sub-stages.

Figure 6 .

 6 Figure6. Top: The activity index against time since calving. Bottom: The progesterone concentration over the same time and cow. Predicted states by the respective univariate models are shown on the horizontal shaded bars. This is an example where an activity-only model has correctly found oestrus events occur within a follicular stage.

Figure 7 .

 7 Figure 7. Left: Kernel density estimate of the time difference between successful AI and the most recently detected oestrus state by the bivariate model. Right: The univariate model against the bivariate model for time to successful AI from detected oestrus. Note the two large values from the univariate model, certainly false negatives.

Table 1 .

 1 Table displaying the ovarian stages of a cow and the related variables. Ovulation occurs within the first day of metoestrus.

	Time	Stage	Metaphase	Activity	Progesterone
	Days 0-13 Days 13-15	dioestrus pro-oestrus	luteal	normal normal	high
	Day 16 Days 17-20	oestrus metoestrus	follicular	high normal	low

Table 2 .

 2 Success and error rates for prediction of a follicular stage by an univariate activity HSMM. TP = True positive, FP = False positive, FN = False negative. The activity-only model failed to detect any oestrus in 7 of the 58 cases as well as the two obvious false negatives that can be seen in the graph, a false negative rate of at least 0.16. The bivariate model essentially uses the activity measure to narrow the window of oestrus with the follicular stage, so it is not surprising the models are comparable in terms of oestrus timing (when the activity model actually finds the oestrus).

		Sensitivity TP TP+FN	Error rate FP TP+FP
	Count	97/137	16/113
	Proportion	0.71	0.14
	lower times indicating better performance.

Table 3 .

 3 Cumulative table of the time between a predicted oestrus event and an actual successful artificial insemination for the bivariate and univariate (activity-only) models. Predictions that occurred less than one day before the actual insemination can be considered successful, predictions between one and two days are borderline. The bivariate model predicts more oestruses within one to two days of an insemination.

	D a y s	1	2	3	4	5	6	7	8
	Time to AI < Days (bivariate)	0.84	0.89	0.93	0.95	0.96	0.96	0.96	1.00
	(activity only)	0.68	0.73	0.77	0.77	0.80	0.84	0.84	0.88

ACKNOWLEDGEMENTS

This work was supported by the BIOSENS project funded by the Danish Ministry of Food, Agriculture and Fisheries and the Danish Cattle Industry via Finance Committee Cattle, and by The Danish Council for Strategic Research through the AUREGAB project.

[Accepted July 2009.]


Moez Rhimi
email: rhimimoez@yahoo.fr

Goran Bajic

Rimeh Ilhammami

Samira Boudebbouze

Emmanuelle Maguin

Richard Haser

Nushin Aghajari
email: n.aghajari@ibcp.fr

The acid-tolerant L-arabinose isomerase from the mesophilic Shewanella sp. ANA-3 is highly active at low temperatures

Keywords: L-arabinose isomerase, acid-tolerant, cold-active, D-tagatose

Background: L-arabinose isomerases catalyse the isomerization of L-arabinose into L-ribulose at insight biological systems. At industrial scale of this enzyme is used for the bioconversion of D-galactose into D-tagatose which has many applications in pharmaceutical and agro-food industries. The isomerization reaction is thermodynamically equilibrated, and therefore the bioconversion rates is shifted towards tagatose when the temperature is increased. Moreover, to prevent secondary reactions it will be of interest to operate at low pH. The profitability of this D-tagatose production process is mainly related to the use of lactose as cheaper raw material. In many dairy products it will be interesting to produce D-tagatose during storage. This requires an efficient L-arabinose isomerase acting at low temperature and pH values.

Results:

The gene encoding the L-arabinose isomerase from Shewanella sp. ANA-3 was cloned and overexpressed in Escherichia coli. The purified protein has a tetrameric arrangement composed by four identical 55 kDa subunits. The biochemical characterization of this enzyme showed that it was distinguishable by its maximal activity at low temperatures comprised between 15-35°C. Interestingly, this biocatalyst preserves more than 85% of its activity in a broad range of temperatures from 4.0 to 45°C. Shewanella sp. ANA-3 L-arabinose isomerase was also optimally active at pH 5.5-6.5 and maintained over 80% of its activity at large pH values from 4.0 to 8.5. Furthermore, this enzyme exhibited a weak requirement for metallic ions for its activity evaluated at 0.6 mM Mn 2+ . Stability studies showed that this protein is highly stable mainly at low temperature and pH values. Remarkably, T268K mutation clearly enhances the enzyme stability at low pH values. Use of this L-arabinose isomerase for D-tagatose production allows the achievement of attractive bioconversion rates of 16% at 4°C and 34% at 35°C. Conclusions: Here we reported the purification and the biochemical characterization of the novel Shewanella sp. ANA-3 L-arabinose isomerase. Determination of the biochemical properties demonstrated that this enzyme was highly active at low temperatures. The generated T268K mutant displays an increase of the enzyme stability essentially at low pH. These features seem to be very attractive for the bioconversion of D-galactose into Dtagatose at low temperature which is very interesting from industrial point of view.

Background

L-Arabinose isomerase, hereafter L-AI, (EC 5.3.1.4) converts L-arabinose into L-ribulose within living systems [START_REF] Kim | Improvement of tagatose conversion rate by genetic evolution of thermostable galactose isomerase[END_REF]. This enzyme is also known as D-galactose isomerase due to its capacity to isomerize D-galactose into D-tagatose in vitro [START_REF] Izumori | Pentose metabolism in Mycobacterium smegmatis: comparison of L-arabinose isomerases induced by Larabinose and D-galactose[END_REF][START_REF] Roh | Bioconversion of D-galactose into Dtagatose by expression of L-arabinose isomerase[END_REF]. Over the past years, rare sugars have made their way to the market and the demands have been increased as in agro-food industry rare natural sugars such as D-tagatose have become of a great interest as sweeteners [START_REF] Levin | Tagatose, the new GRAS sweetener and health product[END_REF]. Indeed this natural sugar has the taste and natural properties of sucrose but it is an antihyperglycemiant factor known to reduce symptoms associated with type 2 diabetes [START_REF] Lu | Tagatose, a new antidiabetic and obesity control drug[END_REF][START_REF] Saunders | Effects of acute and repeated oral doses of D-tagatose on plasma uric acid in normal and diabetic humans[END_REF]. D-Tagatose has not only a very low caloric value of 1.5 kcal/g (only 30% of the energy content of sucrose), but has also been reported to promote weight loss [START_REF] Lu | Tagatose, a new antidiabetic and obesity control drug[END_REF]. All these properties are very interesting for industrials, yet the use of this sugar is quite limited due to its cost [START_REF] Rhimi | Co-expression of Larabinose isomerase and D-glucose isomerase in E. coli and development of an efficient process producing simultaneously Dtagatose and D-fructose[END_REF]. In fact, chemical methods used at present to produce D-tagatose are very expensive. Therefore, biological methods using L- arabinose isomerase as a catalyst are being developed due to its low cost [START_REF] Levin | Sugar substitutes: their energy values, bulk characteristics, and potential health benefits[END_REF].

Many micro-organisms such as Arthrobacter, Lactobacillus, Mycobacterium, Klebsiella and Gluconobacter are able to convert D-galactitol into D-tagatose [START_REF] Manzoni | Biotransformation of D-galactitol to tagatose by acetic acid bacteria[END_REF]. However, this production procedure is not attractive due to the high cost of the raw material D-galactitol. That is why studies are now investigating a new bioprocess, which consists in converting D-galactose into D-tagatose using L-AIs. These studies have shown that the reaction equilibrium shifts towards D-tagatose production with increased temperatures, and that the majority of L-AIs are thermoactive biocatalysts requiring metal ions for their maximal thermostability and activity [START_REF] Roh | Bioconversion of D-galactose into Dtagatose by expression of L-arabinose isomerase[END_REF][START_REF] Kim | Cloning, expression and characterization of L-arabinose isomerase from Thermotoga neapolitana: bioconversion of D-galactose to D-tagatose using the enzyme[END_REF][START_REF] Jorgensen | Enzymatic conversion of Dgalactose to D-tagatose: heterologous expression and characterization of a thermostable L-arabinose isomerase from Thermoanaerobacter mathranii[END_REF][START_REF] Kim | Purification and characterization of an L-arabinose isomerase from an isolated strain of Geobacillus thermodenitrificans producing D-tagatose[END_REF][START_REF] Rhimi | Cloning, purification and biochemical characterization of metallic-ions independent and thermoactive L-arabinose isomerase from the Bacillus stearothermophilus US100 strain[END_REF][START_REF] Prabhu | Cloning and characterization of a novel L-arabinose isomerase from Bacillus licheniformis[END_REF]. Unfortunately, for industrial purposes high temperatures and alkaline pH result in the formation of side reactions generating undesirable by-products that should be removed (together with the metal ions) at the end of the process, thus increasing the production costs [START_REF] Rhimi | Co-expression of Larabinose isomerase and D-glucose isomerase in E. coli and development of an efficient process producing simultaneously Dtagatose and D-fructose[END_REF][START_REF] Lee | Characterization of a thermoacidophilic L-arabinose isomerase from Alicyclobacillus acidocaldarius: role of Lys-269 in pH optimum[END_REF]. For these reasons, highly thermoactive L-AIs deriving from thermophilic strains such as Bacillus stearothermophilus US100, Thermotoga maritima, Thermotoga neapolitana, Geobacillus thermodenitrificans and Thermus sp. were identified [START_REF] Kim | Purification and characterization of an L-arabinose isomerase from an isolated strain of Geobacillus thermodenitrificans producing D-tagatose[END_REF][START_REF] Rhimi | Cloning, purification and biochemical characterization of metallic-ions independent and thermoactive L-arabinose isomerase from the Bacillus stearothermophilus US100 strain[END_REF][START_REF] Lee | Characterization of a thermostable L-arabinose (d-galactose) isomerase from the hyperthermophilic eubacterium Thermotoga maritima[END_REF][START_REF] Hong | Production of D-tagatose at high temperatures using immobilized Escherichia coli cells expressing L-arabinose isomerase from Thermotoga neapolitana[END_REF][START_REF] Kim | Production of tagatose by a recombinant thermostable L-arabinose isomerase from Thermus sp. IM6501[END_REF]. Many industrial applications including bioconversion of D-galactose to D-tagatose during products storage, which must be done at low temperatures, require an L-AI that acts not only at acidic pH but also at low temperatures. To satisfy these requirements, screening of novel L-AIs has been done. As a result, many L-AIs were identified and characterized to be functional at relatively high pH and at high temperatures [START_REF] Kim | Cloning, expression and characterization of L-arabinose isomerase from Thermotoga neapolitana: bioconversion of D-galactose to D-tagatose using the enzyme[END_REF][START_REF] Jorgensen | Enzymatic conversion of Dgalactose to D-tagatose: heterologous expression and characterization of a thermostable L-arabinose isomerase from Thermoanaerobacter mathranii[END_REF][START_REF] Kim | Purification and characterization of an L-arabinose isomerase from an isolated strain of Geobacillus thermodenitrificans producing D-tagatose[END_REF][START_REF] Rhimi | Cloning, purification and biochemical characterization of metallic-ions independent and thermoactive L-arabinose isomerase from the Bacillus stearothermophilus US100 strain[END_REF][START_REF] Prabhu | Cloning and characterization of a novel L-arabinose isomerase from Bacillus licheniformis[END_REF][START_REF] Lee | Characterization of a thermoacidophilic L-arabinose isomerase from Alicyclobacillus acidocaldarius: role of Lys-269 in pH optimum[END_REF][START_REF] Lee | Characterization of a thermostable L-arabinose (d-galactose) isomerase from the hyperthermophilic eubacterium Thermotoga maritima[END_REF][START_REF] Hong | Production of D-tagatose at high temperatures using immobilized Escherichia coli cells expressing L-arabinose isomerase from Thermotoga neapolitana[END_REF][START_REF] Kim | Production of tagatose by a recombinant thermostable L-arabinose isomerase from Thermus sp. IM6501[END_REF]. Lately, four L-AIs from Alicyclobacillus acidocaldarius, Bacillus licheniformis, Lactobacillus plantarum NC8 and Lactobacillus sakei 23 K strains were described to be acid-tolerant. Molecular determinants for the catalytic efficiency were explored by homology modeling followed by site-directed mutagenesis [START_REF] Rhimi | Probing the essential catalytic residues and substrate affinity in the thermoactive Bacillus stearothermophilus US100 L-arabinose isomerase by site-directed mutagenesis[END_REF][START_REF] Prabhu | Probing the molecular determinant for the catalytic efficiency of L-arabinose isomerase from Bacillus licheniformis[END_REF]. Protein engineering was also used to generate enzymes acting at low temperatures, such as the N175H mutant from B. stearothermophilus US100 [START_REF] Lee | Characterization of a thermoacidophilic L-arabinose isomerase from Alicyclobacillus acidocaldarius: role of Lys-269 in pH optimum[END_REF][START_REF] Prabhu | Probing the molecular determinant for the catalytic efficiency of L-arabinose isomerase from Bacillus licheniformis[END_REF][START_REF] Chouayekh | Characterization of an L-arabinose isomerase from the Lactobacillus plantarum NC8 strain showing pronounced stability at acidic pH[END_REF][START_REF] Rhimi | Rational design of Bacillus stearothermophilus US100 L-arabinose isomerase: potential applications for D-tagatose production[END_REF] and at acidic pH (mutant Q268K) [START_REF] Rhimi | Rational design of Bacillus stearothermophilus US100 L-arabinose isomerase: potential applications for D-tagatose production[END_REF]. On the basis of these two mutants a third one was developed (Q268K/ N175K) [START_REF] Rhimi | Rational design of Bacillus stearothermophilus US100 L-arabinose isomerase: potential applications for D-tagatose production[END_REF]. This double mutant is highly active at both acidic pH and low temperatures making it very promising for industrial D-tagatose production. Nevertheless, it remains of interest to set out a stable L-AI active at even lower temperatures.

Here we report the cloning, over-expression, purification and biochemical characterization of an L-AI isolated from Shewanella sp. ANA-3. The activity and stability of the enzyme at low temperatures is discussed in comparison with previously described L-AIs. Site-directed mutagenesis has been performed to generate the T268K mutant and explore its biochemical properties.

Results and discussion

Gene, cloning and over-expression

A DNA fragment of approximately 1.5 kb was obtained by polymerase chain reaction (PCR) amplification using Shewanella sp. ANA-3 chromosomal DNA as template and the two oligonucleotides. DNA was ligated into the pGEMT-Easy vector and then transferred into the Escherichia coli BL21 (DE3) host strain. Numerous white colonies were observed and subsequently analyzed with the appropriate restriction endonucleases; a recombinant clone carrying the araA gene under the control of T7 promoter was obtained (pMR37). Then, the araA gene was sub-cloned in pET-15b under control of the T7 promoter leading to the pMR38 plasmid. Monitoring of liquid culture (BL21/pMR38 strain), followed by an L- AI activity test in the crude extract, showed an activity of 26 U/mg, confirming the molecular cloning and the expression of Shewanella sp. ANA-3 araA under the control of the T7 promoter. This latter gene encodes a 500 amino-acid residues protein with a calculated molecular mass of 55.3 kDa. Inspection of the amino-acid sequence multiple alignment of this protein with several other L-AIs indicated that Shewanella sp. ANA-3 L-AI displays the highest sequence identity (54.5%) with L-AI from the acidophilic strain A. acidocaldarius. Other enzymes such as those from B. licheniformis, L. plantarum NC8 and L. sakei 23 K, share 50.5%, 40.7% and 38.7% sequence identity, with Shewanella sp. ANA-3 L- AI, respectively (Figure 1). Sequence alignments with the thermophilic G. stearothermophilus, B. stearothermophilus US100, B. thermodenitrificans, T. neapolitana and T. maritima showed significative identity of 54.2%, 54.4%, 54.3%, 47.7% and 48.3%, respectively (Figure 1).

Purification of the Shewanella sp. ANA-3 L-AI

After an over-night liquid cell culture of the BL21/ pMR38 strain, the protein crude extract was subjected to fractioned precipitation with ammonium sulfate (Figure 2A). The subsequent step consisted in loading the resolubilized fractions onto a HiTrap™Chelating Immobilized Metal Affinity Chromatography (IMAC) column. A high purity enzyme extract was achieved with size exclusion chromatography (Figure 2A). The single elution peak obtained has an apparent molecular mass of 230 kDa (Figure 2B). Electrophoresis under reducing conditions (12% SDS-PAGE; sodium dodecyl sulfatepolyacrylamide gel electrophoresis) revealed a homogenous band with a molecular mass of about 55 kDa (Figure 2A), suggesting that Shewanella sp. ANA-3 L-AI is a tetramer composed by four monomers of 55-kDa. Previous studies showed that all L-AIs have a similar, homotetrameric quaternary structure [START_REF] Jorgensen | Enzymatic conversion of Dgalactose to D-tagatose: heterologous expression and characterization of a thermostable L-arabinose isomerase from Thermoanaerobacter mathranii[END_REF][START_REF] Rhimi | Cloning, purification and biochemical characterization of metallic-ions independent and thermoactive L-arabinose isomerase from the Bacillus stearothermophilus US100 strain[END_REF][START_REF] Prabhu | Cloning and characterization of a novel L-arabinose isomerase from Bacillus licheniformis[END_REF][START_REF] Lee | Characterization of a thermoacidophilic L-arabinose isomerase from Alicyclobacillus acidocaldarius: role of Lys-269 in pH optimum[END_REF][START_REF] Kim | Production of tagatose by a recombinant thermostable L-arabinose isomerase from Thermus sp. IM6501[END_REF] the exception being E. coli L-AI which forms a homohexamer [START_REF] Manjasetty | Crystal structure of Escherichia coli l-arabinose isomerase (ECAI), the putative target of biological tagatose production[END_REF]. Presently, the only known L-AI 3D-structure is that of E. coli [START_REF] Manjasetty | Crystal structure of Escherichia coli l-arabinose isomerase (ECAI), the putative target of biological tagatose production[END_REF], and further structural studies would be of high interest in order to explore the importance of the quaternary structures in D-galactose binding and/or catalysis. Shewanella sp. ANA-3L-AI presented a relatively high specific activity of 164 U/mg making it more active than most L-AIs from mesophilic bacterial origin (Table 1). Only few of the currently characterized L-AIs exhibit higher specific activity, mainly those from thermophilic bacteria like B. stearothermophilus US100 (185 U/mg) and T. neapolitana (119 U/mg) or lactic acid bacteria as L. sakei (218 U/mg) [START_REF] Kim | Cloning, expression and characterization of L-arabinose isomerase from Thermotoga neapolitana: bioconversion of D-galactose to D-tagatose using the enzyme[END_REF][START_REF] Rhimi | Cloning, purification and biochemical characterization of metallic-ions independent and thermoactive L-arabinose isomerase from the Bacillus stearothermophilus US100 strain[END_REF][START_REF] Rhimi | The acid tolerant L-arabinose isomerase from the food grade Lactobacillus sakei 23 K is an attractive D-tagatose producer[END_REF].

Effects of temperature, pH and metallic ions on L-AI activity

Shewanella sp. ANA-3 L-AI was optimally active within a large range of temperatures from 15 to 35°C, and highly active at low temperatures preserving nearly 90% to 95% of its optimal activity at 4 to 10°C, respectively. At 40 to 45°C, the protein retained 90% to 85% of its activity, respectively, whereas at higher temperatures (> 45°C) the relative activity rapidly decreased. This pronounced activity at low temperatures distinguishes it from all previously described L-AIs (Table 1), and another interesting feature is that more than 85% of its maximal activity is maintained in a wide range of temperatures from 4 to 45°C. This behavior indeed promotes the use of this enzyme in agro-food industry. Study of the effect of pH on this L-AI demonstrates that this protein was highly active at pH 5.5-6.5. Furthermore, the enzyme retains more than 80% of its activity from pH 4.0 to 8.5. Remarkably, this protein preserves more than 60% of its relative activity at pH 3.0 and 3.5, whereas the majority of the earlier reported L-AIs have optimal pH values of 7.0-8.0 (Table 1). For instance, all L-AIs from Bacillus genera possess a rather narrow range of optimal pH from 7.0 to 8.5. Concerning L-AIs from lactic acid bacteria, the optimum pH values are more acidic spanning from 5.0 to 7.5, and the most acid tolerant candidate being L. sakei 23 K L-AI functioned at pH 5.0-7.0 [START_REF] Rhimi | The acid tolerant L-arabinose isomerase from the food grade Lactobacillus sakei 23 K is an attractive D-tagatose producer[END_REF]. In contrast, some L-AIs act at highly alkaline conditions; such is the case of Anoxybacillus flavithermus L-AI that has an optimal pH range from 9.5 to 10.5 [START_REF] Li | Identification and characterization of a novel Larabinose isomerase from Anoxybacillus flavithermus useful in D-tagatose production[END_REF]. However, the acid tolerance is important mainly for industrial purposes. In fact, isomerization of D-galactose under acid pH prevents the formation of unwanted products [START_REF] Lee | Characterization of a thermostable L-arabinose (d-galactose) isomerase from the hyperthermophilic eubacterium Thermotoga maritima[END_REF][START_REF] Lee | Distinct metal dependence for catalytic and structural functions in the L-arabinose isomerases from the mesophilic Bacillus halodurans and the thermophilic Geobacillus stearothermophilus[END_REF].

Treatment of the enzyme with ethylenediaminetetraacetic acid (EDTA) showed a 74% decrease in the activity. Amongst several metal ions only Mn 2+ restored the enzyme activity, which exceeded that observed for the metal-depleted enzyme by 1.6 fold (Figure 3A), at a concentration of only 0.6 mM Mn 2+ (Figure 3B). Analysis of the EDTA-treated enzyme by flame atomic absorption spectrometry demonstrates the presence of only manganese; thus highlighting its tight interaction with protein. This explains the residual activity (26%) obtained after enzyme EDTA treatment. Comparatively to all reported L-AIs, the present enzyme showed a weak requirement for metal ions for its optimal activity (Table 1). This feature is highly advantageous for its potential use in agrofood industry where many metal ions are not authorized, and therefore if used should be discarded from the final products, thereby increasing the process costs [START_REF] Rhimi | Co-expression of Larabinose isomerase and D-glucose isomerase in E. coli and development of an efficient process producing simultaneously Dtagatose and D-fructose[END_REF].

Shewanella sp. ANA-3 thermal and pH stability studies

Enzyme stability is a crucial criterion for industrial applications. Our study showed that without addition of metal ions this enzyme was fully stable at temperatures below 35°C. At 40°C the enzyme upheld more than 80% of its activity after two hours, and at 45 and 50°C the enzyme stability decreased with half-lives of 90 and 40 min, respectively (Figure 4A).

After addition of 0.6 mM Mn 2+ the enzyme was completely stable for 1 hour at 40°C and preserved 90% of its activity after two hours of incubation. Addition of manganese clearly improved the thermal stability at 45 and 50°C, the enzyme retained more than 50% of its initial activity after two hours (Figure 4B). These results give evidence that Mn 2+ ions play a crucial role in the stabilization of the enzyme at high temperatures. This phenomenon is well reported for L-AIs, requiring metal ions for their stability (Table 1). Study of the enzyme stability at different pH values demonstrates that this biocatalyst has high half-lives in a wide pH range from 4.5 to 7.0 (Table 2). Furthermore, at pH 3.0 to 4.0 and 7.5 to 8.5 the half-lives dropped slightly. Variations of temperature from 4 to 30°C were not significant, but highlighting the stability of the Shewanella enzyme in a large range of temperatures and pH. Using sequence alignment, Lee et al. (2005) identified Lys269 in Alicyclobacillus acidocaldarius L-AI as an important determinant for enzyme's acid tolerance [START_REF] Lee | Characterization of a thermoacidophilic L-arabinose isomerase from Alicyclobacillus acidocaldarius: role of Lys-269 in pH optimum[END_REF]. This hypothesis was confirmed by generating mutants that shifted the optimum pH by one unit towards more alkaline pH [START_REF] Lee | Characterization of a thermoacidophilic L-arabinose isomerase from Alicyclobacillus acidocaldarius: role of Lys-269 in pH optimum[END_REF]. Guided by these studies, we have explored the influence of the non conserved lysine residue at position 268 on the biochemical properties of the Shewanella sp. ANA-3 enzyme (Figure 1). By site-directed mutagenesis we replaced the corresponding threonine by a lysine, and our data indicate that this mutation does not affect the enzymatic activity at different temperatures. As for the relative activity profile at different pH values, it seems that the T268K mutation increases the enzyme activity mostly at highly acidic pH. Remarkably, the wild-type Shewanella L-AI upholds more than 60% of its activity at pH 3.0 and the T268K mutant is even more acid-tolerant, retaining 80% of its activity at the same pH. Also, analysis of this mutant stability at different pH and different temperatures demonstrates that the T268 K L-AI mutant displays higher half-lives with low pH (3.0-4.5) when compared to the wild-type enzyme (Table 2). Similar half-lives were obtained for both mutant and wild-type Shewanella sp. ANA-3 L-AIs at pH's from 5.0 to 7.0. Remarkably, T268K L-AI mutant displays a maximal stability at low pH value of 5.5 and high value of 7.0. Thus, this mutant is not only more acid-tolerant than the wildtype enzyme, but also more stable under acidic pH and low temperatures. This takes more importance when we consider that an attractive enzyme, from industrial point of view, should be stable and efficient. Such data confirm the essential role of the K268 residue on both acidtolerance and pH stability of L-AIs.

Kinetic characterization and D-galactose bioconversion

Determination of the kinetic parameters of the Shewanella sp. ANA-3 L-AI showed that it has K m values of 33.7 and 52.1 mM for L-arabinose and D-galactose, respectively. Furthermore the V max was 164 mmole/s/ mg for L-arabinose and 44 mmole/s/mg for D-galactose. To further explore the efficiency of this enzyme to produce D-tagatose, we carried out D-galactose bioconversion at different temperatures. As shown in Table 3 the highest bioconversion rate (34%) was obtained at 35°C. At lower temperatures including 30, 20, 10 and 4°C maximum bioconversion levels of 31%, 25%, 21% and 16% were obtained, respectively. These data underline the catalytic efficiency of this enzyme to produce D-taga- tose at low temperature. In addition, analysis of Table 3 demonstrates that bioconversion rates are improved when temperature is increased, thus confirming previously reported results concerning the relationship between temperature and isomerization reaction equilibrium [START_REF] Rhimi | Cloning, purification and biochemical characterization of metallic-ions independent and thermoactive L-arabinose isomerase from the Bacillus stearothermophilus US100 strain[END_REF][START_REF] Lee | Characterization of a thermoacidophilic L-arabinose isomerase from Alicyclobacillus acidocaldarius: role of Lys-269 in pH optimum[END_REF][START_REF] Lee | Characterization of a thermostable L-arabinose (d-galactose) isomerase from the hyperthermophilic eubacterium Thermotoga maritima[END_REF].

Such enzyme features are very attractive from an industrial point of view. In fact, until now the majority of characterized L-AIs are thermoactive and not acid-tolerant. Generating a new acid-tolerant L-AI being highly active at low temperatures would allow the efficient bioconversion of D-galactose into D-tagatose in dairy products during fermentation and/or storage.

Conclusions

The L-AI from psychrotolerant Shewanella sp. ANA-3 is highly active at a wide pH range and converts efficiently D- galactose into D-tagatose even at low temperatures. These biochemical properties are of interest for industrial purposes. Indeed, isomerization of D-galactose at acid pH prevents the formation of unwanted side-products and high stability at low temperatures is convenient for food industry applications. The T268K mutation did not change enzyme behavior regarding temperature but clearly enhanced the enzyme activity and stability at low pH values.

Further work will target the overexpression of this particular enzyme in a food-grade system such as lactic acid bacteria. Given the attractive properties of the enzyme it would be of high interest to employ it as a model for optimization of the bioconversion of D-galac- tose into D-tagatose in milk and dairy products such as yoghurts. In order to improve the suitability of this biocatalyst for industrial applications it is important to increase its affinity to D-galactose. For this purpose crystallization followed by the determination of the threedimensional structure of this protein in complex with substrates by X-ray crystallography in combination with protein engineering are underway. Such data will allow improving the efficiency of other food-grade L-AIs.

Methods

Bacterial strains, plasmids and media Shewanella sp. ANA-3 strain was routinely grown as reported by Saltikov et al. [START_REF] Saltikov | The ars detoxification system is advantageous but not required for As(V) respiration by the genetically tractable Shewanella species strain ANA-3[END_REF]. Escherichia coli BL21 (DE3) were used in this study as host strains. Culture of different E. coli strains was done in Luria Bertani (LB) medium. These media were supplemented, when necessary, with ampicillin (100 μg/ml) and IPTG (Isopropyl β-D-thiogalactopyranoside) at 160 μg/ml. The pGEMT-Easy (Promega) and pET-15b plasmids (Novagen) were used according to the manufacturer's instructions.

DNA manipulation and PCR

Genomic DNA was prepared as reported by Johnson [START_REF] Johnson | Genetic characterization[END_REF]. Preparation of plasmid DNA, digestion with restriction endonucleases and separation of fragments by agarose gel electrophoresis were performed as described by Sambrook et al. [START_REF] Sambrook | Molecular cloning: a laboratory manual[END_REF]. PCRs were carried out using Gene Amp ® PCR System 9700 (Applied Biosystems). The amplification reaction mixtures (100 μl) contained High Fidelity Taq amplification buffer, 20 pmol of each primer, 100 ng of DNA template, and 10 units of High-Fidelity Taq DNA polymerase (Invitrogen). The cycling parameters were 94°C for 5 min, followed by 40 cycles of 94°C for 30 s, 55°C for 60 s, and 72°C for 120 s. The amplified fragment containing the Shewanella sp. ANA-3 araA gene was sequenced using an automated DNA sequencer (Applied Biosystems).

Cloning and over-expression of the araA Shewanella sp.

ANA-3 gene

Aiming the amplification of the Shewanella sp. ANA-3 araA gene, we designed two primers based on the complete genome available in the NCBI databank. Oligonucleotide sequences were F-araA 5' GATAAGCACTA TTTGCGTAAGCATG 3' and R-araA 5' TTATTGCTGC GCAGAATTATAAACCTCTTGCTAAACCG 3' . Chromosomal DNA isolated from Shewanella sp. ANA-3 strain was used as a template and High-Fidelity Taq DNA polymerase for amplification. The PCR product was purified using QIAquick Gel Extraction Kit (QIAGEN ® ) by following the manufacturer's instructions. Subsequently, the resulting fragments were cloned in pGEMT-Easy Vector (Promega) and transformed into E. coli BL21 (Invitrogen) competent cells. Recombinant white clones were selected on LB agar medium added with ampicillin (100 μg/ml), IPTG (160 μg/ml) and X-gal (100 μg/ml). Then, the araA gene was subcloned in the pET-15b (Novagen) under the control of the T7 promoter with six histidines at the N-terminal end of the protein. The T268K L-AI mutant was generated by using a site-directed mutagenesis Kit purchased from Stratagene (used as recommendation by supplier) and oligonucleotides T268Kd 5' AACCTCA CTGGCATGAAAGGATTACCCGGACTG 3' and T268Kr 5' CAGTCCGGGTAATCCTTTCATGCCAGTGAGGTT 3' .

Amino acid sequence alignment

The multiple alignment of the L-AI amino-acid sequence was done using the program ClustalW [START_REF] Thompson | CLUSTAL W: improving the sensitivity of progressive multiple sequence alignment through sequence weighting, position-specific gap penalties and weight matrix choice[END_REF] and the figure rendering was done using the ESPript sequence analysis server [START_REF] Gouet | ESPript/ENDscript: Extracting and rendering sequence and 3D information from atomic structures of proteins[END_REF].

Crude cell-lysate preparation and enzyme purification

The recombinant Escherichia coli strains, overexpressing Shewanella sp. ANA-3 L-AI, were grown in LB medium with ampicillin (100 μg/ml). From this culture a second one was launched under the same conditions at 200 rpm and 37°C. The culture was supplemented with IPTG (0.5 mM) at OD 600 nm of 1.0 and subsequently incubated at 37°C and 220 rpm for over-night. Cells were harvested by centrifugation (6500 × g for 10 min at 4°C) and the pellets were re-suspended in 100 mM MES buffer at pH 5.5 supplemented with one Complete™ protease inhibitor cocktail tablet (Roche ® ). Cell disruption was carried out, twice, using a French Press at 1000 psi. For the purification, crude cell extract from the E. coli BL21 (DE3) strain was centrifuged (30 000 × g, for 30 min at 4°C). Proteins precipitated between 60 and 70% ammonium sulphate saturation, suspended in 100 mM MES buffer pH 5.5, concentrated and desalted in centrifugal micro-concentrators (Amicon, Inc) with a 30 kDa cut-off membrane. Then, an affinity chromatography was done using IMAC column on the ÄKTA Purifier FPLC system (Amersham Pharmacia Biotech) equilibrated with 100 mM MES buffer pH 5.5. Proteins were eluted at a flow rate of 3 ml/min using a linear imidazole gradient ranging from 0 to 500 mM in the same MES buffer with imidazole (500 mM). Fractions containing L-AI activity were pooled and the final purification phase was performed using size exclusion chromatography on a Sephacryl S-200 column (GE-Healthcare).

Protein quantification, electrophoresis and molecular mass determination

Protein concentration was determined using Bradford's method with bovine serum albumin as standard [START_REF] Bradford | A rapid and sensitive method for the quantification of microgram quantities of protein utilizing the principle of protein-dye binding[END_REF].

The protein samples were separated in 12% SDS-PAGE gel according to the Laemmli method, and bands were visualized by Coomassie brilliant blue R-250 (Biorad) staining. The estimated molecular mass of the purified L-AI was determined by size exclusion chromatography.

Enzyme activity assays

L-AI activity was measured by determining the amount of formed L-ribulose or D-tagatose. Under standard conditions, the reaction mixture contained 100 μl of enzyme (L-AI) preparation at a 1 mg/ml, 50 mM of L-arabinose (D-galactose) and 100 mM MES buffer (pH 5.5) in a final volume of 1 ml. The reaction was stopped by heating the mixture at 99°C for 2 min. The amount of L-ribulose (or D-tagatose) was colorimetrically determined using the cystein-carbazol-sulphuric-acid method and the absorbance was measured at 560 nm [START_REF] Laemmli | Cleavage of structural proteins during the assembly of the head of bacteriophage T4[END_REF]. D-Tagatose production was also detected with high-performance ionic chromatography (HPLC) using a Polypore CA column (250 × 4.6 mm). The products were separated by isocratic elution with water at a flow rate of 0.3 ml/min and detected with a refractive index detector (SHIMADZU, Refractive Index Detector). Solutions of 5 g/l D-galactose and 5 g/l D-tagatose were used as standards.

One unit of L-AI activity was defined as the amount of enzyme catalyzing the formation of 1 μmol keto-sugar per min under the above-specified conditions.

Temperature, pH and thermostability profiles

The effect of temperature on the activity was studied by incubating the purified protein at temperatures from 4 to 55°C. The enzyme pH profile was obtained at pH values between 3.0 and 8.5 (3.0 to 5.0 with sodium acetate buffer, 5.5 to 7.0 with MES buffer and 7.5 to 9.0 with Bicine buffer). Determination of the enzyme stability as a function of temperature and pH was carried out by measuring the residual activity at periodic intervals.

Effect of metal ions on the enzymatic activity

Shewanella sp. ANA-3 L-AI was purified as described in the enzyme purification section but without addition of metal ions to the crude extract and dialyzed against 100 mM MES buffer (pH 5.5) containing 10 mM EDTA for 48 h at 4°C (buffer renewed twice). Subsequently, we twice dialyzed the EDTA-treated enzyme against 100 mM MES buffer (pH 5.5). Thereafter, the enzyme was pre-incubated during 10 min in solutions containing metal ions (CuSO 4 , CaSO 4 , MnSO 4 , FeSO 4 , BaSO 4 , CoSO 4 , FeSO 4 , MgSO 4 ) followed by L-AI activity assays under standard conditions. The optimal metal ion concentrations were monitored by incubating the enzyme samples at various concentrations followed by activity assays.

Flame atomic absorption spectrometry analysis

For sample preparation, each mineralized sample was employed at a final concentration evaluated to 10 mg/ ml. A Perkin-Elmer Analyst 200 atomic absorption (Norwalk, USA), equipped with a deuterium lamp background correction system, was used for metal binding quantification. Hollow cathode lamps (Perkin-Elmer) were used as primary radiation source. Analytical measurements were based on time average absorbance, under conditions recommended by the manufacturer.

Determination of kinetic parameters

Kinetic parameters were determined on the basis of Lineweaver-Burk plots. Assays were carried out in 100 mM MES buffer (pH 5.5), 0.6 mM Mn 2+ and 1 to 800 mM substrate (L-arabinose or D-galactose). Samples were incubated at 35°C and the amount of keto-sugar generated (L-ribulose or D-tagatose) was determined by the cysteine-carbazole-sulfuric acid method.

Bioconversion assays

Bioconversion of D-galactose into D-tagatose was monitored in a solution containing 100 mM MES buffer (pH 5.5), 0.6 mM Mn 2+ and 50 mM D-galactose. Kinetic conversion of D-galactose was investigated for 7 h at 4, 10, 20, 30 and 35°C. Samples were taken periodically and the concentration of generated D-tagatose was determined as described earlier and confirmed by HPLC as indicated in the enzyme assay section.

List of abbreviations

L-AI: L-arabinose isomerase; MES: 2(-N-morpholino)ethanesulfonic acid; EDTA: Ethylenediaminetetraacetic acid; HPLC: High-performance liquid chromatography; IPTG: Isopropyl β-D-thiogalactopyranoside; PCR: Polymerase Chain Reaction; RT: Retention Time; SDS-PAGE: Sodium Dodecyl Sulfate Polyacrylamide Gel Electrophoresis.

Figure 1

 1 Figure 1 Multiple sequence alignment. Invariant residues between sequences are typed red on a white background and residues conserved within each group are displayed as white letters on a red background. ECOL, E. coli (accession number PO8202); SHEW, Shewanella sp. ANA-3 (accession number A0KWX7); ACCE, A. cellulolytics ATCC 43068 (accession number GU188440); ALAC., Alicyclobacillus acidocaldarius (accession number AAY68209); BAST, B. stearothermophilus US100 (accession number CAI29261); BAHA, B. halodurans (accession number NP_242739); BALI, Bacillus licheniformis (accession number YP_080173.1); BASU, B. subtilis (accession number ACT82395); GEST, G. stearothermophilus (accession number ABY84698); GETH, G. thermodenitrificans (accession number Q5QKM7); LAFE L. fermentum (accession number ADJ94948); LAPL, L. plantarum (accession number Q88S84); LASA, Lactobacillus sakei 23 K (accession number YP_396468.1); THMA, T. maritima (accession number Q9WYB3); THEN, T. neapolitana (accession number Q8RMB9); THSP, Thermus sp. (accession number AY225311); THMH T. mathranii (accession number ADH61392).

Figure 2

 2 Figure 2 Electrophoretic and size exclusion FPLC analysis of the purified Shewanella enzyme. (A) SDS-PAGE showing the various steps in the purification procedure under reducing conditions. Lane M, protein markers (molecular masses in kDa); lane 1, total cell extract; lane 2, sample from lane 1 after ammonium sulphate precipitation; lane 3, sample from lane 2 after affinity chromatography; lane 4, purified recombinant L-AI after size exclusion chromatography. (B) Size exclusion HPLC analysis of the purified recombinant Shewanella sp. ANA-3 L-AI shows a single peak of 230 kDa (RT = 14.1 min), using protein markers of 669 kDa (RT = 8.755 min), 440 (RT = 13.165 min), 232 (RT = 14.02 min; RT, retention time), 140 (RT = 17.163 min) and 66 (RT = 21.438 min).

Figure 3

 3 Figure 3 Shewanella sp. ANA-3 L-AI activity. (A): Effect of divalent metal ions. (B): Effect of various concentrations of Mn 2+ (♦). "MD" indicates enzyme activity measured in a metal depleted solution. Activities at optimal metal ion concentrations were defined as 100%. Error bars represent the standard deviation from three independent experiments.

Figure 4

 4 Figure 4 Thermostability profiles of (A) the metal depleted wild-type Shewanella sp. ANA-3 L-AI and (B) in the presence of 0.6 mM Mn 2+ : (■) 35°C, (♦) 40°C, (▲) 45°C and (•) 50°C. The initial activity was defined as 100%. Error bars represent the standard deviation from three independent experiments.

Table 1

 1 Biochemical properties of previously reported L-AIs

	Organism	Topt (°C) pH	Metal ion	Km (mM) Vmax (mmole/s/mg) References
				Requirement (mM)
	Shewanella sp. ANA-3	15-35	5.5-6.5	Mn 2+ (0.6)	33.7*	164*	This study
	Acidothermus cellulolytics	75	7.5	Mn 2+ (1.0), Co 2+ (0.5)	NR	NR	[35]
	Alicyclobacter acidocaldarius	50	6.0	Mn 2+ (1.0), Co 2+ , (1.0) Mg 2+ (1.0) 48	35.5	[15]
	Anoxybacillus flavithermus	95	9.5-10.5 Ni 2+ (1.0)	78.5	94.3	[25]
	Bacillus stearothermophilus US100	80	7.5-8.0	Mn 2+ (1.0), Co 2+ (0.2)	28.57	185	[13]
	Bacillus halodurans	50	7.5-8.0	Mn 2+ (1.0)	36	33.1	[27]
	Bacillus licheniformis	50	7.5	Mn 2+ (1.0), Co 2+ (1.0)	369	NR	[14]
	Bacillus subtilis	32	7.5	Mn 2+ (1.0)	120	NR	[36]
	Escherichia coli	30	8.0	Mn 2+ (0.5)	60	NR	[37]
	Geobacillus stearothermophilus T6	70	7.0-7.5	Mn 2+ (1.0)	63	36.5	[1]
	Geobacillus thermodenitrificans	60	7.5	Mn 2+ (5.0), Co 2+ (3.0)	142	86	[12]
	Lactobacillus fermentum	65	6.5	Mn 2+ (1.0), Co 2+ (2.0)	29.9	24.3	[38]
	Lactobacillus plantarum NC8	60	7.5	Mn 2+ (1.0), Co 2+ (0.5)	43.42	40	[21]
	Lactobacillus sakei 23 K	30-40	5.0-7.0	Mg 2+ (0.8), Mn 2+ (0.8)	31.6	264	[24]
	Thermotoga maritima	90	7.0	Mn 2+ (5.0), Co 2+ (1.0)	63	41.3	[26]
	Thermotoga neapolitana	85	7.0-7.5	Mn 2+ (1.0), Co 2+ (1.0)	60	119	[10]
	Thermus sp. IM6501	60	8.5	Mn 2+ (5.0)	NR	NR	[18]
	Bacillus stearothermophilus IAM1101 65	6.0	MnP 2+P (1.0)	NR	NR	[39]

NR: Not Reported Top: Optimal temperature *: Substrate is L-arabinose

Table 2 pH

 2 -Stability of wild-type L-AI and its T268K mutant (values within brackets) at different temperatures without addition of metal ions

	Temperature (°C)

Table 3

 3 Kinetics of D-galactose bioconversion rates using the purified enzyme (1 mg/ml) at different temperatures

	Temperature (°C)	4	10	20	30	35
	Time (h)		Bioconversion rate (%)
	1.0	6.0 ± 0.1	6.0 ± 0.2	9.0 ± 0.1	10 ± 0.5	11 ± 0.1
	2.0	9.0 ± 0.2	8.0 ± 0.5	12 ± 0.5	14 ± 0.2	16 ± 0.3
	3.0	12 ± 0.1	12 ± 0.3	17 ± 0.8	19 ± 0.7	22 ± 0.4
	4.0	15 ± 0.4	18 ± 0.7	20 ± 0.4	25 ± 0.9	28 ± 0.7
	5.0	16 ± 0.8	20 ± 0.6	24 ± 0.6	29 ± 0.6	32 ± 0.6
	6.0	16 ± 0.2	21 ± 0.2	25 ± 0.5	30 ± 0.4	34 ± 0.1
	7.0	16 ± 0.3	21 ± 0.5	25 ± 0.2	31 ± 0.3	34 ± 0.8

Acknowledgements

We acknowledge Dr. Chad SALTIKOV from California Institute of Technology for generously providing us the Shewanella sp. ANA-3 strain. The authors would like to express their gratitude to Dr Michel JUY and Magali DEJOB in the Laboratory for BioCrystallography and Structural Biology of Therapeutic Targets and to Dr Annie CHABOUD for protein Production and Analysis facility, SFR BioSciences Lyon-Gerland UMS3444/US8 (France).

Author details

Competing interests

The authors declare that they have no competing interests.

Lifang C, Wanmeng M, Bo J: Thermostable L-arabinose isomerase from

Bacillus stearothermophilus IAM 11001 for


Mickaël Henry1
email: mickael.henry@avignon.inra.fr

Maxime Beguin2

Fabrice Requier3

Orianne Rollin1

Jean-François Odoux4

Pierrick Aupinel4

Jean Aptel1

Sylvie Tchamitchian1

Axel Decourtye5

A common Pesticide Decreases Foraging Success and Survival in Honey Bees

Non-lethal exposure of honey bees to thiamethoxam (neonicotinoid systemic pesticide) causes high mortality due to homing failure at levels that could put a colony at risk of collapse. Simulated exposure events on free-ranging foragers labeled with an RFID tag suggest that homing is impaired by thiamethoxam intoxication. These experiments offer new insights into the consequences of common neonicotinoid pesticides used worldwide.

One Sentence Summary: Non-lethal exposure to thiamethoxam leads to high mortality due to homing failure in honey bees.

Main Text:

Colony collapse disorder (CCD) is a recent, pervasive syndrome affecting honey bee (Apis mellifera) colonies in the Northern hemisphere, which is characterized by a sudden disappearance of honey bees from the hive (1). Multiple causes of CCD have been proposed, such as pesticides, pathogens, parasites, and natural habitat degradation (2-4). However, the relative contribution of those stressors in CCD events remains unknown. Some scientists and beekeepers suspect pesticides to hold a central place in colony weakening processes (1) or at least in interaction with other stressors (5, 6). In modern cereal farming systems, honey bees are readily exposed to pesticides because they rely heavily on common blooming crops, like oilseed rape (Brassica napus), maize (Zea mays) or sunflower (Helianthus annuus), that are now routinely treated against insect pests (3). Systemic pesticides, in particular, diffuse throughout all the tissues as plants grow-up, and eventually contaminate nectar and pollen (7). Foraging honey bees are

Manuscrit d'auteur / Author manuscript therefore directly exposed, but also the rest of the colony as returning foragers store or exchange contaminated material with hive conspecifics (7,8). Those exposure pathways are of important concern and pesticide manufacturers pay special attention to reduce non-intentional intoxications in field conditions. Pesticide authorization procedures now require running mortality surveys to ensure doses encountered in the field remain below lethal levels for honey bees.

However, a growing body of evidence shows that sublethal doses, i.e. doses that do not entail direct mortality, still have the potential to induce a variety of behavioral difficulties in foraging honey bees, such as memory and learning dysfunctions and alteration of navigational skills (9). Neonicotinoid pesticides used to protect crops against aphids and other sap-sucking insects are especially liable to provoke such behavioral troubles. They are highly potent and selective agonists of nicotinic acetylcholine receptors, which are important excitatory neurotransmitter receptors in insects (10,11). Effects of sublethal neonicotinoid exposures in honey bees may include abnormal foraging activity (12)(13)(14), reduced olfactory memory and learning performance (15)(16)(17) and possibly impaired orientation skills (18). Yet, the consequences of such behavioral difficulties on the fate of freeranging foragers and on colony dynamics are extremely difficult to assess in the field and remains poorly investigated.

In this study, we tested the hypothesis that a sublethal exposure to a neonicotinoid indirectly increases hive death rate through homing failure in foraging honey bees. We focused our attention on thiamethoxam, a recently marketed neonicotinoid substance (19), currently being authorized in an increasing number of countries worldwide for the protection of oilseed rape, maize and other blooming crops foraged by honey bees. We proceeded in two steps. First, we assessed mortality induced by homing failure (mhf) in exposed foragers. This was achieved by monitoring free-ranging honey bees using RFID tagging technology (20). Second, we assessed the extent to which mhf, in combination with natural forager mortality, may upset colony dynamics. For that purpose, mhf was introduced into a model of honey bee population dynamics (21).

We used a custom-made RFID device (20) to monitor the fate of 653 individual free-ranging foragers in the course of four separate treatment-vs.-control homing experiments [START_REF]Materials and methods are available as supporting material on Science Online[END_REF]. The study was conducted in an intensive cereal farming system of western France (Zone Atelier Plaine et Val de Sèvre research facility, CEBC) and in a suburban area in Avignon, southern France. To simulate intoxication events, foragers received a field-realistic, sublethal dose of thiamtethoxam (a real dose of 1.34 ng in a 20-µl sucrose solution) and were released away from their colony with a microchip glued on their thorax (Fig. 1A). RFID readers placed at the hive entrance (Fig. 1B) were set to detect on a continual basis tagged honey bees going through the entrance. Mortality due to post-exposure homing failure, mhf, was then derived from the proportion of non-returning foragers. To further discriminate mhf from other causes of homing failure in treated foragers, e.g. natural mortality, predation or handling stress, we simultaneously released equal numbers of control foragers -fed with an untreated sucrose solution. Hence, mhf was calculated as the proportion of non-retuning treated foragers relative to expectations given by the proportion of returning control foragers. Depending on the experiment, tagged honey bees where released up to 1 km away from their respective colony, i.e. at a distance usually covered by foragers during normal foraging flights [START_REF] Steffan---Dewenter | [END_REF]. Experiments were conducted on individuals from three different colonies [START_REF]Materials and methods are available as supporting material on Science Online[END_REF].

Our strategy was not to get an estimate of mhf per se. Instead we assessed its upper and lower bounds, depending on whether foragers were familiar or not with the foraging site they might get Manuscrit d'auteur / Author manuscript Manuscrit d'auteur / Author manuscript Manuscrit d'auteur / Author manuscript intoxicated in. Indeed, one might expect that foragers familiar with the pathway back to the colony are less prone to homing failure than unfamiliar foragers. Under field conditions, many foragers are probably familiar with the pathway back to the colony because they repeatedly forage on the same site (24). However, many others are unfamiliar too. Those include young honey bees at the onset of foraging, scouting honey bees that look for new food sources, and foragers newly recruited by scouting bees on the basis of the dance information (25). Most importantly, systemic pesticides like thiamethoxam are readily present in the nectar and pollen when flowering starts and receive the first visits of honey bees, hitherto unfamiliar with this newly available food source.

To account for individuals' past foraging experience, we conducted two distinct homing experiments. Experiment 1 simulated intoxication at a familiar foraging site and experiment 2 at a random site regarding past foraging experience. These experiments were assumed to return the lower and upper bounds of mhf, respectively. In experiment 1, we referred to as familiar foragers those foragers for which we could make sure they covered at least once the pathway from the release site back to the colony. For that purpose, we selectively captured foragers returning to the colony with pollen loads from a known location, and subsequently released them at that location. To ascertain pollen origin, we sowed beforehand a 1-ha field with scorpion weed Phacelia tanacetifolia, a highly attractive floral resource with bright blue pollen that is easily recognizable (26). Given that no other phacelia fields occurred in the area, we could ensure that phacelia-carrying honey bees came back from our experimental field. The colony was specifically placed 1 km away from the field for subsequent forager release (Fig. 2). In experiment 2, we used the nonphacelia pollen foragers. They were released in equal groups at six sites equally spaced on a 1-km circle around the colony (Fig. 2). Following that design, release sites were considered as random locations regarding the past experience of foragers.

Both experiment 1 and 2 evidenced substantial mortality due to post-exposure homing failure, mhf, with the proportion of treated foragers returning to the colony being significantly lower than that of control foragers (exact binomial tests, P=0.033 and P<0.001, respectively; Fig. 3, Table S1). Additionally, mhf was greater in treated foragers that tended to be unfamiliar with the foraging site, as indicated by their significantly lower homing proportions compared to familiar foragers (exact binomial tests, P<0.001). Experiments 1 and 2 returned mhf estimates of 0.102 and 0.316, respectively, potentially setting the lower and upper bounds for real mhf values. In other words, 10.2% to 31.6% of exposed honey bees would fail to return to their colony after foraging in a treated crop. For the sake of comparison, foragers live about 6.5 days, and therefore die at an average rate of 1/6.5 = 0.154 individual.day-1 (27). Therefore, the probability that a forager would die due to homing failure after visiting a treated crop (up to 0.316) may attain twice the probability this same forager has to die naturally that day (about 0.154).

Such an additional mortality might represent a heavy burden to bear for colonies exposed to treated crops in their environment. When implementing mhf into a honey bee population dynamics model (21), all the tested scenarios predicted a major deviation from the expected dynamic (Fig. 4). In our simulations, we considered the evolution of a typical colony during the first three months of a beekeeping season, encompassing the oilseed rape blooming period -April-May in our study area [START_REF]Materials and methods are available as supporting material on Science Online[END_REF]. At this time of the year, colonies emerge from the wintering period. Population size is rather low (<20,000 individuals) and gradually expands in order to rapidly increase food storage and ensure colony sustainability. The daily egg-laying rate of the queen is a critical parameter in this colony dynamic because it determines the daily egg-hatching rate, and in turn the rate at which honey bees working in the hive will be replaced as they become themselves foragers. We simulated three scenarios with realistic levels of egg-laying rate Manuscrit d'auteur / Author manuscript (28), namely a rate allowing for a normal colony development (Fig. 4A), a rate ensuring equilibrium population (Fig 4B ) and a slightly deficient rate forcing the population to stabilize at a lower size (Fig. 4C). In each case, we also computed the expected trends if most foragers (90%) were exposed to nectar of treated oilseed rape each day, and therefore had a natural mortality increased by a homing failure probability mhf. Regardless of the queens' egg-laying rate, populations from colonies exposed to the treated nectar would follow a marked decline during the blooming period, and would hardly recover afterwards ( Figs 4A-4C). When combined with natural forager mortality, mhf raised total forager death rate up to a point that could hardly be compensated for by the rate at which new foragers are recruited. In the worse scenarios, populations would fall down to 5,000 individuals, which is the lowest level one can usually observe in current beekeeping practices. With an exposure rate reduced to 50% of foragers exposed to treated nectar each day (Figs. 4D-F), the model still predicts a major deviation from normal conditions.

In an attempt to verify the applicability of these results to other contexts, we repeated two additional experiments with two different colonies (Fig S2 ,Tale S1). In experiment 3, we tested whether mhf was still significant when exposure occurred in the least challenging situation, i.e. in the direct vicinity of the colony and with honey bees familiar with the foraging site. Herein, we repeated experiment 1 with phacelia foragers captured from a beehive placed at the phacelia field margin, and released from inside the phacelia field, only 70 m away. Homing failure (mhf = 0.061, Fig. S2A, Table S1) was much reduced compared to experiment 1 (mhf = 0.102), but was still significant (exact binomial test, P=0.003). In experiment 4, we transposed experiment 2 into a different landscape. A beehive was placed in a suburban area in southern France, including a mosaic of mixed farming fields and orchards of moderate size. Foragers were released 1 km away at six equidistant sites. Homing failure (mhf = 0.098, Fig. S2B, Table S1) was significant as well (exact binomial test, P=0.029), but much smaller than in experiment 2 (mhf = 0.316).

Our study clearly demonstrates that exposure of foragers to non-lethal but commonly encountered concentrations of thiamethoxam can impact forager survival, with potential contributions to collapse risk. Furthermore, the extent to which exposures affect forager survival appears dependant on the landscape context and the prior knowledge of foragers about this landscape. Higher risks are observed when the homing task is more challenging. As a consequence, impact studies are likely to severely underestimate sublethal pesticide effects when they are conducted on honey bee colonies placed in the immediate proximity of treated crops. Finally, this study raises important issues concerning exposed solitary bee species, whose population dynamics are probably less resilient to forager disappearance than honey bee colonies. Fig. 1. Honey bee RFID monitoring equipment. (A) A pollen-forager honey bee fitted with a 3-mg RFID tag. (B) A hive entrance equipped with RFID readers for detecting returning marked foragers. Fig. 2. Study area and location of honey bee release sites relative to the colony hive in experiments 1 and 2. Fig. 3. Cumulative homing probability of foragers released 1km away from the hive. Temporal gaps denote the nighttime between the first and second days of release. Homing experiment 1 was carried out with foragers familiar with the release site (A) and experiment 2 with foragers released at random sites regarding their past experience (B). In both cases, treated honey bees that received a non-lethal dose of thiamethoxam returned to the hive in significantly lower proportions than control honey bees (Table S1). Fig. 4. Comparison of honey bee population dynamics between simulated colonies exposed to thiamethoxam (red lines) or not exposed (blue lines), following six demographic scenarios. L: queens daily laying rate (nb. of eggs per day). Exp: proportion of foragers exposed to treated crops during the day. The non-exposed colony follows either (A and D) a normal development trajectory (at L=2000), (B and E) an equilibrium dynamic (L=1800), or (C and F) a slightly declining trajectory (L=1600). Shaded areas delineate the exposure period (e.g. oilseed rape). Pairs of trajectories in exposed colonies were obtains with the lower and upper bounds of homing failure mortality (0.102 and 0.316), in order to delineate the best and worse estimates for population dynamics, respectively. Dotted lines extend the declining trajectory expected for a sustained exposure. [simulations derive from demographic models in (21)].

Supplementary Materials: Materials and Methods Figures S1-S2

Fig. 3

 3 FIG.2

Manuscrit d'auteur / Author manuscript Manuscrit d'auteur / Author manuscript Fig.4

Acknowledgments: This study was funded by the European Community program (797/2004) for French beekeeping coordinated by French Ministry of Agriculture (convention FranceAgriMer 11-45R). Special thanks go to L. Belzunces, J.L. Brunet, B. Vaissière, A. Maisonnasse, D. Fortini, Y. Le Conte and C. McDonnell for valuable help and corrections, as well as three anonymous reviewers for useful comments on the manuscript. We are grateful to M. Grijolot for allowing experiments to take place in his fields. RFID devices were designed by Tag Tracing Solutions Inc., Valence, France. Data are in the supporting online material.