The Insecticide Imidacloprid Causes Mortality of the Freshwater Amphipod Gammarus pulex by Interfering with Feeding Behavior

If an organism does not feed, it dies of starvation. Even though some insecticides which are used to control pests in agriculture can interfere with feeding behavior of insects and other invertebrates, the link from chemical exposure via affected feeding activity to impaired life history traits, such as survival, has not received much attention in ecotoxicology. One of these insecticides is the neonicotinoid imidacloprid, a neurotoxic substance acting specifically on the insect nervous system. We show that imidacloprid has the potential to indirectly cause lethality in aquatic invertebrate populations at low, sublethal concentrations by impairing movements and thus feeding. We investigated feeding activity, lipid content, immobility, and survival of the aquatic arthropod Gammarus pulex under exposure to imidacloprid. We performed experiments with 14 and 21 days duration, both including two treatments with two high, one day pulses of imidacloprid and one treatment with a low, constant concentration. Feeding of G. pulex as well as lipid content were significantly reduced under exposure to the low, constant imidacloprid concentration (15 µg/L). Organisms were not able to move and feed – and this caused high mortality after 14 days of constant exposure. In contrast, feeding and lipid content were not affected by repeated imidacloprid pulses. In these treatments, animals were mostly immobilized during the chemical pulses but did recover relatively fast after transfer to clean water. We also performed a starvation experiment without exposure to imidacloprid which showed that starvation alone does not explain the mortality in the constant imidacloprid exposure. Using a multiple stressor toxicokinetic-toxicodynamic modeling approach, we showed that both starvation and other toxic effects of imidacloprid play a role for determining mortality in constant exposure to the insecticide.


Introduction
To protect crops and seeds from pests, about 3 billion tons of pesticides are applied annually to fields worldwide [1]. A fraction of this reaches other environmental compartments such as surface waters via runoff, spray drift and leaching. One of the world's bestselling insecticide is imidacloprid, 1-(6-chloro-3-pyridylmethyl)-Nnitroimidazolidin-2-ylideneamine, which belongs to the chemical group of neonicotinoid insecticides [2]. Neonicotinoids have selective toxicity for insects and act by binding to the nicotinic acetylcholine (ACh) receptors in the receiving nerve cells of the central nervous system [3,4]. Mammals have lower numbers of nicotinic receptors with high affinity to neonicotinoids, which is why the toxicity of these insecticides is low in mammals [5]. Imidacloprid has a relatively high water solubility (610 mg/L in 20uC H 2 O; log K ow = 0.57) and therefore, a great potential to reach water bodies. Accordingly, several studies have reported the occurrence of imidacloprid in surface waters [6,7] where it may affect non-target organisms such as Gammarus pulex (Crustacea, Amphipoda, Gammaridae). The concentrations of imidacloprid in surface waters in Sweden reported by Kreuger and coworkers (max. 15 mg/L) [7] are below lethal acute toxicity levels in G. pulex (50% of the test individuals die after constant exposure to 270 mg/ L for 4 days [8]). However, the lower concentrations found in water bodies might cause sublethal effects.
In aquatic environments pesticide contamination generally occurs in pulses due to fluctuation in rainfall, seasonal application of pesticides, and accidents [9,10]. Because neonicotinoids lack ester bonds and thus cannot be hydrolyzed by ACh esterase, also temporary exposure to these insecticides can generate sustained activation in receptors and cause long lasting effects. However, it has been shown that imidacloprid can be dissociated (dissociation constant 0.419 min 21 ) and removed from ACh receptors by ACh and other ligands [11]. Therefore, it is possible that organisms recover between imidacloprid pulses. On the other hand the elimination of imidacloprid in G. pulex is very slow [12] and the substance is not biotransformed [13]. Thus, one could also expect cumulative effects from subsequent exposure events. Imidacloprid has not shown cumulative effects on Gammarus roeseli survival after repeated pulses of the insecticide [14]. However, a cumulative sublethal effect (increased drifting) has been reported [15].
By binding to the ACh receptors and interfering with nerve impulses, imidacloprid causes twitching, cramps and muscle weakness. Therefore, it impairs invertebrate movements and can lead to starvation and death via dysfunctional feeding behavior. It has been shown previously that imidacloprid inhibits feeding of many non-target aquatic species [16][17][18]. However, the connection between impaired feeding and mortality via starvation has not been further investigated, in spite of the importance of conserving populations of aquatic shredding invertebrates like Gammarus pulex.
We studied the effects of imidacloprid on feeding rate, lipid content, immobility and survival of G. pulex in 14-day and 21-day long experiments. As exposures in aquatic environments generally occurs in pulses, we exposed the animals to two high imidacloprid concentrations. To further study the impact of imidacloprid under low constant exposure, we also exposed G. pulex to the time weighted average concentration, which was 15 mg/L. Using the time weighted average concentration allowed us to compare effects of different exposure patterns while still employing the same overall dose (time 6 concentration). We chose the low concentrations for the constant treatments [8], because we hypothesized that starvation causes death in constant treatments due to impaired movements. In the pulsed treatments, however, feeding activity might recover between chemical pulses and thus starvation would not play a big role in determining mortality. In contrast, we can observe other toxic effects of imidacloprid, e.g. direct mortality, in response to high imidacloprid peaks. Thus, there might be different mechanisms behind mortality in pulsed and constant treatments. To investigate this, we used toxicokinetic-toxicodynamic (TKTD) modeling to analyze the survival data and tested whether fitted model parameters indicate different effect mechanisms in the pulsed and constant exposures. We performed also a starvation experiment, without adding imidacloprid, to further investigate the effect of starvation on survival and to test if using a calibrated starvation model would predict survival in our constant imidacloprid treatments. The chemical effect and the starvation models were also combined to develop a multiple stressor model which again was tested by simulating survival in constant exposure treatments.

Test Animals and Chemicals
Gammarus pulex is an important invertebrate species in lentic waters for e.g. decomposition of organic material and nutrient cycling [19]. The G. pulex test individuals in our study were collected from a small headwater stream in the Itziker Ried, Switzerland (E 702150, N 2360850). No permission to collecting was required as G. pulex is not an endangered species and the site is located on public land. The test animals were maintained for 5-7 days prior to the experiments in a large aquarium in a temperature controlled room (13uC, 12:12 light:dark photoperiod) and were fed with horse chest-nut (Aesculus hippocastanum) leaves which were inoculated with the fungi Cladosporium herbarum for at least 10 days [20]. The water in the aquarium was preaerated artificial pond water (APW, Table S1 in File S1). 14 C-labelled imidacloprid (radiochemical purity 96.97%) was purchased from the Institute of Isotopes Co., Ltd. Budapest, Hungary and unlabeled material (chemical purity 99.9%) from Sigma-Aldrich. A mixture of both was dissolved in acetone and used for dosing.

Imidacloprid Experiments
A 14-d and a 21-d toxicity experiment including three treatments plus controls in each were conducted. Two of the treatments (A, B) included two 1-day imidacloprid pulses with differing recovery time between pulses in uncontaminated APW. In another treatment (C), the concentration was maintained constant (15 mg/L, 0.06 mmol/L in both experiments) but the overall dose was the same over time as in the pulsed treatments (i.e. time-weighted average concentration). All treatments included 7 replicate beakers, one plain and one solvent control beaker. Each 600 mL Pyrex beaker contained 500 ml of APW and 5 leaf discs (diameter of 20 mm, Aesculus hippocastanum leaves inoculated with Cladosporium herbarum). Ten G. pulex were placed in each beaker a day prior to the experiments. The beakers were covered with parafilm and kept in a climate chamber (13uC, 12:12 light:dark photoperiod). The beakers were spiked individually and after spiking, test solutions were stirred with a glass rod and 1 mL samples were taken from the solution to quantify the initial chemical concentration in medium. Ten mL of Ecoscint A scintillation cocktail (Chemie Brunschwig, Switzerland) were added to the samples and activities were counted using a liquid scintillation counter (LSC, Tri-Carb 2200CA, Packard, USA). Samples to determine imidacloprid concentrations in water were taken throughout the experiments (see time points and concentrations in Tables S2 and 3 in File S1). Only the total radioactivity in the aqueous samples was measured and therefore imidacloprid could not be differentiated from its possible breakdown products. For example, organisms might be exposed to the breakdown products of imidacloprid rather than the parent compound due to fast photolysis of imidacloprid in aqueous medium with a half-life of 1.2 h at 290 nm irradiation [21]. However, it is also shown that wavelength has a great impact on the photolysis and already in 365 nm, the half-life is extended to 18 h [22]. The wavelengths in our experiments resemble those of day light ranging from 380-730 nm (relative intensity being the highest in 580 nm). Thus under the conditions of our experiments, imidacloprid most likely is less susceptible to photolysis and in an earlier study no breakdown products of imidacloprid were observed in G.pulex samples [13]. The test solution was changed at least every 5 days. Always during water change and every time when eaten leaf discs were observed, they were replaced by new ones. Water pH, conductivity and oxygen concentration were measured in exposed and non-exposed conditions during experiments (see in more detail in Tables S4 and S5 in File S1).
In the 14-day experiment mortality, immobility and consumption of leaf discs were observed. In addition, internal concentrations of imidacloprid in G. pulex were measured. The pulsed treatments (A, B) had 4 (A) and 8 (B) days between 1-day pulses (concentration 90 mg/L = 0.35 mmol/L) and individuals in treatment C were exposed constantly to a concentration of 15 mg/L (0.06 mmol/L). Immobility was defined as incapability of moving after ten gentle prods with a glass rod. Immobile individuals were taken out of the beakers and frozen until analysis of internal concentrations. In addition, mobile individuals were sampled for analysis of internal concentrations at the end of the experiment (A and B: n = 21 per treatment, C: n = 10) and during the experiment from additional beakers of the treatments A and B. These supplemental beakers were not used for observation of survival, immobility and consumption of leaf discs. See detailed sampling times and internal concentrations in Tables S6, S7, S8 in File S1.
Sample processing and quantification of radioactivity in G. pulex were measured similarly to Ashauer and co-workers [12]. In short, individuals were plotted dry with tissue paper, weighed in preweighed glass vials and frozen in 220uC. For analysis, 3 mL of the tissue solubilizer Soluene-350 (Perkin Elmer, USA) was added to vials. Vials were placed in a water bath (60uC, 24 h) and after cooling down, 15 mL of scintillation liquid Hionic Fluor (Perkin Elmer, USA) was added. Radioactivity was counted using a liquid scintillation counter; color quenching and efficiency were corrected using external standards and background activity (i.e. activity in control samples was subtracted from counts of the samples).
The second, 21-day long experiment included also observations of survival, immobility and food consumption. In addition, lipid content was measured from immobile individuals, which were sampled and frozen any time they were observed as well as from mobile individuals sampled at the end of the experiment. The pulsed treatments (A, B) had 4 (A) and 11 (B) days between 1-day pulses (concentration 140 mg/L = 0.59 mmol/L) and individuals in treatment C were exposed constantly to a concentration of 0.06 mmol/L.

Feeding Rate
In the 14-day experiment, food consumption was measured as the number of leaf discs consumed by G. pulex individuals in each beaker. Every time when a leaf disc was fully eaten, it was replaced by a new leaf disc and the exchange was noted. This way of measuring was not based on the mass of leaf discs but only on the   Pulsed treatment with a short interval in uncontaminated water between imidacloprid pulses (4 days). 2 Pulsed treatment with a long interval in uncontaminated water between imidacloprid pulses (8 days). 3 Pulsed treatment with a long interval in uncontaminated water between imidacloprid pulses (11 days). 4 Between control and treatment. 5 Among all treatments within one experiment. doi:10.1371/journal.pone.0062472.t001 number of same sized discs. During the 21-day experiment food consumption was measured as the mass of leaf material. Wet weight of the leaf discs of each beaker was measured before providing them to the test organisms and when removing the rest of them from the beakers. Food consumption is given as cumulative amount consumed over time (either amount of leaf discs (14-day experiment) or mg (21-day experiment)). The food consumption was divided by the number of mobile organisms in the respective beaker (amount consumed/G. pulex). Statistical testing to compare feeding among treatments was performed for cumulative food consumption at the end of the experiments using the non-parametric Kruskal-Wallis test and further pairwise testing with the Wilcoxon rank sum test. The assumption of normality could not be tested due to the too small sample size (i.e. number of replicate beakers, 7 per treatment) and thus a normal distribution of the data could not been assumed and one-way analysis of variance could not been used. Analyses were performed using the software R (www.r-project.org).

Lipid Content
The lipid content was analysed from immobile G. pulex sampled during the 21-day experiment (treatment A: 33 samples, treatment B: 28 samples, treatment C: 1 sample) and mobile individuals sampled at the end of the experiment (n = 21 for each treatment, total of 30 for controls: 5 from plain and 5 from solvent control beaker of each treatment). A gravimetric method was used to determine the lipid content according to Kretschmann and coworkers [23]. In short: Extraction was done using H 2 O, i-PrOH, and cyclohexane (11:8:10) in 2 mL Eppendorf tubes. 300 mg of zirconia/silica beads (Ø 0.5 mm, BioSpec Products, Bartlesville, OK, USA) was added to iPrOH/cyclohexane solution together with the sample and FastPrepH FP120 Bio 101 (Savant Instruments, Inc., NY, USA) was used to break down the tissues of G. pulex. Nanopure water was added to samples. A water content of 77% of G. pulex wet weight was assumed to achieve a 11:8:10 ratio of H 2 O, i-PrOH, and cyclohexane. Then, samples were vortexed, centrifuged (20 min, 450 g, 20uC) (Centrifuge 5417R, Vaudaux-Eppendorf AG, Schönenbuch/Basel, Switzerland) and the organic phases were separated. Volumes of 435 mL of cyclohexane and 65 mL of iPrOH were added once more, and after vortexing, centrifugation and separation of the organic phase, the solvents of the combined organic phase aliquots were evaporated under nitrogen flow and extracts were dried at 60uC for 14 hours. The remaining phases, i.e. the lipids, were weighted. The weight of lipids was divided by total wet weight to obtain the lipid content as a percentage. Because the lipid content in treatment C did not follow a normal distribution (see Figure 1 and Figure S1, raw data are provided in Table S11 in File S1), the differences among treatments at the end of experiment were compared by the Kruskal-Wallis rank sum test and further pairwise testing (controltreatment) using the Wilcoxon rank sum test. The software R was used for the analysis.

Chemical Stress Modeling
The survival/mobility data of the imidacloprid experiments was analyzed using the toxicokinetic-toxicodynamic (TKTD) model GUTS (General Unified Threshold model for Survival) published by Jager and co-workers [24]. The model was implemented in the software ModelMaker 4 (Cherwell Scientific Ltd., Oxford, UK). The fraction of mobile animals over time was used to calibrate the model. Note that immobile animals were removed from the experiment, allowing us to apply the same assumptions about error structure to our immobility data as to survival data and to use GUTS for modeling the mobility/survival data. As we suspect that constant low concentration of imidacloprid (treatments C) might have a different mechanism for survival/mobility than the pulsed treatments (A and B), the TKTD model was calibrated separately using the constant treatments C of both experiments together and the pulsed treatments (A, B) of both experiments together. In addition, the model was fitted to all data in order to compare the parameter estimates with those of the separately fitted data sets. The parameter estimates from the fit to all data were used as initial values for fitting the model separately to pulsed and constant treatments.
Imidacloprid is not biotransformed in G. pulex [13]. Therefore, a one-compartment toxicokinetic model (Eqn 1) was used to simulate the internal concentration of imidacloprid.  [12]. We validated this TK model by comparing its predicted internal concentrations with measured internal concentrations in our first experiment.
Survival modeling was based on the GUTS model and the modeling is similar to our previous study with G. pulex and propiconazole [25]. However, here least squares optimization together with the Marquardt algorithm was used to fit models to the data. Confidence intervals (95%) were calculated from standard errors as described in Motulsky & Christopoulos 2003 [26]. Fitting was also performed by maximizing the log-likelihood function, which maximizes the likelihood of yielding the parameter set which best describes the number of death events between time intervals (see Table S14 in File S1). In this study, the cumulative fraction of survivors over time was better described by the parameters found via least squares optimization.
Two models, either assuming stochastic death (SD) or individual tolerance distribution (IT) were used separately to test which hypothesis of death applies for imidacloprid. SD models have one value for the threshold of survival and after exceeding it, an organism has an increased probability to die. In contrast, according to IT models the threshold is distributed within the population and death is instantaneous after exceeding the individual threshold. To date, it is not known which of the two hypotheses describes our data better (see also discussion in [24] and [25]). Therefore, both models, SD and IT, were calibrated and the goodness of fit values were compared. Equation 1 was used to simulate the internal concentrations (C int ) in the survival model for both, GUTS-SD and GUTS-IT. The implementation of the stochastic death model (GUTS-SD) is given in Eqns 2 to 5. Eqns (2) and (3) were used to calculate the cumulative hazard at time t (H(t)). where where Sb is the background survival probability [2] describing survival in unexposed conditions.   Once the cumulative hazard H(t) is obtained, the survival probability, S (t) [2], was calculated using Eqn 5.
The model that assumes the threshold for death to be drawn from an individual tolerance distribution (GUTS-IT model) is presented in Eqns 6 and 7. The IT model uses the same dose metric, scaled damage D * , as the SD model (Eqn 2). Cumulative threshold distributions are based on a log-logistic cumulative distribution function (Eqn 6). The resulting survival probability is given by Eqn 7.
where F(t) is the log-logistic cumulative distribution function for the threshold [2], a is the median of the distribution [mmol/kg], b determines the width of the distribution [2] and the 'max' function selects the largest value of the dose metric D* that occurred until time t.
To compare the goodness of fit among models and calibration data sets, the mean percentage error (MPE) was calculated (Eqn 8) [25]. The MPE was calculated for each treatment separately and therefore for each treatment the goodness of fit could be compared between a) stochastic death and individual tolerance models and b) models calibrated with different data sets (i.e. pulsed or constant exposures).   . In the starvation model (II), lack of food (LF) for the 14-day experiment was set to 1.0 and for the 21-day experiment 0.5 due to differences in feeding activity (no feeding in 14-day experiment, ca. 50% reduced feeding in the 21-day experiment). The chemical stress model was GUTS calibrated with pulsed toxicity data sets. doi:10.1371/journal.pone.0062472.g006 survivors, S obs is the observed fraction of survivors [2], S model is the model prediction of the fraction of survivors [2] and n is the number of data points used in the calculation.

Starvation Experiment and Modeling
As we hypothesise that organisms in the constant imidacloprid exposure die due to impaired movements leading to starvation, an experiment studying the effect of starvation, without exposure to imidacloprid, on survival of G. pulex was conducted. There were 30 replicates for both, the control and the starvation treatment. In the control group, one leaf disc was provided for food and more was given when the disc was eaten. In the starvation treatment, no food was given. G. pulex were placed individually in 100 mL beakers in order to prevent cannibalistic behavior of the test animals, which is more likely without leaf discs. Mortality was monitored for 34 days, i.e. long enough to observe mortality of at least half of the animals. Experimental water was APW (Table S1 in File S1) which was renewed weekly.
The mortality observed in this experiment was compared with mortality under the low, constant exposure to imidacloprid. For easier comparison, modified versions of the TKTD models described above (SD and IT) were calibrated using the starvation data and the survival in the constant imidacloprid treatments was predicted using this new starvation model. Instead of using chemical internal concentration causing the scaled damage D * (t), a new concept, lack of food (LF), leading to the damage (Eqn 9) was introduced. In other words, the dose metric is LF, defined as the relative lack of food compared to control conditions (LF = 1 -(available food/food available in control)). Thus the survival model for starvation consists of equations identical to Eqns 3-7, Eqn 2 being replaced by Eqn 9 (LF replaces C int (t)), which also leads to different dimensions of the model parameters. The LF was set to 1 when calibrating the model with starvation data as well as simulating the survival in the constant imidacloprid treatment in the 14-day experiment where hardly any food consumption was observed. When simulating the survival in 21-day imidacloprid experiment, the LF was set to 0.5 as food consumption in this experiment appeared to be only partially inhibited (the feeding activity was approximately half of control levels, see Figure 1).
where D (t) describes the damage caused by lack of food [2], k d is damage recovery rate [1/d] and LF (t) is lack of food [2]. Units of the following parameters in Eqns 3-7 were therefore different from described above: and a [2]. The background hazard rate was calibrated for each experiment separately.

Multiple Stressor Modeling
The starvation model was combined with the chemical stress model and the survival in constant imidacloprid exposure was simulated in order to test whether survival is determined by both, the effect of starvation and other toxic effects of imidacloprid. The chemical stress models were the GUTS models described above (Eqns 1-7) calibrated with the data from the pulsed toxicity treatments alone. In this model the effect of starvation is excluded from other chemical effects because we can assume that in the pulsed toxicity treatments starvation does not play a role due to possible recovery of the movements and feeding between chemical pulses. To implement the multiple stressor model, equations for both, chemical stress (SD: Eqns 1-4; IT: Eqns 1-2, 6) and starvation (SD: Eqns 9, 3-4; IT Eqns 9, 6) were followed. Then the cumulative hazards H (t) of both chemical and starvation stress were added (SD model, Eqn 10) or both cumulative distribution functions for the threshold F (t) were subtracted (IT model, Eqn 11) to predict survival in the constant imidacloprid treatments where both processes likely play role.

Feeding Rate
Inhibition of feeding by imidacloprid has been observed in many invertebrate species [16][17][18][27][28][29]. We observed that feeding of Gammarus pulex was heavily inhibited by imidacloprid in the constant treatment of the 14-day experiment ( Figure 1, Table 1, Table S9 in File S1) while in the pulsed treatments the effect was not as strong ( Figure 2, Table 1, Table S9 in File S1). However, in this experiment, the method of measuring feeding activity was only semi-quantitative because it was based on the number of same-sized, however undefined by their mass, leaf discs consumed by G. pulex individuals. Thus, we focused on the 21-day experiment to draw conclusions with regards to the effects of imidacloprid on feeding activity. In the 21-day experiment, feeding activity was measured as the mass of leaf material eaten [mg]. There the effect on feeding in the constant treatment was again evident, however, not as strong as in the 14-day experiment ( Figure 1, Table 1, Table S10 in File S1). In the pulsed treatments, no effects were observed on feeding activity -the organisms started to feed roughly 2 days after they were transferred to uncontaminated media ( Figure 2, Table 1, Table S10 in File S1). This finding is in agreement with the fact that imidacloprid binding to ACh receptors in insect membranes is reversible, i.e. it can be dissociated as well as removed by ACh and other ligands [11,30,31]. However, in spite of the fairly fast dissociation of imidacloprid from the ACh receptors (0.419 min 21 [11]), the elimination of imidacloprid has been shown to be slow in G. pulex [12] and can be associated with the receptors again. The internal concentrations in our experiments decrease relatively fast ( Figure 3) and 2 days after a pulse around 60% of imidacloprid is left inside the organisms. Once feeding activity is recovered in 2 days, it seems that the internal concentrations during these days fall below the threshold of preventing animals from feeding.
The ability to recover from imidacloprid pulses has been observed also by other authors [17,18,32]. Alexander and coworkers [17] observed that both mayfly (Epeorus longimanus) larvae and the oligochaete worm (Lumbriculus variegatus) could recover from 1-day exposure to imidacloprid in 4 days and the recovery potential was concentration dependent [17]. Concentration dependency was not observed in this study -in fact in the 21day experiment where we had higher imidacloprid pulses, the feeding rate was not different from the controls while in the 14-day experiment we observed an effect. However, this observation could be caused by differences in the measurement method (number of leaf discs versus mg) and variation in organism fitness between the experiments, which can be seen for example in the control mortality (Figure 1).

Lipid Content
We observed reduced lipid content in G. pulex after exposing them constantly to imidacloprid for 21 days (21-day experiment, Figure 1). In the pulsed treatments, lipid content was not different from that of the controls (medians of A: 1.36%, B: 1.45%, Figure  S1). When comparing among all treatments using the Kruskal-Wallis rank sum test, the p-value of 0.017 indicated significant differences and the pairwise comparisons showed that the constant treatment was responsible for the difference (Wilcoxon rank sum test between control and C: p = 0.0053). The differences in lipid content between pulsed treatments and controls were not significant (A: p = 0.4646; B: p = 0.6834). The decreased lipid content was observed in the same treatment where the feeding was inhibited (C). Therefore, we can conclude that lipid content is affected by a decreased feeding rate and might imply starvation of G. pulex in the presence of imidacloprid.

Mortality and Immobility due to Imidacloprid
We observed a sudden drop in survival in the end of the constant treatments (C), especially in the 14-day experiment ( Figure 1, part III, Tables S12 and S13 in File S1). The percentage of dead out of all immobile individuals (dead+only immobile) was high in the constant treatments. In fact, we did not observe many individuals classified as immobile in the constant treatment (i.e. 14day experiment: 3 out of 70, 21-day experiment: 1 out of 70). However, almost all the ''mobile'' animals were close to the limit of immobility, they were passive and could not move in a normal way. Similar behavior has been observed before in Gammarus roeseli exposed to 12 mg/L imidacloprid [15] which is close to the concentrations in our constant treatments. This effect of imidacloprid decreased the feeding rate and might have caused starvation in our experiments. In our treatments with high imidacloprid pulses, the organisms were mostly immobile (Figures 1 and 2), thus even in the highest concentration that we used (0.59 mmol/L), the acute lethal toxicity was not reached within one day. There were no differences in the mobile fraction at the end of both experiments between treatments with short (A) and long (B) recovery time between pulses, also not in feeding activity and lipid content in the 21-day experiment (Figure 2, Figure S1). This indicates that organisms recovered fast from imidacloprid exposure between the pulses, even in the treatments A which had short intervals between pulses. Calculated 95% organism recovery times were 12.7 and 12.3 days according to IT and SD models. Organism recovery times were calculated as the time when the modeled internal damage has dropped to 5% from the maximum in a pre-defined exposure scenario [25,33,34]. Note that organism recovery refers to the recovery of the underlying damage that causes effects on survival until it reaches levels far below those causing mortality. Thus, the observed fast recovery with regards to mortality after the pulses does not necessarily conflict with organism recovery times of 12 to 13 days.
We observed more mortality in the constant treatment of the 14-day experiment than in the constant treatment of the 21-day experiment, even though the exposure concentration was the same and the animals were exposed longer in the 21-day experiment. This variation in our results was also seen in the background mortality: during the 14-day experiment, mortality in the controls was much higher (Figure 1). This observation can be explained by organism fitness, which varies when we collect animals from the field for each experiment. For instance, season has been shown to influence the condition of gammarids [35][36][37]. One important seasonal factor is food availability. After leaf fall during autumn, organisms have shown to have better lipid reserves during winter (maximum for males in November and for females in January) while during summer lipid reserves are the lowest due to scarcity of food [35]. This can make summer populations more sensitive to toxicants [37], especially to imidacloprid which interferes with feeding behavior.
We conducted the experiments in November (14-day experiment) and February (21-day experiment) -both are months of the season where we expect high lipid contents according to Stroda and Cossu-Leguille [35]. We can hypothesize that in November, when we conducted the 14-day experiment resulting in high mortality in the constant treatment, the organisms had not had enough time to build up and store lipids after the leaf fall. Thus, in our 21-day experiment, organisms possibly had better lipid reserves, which might be why less individuals died than in the 14-day experiment. An alternative explanation could be preexposure of the animals to pollutants in the field before collection, which is more likely in November and the weeks before. Although the collection site is in a headwater stream with low probability for pollution, it cannot be ruled out that some exposure occurred, for example to pesticides applied in autumn and transported via runoff in autumn rains. Independent of the possible causes, the variance in background mortality among experiments was corrected in our survival models by using experiment specific background hazard rates (Eqn 4, 14-day experiment: 0.0144 [1/ d]; 21-day experiment: 0.0079 [1/d]).
From our survival models, stochastic death (SD) and individual tolerance distribution (IT), the IT model fitted better to the data ( Table 2). The mean percentage error (MPE) in the treatments was always below 5.4% in the IT models (mean 4.361% when all data was used for calibration) while in the SD models the maximum error was as high as 22.7%, with a mean of 10.366.8% when all data was used for calibration. This would imply that the data provided here supports the individual tolerance distribution hypothesis, which assumes that organisms have individual effect doses which are distributed within a population. A possible interpretation is the variability of lipids and other energy reserves within the population. Similar and opposing findings, i.e. studies where the SD hypothesis described the data better, have been published earlier [34,[38][39][40]. The applicability of either extreme hypothesis seems to be chemical and species specific. Because our data supports the IT theory, we use the IT model results to compare survival among pulsed and constant treatments. The toxicokinetic sub-model calibrated by Ashauer and coworkers [12] predicted well the internal concentrations, which were measured in the 14-day experiment (i.e. model validation, see Figure 3).
Parameters differed when the models were calibrated with either pulsed or constant survival/immobility data. In Figure 4, parameter estimates of the damage recovery rate (k r ), the median of the threshold distribution (alpha), and the width of the distribution (beta) for pulsed treatments, constant treatments and all data are illustrated. The biggest differences between constant and pulsed calibration data can be seen in the values of the parameter beta, which was extremely high for the constant data compared to the pulsed data. The value of beta determines the width of the threshold distribution -the higher the value is, the narrower is the distribution and the steeper is the drop in the fraction of survival/mobility close to the median value alpha. This would imply that the individuals in a population react fairly similarly to starvation, or the damage describing decreasing lipid content.
However, parameter estimates should not be over-interpreted because they are linked to each other (parameter co-variation). To compare among calibration data sets in a more reliable manner, the whole set of parameters can be used to simulate the survival in another exposure type (i.e. from calibration to the constant data set to simulate survival in the pulsed scenario and vice versa). When the pulsed calibration data was used, mobility in the constant scenario was rather well predicted with MPE of 5.569.5% (14-day experiment) and 4.462.4% (21-day experiment) (Figure 4, II). However, the high mortality in the end of constant treatments was not captured (i.e. the effect of starvation). When the constant calibration data was used, the mobility in the pulsed scenario could not be predicted well; the mobile fraction went directly to zero in all pulsed treatments (Figure 4, III).
This simulation result indicates that there are different processes, described by very different model parameters, governing the mobility under long, constant exposure to imidacloprid vs. the pulsed exposure scenario. The toxic processes in the constant treatments can be related to a low degree of imidacloprid binding in nicotinic receptors making organisms passive and causing death via impaired movements and starvation. In the pulsed treatments other toxic effects from higher degree of binding to the target sites play also a role. This can be seen also as a difference between the adverse effects of imidacloprid in pulsed treatments, where imidacloprid immobilizes the organisms immediately, and more chronic effects like starvation which appears as a result from impaired movements in long constant exposure. However, we cannot rule out the possibility that the differences in model parameters can also be caused by differing proportions of dead and immobile individuals among treatment types (the model treats them equally, because immobile individuals were removed). Almost only mortality was observed in the constant treatments while there were much more immobile individuals than dead ones in the pulsed treatments (pie charts in Figures 1 and 2, part III).

Mortality in the Starvation Experiment
In one of the target organisms, the tobacco whitefly (Bemisia tabaci), imidacloprid has been shown to cause starvation [41,42]. Even though the exposure route is partly different for aquatic species (oral uptake alone for the tobacco whitefly versus oral uptake and diffusion from water for G. pulex), low concentrations of imidacloprid can cause death of non-target amphipods in aquatic environments via behavioral changes which prevent the organisms from feeding. To investigate whether other effects of imidacloprid than starvation play a role in the constant exposures, we performed a pure starvation experiment where no imidacloprid was added. We calibrated the TKTD models with food limitation as dose metric and the models were then used to simulate mortality in the constant imidacloprid treatments. Results of the starvation experiment and calibration of the models are shown in Figure 5. The simulation of survival in the constant lack of food conditions, representing our constant imidacloprid exposure, showed rather poor agreement with the measured values, especially for the 14day experiment ( Figure 6, part II). However it must be noted, that this prediction is based only on starvation stress. As cholinergic neurotransmission has been suspected to have a central role in neurotransmission in invertebrate central nervous system [43], likely also other pathways of toxicity, e.g. failure of respiration, than starvation through impaired movements occur in imidacloprid exposure.

Survival in Multiple Stress Conditions
We developed a multiple stress model, combining the effects of starvation and the other toxic pathways of the chemical, for simulating survival in the constant imidacloprid exposure. The SD model for multiple stressors did not predict survival well when we compare mean percentage error (MPE) values between three different model types: The 14-day experiment was better predicted by the chemical stress model and the 21-day experiment by the starvation model ( Figure 6). However, as discussed above, the IT model had a better fit to the data ( Figure 2, Table 2) and therefore we should focus on comparing results given by this model (see solid line in Figure 6). The agreement of the simulation with the data increased with the multiple stress model, when compared to simulations using either of the individual stressor models (i.e. chemical stress and starvation) alone especially in the case of the 14-day experiment ( Figure 6). Using the multiple stress model did improve the predictive power, however, one pattern in the data, i.e. the sudden drop in survival at the end of the experiments, was not predicted quantitatively.

Conclusions
Organisms in nature are facing multiple stress conditions, natural and anthropogenic. We have investigated the effects of the insecticide imidacloprid in Gammarus pulex and showed that multiple stress pathways might influence organism survival, rather than just a single mechanism of toxicity or stress pathway alone. By binding to acetylcholine receptors, imidacloprid impairs invertebrate movements and feeding behavior and therefore, in addition to other pathways of toxicity, starvation can influence organism survival if the exposure lasts longer than the duration of standard toxicity tests. Another noteworthy point of our findings is that the effect of reduced feeding on invertebrate survival is affected by food availability and initial nutritional status of organisms which vary along with the season. Nevertheless, in aquatic systems the exposure times are usually short due to water flow and dilution. We showed that in these conditions, organism movements and feeding can recover fast after the exposure and therefore the stress from starvation is reduced, however, other toxic pathways can still affect organism survival. File S1 Supporting Tables. Table S1. Composition of artificial pond water (APW) and stock solutions. Table S2. Imidacloprid concentrations in water in 14-day experiment. Table S3. Imidacloprid concentrations in water in 21-day experiment. Table S4. Water characteristics in 14-day experiment. Table S5. Water characteristics in 21-day experiment. Table S6. Data on internal concentrations in 14-day experiment: Tr A. Table S7. Data on internal concentrations in 14-day experiment: Tr B. Table S8. Data on internal concentrations in 14-day experiment: Tr C. Table S9. Cumulative food consumption (leaf discs/G. pulex individual) in 14-day experiment. Value is corrected with number of mobile individuals in the beaker. Table  S10. Cumulative food consumption (mg/G. pulex individual) in 21day experiment. Nm denotes ''not measured'' in that time point. This is taken into account in the value in the next time point. Table S11. Lipid content of immobile Gammarus pulex individuals in the 21-day experiment and mobile individuals in the end of experiment. Table S12. Number of mobile and immobile individuals in the 14-day experiment. Cint A and Cint B denotes beakers which were used to sample mobile individuals and were not used for survival modeling. Table S13. Number of mobile and immobile individuals in the 21-day experiment. Table S14. Parameter estimates for different calibration data sets. Parameter estimates of stochastic death (SD) and individual tolerance distribution (IT) models based on fit using log-likelihood function or ordinary least squares fit. (DOCX)