Virgile Baudrot
email: virgile.baudrot@univ-fcomte.fr

Clémentine Fritsch
email: clementine.fritsch@univ-fcomte.fr

Antoine Perasso
email: antoine.perasso@univ-fcomte.fr

Malay Banerjee
email: malayb@iitk.ac.in

Francis Raoul
email: francis.raoul@univ-fcomte.fr

Effects of contaminants and trophic cascade regulation on food chain stability: application to cadmium soil pollution on small mammals -raptor systems

Keywords: critical transitions, contaminant transfer, food chain, ordinary differential equations

Within the current context of global change, ecosystems face simultaneously multiple stressors such as environmental pollution, excessive enrichment, and additional disturbances affecting population dynamics. Numerous theoretical and empirical studies explored the transfer of contaminants by food ingestion between trophic levels, and how contaminant effects on survival and fecundity may change regime shifts of classical predator-prey dynamics. However, the extent to which those contaminants may influence the cascade effect of external stressors, as a change in resources (bottom-up cascade) or a variation in predator abundance (top-down cascade) is still poorly understood. We develop a data-driven model to explore how soil contamination modulate the food chain stability and resilience to changes in prey nutrient and in the apex predator mortality rate. We particularly focus on the ecotoxicological impact of the trace metal, cadmium, on a widespread raptor, the barn owl (Tyto alba), feeding on several prey distinct by their trophic positions: herbivores (Microtus spp), omnivores (wood mouse Apodemus sylvaticus and bank vole Myodes glareolus) and insectivores (shrews: Crocidura spp and Sorex spp). Our model reveals the alternative steady states in population dynamics and the occurrence and position of regime shifts where a subtle change in conditions causes a sudden shift in the ecological system. Based on mathematical modelling and bifurcation analysis, the results show for instance that under toxicity threshold, where no population decline is observed, the contaminant weakens food-chain resilience. Then at higher contamination, the toxic effects on predator releases the top-down control over prey that may increase. This range of chemical stress overturns the paradox of enrichment, a central concept in trophic cascade theory. The transition phase at the highest contamination, where the whole community collapse, exhibits multiple patterns, from smooth to abrupt, depending on external stressors and the prey population. Thus, this work provides a methodology to identify ecological traits of preys that are critical for transferring adverse effects of contaminants across the whole community.

Introduction

Across the globe, anthropogenic activities increase the input of contaminants in ecosystems, disturbing their regime and potentially threaten their health. One of the most widely observable impact is the lost of apex predators [START_REF] Estes | Trophic downgrading of planet Earth[END_REF][START_REF] Heath | Understanding patterns and processes in models of trophic cascades[END_REF]. As a consequence, there is a surge in attention concerning the effect of pollution at the food web scale. The response of a food-chain to external disturbances may be gradual (i.e., from a fixed state to another or to an oscillating pattern), or remains inert until a sudden regime shift, a critical transition, where the food-chain structure and composition are drastically altered [START_REF] Scheffer | Critical transitions in nature and society[END_REF]. These sudden changes often correspond to abrupt degradation of ecosystem services which are difficult, if not impossible, to recover [START_REF] Rockström | A safe operating space for humanity[END_REF]. Due to their socio-economical importances, analysis of discontinuous transitions (i.e., regime shift) have been widely studied [START_REF] Scheffer | Generic Indicators of Ecological Resilience: Inferring the Chance of a Critical Transition[END_REF], but very few theoretical works have considered terrestrial food-chain impacted by the combination of a chemical contaminant and external disturbances (i.e., trophic cascades). This lack of a mathematical description on how polluted ecosystems respond to external stressors (e.g., eutrophication, climate change, diseases) is due to the combination of many direct and indirect effects [START_REF] Fleeger | Indirect effects of contaminants in aquatic ecosystems[END_REF][START_REF] Rohr | Community ecology as a framework for predicting contaminant effects[END_REF][START_REF] Relyea | Assessing the ecology in ecotoxicology: a review and synthesis in freshwater systems[END_REF][START_REF] Laender | Theoretically exploring direct and indirect chemical effects across ecological and exposure scenarios using mechanistic fate and effects modelling[END_REF] Direct effects of severe contaminant exposure as population declines (by increasing mortality or reproductive failure) are possibly predictable from controlled laboratory toxicity tests. Those tests are the classical approach to assess the potential risk of contaminants in ecosystems, and most of them use one single species and one chemical compound, and quantify the toxicological effect through several threshold of exposure concentration (i.e., NOAEL, LOAEL, LC 50 ) [START_REF] Villeneuve | Vision & strategy: Predictive ecotoxicology in the 21st century[END_REF]. However, there is a difficulty in predicting the ecological effects of an environmental pollution from those laboratory results [START_REF] Mineau | A review and analysis of study endpoints relevant to the assessment of "long term" pesticide toxicity in avian and mammalian wildlife[END_REF][START_REF] Villeneuve | Vision & strategy: Predictive ecotoxicology in the 21st century[END_REF].

Indeed, the exposure of vertebrate wildlife to environmental contaminants integrates complex interactions between the routes of exposure (inhalation, ingestion, dermal absorption and maternal transfer), the intensity (concentration in the vector agent and bioavailability), the frequency and duration of the exposure, and the life history traits of exposed organisms [START_REF] Smith | Contaminant exposure in terrestrial vertebrates[END_REF]. For instance, among factors affecting metals body burden in small mammals, several field studies enlightened the critical role of the feeding behavior and the habitat use (Veltman et al., 2007b;Fritsch et al., 2010a;[START_REF] Van Den | Cadmium accumulation in small mammals: species traits, soil properties, and spatial habitat use[END_REF]. In addition to this great variability in the exposure, the fate of a contaminant within the organism (i.e., toxicokinetics) and its toxicological effect (i.e., toxicodynamics) differ greatly between and within species [START_REF] Rainbow | Trace metal concentrations in aquatic invertebrates: why and so what?[END_REF][START_REF] Burger | Assessment and management of risk to wildlife from cadmium[END_REF]Fritsch et al., 2010a). This complexity in exposure and in the response of organisms induces a disproportional impact on the constituent species of a community and therefore many indirect effects (e.g., changes in the diversity of resources and/or consumers) [START_REF] Rohr | Community ecology as a framework for predicting contaminant effects[END_REF]. For instance, a non-sensitive species may still be affected when one of its resource, predator, mutualistic species or competitor is affected by the contaminant [START_REF] Fleeger | Indirect effects of contaminants in aquatic ecosystems[END_REF][START_REF] Clements | Community responses to contaminants: using basic ecological principles to predict ecotoxicological effects[END_REF].

One of the most known type of propagation of indirect effects in food web is the "trophic cascade" where the food chain is disturbed by either a change in the abundance of the highest trophic level (i.e., top-down cascade) or by a change in the resources (i.e., bottom-up cascade) [START_REF] Carpenter | Cascading trophic interactions and lake productivity[END_REF][START_REF] Heath | Understanding patterns and processes in models of trophic cascades[END_REF]. For instance, a top-down cascade initiated by a decrease in predator abundance releases the pressure on its prey which consequently increase in density and subsequently induce more pressure on prey resources. The reverse process of a bottom-up cascade in a three-level food chain, may start by an increase/decrease of prey resources favoring/constraining in first prey density and then its predator [START_REF] Heath | Understanding patterns and processes in models of trophic cascades[END_REF]. Another well-known bottom-up effect is the "paradox of enrichment" where the increase of prey resources (i.e., the carrying capacity) has a destabilizing effect where population dynamics of prey and predators shift from a fixed stable state to an oscillating behavior [START_REF] Rosenzweig | Paradox of enrichment: destabilization of exploitation ecosystems in ecological time[END_REF]. However, recent studies in aquatic systems showed that exposure to contaminants may inhibit the paradox of enrichment by reducing oscillations and driving the system back to a fixed equilibrium [START_REF] Prosnier | Modeling the direct and indirect effects of copper on phytoplankton-zooplankton interactions[END_REF][START_REF] Huang | The impact of environmental toxins on predator-prey dynamics[END_REF].

Many individual-based models for exposure of terrestrial wildlife to contaminants have been proposed (Loos et al., 2010b;[START_REF] Schipper | Contaminant exposure in relation to spatio-temporal variation in diet composition: A case study of the little owl Athene noctua[END_REF][START_REF] Schmitt | An example of population-level risk assessments for small mammals using individual-based population models[END_REF]. Those highly detailled models are very efficient for the management of a specific ecosystem [START_REF] Schmitt | An example of population-level risk assessments for small mammals using individual-based population models[END_REF][START_REF] Deangelis | Individual-based models in ecology after four decades[END_REF], but the wide variety of components and mechanisms induces a black box effect hiding more or less the possibility of an analytical description of critical transitions [START_REF] Gómez-Mourelo | The differential equation counterpart of an individual-based model for yeast population growth[END_REF]. Therefore, in order to achieve a stability analysis of steady state(s) (i.e., namely a bifurcation analysis) [START_REF] Kooi | Sublethal toxic effects in a simple aquatic food chain[END_REF][START_REF] Huang | The impact of environmental toxins on predator-prey dynamics[END_REF], we develop an ODE model including the ecological dynamics of a two-level food chain with a soil compartment, small mammals (rodents and shrews) as preys of the raptor barn owl (Tyto alba) contaminated by the trace metal cadmium (Cd). Presence of Cd occurs in the Earth's crust at relatively low concentration (0.1-0.5 ppm), and Cd is well-known to be highly toxic since non-biodegradable and because of its ability to bind with many organic molecules that distribute it in all or part of the organism [START_REF] Hopkin | Principles of Ecotoxicology, Fourth Edition[END_REF]. The main anthropogenic sources are non-ferrous metal mining and refining, application of phosphate fertilizers and byproducts of battery manufacturer [START_REF] Burger | Assessment and management of risk to wildlife from cadmium[END_REF][START_REF] Faroon | Toxicological profile for cadmium[END_REF].

Mammals are mainly contaminated through oral consumption of contaminated water and food items [START_REF] Smith | Contaminant exposure in terrestrial vertebrates[END_REF]. Moreover, there is a great variety of small mammal diet composition with a gradient from herbivores to carnivores by including omnivores species which can induce various patterns of exposure and responses to environmental contaminants within this mammal group and shape the exposure of their predators. Similarly, while birds may be contaminated by direct unintentional ingestion of the contaminant (human poisoning, chemically treated seeds, polluted water), secondary poisoning through the consumption of contaminated food items is more common [START_REF] Smith | Contaminant exposure in terrestrial vertebrates[END_REF]. Due to their high trophic positions, their scavenging activities, their large spatial living area and long lifespan over which to accumulate contaminants raptor species are highly exposed to persistent and/or bioacumulative contaminants [START_REF] Gómez-Mourelo | The differential equation counterpart of an individual-based model for yeast population growth[END_REF][START_REF] Bustnes | Temporal trends (1986-2005) of essential and non-essential elements in a terrestrial raptor in northern Europe[END_REF][START_REF] Espín | Tracking pan-continental trends in environmental contamination using sentinel raptors-what types of samples should we use?[END_REF].

Based on empirical scenarios and data, our aims are to provide a theoretical model on how external stressors distributed in food chains by trophic cascades modulate the occurrence and strength of regime shifts emerging in polluted ecosystems. We first explore different exposure scenarios of the predator according to trophic position of preys (herbivores, carnivores or omnivores).

Then, we analyze the occurrence and position of regime shifts in the predator-prey dynamics with variations in Cd concentration in soils. Finally, we test the sensitivity of the system to external disturbance inducing trophic cascade effects: variations in prey resources and in intrinsic predator mortality rate.

Methods

Formulation of the general model

Contaminants are commonly measured in concentration per biomass of an individual, [µg.g -1 ] denoted ppm. In the terrestrial ecotoxicological system, we consider a soil with contaminant concentration denoted by C s in [µg.g -1 ]. Cadmium being a persistent contaminant (neither degradation nor metabolization), we assumed C s to be constant since the transfer of contaminant to prey species cannot significantly change the soil concentration (i.e., negligible over the considered period owing to the quantities in respectively the soil and biomass compartments, and included in a balanced uptake/release cycle), and no external inputs of Cd in soils can be expected since the main point source is no longer emitting cadmium. Also, in the system under consideration, several publications showed the same order of magnitude of soil contamination in samplings from 1995 to 2006 (while not exactly on the same locations) [START_REF] Sterckeman | Vertical distribution of Cd, Pb and Zn in soils near smelters in the North of France[END_REF][START_REF] Sterckeman | Assessment of the contamination of cultivated soils by eighteen trace elements around smelters in the North of France[END_REF][START_REF] Pruvot | Heavy metals in soil, crops and grass as a source of human exposure in the former mining areas (6 pp)[END_REF][START_REF] Douay | Contamination of woody habitat soils around a former lead smelter in the North of France[END_REF].

Then, the contaminant is transferred to a barn owl's preys (i.e., small mammal species), which are exposed to the environmental contaminant mainly through ingestion of contaminated foodstuffs (i.e., trophic transfer), and Cd accumulates in their tissues. We denote x the density of the prey species and B x the mean biomass of an individual of that species. The notation C x holds for the mean concentration of the contaminant in one individual, commonly called the body burden of the prey. The growth function of the prey population is a function g x (x, C x ) in [day -1 ] depending on the population density x and the concentration of the contaminant C x . The additional mortality due to pollution is a function µ x (C x ) in [day -1 ]. The dose-response curve µ x (C x ) is defined by a lognormal cumulative distribution function as in Loos et al. (2010a) (see detailled parameterization hereafter, section 2.2.2 and Supporting Information). We denote κ(C s ) the rate of transfer of the contaminant from soil to prey which follow a log-log linear regression (i.e., ln(C x ) = a + b ln(C s ) defined from data, see section 2.2.1). The rate of excretion of the contaminant in the prey population is a constant denoted by k out,x .

The highest level of the food chain is the predator with density z. The predator feeds on prey at a rate Φ(x), which is called the functional response. We assume that predator is only contaminated through this trophic interaction with an absorption coefficient η (i.e., bioavailability). The excretion rate of the contaminant by the predator is denoted k out,z . Previous ecotoxicological model based on ODEs were considered in aquatic system, where every prey species, whatever their trophic level, are exposed to the contaminant present in the water compartment [START_REF] Kooi | Sublethal toxic effects in a simple aquatic food chain[END_REF][START_REF] Bontje | Modelling long-term ecotoxicological effects on an algal population under dynamic nutrient stress[END_REF][START_REF] Huang | The impact of environmental toxins on predator-prey dynamics[END_REF]. In our terrestrial ecotoxicological system, we assume the predator not to be exposed directly to soil contaminant. We call ψ(C z ) the stepwise linear function, equation (4), of reproduction efficiency of the predator depending on contaminant concentration within the predator denoted C z , m the natural death rate of the predator, and µ z (C z ) the linear function, equation ( 5), standing for the additional mortality rate due to the contaminant.

With all the previous assumptions, the model for the transfer of contaminant is given in Supporting Informations S1.1. Also, a classical growth function for rodents is the logistic type growth rate. We then assume that the contaminant affects the intrinsic growth rate: g x (x, C x ) = π(C x )(1 -x/K), with π(C x ) following a log-normal cumulative distribution function (Loos et al., 2010a) (see parameterization in section 2.2.2 and Supporting Information). The carrying capacity K represents the resource of prey and is used as the parameter for bottom-up cascade in the bifurcation analysis. For top-down cascade effect, we used the natural mortality rate of the predator, m. Prey-dependent functional response satisfies the properties: Φ(0) = 0 and Φ (x) > 0, x > 0.

Thus, it is straightforward to express the dynamics as:

                                               dx dt = growth rate xπ(C x ) 1 - x K - additional mortality µ x (C x )x - predation Φ(x)z dC x dt = exposure κ(C s ) -C x     release rate k out,x excretion + π(C x ) 1 - x K populational     dz dt = growth rate zψ(C z )Φ(x) - mortality rate z(µ z (C z ) + m) dC z dt = intake rate η B x B z Φ(x)C x -C z    release rate k out,z excretion + ψ(C z )Φ(x) populational    (1)

Parameterization of transfers and effects of cadmium in a soil-small mammals-raptor food chain

For parameterization, we consider the case of cadmium (Cd) transfer in a trophic chain from soil to small mammals and finally to a raptor feeding on these preys (see Fig. 1). Cadmium does not appear to have biological function in small mammals and birds, but is well-known to be one of the most mobile toxic metals [START_REF] Cooke | Cadmium in small mammals[END_REF][START_REF] Wayland | Cadmium in birds[END_REF]. In this section, we parameterize model (1) to obtain five food chains corresponding to five different small mammal prey species (herbivores: Microtus spp ; omnivores: Apodemus sylvaticus and Myodes glareolus ; carnivores: Sorex spp and Crocidura spp) and one raptor predator, the barn owl (Tyto alba).

Prey exposure to the contaminant

When Cd is ingested through the diet, the absorption is around 1 to 5 % [START_REF] Andersen | Experimental localization of intestinal uptake sites for metals (Cd, Hg, Zn, Se) in vivo in mice[END_REF][START_REF] Schipper | Modeling the influence of environmental heterogeneity on heavy metal exposure concentrations for terrestrial vertebrates in river floodplains[END_REF][START_REF] Cooke | Cadmium in small mammals[END_REF]. Many study use a fixed contaminant absorption efficiency at 2% which is also used in the present study [START_REF] Schipper | Modeling the influence of environmental heterogeneity on heavy metal exposure concentrations for terrestrial vertebrates in river floodplains[END_REF]Loos et al., 2010a). Exposure to cadmium is mainly chronic (a gradual accumulation from diet). The excretion rate of cadmium is long because of the binding to metallothionein (MT) (Fritsch et al., 2010a;[START_REF] Cooke | Cadmium in small mammals[END_REF].

The half-lifes of Cd is in the range of [100 -300] days in rats [START_REF] Cooke | Cadmium in small mammals[END_REF] what would give an excretion rate of: -ln(1/2)/[100 -300] ≈ [0.0023, 0.0069] day -1 . Most of the studies [START_REF] Gorree | BIOMAG: a model for biomagnification in terrestrial food chains. The case of cadmium in the Kempen, The Netherlands[END_REF]Veltman et al., 2007a;[START_REF] Van Den | Effects of soil properties on food web accumulation of heavy metals to the wood mouse Apodemus sylvaticus[END_REF] et al. (2010a)). Soil Cd level ranges from 0.1 ppm to 2402.1 ppm of dry matter based on more than 260 samples of composite 0 -25 cm soil [START_REF] Douay | Contamination of woody habitat soils around a former lead smelter in the North of France[END_REF]Fritsch et al., 2010a). Small mammals were captured in sub-samples of 500×500 squares selected along a soil pollution gradient.

Small mammals were identified at the species level using morphometric, skull and teeth criteria.

Total body concentration in the different small mammal species were calculated from kidney and liver concentrations measured in the animals Veltman et al. (2007a). From those data, we define the function κ(C s ) by a log-log linear regression (i.e., ln(C x ) = a + b ln(C s )) as previously used (Veltman et al., 2007a;Fritsch et al., 2010a) (see Table S1 in Supplementary Information and Fig. 2.A).

Toxicity effect of Cd on reproduction and mortality of prey

The effect of Cd on the small mammals organism are, among others, a reduction of food and water intake, growth depression, renal dysfunction, hypertension, anemia and cancer [START_REF] Cooke | Cadmium in small mammals[END_REF].

According to [START_REF] Świergosz-Kowalewska | Cadmium distribution and toxicity in tissues of small rodents[END_REF], the reproductive organs are more sensitive than other organs and a low doses may affect spermatogenesis and ovulation. In the report of US-EPA (2005), the NOAEL for reproduction and growth is a diet with a concentration of 1.86µg.g -1 .day -1 , that is 10.63 µg.g -1 in the whole body. Using data and methodology of the same report, the LOAEL (Lowest Observed Adverse Effect Level) for reproduction is estimated at 6.27µg.g -1 .day -1 , that is 35.83µg.g -1 which is in the range estimated previously from [START_REF] Shore | The ecotoxicological significance of cadmium intake and residues in terrestrial small mammals[END_REF]. Again, with the data provided by US-EPA ( 2005) we obtain a NOAEL (Non Observed Adverse Effect Level) for mortality of 6.72 µg.g -1 .day -1 (i.e., 38.4 µg.g -1 ), and from the same report, a LOAEL of 9.72 µg.g -1 .day -1 (i.e., 55.54 µg.g -1 ).

We derive dose-response curves using log-normal law (Loos et al., 2010a) defined by:

π(C x ) = r 1 - 1 2 1 + erf ln(C x ) -µ r σ r √ 2 (2)
where the function erf() is the Gauss error function defined as: erf

(x) = 2 √ π x 0 e -t 2 dt.
For parameterization, we use the equivalence NOAEL= EC 10 and LOAEL= EC 30 [START_REF] Greenberg | Regulatory considerations for the potential development and application of metal cleanup values[END_REF], and then a least square regression to estimate parameters of the log-normal curve.

For reproduction, we estimate EC 50 = 82.91µg.g -1 , the parameter of the log-normal curve, equation (3) are µ r = 4.418 and σ r = 1.601 (see Fig. 2.A in Supporting Information). Combination of soil-prey contamination (i.e., κ(C s ) described in Table S1 in Supporting Information) with the dose-response function (i.e., π(C x )) gives the probability of prey to reproduce in function of concentration in soil (see curves plotted in Fig. 2.B in Supporting Information). The equation of the dose-response curve for prey reproduction is given by the cumulative distribution function of a log-normal law: The reproduction rate of small mammals are well-known to be high. From Quéré and Le Louarn (2011), we estimate the maximal growth rate of preys r around [1/365 -25/365] per day, and fixed it at r = 20/365 for every prey species.

For mortality, we infer on a log-normal cumulative function, equation ( 3), with µ mx = 4.272 and σ mx = 0.486 (see Fig. 2.A in Supporting Information).

µ x (C x ) = m x 1 2 1 + erf ln(C x ) -µ mx σ mx √ 2 (3)
The combination of soil-prey transfer with dose-response curve for mortality, µ x (C x ) gives the probability of prey mortality as a function of concentration in soil (see curves plotted in Fig. 2.B in Supporting Information when assuming the maximum of mortality rate due to Cd at m x = 1).

Compared to rodents, the insectivores shrews have often highest exposure to cadmium but they are able to tolerate high concentrations of Cd in their kidney and liver without adverse effect [START_REF] Cooke | Cadmium in small mammals[END_REF], what is confirmed with the data of US-EPA ( 2005) where shrews have highest NOAEL and/or LOAEL for growth and survival than other species. However the set of data was too small to give specifics NOAEL and LOAEL for shrews. Therefore, we gave different maximum rate of mortality due to contamination: m x = 2/365 for rodents (Microtus, Apodemus and Myodes), and m x = 0.5/365 for shrews (Crocidura and Sorex).

Exposure of raptors to cadmium

As for small mammals, birds accumulate Cd mainly through their diet [START_REF] Wayland | Cadmium in birds[END_REF][START_REF] Burger | Assessment and management of risk to wildlife from cadmium[END_REF]. While there are some works on avian toxicity of Cd, few concern raptors. In the review of US-EPA ( 2005), the species are mainly herbivores and insectivores. From these species, US-EPA (2005) estimates the NOAEL for reproduction and growth at 1.47µg.g -1 .day -1 and the LOAEL at 6.35 µg.g -1 .day -1 . For mortality, the NOAEL is estimated at 10.33µg.g -1 .day -1 and the LOAEL around 22.97 µg.g -1 .day -1 .

Also, from [START_REF] Gorree | BIOMAG: a model for biomagnification in terrestrial food chains. The case of cadmium in the Kempen, The Netherlands[END_REF] (on kestrel and barn owl) we assume an excretion rate for bird k out,z = 0.0035 day -1 . From [START_REF] Wayland | Cadmium in birds[END_REF], the absorption rate, denoted η in the model, is in the range of [0.01-0.07], that we fixed at 2% [START_REF] Schipper | Modeling the influence of environmental heterogeneity on heavy metal exposure concentrations for terrestrial vertebrates in river floodplains[END_REF]Loos et al., 2010a). Therefore, to estimate the whole body concentration, we use the same equation as for small mammals, and we obtain for reproduction: NOAEL= 8.4 µg.g -1 and LOAEL= 36.28 µg.g -1 , and for mortality, we have: NOAEL= 59.03 µg.g -1 and LOAEL= 131.26 µg.g -1 .

Toxicity effect of Cd on reproduction and mortality in raptors

The natural mortality rate of the barn owl, denoted m, is based on the usual life span which is about 4 years, so m = 1/(365 × 4) in day -1 . For barn owl dynamics, the bifurcation analysis constraints the use of a log-normal distribution for the energy conversion and the predator's mortality rate.

As a consequence, for the function of energy conversion in predator, we therefore use the following stepwise linear function:

ψ(C z ) = ε max{αC z + 1 , 0} (4)
We estimate α in function of NOAEL and LOAEL. Using a classical linear regression model, we get α = -0.008, and therefore, reproduction is not possible for the threshold concentration

C z > T R = 118.287µg.g -1
, where T R is the threshold concentration such as αT R + 1 = 0. We set the conversion efficiency to ε = 0.001. Indeed, if we assume that 4 to 8 prey are eaten every day (around 70 to 100 grams of prey see Table S1 in Supplementary Information), we obtain a reproduction rate of 0.001 × [4 -10] × 365 = [1.46 -3.65] on an annual basis.

For mortality function, we compute a linear regression based on NOAEL and LOAEL giving:

µ z (C z ) = m z γC z , with γ = 0.002 (5)
3 Results

Bifurcation diagram on soil concentration

A classical summary of the long-term behavior of populations dynamics is with bifurcation diagram.

Fig. 3 represents the steady points of the prey populations for Cd concentration in soil ranging from 0.4 to 2402.1 ppm with carrying capacity fixed at K = 1121.8 (see also Fig. S2.A in Supplementary Information), and the natural mortality rate of predator at m = 10/(4 × 365). When external changes occur, here Cd concentration in soil, C s , ecosystems changes can be in a continuous way or exhibit a sudden discontinuity. On a certain domain of parameters, ecosystems may have several stable points (i.e., bistability for 2 stables points) where stochastic disturbances may drive the system toward a different steady state.

For prey densities (see Fig. 3), we observe three situations: at low concentration, depending on initial condition, prey tends to a prey-only equilibrium E 1 (predator is extirpated), a prey-predator coexistence equilibrium E * (with limit cycle) or a no-species steady point E 0 . Prey-only steady point starts at the carrying capacity and decreases faster and faster with the increase in C s . For the prey-predator coexistence, prey point starts lower and increases slightly with the increase of soil contamination. This point of coexistence is unstable but with a limit cycle, and the amplitude of the limit cycle decreases until the Hopf-bifurcation point, denoted HB. Simultaneously, predator population is either extirpated or decreases slightly (see Fig. 3.A-B). Again, predator oscillates around the coexistence point decreasing with the increase in contaminant. At the coexistence point, oscillation stops at HB but prey still increases and predator decreases until the transcritical bifurcation (TB) at concentration C s,2 where predator are extirpated and only prey-only steady point is feasible. Finally, prey is also extirpated by a TB at concentration C s,1 (see Fig. 3.C)

where no species can survive. For some range of parameter values, we may observe a homoclinic bifurcation when the periodic orbit collides with the prey-only saddle steady point.

To better understand the resilience of the system, the "marble-in-a-cup" diagram in Fig. 3.C represents the stability landscape where valleys are stable zone and top of hills are unstable points.

In this graph, the concentration of Hopf bifurcation, C s,H , and transcritical bifurcation C s,1 and C s,2 are represented.

Cascade effects with different small mammals species

In order to study at the bottom-up and top-down cascade effects on the steady points for the five groups of small mammals, we analyse and compute 2D bifurcation digrams in Fig. 4 and Fig. 5.

Bottom-up cascade effects

The null model of Fig. 4 shows that at a very low carrying capacity, the increase in contaminant changes the system from prey-only to no-species steady point through a transcritical bifurcation.

The resources of preys are sufficient to maintain the prey population but without predation. In this situation, there is no limit cycle since the predator is extirpated. Limit cycle can only appears with coexistence. When the carrying capacity increases, coexistence is possible (prey-only and nospecies are unstable points). Then, prey-only equilibrium is stable, and the predator is extirpated.

Finally, after ≈ 80 ppm, no-species can survive, the soil contamination is too high. The two changes in steady states stability are called transcritical bifurcations (green lines). When the carrying capacity increases, the coexistence point becomes unstable with a limit cycle. This represents the famous "paradox of enrichment" showed by [START_REF] Rosenzweig | Paradox of enrichment: destabilization of exploitation ecosystems in ecological time[END_REF]. However, the increase of contaminant stabilizes the system since the oscillation of the coexistence points disappears at the Hopf bifurcation (red lines). And then, as previsouly, the system undergoes two transcritical bifurcations. The fourth dynamic is the one described previously in Fig. 3 with bistability of coexistence and prey-only steady point. When carrying capacity is low, the ecosystem changes are continuous between steady points. But at high carrying capacity bistability occurs, where discontinuous changes happen at the tipping points represented by fold bifurcations (blue lines).

For other panels of Fig. 4, for every species, the increase in K has a destabilizing effect. At low carrying capacity, there is no limit cycle, and a limit cycle appears in every chain when K increases.

Also, the increase in Cd has a stabilizing effect since we observe a change from oscillations to stable coexistence, prey-only steady point and finally extinction of species. At high carrying capacity, we observe a zone of alternative stable (or limit cycle) states of coexistence and prey-only. This illustrates that the enrichment in prey resources favors the occurrence of regime shifts.

Top-down cascade effects

In Fig. 5, we observe that m has a stabilizing effect. At low predator mortality, every system exhibits oscillation when soil concentration is low. The oscillations disappears with the increase in m. As for bottom-up bifurcation diagrams (Fig. 4), the increase in Cd in soil has a stabilizing effect with the decline of oscillations until their disappearances. When the predator mortality rate is high, the transition from coexistence to prey-only are smoother. Under this condition, the predator fitness is weak and a slight increase in contaminant reducing the fitness drives the predator to extinction. Then at intermediate predator mortality rate, there are tipping points due to bistable zone (see Fig. 5.A-B,F). Finally, at low predator mortality rate, there is no transition from coexistence to prey-only but directly from coexistence to no-species. This last pattern means that the predator survival is not driven by the contaminant but by the presence of prey.

Predator and prey extirpation

The transcritical bifurcation from prey-only to no-species happens at the threshold value C s,1

where the mortality rate of a prey population is greater than its reproduction rate leading to the extinction of prey. Since the prey is extincted, predator cannot survive. For the null model where concentration in prey is the same that the Cd concentration in the soil C x = C s , the threshold is C s,N ull M odel,1 ≈ 74.1 ppm. For the herbivores Microtus spp, C s,M i,1 ≈ 79.5 ppm.

This transcritical bifurcation for omnivores and for Sorex spp is not in the range considered for Cd soil concentration [0.4 -2402.1] ppm. For Crocidura spp C s,Cr,1 ≈ 62.9 ppm.

Discussion

This work provides a set of conditions on stressors (e.g., environmental pollution, excessive enrichment, predator depletion) that may introduce regime shifts in small mammals-raptor dynamics taking place in metal-polluted ecosystems. Along the soil contaminant concentration gradient, the food chain exhibits four types of outcomes: (i) predator-prey coexistence, (ii) prey-only, (iii) population extinction and (iv) a bistable state where a coexistence point and a prey-only steady point form two stable valleys separated by an unstable coexistence area. This last case means that a slight change in populations densities, due to stochastic variations, may lead the system to a totally different state (i.e., from coexisting predator-prey state to a prey-only state where predators are extirpated). Also, we observed oscillations when prey and predator coexist even in a bistable context. These results provide a conceptual framework, with its limitations discussed hereinafter, to modelize and analyze the effect of chemical pollution in the classical trophic cascade theory, in term of stability (i.e., resilience and oscillations) and strength of transition between states (i.e., smooth transition or abrupt change due to tipping point). For environmental risk management, predicting alternative stable states and system sensitivity to external variations is essential to anticipate or prevent abrupt changes in ecosystem [START_REF] Scheffer | Critical transitions in nature and society[END_REF].

Gradient of soil contamination

Within the range of low contamination of soil, the observation of population densities reveals no impact of Cd pollution at the ecosystem level, but the bifurcation diagrams shows that when the contamination increases, the resilience of the system decreases as the system moves closer to preyonly or no-species stable states. Also, when the Cd concentration increases, there is an increase in prey densities simultaneous with a decrease in predator population. The reason is that the decrease in predator due to the contaminant releases the pressure on prey population. The prey population is also impacted by the contaminant but the balance between the release of predation and the stress of contamination is beneficial for the prey. Then, at a higher level of contaminant concentration in the soil, predator population goes extinct. This leads the system to a prey-only steady point where the contaminant concentration is the only external factor affecting prey survival until its extinction. If prey is more sensitive and/or more exposed to the contaminant than the predator, then, a transition from coexisting predator-prey steady state to a total extinction without the transition at a prey-only equilibrium may occur.

Combined effects of trophic cascades and contamination

The transition from predator-prey coexistence to prey-only state is different according to bottomup and top-down pressures: it may change smoothly, strongly around a threshold concentration or exhibits tipping points with two alternative stable states.

Several theoretical and empirical studies predicted that enrichment (i.e., increase in resources) leads to oscillation of population dynamics, what is known as the "paradox of enrichment" [START_REF] Rosenzweig | Paradox of enrichment: destabilization of exploitation ecosystems in ecological time[END_REF][START_REF] Fussmann | Crossing the Hopf bifurcation in a live predator-prey system[END_REF]Roy and Chattopadhyay, 2007b). The system we consider

shows also an increase in oscillation amplitude with the increase in resource enrichment (i.e., carrying capacity of prey). But we also show that the increase in soil contamination overturns oscillating dynamics. Indeed, the transfer of contaminant in the food-chain, affecting the reproduction and mortality of the species involved a stabilizing effect. In the same way, considering a one predator-two preys model, where one prey is toxic while the other is palatable, Roy and Chat-topadhyay (2007a) showed that the introduction of a toxic prey reduces oscillations. In aquatic ecosystem, [START_REF] Kooi | Sublethal toxic effects in a simple aquatic food chain[END_REF]; [START_REF] Prosnier | Modeling the direct and indirect effects of copper on phytoplankton-zooplankton interactions[END_REF]; [START_REF] Huang | The impact of environmental toxins on predator-prey dynamics[END_REF] made the same observations for respectively persistent toxicant-bacterium-ciliate, copper-Scenedesmus-Daphnia and methylmercury-aquatic insects-rainbow trout systems. These models, showing that contaminant may contribute to reduce oscillations in the community dynamics, concern chronic exposure to contaminant in bi-trophic food chain. However, the theoretical work of Garay-Narváez et al. (2013) showed that more complex systems (i.e., with a greater species richness) are destabilized in terms of reduction of food-web connectance, by the pulsed introduction of a contaminant into the environment. Therefore, further works are necessary to understand the different effects between pulsed and chronic exposure, from low to high species richness including variations in both structure and composition of the food web.

We show that a decrease in predator mortality rate, that is a decrease in top-down regulation, has also a destabilizing effect on the predator-prey dynamics. This result means that when topdown pressure effect is high (i.e., high intrinsic predator mortality rate), the predator survival is very sensitive to contaminant variation until a point where contaminant concentration drives the predator population toward extinction. While with a low intrinsic mortality rate of the predator, it is the density of prey that drives predator to extinction, and no more the contaminant. Also, similarily to the paradox of enrichment, at low top-down regulation, contaminant has a stabilizing effect. To our knowledge, the observation of this pattern is new and could be important for environmental management. Indeed, if empirical studies confirm our prediction, we could argue that to prevent the extinction of apex predators, when the ecosystem contamination is high, manager should focus on reducing exposure to contamination first, and then, when contaminant level is lower, care has to be taken on prey population. Our model also shows that the strength of steady state transitions increases with the decrease of the top-down regulation. Again, this illustrates that predator survival is more dependent on prey density and less on its intrinsic mortality rate.

Trophic level of prey

In this study, we consider five groups of prey pooled according to their diet (i.e., herbivores, omnivores and insectivores). We compare these five models to a null model using a null-prey having the same body burden as the concentration in the soil. From field data, the five groups can be clearly distinguished by their responses along the pollution gradient in terms of bioaccumulation.

At low concentration of the contaminant in soil, the herbivores (i.e. Microtus spp) are the least contaminated, and the insectivores (i.e., shrews) are the most contaminated. Cd body burden in the omnivores are in between. This ranking follows previous studies (Veltman et al., 2007a;Fritsch et al., 2010a) Those results reflects previous observations that the body burden of Cd in more related to diet preferences for mobile species as A. sylvaticus rather than soil concentration, and is the opposite for a nonmobile species as Microtus spp [START_REF] Van Den | Cadmium accumulation in small mammals: species traits, soil properties, and spatial habitat use[END_REF]. Finally, as expected,

insectivore have the highest level of contamination but the response to contaminant level variation are not the same between the two groups. While Sorex spp present a high Cd body burden at low concentration in soil, field data did not exhibit an increase of Cd body burden along the pollution gradient. So, as for omnivores, this buffers the impact of contamination when its level increases.

For the other insectivores, Crocidura spp, the response to contaminant level variations is similar to that of Microtus spp species, but the body burden is higher. Therefore, Crocidura spp are the most risky species to consume for a predator. This knowledge about critical prey species is important to assess the exposure of apex predators to environmental contaminants [START_REF] Schipper | Modeling the influence of environmental heterogeneity on heavy metal exposure concentrations for terrestrial vertebrates in river floodplains[END_REF]Loos et al., 2010a). Also, the position of raptor at the top of food webs informs on the bioavailability of a wide range of contaminants and on the exposure of species at lower trophic levels and is therefore an important challenge for environmental risk assessment [START_REF] Bustnes | Temporal trends (1986-2005) of essential and non-essential elements in a terrestrial raptor in northern Europe[END_REF].

Limitations and perspectives

Our model explore the changes in the balance between the necessity for a predator to feed on prey and the damaging consequences of predation in terms of exposure to contaminants. This modeling of prey-predator relationships under a panel of pollution stress scenarios allows not only predicting long-term direct and indirect effects of contaminants on populations, but also, it shows how extinctions may occur within communities free-living in polluted sites. The model reveals the occurrences and the positions of bifurcation points where a small change in parameters (i.e., soil contaminant concentration, prey carrying capacity, predator mortality rate) induces a sudden change in the food chain dynamics. The analysis of these bifurcations helps to identify which prey ecological traits are critical for propagating adverse effects across trophic levels, how a whole community is sensitive to contaminant. It also provides inputs for concepts about how contaminants modulate trophic cascades, and further how multiple stressors can shape population dynamics and tipping points in ecosystem functioning.

However, in this simplified model, the complexity of trophic interactions is reduced to only one prey-one predator system even if different functional responses have been considered. Yet, a variation in the diet of a species is likely to change the exposure of that species to contaminants.

The direct next step is therefore the inclusion of multi-prey predation to explore the influence of the adaptive foraging behavior of the predator [START_REF] Baudrot | The adaptation of generalist predators diet in a multi-prey context: insights from new functional responses[END_REF] on similar ecotoxicological systems, and the role of biodiversity in prey community on shaping the responses of food-webs to environmental disturbances such as environmental pollution. Among many complexification of the models, since collection of data have been thought into a landscape approach, we suggest to take this spatial heterogeneity of habitats into account [START_REF] Fritsch | Influence of landscape composition and diversity on contaminant flux in terrestrial food webs: A case study of trace metal transfer to European blackbirds Turdus merula[END_REF].

dC x dt = d w x /B T x dt = 1 B 2 T x dw x dt B T x -w x dB T x dt = 1 B 2 T z (B T x (κ(C s )B T x + w x (w x γ x g x (x, C x ) -k out,x - µ x (C x )) -C x Φ(x)z) -w z B T z (ψ(C z )Φ(x) -m -µ z (C z ))) = κ(C s ) -k out,x C x -C x µ x (C x ) - C x B x B T x Φ(x)z- C x g x (x, C x ) + C x µ x (C x ) + C x B x B T x Φ(x)z = κ(C s ) -C x (k out,x + g x (x, C x )) (S3)
and for the predator:

dC z dt = 1 B 2 T z dw z dt B T z -w z dB T z dt = 1 B 2 T z (B T z (-k out,z w z + w z (ψ(C z )Φ(x)γ z -m -µ z (C z ))+ ηC x B x Φ(x)z) -w z B T z (ψ(C z )Φ(x) -m -µ z (C z ))) = -k out,z C z -C z (m + µ z (C z )) + ηC x Φ(x) B x B z - C z (ψ(C z )Φ(x) -m -µ z (C z )) = ηC x Φ(x) Bx Bz -C z (k out,z ψ(C z )Φ(x)) (S4)

S1.3 Reduction of system using slow-fast dynamics

With ξ the low parameter, system (1) becomes:

                         dx dt = x π(C x ) 1 - x K -µ x (C x ) -Φ(x)z ξ dC x dt = κ(C s ) -C x k out,x + ξπ(C x ) 1 - x K dz dt = z (ψ(C z )Φ(x) -µ z (C z ) -m) ξ dC z dt = η B x B z Φ(x)C x -C z (k out,z + ξψ(C z )Φ(x)) (S5)
and letting ξ → 0, we obtain:

⇒                            dx dt = x π(C x ) 1 - x K -µ x (C x ) -Φ(x)z C x = κ(C s ) k out,x dz dt = z (ψ(C z )Φ(x) -µ z (C z ) -m) C z = ηB x k out,z B z Φ(x)C x (S6)

S1.4 Elasticity of the regression

Note that we have:

ln(C x ) = b ln(C s )+a ⇔ C x = e a C b s ⇔ ∂C x ∂C s = e a bC b-1 s ⇔ b = C s C x × ∂C x ∂C s (S7)
Therefore, the parameter b is the elasticity of our model, that is a measure of the proportional change in C x in response to a proportional change in C s .

Biomagnification:

The reduction formulation of system (S9) gives the following coefficient of biomagnification (BMF ): S1: Log-log linear model of internal concentration in small mammals depending on soil concentrations (see Figure 2). In the log-log linear regression, ln(C x ) = a + b ln(C s ), a is the intercept and b is the elasticity of the model (see Supporting Informations S1.4), which is a measure of proportional change in prey concentration, C x , in response to a proportional change in soil concentration C s .

BM F = η B x k out,z B z Φ(x) = 0.

S1.6 Description of a 2D bifurcation diagram

S1.7 Mathematical analysis

The ecotoxicological model ( 1) is a 4-dimensional system resulting from the coupling of contaminant and population dynamics. This complexity involving a great number of variables and parameters challenges the mathematical stability analysis. We therefore take advantage of the two time scales included in the system, to reduce the model (1) to a 2-dimensional system where fast events are neglected [START_REF] Auger | Methods of aggregation of variables in population dynamics[END_REF]. In our context, we consider the common assumption [START_REF] Kooi | Sublethal toxic effects in a simple aquatic food chain[END_REF][START_REF] Huang | The impact of environmental toxins on predator-prey dynamics[END_REF] that the contaminant dynamics are faster than the dynamic of population growth. This assumption means that the time scale of population growth and the efficiency to convert ingested prey into new-born are much slower than the intake and elimination of contaminants and which approaches a quasi-steady state. Due to this slow-fast dynamics process, we assume a constant concentration rate in populations [START_REF] Kooi | Sublethal toxic effects in a simple aquatic food chain[END_REF][START_REF] Huang | The impact of environmental toxins on predator-prey dynamics[END_REF]. C s = C s,2 where , C s,2 is a positive root of the equation K 1 -µx (Cs) π(Cs)

= x * , that is when

x 1 = x * (see appendix S1.10 for details).

Steady points

Conditions of existence and stability E 0 = (0, 0) Always feasible. Stable when C s > C s,1 , that is when only E 0 is feasible. Saddle for C s < C s,1 , and transcritical bifurcation with E 1 at C s,1 .

E 1 = (x 1 , 0) Feasible for C s < C s,1 . Unstable when E * is feasible, that is for C s < C s,2 . Stable for C s,2 < C s < C s,1 . E * = (x * , z * ) Feasible for C s < C s,2 . See text.
Table S2: Steady points with conditions of existence and stability.

S1.8.2 Stability of the steady states

The analysis of stability of the steady points E 0 , E 1 and E * is totally given in the appendix S1.11.

Here, we describe the mains results.

From the Jacobian matrix for the system (S10) around the trivial equilibrium point E 0 , one eigenvalue is always negative and

E 0 is stable if π(C s ) < µ x (C s ) and is a saddle-point if π(C s ) > µ x (C s ).
The first eigenvalue becomes zero when C s = C s,1 and the model under consideration undergoes first transcritical bifurcation. The transversality condition has to be satisfied, that is the real part of the eigenvalues have to be non-zero.

For the second steady point E 1 where only the predator is extirpated, we know this point is feasible when µ x (C s ) < π(C s ) and hence one eigenvalue is always negative. Then, using the condition of feasibility of E * , the point E 1 is locally asymptotically stable whenever E * is not feasible, that is when C s > C s,2 . Conversely, E 1 is a saddle-point when E * is feasible.

For the coexisting point E * , since the determinant of the Jacobian matrix of system (S10), denoted J * , is positive at this point, the stability of E * is determined by the trace of the matrix as whenever E * is feasible. Consequently, E * is stable if T r(J * ) < 0, is unstable if T r(J * ) > 0 and undergoes Hopf-bifurcation if T r(J * ) = 0. We call C s,H the Hopf-bifurcation threshold. With

Holling type II functional response, the sign of T r(J * ) is the same as:

-2ahx * π(C s ) + ahK(π(C s ) -µ x (C s )) -π(C s ) (S12)

S1.9 Explicit formulation of threshold concentration for predator-prey coexistence

For the predator dynamics, substituting equations ( 4) and (5) into system (S10), we obtain:

dz dt = z εα ηB x k out,z B z Φ(x) 2 C x + εΦ(x) -m z γ ηB x k out,z B z Φ(x)C x -m (S13)
Then, as demonstrated in appendix S1.12, E * has two solutions for C s < C s,3 , denoted E * ,a and E * ,b , one solution for C s = C s,3 and no solution otherwise. The threshold contamination C s,3 is defined by: the point E 1 at different level of contamination (case of the Fig. 3). When C s,3 ≤ C s,2 , the point coexisting E * appears at C s,3 and then is divided into the two possible points.

C s,3 = e -a/b × k out,z B z ε ηB x (m z γ) 2 m z γ -2αm ± 2 αm(αm -m z γ) 1/b (S14) If C s > C s,
However, the last condition required for the feasibility of E * ,{a,b} is that x i, * ,{a,b} < x 1 . Therefore, under the previous necessary conditions for coexistence, there are three possible situations: (1)

ψ(C s , x * )Φ(x * ) -µ z (C s , x * ) -m = 0 (S15)
and z * is:

z * = x * Φ(x * ) π(C s ) 1 - x * K -µ x (C s ) (S16)
Therefore, if x * is a positive root of (S15), then z * > 0 is positive when x * < K 1 -µx(Cs) π(Cs) .

Feasibility condition for E * is given by 0

< x * < K 1 - µ x (C s ) π(C s ) := x 1 with µ x (C x ) < π(C s ) (S17)
The interior equilibrium point E * bifurcates from the axial equilibrium point E 1 at the bifur-

cation threshold C s = C s,2 where , C s,2 is a positive root of the equation K 1 -µx(Cs) π(Cs) = x * , that is when x 1 = x * . Note also that x * < x 1 when C s < C s,2 , otherwise, x * > x 1 but in this
situation E * is no more feasible. One requires x * > 0 and functions π(C s ) is a decreasing function and µ x (C s ) is an increasing function. Therefore, we have C s,2 < C s,1 . Now we can find the local asymptotic stability conditions for various equilibrium of the system (S10).

S1.11 Stability of equilibria

Stability conditions for point E 0 = (0, 0)

First we calculate the Jacobian matrix for the system (S10) around the trivial equilibrium point E 0 and find:

J 0 =   π(C s ) -µ x (C s ) 0 0 -m   (S18)
One eigenvalue is always negative and

E 0 is stable if π(C s ) < µ x (C s ) and is a saddle-point if π(C s ) > µ x (C s ).
The first eigenvalue becomes zero when C s = C s,1 and the model under consideration undergoes first transcritical bifurcation. The transversality condition has to be satisfied, that is the real part of the eigenvalues have to be non-zero.

Stability conditions for point E

1 = (x 1 , 0)
Now we evaluate the Jacobian matrix for the system (S10) around E 1 to obtain:

J 1 =   -π(C s ) + µ x (C s ) -Φ(x 1 ) = -Φ K 1 -µx(Cs) π(Cs) 0 ψ(C s , x 1 )Φ(x 1 ) -µ z (C s , x 1 ) -m   (S19)
The eigenvalues of this matrix are given by the diagonal elements of the above matrix. For feasible existence of E 1 , we need µ x (C s ) < π(C s ) and hence the first eigenvalue is negative.

Next we know that E * is not feasible for x * > K 1 -µx(Cs) π(Cs) . Then using the relation

that x * is the positive root of the equation ψ(C s , x * )Φ(x * ) -µ z (C s , x * ) -m = 0, then we have x * = Φ -1 µ z (C s , x * ) + d ψ(C s , x * )
, and we can derive that E * is not feasible for:

Φ -1 µ z (C s , x * ) + d ψ(C s , x * ) > K 1 - µ x (C s ) π(C s ) (S20)
Since Φ(.) is an increasing function and ψ(.) > 0, the last inequality is equivalent to:

µ z (C s , x * ) + m > ψ(C s , x * ) × Φ K 1 - µ x (C s ) π(C s ) (S21)
And therefore, we can derive that ψ(C s , x * )Φ(x * ) -µ z (C s , x 1 ) -m < 0 whenever E * is not feasible.

Hence E 1 is locally asymptotically stable whenever E * is not feasible, when

C s > C s,2 . If the condition ψ(C s , x 1 )Φ(x 1 ) = µ z (C s , x * ) + m is satisfied then one eigenvalue of J 1 is
zero and E 1 loses stability through another Transcritical bifurcation. Proceeding with the similar argument we can prove that E 1 is a saddle-point if E * is feasible.

Stability conditions for point E

* = (x * , z * )
Finally we consider the local asymptotic stability of the coexisting equilibrium point E * . Evaluating the Jacobian matrix around E * we find:

J * =    π(Cs) - 2π(Cs)x * K -µx(Cs) -z * ∂Φ(x) ∂x x * -Φ(x * ) z * - ∂µz(Cs, x) ∂x x=x * + ψ(Cs, x * ) × ∂Φ(x) ∂x x=x * + Φ(x * ) ∂ψ(Cs, x) ∂x x=x * 0    (S22)
Since the determinant is positive (i.e., Det(J * ) > 0 ), the stability of E * is determined by the trace of the matrix J * as whenever E * is feasible. Consequently, E * is stable if T r(J * ) < 0, is unstable if T r(J * ) > 0 and undergoes Hopf-bifurcation if T r(J * ) = 0. Solving the equation T r(J * ) = 0 for C s , we find the Hopf-bifurcation threshold C s,H . Then satisfaction of the condition . Therefore, the trace is equivalent to:

T r(J * ) = ∂F i ∂x -z i * ∂Φ(x) ∂x = ∂F i ∂x - F i (x * ) Φ(x * ) × ∂Φ(x) ∂x = Φ(x * ) × ∂ ∂x F i (x) Φ(x) (S23)
Since Φ(x * ) > 0 when E * = (x * , z * ) is feasible, then the sign of T r(J * ) is the sign of

∂ ∂x F i (x) Φ(x) (S24)
The case where ∂ ∂x F i (x) Φ(x) < 0 means that if we move a little from the point E * , then the number of prey eaten is greater than the number of prey produced, and therefore, the system goes back to the steady point E * . However, after a little change, if ∂ ∂x F i (x) Φ(x) > 0, then the number of prey produce will be greater than the number of prey eaten, and the system moves away from the steady point E * . With Holling's type functional response, the sign of T r(J * ) is the sign of: We define φ the roots of the polynomial G i (x, C x ) giving:

εα ηB x k out,z B z C x φ 2 + ε -m z γ ηB x k out,z B z C x φ -m = 0 (S27)
Since εα ηB x k out,z B z C x < 0, to have solution for the functional response, φ, we have to satisfy the positivity of the discriminant we denote ∆ ≥ 0 with:

∆ := ε -m z γ ηB x k out,z B z C x 2 + 4εα ηB x k out,z B z C x m > 0 (S28)
The roots of ∆ = 0 according to C x are:

C * i,
φ i,1 = k out,z B z 2εαC x ηB x γm z C x ηB x k out,z B z -ε - √ ∆ and φ i,2 = k out,z B z 2εαC x ηB x γm z C x ηB x k out,z B z -ε + √ ∆ (S31)
In the expression of ∆, we have εαD i C x m < 0, so there are always two positive solutions, and we have:

ε - γm z C x ηB x k out,z B z - ε -m z γ ηB x k out,z B z C x 2 + 4εα ηB x k out,z B z C x m > 0 ⇔ ε - γm z C x ηB x k out,
C * i,c = k out,z B z ε (m z γ) 2 ηB x γm z -2αm -2 αd(αm -γm z ) (S35)
Also, we previously found the threshold values:

C * i,a = ε D(m z γ) 2 m z γ -2αm -2 αm(αm -m z γ) (S36)
and

C * i,b = ε D(m z γ) 2 m z γ -2αm + 2 αm(αm -m z γ) (S37)
As a consequence, we define C i,3 = min{C i,a , C i,b , C i,c } = C i,c = C i,a , which is the threshold value of prey concentration for the existence of E * .

We defined the soil concentration threshold C s,3 the point: As a consequence, if C s < C s,3 , then G i (x, C x ) > 0 has two positive solutions on Φ(x) denoted φ i,1 and φ i,2 . And if C s > C s,3 , then G i (x, C x ) has no solution.

C * i,

S1.13 Positivity of the determinant

Fig. 4

 4 Fig.4represents bottom-up cascade effect by an increase in prey resources parameterized with the carrying capacity K ranging from 0 to 5000. The predator mortality rate is fixed at m = 1/(4×365).

Fig. 5

 5 Fig. 5 represents top-down cascade effect by an increase in the predator mortality rates, m, ranging from 0 to 20/(4×365). The carrying capacity is fixed at K = 2500. In both Fig. 4 and Fig. 5, panels A are the null models where the concentration in the prey is the concentration is the soil C x = C s with a maximal intrinsic prey mortality rate fixed at m x = 2/365. For a better distinction of bifurcation zones, for each panel of Fig. 4 and Fig. 5 the y-axes (i.e., carrying capacity or predator death rate) are square-root scaled and the x-axis (i.e., concentration in soil) are log-transformed in the range 0.4 to 2402.1 ppm as for data representation in Fig. 2.A.

Figure 1 :

 1 Figure 1: Flow diagram of the transmission of a contaminant along a trophic chain. The plain line is the trophic interaction, dashed lines are transfers of the contaminant and dotted lines are population intrinsic dynamics. Grey elements are bottom-up and top-down regulation parameterized with respectively the carrying capacity K and the predator mortality rate m.

Figure 2 :

 2 Figure 2: (A) Log-log linear model with 95% credible interval from data from Metaleurop. Whole body concentration in fresh weight and concentration in soil (dry weight). See equation in Table S1 in Supplementary Information. (B) Rate of reproduction with r = 2/365 day -1 , and rate of mortality with m M i,Ap,M y = 2/365 (maximal reproduction rate r) and m Cr,So = 0.5/365 (1/4 of maximal reproduction rate to have reproduction greater than mortality at C s = 0 for Sorex spp).

  Figure 3: 1D-Bifurcation diagram illustrating smoothly changes and tipping points in the null model (prey Cd-concentration is equal to soil Cd-concentration C x = C s ). TB: transcritical bifurcation. HB: Hopf Bifurcation. FB: Fold bifurcation. (A-B) Bifurcation diagram of predator and prey densities with respect to contaminant level in soil C s . (C) "Marble-in-a-cup" representation of bifurcation diagram in the range [45 -80] ppm. The marble rolls towards the bottom of valleys (stable steady points) moving away from top of hills (unstable steady points).

Figure 4 :Figure 5 :

 45 Figure 4: 2D-Bifurcation diagram of bottom-up forcing for the five preys. (A) Null model (B) Microtus spp, (C) Apodemus sylvaticus, (D) Myodes glareolus, (E) Sorex spp and (F) Crocidura spp.

Figure

  Figure S1: (A) Log-normal probabilities of reproduction and mortality for prey depending on the concentration in the species individuals. (B) Log-normal probabilities of reproduction and mortality depending on the concentration in the soil. We use same parametirization than in graph A.

Figure

  Figure S2: (A) 2D-Bifurcation diagram illustrating smoothly changes and tipping points in the null model (prey concentration is equal to soil concentration C x = C s ). This diagram results from 500 1D-bifurcation diagram along a range of carrying capacity K ∈ [0, 5000]. (B) 1D-Bifurcation diagram showing a continuous change in ecosystem when the concentration in the soil increase.

  that system undergoes Hopf-bifurcation at C s = C s,H . We need to keep in mind that the occurrence of Hopf-bifurcation is related to the existence of a positive root of the equation T r(J * ) = 0 in terms of C s , apart from the satisfaction of transversality condition.Threshold for Hopf-Bifurcation: C i,H Let F i (x) = x π(C s ) 1 -x K -µ x (C s ) . At E * , we have z i * = F i (x * ) Φ(x * )

  Holling type I:= -π(C s ) ≤ 0Holling type II:= -2ahx * π(C s ) + ahK(π(C s ) -µ x (C s )) -π(C s ) Holling type III: = -2ahx 3 * π(C s ) + (K -ahKx 2 * )(µ x (C s ) -π(C s )) ψ(C s )Φ(x) -µ z (C z ) -m) = z (ε(αC z + 1)Φ(x) -m z γC z -m) = z ε α ηB x k out,z B z Φ(x)C x + 1 Φ(x) -m z γ ηB x k out,z B z Φ(x)C x -m = z εα ηB x k out,z B z Φ(x) 2 C x + εΦ(x) -m z γ ηB x k out,z B z Φ(x)C x -m := zG i (x, C x ) (S26)

  a = εk out,zBz ηB x (m z γ) 2 m z γ -2αm -2 αm(αm -m z γ) ηB x (m z γ) 2 m z γ -2αm + 2 αm(αm -m z γ) (S30)Note that the part under the square root is positive since α < 0.Now, since mγ ηBx k out,z B z 2 > 0, then ∆ > 0 for C x ∈]C * i,a , C * i,b [.The solution of the equation (S27) are:

  3 = e a × C b s,3 ⇔ C s,3 = e -a/b (C * i,3 ) 1/b (S38)

  use an excretion rate based on[START_REF] Friberg | Cadmium and health: a toxicological and epidemiological appraisal[END_REF] which estimates the half-time for the whole body to be around 200 days, which gives an excretion rate of k out,x = 0.0035 day -1 .

	To parameterize the function of transfer from soil to prey, denoted κ(C s ), we use data collected
	from 2006 to 2010 in a 40 km 2 area surrounding the former smelter "Metaleurop-Nord", north of
	France (see publications Fritsch et al. (2010b); Fritsch et al. (2011); Tête et al. (2013); Fritsch

  cascades, the bifurcation diagram of Microtus spp is similar with the null model. This similarity is due to the elasticity of exposure close to one (1.56, see TableS1 in SupplementaryInformation), meaning that the Cd body burden in Microtus spp varies with the similar amplitude that soil concentration. Omnivores rodents, A. sylvaticus and M. glareolus, are the least sensitive to soil contamination since they have also a positive response, but this response diminishes with the increase of Cd in soil. Therefore, at the food-chain level, omnivores preys have a buffer effect.

suggesting a biomagnification of Cd in food chain, where the species at the highest trophic level is the most contaminated. However, when Cd increases in soil, responses of the various groups are different: Microtus spp are the most sensitive, with a positive response which is accentuated with an increase in soil contamination. As a consequence, for bottom-up and top-down

Table 1 :

 1 Definition of parameters and their values. References are[a]

	Symbols Definitions	Values [References]
	π(C x )	Contaminant effect on prey reproduction	equation (2)
	r	Maximal reproduction rate of prey	2/365 [a]
	µ x (C x )	Mortality in prey due to the contaminant	equation (3)
	m x	Maximal prey mortality due to contaminant	2/365 rodents, 0.5/365 shrews
	K	Carrying capacity in prey	[0 -5000]
	κ(C s )	Function of prey exposure	data in Table S1 in Sup. Inf. [b]
	k out,x	Excretion rate of contaminant in prey	0.0035 [c]
	B x	Mean biomass of prey individual	data in Table S1 in Sup. Inf.
	Φ(x)	Ingestion rate of prey: Holling III	a = 0.003 and h = 0.3
	ψ(C z )	Contaminant effect on predator reproduction see text
	ε	Maximal predator reproduction efficiency	0.001
	µ z (C z )	Contaminant effect on predator mortality	1/100
	η	Percentage of absorption of ingested Cd	0.02 [d]
	k out,z	Excretion rate of contaminant in predator	0.0035 [c]
	m	Natural mortality rate of predator	1/(4 × 365)
	B z	Mean biomass of individual predator	350 grams

Parameterization for probability of reproduction and mortality

		02 × [70 -100] 0.0035 × 350	= [1.14 -1.63] , for 70 to 100 grams ingested daily
				(S8)
	S1.5 Species	Mean individual biomass	Function of ln(C x ):
		B x (Std. Error)	Estimate (Std. Error)
	Microtus spp	16.33 (1.03)		1.56 (0.33) × ln(C s ) -2.52 (0.72)
	Apodemus sylvaticus 17.28 (0.13)		0.23 (0.03) × ln(C s ) + 0.32 (0.08)
	Myodes glareolus	19.29 (0.27)		0.25 (0.04) × ln(C s ) + 1.83 (0.10)
	Sorex spp	7.06 (0.25)		-0.03 (0.09)×ln(C s )+5.11 (0.14)
	Crocidura spp	9.81 (0.12)		0.65 (0.10) × ln(C s ) + 2.82 (0.23)
	Table

  3 , there is no equilibrium. If C s ≤ C s,3 , then there is two (or one if equal) possible steady points E * appearing at C s,3 . If C s,3 > C s,2 then the two points E * appears when crossing

  If x i, * ,{a,b} < x 1 , then there is two coexistence points (i.e., E * ,a and E * ,b ). (2) If x i, * ,a < x 1 < x i, * ,b ,then there is only one coexistence point (i.e., E * ,a or E * ,b ). (3) Finally when x i, * ,{a,b} > x 1 , then there is no coexistence point.

S1.10 The coexisting point E *

Since E * (x * , z * ) is the coexisting equilibrium point, then x * is the root of:

  γm z C x ηB x k out,z B z > 0 ⇔ C x < k out,z B z m z γηB x ε (S33)Also, let us define the equationP (C x ) := ε-γm z C x ηB x k out,z B z -2 -εα ηB x k out,z B z C x m.The function P is decreasing with C x . Then, solving P (C * i,c ) = 0 we get:

		ε -
	C * i,c =	k out,z B z ε (m z γ) 2 ηB x	γm z -2αm ± 2 αd(αm -γm z )	(S34)
	Therefore, for C x <	k out,z B z ε m z γηB x	, we get:
						(S32)
	z B z	-2 -εα	ηB x k out,z B z	C x m	> 0
	So the previous inequality (S32) implies that:
						S9

  Determinant of the Jacobian matrix at the coexistence steady point E * is:det(J * ) = Φ(x * )z * -∂µ z (C s , x) ∂x x=x * + ψ(C s , x * ) ×

							∂Φ(x) ∂x x=x * + Φ(x * )	∂ψ(C s , x) ∂ x	x=x *
											(S39)
	We have:
					∂Φ(x) ∂x	> 0					(S40)
	and since ψ(.) is an increasing function:
	∂ψ(C s , x) ∂x	=	∂ψ	ηB x k out,z B z ∂x	Φ(x)C x	=	ηB x k out,z B z	C x	∂ψ (Φ(x)) ∂x	> 0	(S41)

Acknowledgment

We are particularly grateful to M. Coeurdassier and R. Scheifler for data collection and constructive discussions. We also express our gratitude to S. Charles, J.-C. Poggiale and B. Laroche for constructive comments on the modelling approach. This work was funded by the Agence Nationale de la Recherche (ANR) and the ADEME within the STARTT program, and by the SNOWMAN network within the framework of the INSPECT program.

https://dataosu.obs-besancon.

S1 Supporting Information

Supporting Information provides all mathematical developments done in the present study.

S1.1 General model

With all the assumptions made in the manuscript (see also description of notations in the manuscript), especially the assumption that individuals of a same population have the same biomass and the same amount of contaminant in their bodies, the model for the transfer of a contaminant is given by:

In the general model (S1), the equations (S1a) and (S1c) present the dynamic of the total biomass of prey and the predator respectively. The equations (S1b) and (S1d) stand for the dynamics of contaminant concentrations within individuals for respectively the prey and the predator, where the death of individuals (prey and predator) implies a decrease of amounts of the contaminant, respectively -w x µ x (C x ) and -w z (m + µ z (C z )). The transfer of the contaminant through the predation is express by the term C x B x Φ(x)z, which is removed in prey (S1b), and added in the predator (with a factor of transfer efficiency η) (S1d).

While we could consider the population dynamic in term of biomasses, it is more convenient to use density of population as suggested by the change of variable presented in equations (S1f) which give directly the following equalities:

S1.2 Dimensional transformation of variables

From the change of variables in equations (S1e), it is straightforward to give the contaminant dynamics for the prey:

Therefore, with ξ the low parameter, and letting ξ → 0, system (1) becomes (see appendix S1.3

for details):

Cadmium excretion is very low compared to other pollutant (k out,{x,z} = 0.0035). Also, from the parameterization done in the previous section (see Table 1), we have π(C x ) ≤ r = 2/365 ≈ 0.0055 and ψ(C z ) ≤ ε = 0.001, and since reproduction rates are decreasing with the increase in contaminant, there is a range of contamination where the slow-fast process is relevant. Also, the prey reproduction is reduced with the saturating part of the logistic growth rate, (1 -x/K).

S1.8 Steady states and their stability

The contamination of prey population and the transmission of contaminating substance through the predation is described in model (S9). It is important to mention here that the properties of π, ψ, µ x and µ z can not be prescribed beforehand since it varies from case to case. From now and onwards, we use the notation π(C s ), µ x (C s ), ψ(C s , x) and µ z (C s , x) instead of respectively π(C x ), µ x (C x ), ψ(C z ) and µ z (C z ). Hence, the change in the system dynamics can be studied with C s as bifurcation parameter (see part S1.1):

S1.8.1 Existence of steady states

Above system possesses two axial equilibrium points, E 0 (0, 0) and

, system undergoes a transcritical bifurcation and the threshold of this bifurcation is given by:

Graphical representation of C s,1 is given by Fig.