Y Reyjol

J.-P Léna

F Hervant

D Pont

Effects of temperature on biological and biochemical indicators of the life-history strategy of bullhead Cottus gobio

Keywords: experimental ecology, freshwater fish, phenotypic plasticity, triglycerides

The biological and biochemical effects of temperature on life-history strategy of female bullhead Cottus gobio were investigated. Fish from two populations (Bez Basin, south-east France) experiencing contrasted thermal environments (i.e. more or less stable) were reared during 4 months at three distinct temperatures (7, 9 or 12 • C). Both somatic (soma fresh mass and muscle triglyceride content) and reproductive (gonad fresh mass, fecundity, mean diameter of eggs and gonad triglyceride content) indicators were examined. Mixed models indicated that an increasing temperature had significant negative effects on all life-history indicators except for soma fresh mass. Differences in life-history strategy with regard to muscle and gonad triglyceride contents, however, suggest that populations experiencing more variable thermal environments may be better adapted than others to cope with an increasing temperature. These findings may have important implications for C. gobio populations, within the context of climate warming.

INTRODUCTION

Because of its large influence on metabolic reactions and physiological processes (e.g. digestive and assimilation efficiencies, oxygen consumption, growth rate and reproduction), temperature is a keystone environmental variable for ectotherms such as fishes [START_REF] Hochachka | Biochemical Adaptation: Mechanism and Process in Physiological Evolution[END_REF]. Considering reproduction, the most common effects of temperature is a modulation in the timing of spawning and hatching [START_REF] Pepin | Time to hatch and larval size in relation to temperature and egg size in Atlantic cod (Gadus morhua)[END_REF]. It has also been shown to influence a variety of reproductive life-history features, among which fecundity [START_REF] Portner | Climate induced temperature effects on growth performance, fecundity and recruitment in marine fish: developing a hypothesis for cause and effect relationships in Atlantic cod (Gadus morhua) and common eelpout (Zoarces viviparus)[END_REF], egg survival [START_REF] Sandström | Temperature effects on spawning and egg development in Eurasian perch[END_REF] and egg diameter [START_REF] King | Effect of elevated summer temperatures on gonadal steroid production, vitellogenesis and egg quality in female Atlantic salmon[END_REF] are the most documented. Within the context of climate warming, it has recently been shown that †Author to whom correspondence should be addressed. Tel.: +1 8193765011; fax: +1 8193765084; email: yorick.reyjol@onema.fr an increase in water temperature of 1 • C could modify the life-history strategy of the species, by accelerating the development of gonads and affecting spawning date for roach Rutilus rutilus (L.) in Lake Geneva [START_REF] Gillet | Effect of temperature changes on the reproductive cycle of roach in Lake Geneva from 1983 to 2001[END_REF]. In rivers, a variety of studies have shown that both recruitment and year-class strength are affected by density-independent mechanisms involving temperature and hydrology [START_REF] Lobón-Cerviá | Population size in stream-living juveniles of lakemigratory brown trout Salmo trutta L.: the importance of stream discharge and temperature[END_REF], which could in turn have significant effects on fish population dynamics and life-history strategies [START_REF] Cattanéo | Relationship between hydrology and cyprinid reproductive success in the Lower Rhône at Montélimar, France[END_REF][START_REF] Daufresne | Population fluctuations, regulation and limitation in stream-living brown trout[END_REF]. Climate warming could have significant effects on fish populations and communities in rivers by affecting both thermal and hydrological regimes [START_REF] Parry | Climate Change 2007: Impacts, Adaptation and Vulnerability. Contribution of Working Group II to the Fourth Assessment Report of the Intergovernmental Panel on Climate Change[END_REF]. An increase in warmwater species abundances parallel to a decrease in coldwater species abundances has been predicted [START_REF] Mohseni | Global warming and potential changes in fish habitat in U.S. streams[END_REF] and recently documented (Daufresne & Boët, 2007), indicating that the consequences of climate warming on riverine fish communities are becoming increasingly obvious [START_REF] Parry | Climate Change 2007: Impacts, Adaptation and Vulnerability. Contribution of Working Group II to the Fourth Assessment Report of the Intergovernmental Panel on Climate Change[END_REF]. Besides, recent studies emphasized that the viability of some fish populations with temperature-dependent sex determination can be compromised through alterations in their sex ratios, as a response to temperature fluctuations of the magnitude predicted by climate change [START_REF] Ospina-Álvarez | Temperature-dependent sex determination in fish revisited: prevalence, a single sex ratio response pattern, and possible effects of climate change[END_REF].

Lipid contents are known to be dependent on temperature conditions [START_REF] Hochachka | Biochemical Adaptation: Mechanism and Process in Physiological Evolution[END_REF][START_REF] Kamler | Resource allocation in yolk-feeding fish[END_REF], and the absorption rate of yolk lipids typically increases with temperature [START_REF] Jaworski | Development of a bioenergetics model for fish embryos and larvae during the yolk feeding period[END_REF]. More precisely, yolk lipids act as the principal source of metabolic organic energy for egg development and participate in embryonic tissue formation in freshwater fishes [START_REF] Kamler | Resource allocation in yolk-feeding fish[END_REF]. The importance of free amino acids as a significant energy substrate in metabolism and during embryogenesis in marine fishes has been suggested by many studies [START_REF] Kamler | Resource allocation in yolk-feeding fish[END_REF], but neutral lipids (mainly triglycerides) are the dominant endogenous energy source during embryonic development in freshwater species, and in numerous marine ones (e.g. 80-90% in gilthead seabream Sparus aurata L. yolk-feeding larvae: [START_REF] Ronnestad | Energy metabolism during development of eggs and larvae of gilthead sea bream (Sparus auratus)[END_REF]. In addition, triglycerides are the major lipid class for energy storage in freshwater fishes [START_REF] Weber | Development and validation of methods for measuring multiple biochemical indices of condition in juvenile fishes[END_REF], and their dynamics are particularly important for fish condition and future population success [START_REF] Adams | Ecological role of lipids in the health and success of fish populations[END_REF]. They are known to play a role as energy reserves (especially during non-feeding and reproductive periods) and as regulators of body density, cellular metabolism, detoxication, behaviour and reproduction [START_REF] Morris | Fish. In Marine Biogenic Lipids, Fats and Oils[END_REF]. Triglycerides deficiency has previously been shown to negatively affect gonad development and fecundity, and inadequate lipid reserves have been related to reduced reproductive effort for several fish species, primarily through reduced fecundity and in altered quality of eggs and larvae [START_REF] Weber | Development and validation of methods for measuring multiple biochemical indices of condition in juvenile fishes[END_REF].

Bullheads (Cottidae, Scorpaeniformes) comprise a family of small benthic fishes found throughout Europe, in high-order streams with cold oxygen-rich water [START_REF] Reyjol | Patterns in species richness and endemism of European freshwater fish[END_REF], where it is usually associated with brown trout Salmo trutta L. [START_REF] Maitland | Guide to Freshwater Fish of Britain and Europe[END_REF]. In France, bullhead Cottus gobio L. mainly dwell in areas with mean annual air temperatures of c.11 • C [START_REF] Pont | Modelling habitat requirement of European fishes: do species have similar responses to local and regional environmental constraints[END_REF][START_REF] Pont | Conséquences du changement global sur la distribution et la démographie du chabot (Cottus gobio). Une approche intégrative[END_REF]). Because of its low mobility and high fidelity to sampling site and its coldwater preferences, and also because it does not benefit from massive restocking operations in contrast to S. trutta, C. gobio is a very reliable target species for investigating the potential effects of climate warming on French coldwater fish populations. The present study was conducted on female C. gobio to (1) document the effects of temperature on biological (soma fresh mass, gonad fresh mass and egg number and diameter) and biochemical (muscle and gonad triglyceride content) indicators of life-history strategy, (2) determine whether related changes in these indicators were affected by natal-stream thermal regimes and (3) investigate whether or not trade-offs were apparent among life-history indicators. Indeed, a trade-off between growth and reproduction, or between fecundity and egg diameter, has been repeatedly documented in the study of life-history theory [START_REF] Stearns | The Evolution of Life Histories[END_REF] and it was therefore interesting to test if these classical forms of trade-offs were identifiable in this experimental work. Moreover, it was also interesting to investigate whether or not a trade-off could be identified between biochemical life-history indicators (i.e. muscle and gonad triglyceride content), as, as far as it is known, it has never documented before for freshwater fishes and could provide a finer resolution of the analysis compared with classical biological indicators. The first hypothesis was that growth and reproduction would be optimal at natal-stream temperatures. The second hypothesis was that a population originating from a more variable thermal environment would be better adapted to an increase in water temperature than a population originating from a relatively stable thermal environment. The third hypothesis was that a trade-off between reproductive and growth features, either biological or biochemical, or between different facets of reproduction (i.e. fecundity and egg diameter) should be apparent here.

MATERIALS AND METHODS

TARGET SPECIES, STUDY AREA AND SAMPLING SITES

The study was carried out in the Bez River basin (France) [Fig. 1(a)]. This basin drains 275 km 2 of relatively undisturbed forest and brush land, ranging in elevation from 479 to 1781 m (mean annual flow 4•5m 3 s -1 ). Two brooks were considered for the present work: Borne and Nonnières [Fig. 1(a)]. These brooks were chosen because they are known to shelter abundant densities of C. gobio (Borne: 2•2-3•8 individuals m 2 and Nonnieres 1•0-2•7 individuals m 2 ; [START_REF] Abdoli | Influence of female age, body size and environmental conditions on annual egg production of the bullhead[END_REF]. Habitat conditions are very close in Borne and Nonnières; the microhabitat characteristics (in terms of water depth, water velocity and substratum composition) are typical of the preferential habitat of C. gobio; reaches with relatively steep slopes and coarse substratum, mainly cobbles and large stones [START_REF] Maitland | Guide to Freshwater Fish of Britain and Europe[END_REF] (Table I). The mean temperature was similar for Borne (8•78 • C) and Nonnières (8•46 • C), but due to the complex hydrology of the Bez Basin (i.e. karstic water input for Nonnières), water temperature was more stable in Nonnières (mean daily temperature: 4

•5-11•5 • C, c.v., = 0•21) than in Borne (mean daily temperature: 0•5-19•2 • C, c.v. = 0•48) [Fig. 1(b)].
The levels of nutritional resources (macro-invertebrates) were close for the two brooks (Table I). Based on the whole of these habitat characteristics, the assumption was made that the variability in thermal regimes constituted the main difference between the Borne and Nonnières Brooks.

EXPERIMENTAL DESIGN

Cottus gobio were collected at two sampling sites [one in each brook; Fig. 1(a)] using electrofishing (DEKA 3000 Lord apparatus; www.efko-gmbh.de), from 17 to 28 November, 2003. Only sexually mature fish were collected (total length, L T , ≥71 mm; [START_REF] Abdoli | Influence of female age, body size and environmental conditions on annual egg production of the bullhead[END_REF]. During the transport to the laboratory, fish were put into two plastic tanks (one for each population) filled with natural water and kept artificially oxygenated. In the laboratory, fish were dispatched into three rows of six plastic tanks (80 × 50 × 50 cm). Within each row, each tank was alternately filled with fish from one of the two populations (Borne, PopB; Nonnières, PopN). Eleven fish were placed in each tank after an acclimatization period of 10 days. A total of 198 fish (11 fish, six tanks and three rearing temperatures) were thus used for the experiments. Tank densities were higher than the densities observed in the brooks, but ensured a reasonable number of samples for statistical analyses. Some plastic half-cylinders (length = 12 cm; diameter = 8 cm) were added in each tank to mimic shelters. More shelters than fish were provided (13 v.11), so that competition for space was negligible. This was verified by daily observation at different times of day, during which no aggressive behaviour among C. gobio was observed. Ad libitum feeding using defrosted chironomids Chironomus sp. was provided daily, following recommendations made by Jeevanthi De [START_REF] De Silva | Production of the common sculpin, Cottus gobio L. in a Belgian chalk stream, the Samson, and the contribution of benthic macroinvertebrate fauna to its diet[END_REF] for the species.

Both the L T and the mass (M, g) were randomly distributed among tanks (one-way ANOVA, both L T and M:F 17,180 , P>0•05). A cooling machine was connected to each row to keep water temperature constant. One row was maintained at 7 • C( T 1 ), one at 9 • C( T 2 )a n d one at 12 • C( T 3 ). The lowest temperature (T 1 ) was chosen to be as close as possible to the mean temperature C. gobio would have experienced in the brooks during the rearing ; www.vemco.com). Each row of tanks worked as a water recirculation system, so that water did not need to be replaced during the whole study duration. Tests of physical and chemical characteristics of the water (conductivity, pH, DOC, nitrates and phosphates) were made weekly to ensure that no significant variation in quality occured. JBL ® tests (http://www.jbl.de/factmanger/index.php) were used for nitrates and phosphates concentrations. As tanks were situated in a greenhouse, a natural cycle of illumination was applied.

QUANTIFICATION OF BIOLOGICAL AND BIOCHEMICAL INDICATORS

At the end of the rearing period (20 April), all fish were killed. The killing procedure took place tank by tank: all fish from a given tank were sampled and placed for 5 min in a0 •5gl -1 tricaine methane sulphonate solution (Sandoz MS-222; www.sandox.com) before dissection. Fish were then weighed, measured and dissected. The two gonads as well as a percutaneous biopsy (50-100 mg) of dorsal muscles (twitch fibres) were taken from each female. One entire gonad (always the left one) was used for counting the eggs, while the other one (always the right one) was used for assays of triglycerides. Both the gonads and muscles tissues used for triglyceride assays were immediately stored in capped vials at -75 • C until the assay. The second gonad was placed in a vial with Gilson fluid in order to harden eggs and dissolve ovarian membranes [START_REF] Abdoli | Influence of female age, body size and environmental conditions on annual egg production of the bullhead[END_REF]. After 1 month in the Gilson fluid, eggs were poured into a Petri dish and maintained in a small volume of water to prevent desiccation. Diameter was determined using a binocular microscope fitted with an integrated ruler. Triglycerides were extracted according to [START_REF] Hervant | Behavioural, physiological and metabolic responses to long-term starvation and refeeding in a blind cave-dwelling Proteus anguinus and a surface-dwelling Euproctus asper salamander[END_REF] and assayed using specific test combinations (Boehringer-Mannheim, Roche Applied science; www.rocheapplied-science.com). Assays were performed in a recording spectrophotometer (Beckman DU-6; www.beckman.com) at 25 • C. The accuracy of each analysis was tested by assaying samples with and without an added internal standard. Triglyceride content was expressed in mg of triglycerides g -1 of fresh mass (mg g -1 ).

S TAT I S T I C A L T R E AT M E N T

ANCOVA was used to examine the effects of temperature (three modalities: 7, 9 or 12 • C) and population (two modalities: PopB or PopN) on allometric indicators of growth and reproduction for female C. gobio. For growth, either the soma fresh mass (M S ; the fish mass before dissection) or the muscle triglyceride content was used as dependent variables, while L T was used as adjustment covariate. This procedure enabled the removal of the effect of fish size on allometric indicators [START_REF] Glazier | Trade-offs between reproductive and somatic (storage) investments in animals: a comparative test of the Van Noordwijk and De Jong model[END_REF]. For reproduction, the gonad fresh mass (M G ), the number (N E ) and mean diameter (D E ) of eggs and the gonad triglyceride content were used as dependent variables. Following [START_REF] Glazier | Trade-offs between reproductive and somatic (storage) investments in animals: a comparative test of the Van Noordwijk and De Jong model[END_REF], the M S was used as adjustment covariate when the M G and N E were predicted, while the muscle triglyceride content was used as adjustment covariate when the relative triglyceride allocation to gonads was examined. The N E was used as adjustment covariate when mean D E was predicted.

Mixed models were applied. For this purpose, the residual variance was allowed to vary according to both temperature and population. The possible correlation of the dependent variable among animals from the same tank was taken into account by introducing the tank as a random effect on the mean value of the dependant variables. All models were structured as follows:

y ij = β 0 + p f =1 β f x fij + (µ 0j + e kij ), (1)
where y ij was the value of the dependent variable for the i-th fish in the j -th tank, x the explicative fixed effects (i.e. temperature, population, adjustment covariate and all possible interactions), β 0 the intercept, β f the fixed effects associated to the i-th fish in the j -th tank, µ 0j the random effect of the tank j on the mean value of y ij conditioned to the fixed effects and e kij the residual error term specific to the group k (one population at one rearing temperature). To improve variance homogeneity, all continuous variables were ln-transformed prior to analysis. For easier interpretation, they were also centred on their grand mean.

For each dependent variable, the highest order model was built, i.e. with all interaction terms between fixed effects and all variance-covariance statistics. Heteroscedasticity according to temperature and population was tested using likelihood ratio tests (LRT), based on the likelihood estimation of variance statistics [START_REF] Goldstein | Multilevel Statistical Models[END_REF]. When non-significant, variance statistics were removed from the model. Restricted maximum likelihood estimates were used to test the significance of the fixed effects. All terms were successively removed from the model, and the final model was selected based on the Akaike information criterion (AIC) and Bayesian information criterion (BIC). Because AIC and BIC gave congruent results in the present case, only models selected using BIC are presented here [START_REF] Burnham | Model Selection and Multimodel Inference: a Pratical Information-theoretic Approach[END_REF]. The significance of each explicative term that remained in the final model was examined using non-sequential F -tests. The denominator d.f. was computed using the Kenward-Roger method. This method is similar to the Satterthwaite estimation, and it uses the covariance statistic to adjust the d.f. according to the level of dependency among observations made on same tank [START_REF] Littell | SAS for Mixed Models[END_REF]. Responses observed in one tank were often found to be uncorrelated (i.e. null or not positive covariance). In such cases, it follows that the denominator d.f. simplify to the residual d.f. When an interaction between the adjustment covariate and temperature or population was detected, local tests (i.e. at the first, second and third quartile of the adjustment covariate) were performed, to examine the effect of the factor with regard to the level of the adjustment covariate. Tests for linear trends in the response of the dependent variables to increasing temperature were also conducted, using orthogonal polynomials to describe the linear and quadratic components in the responses.

Statistical analyses were carried out using SAS (www.sas.com).

RESULTS

At the end of the rearing phase, 169 C. gobio of the initial total of 198 were alive. The post mortem determination of sex indicated that, among these 169 fish, 84 were females and 85 were males. Because mortality during experiment was high for two tanks of PopB at T 2 due to an unidentified reason (no obvious diseases, problem with water quality or problem with the water recirculation system), only females from T 1 and T 3 were used to compare the effects of temperature between populations (n = 53). For the same reason, the analyses investigating the effect of temperature alone were performed on fish from PopN only (n = 49). The total number of females considered for statistical treatment was 76. It was not possible to use all 84 females alive at the end of the experiment because some samples (eight) were unusable for triglycerides assays. The total number of females used (76) does not match the sum of the number of females used in the two separate analyses (53 + 49 = 102), simply because some females were used in both analysis (those belonging to PopN and reared at T 1 and T 3 ).

EFFECTS OF TEMPERATURE ALONE (PO PN)

Soma fresh mass

Temperature did not alter the mean value of M S but significantly altered the allometric relationship between M S and L T , indicating a different effect of temperature on body condition according to size [Table II andFig

Muscle triglyceride content

Muscle triglyceride content varied significantly more at T 2 and T 3 than at T 1 (LRT for a model allowing respectively variance heterogeneity among T 1 ,T 2 and T 3 ,o r variance heterogeneity between T 1 and both T 2 and T 3 : χ 2 ,d .f. = 2, P>0•05; χ 2 , d.f. = 1, P<0•05).

The muscle triglyceride content did not vary with L T , but was significantly altered by temperature (Table II). Contrast analysis indicated a non-linear reduction of muscle triglyceride content with increasing temperature [Table II and Fig. 2(b)]: there was a significant difference between T 3 and T 2 (F 1,9.19 , P<0•001), but not between T 1 and T 2 (F 1,14.4 , P>0•05).

Gonad fresh mass

The relative mass allocation to gonads varied significantly more at T 2 and T 3 than at T 1 (LRT for a model allowing respectively variance heterogeneity among T 1 ,T 2 and T 3 , or variance heterogeneity between T 1 and both T 2 and T 3 :

χ 2 ,d .f. = 2, P<0•01; χ 2 ,d .f. = 1, P<0•01).
Temperature significantly altered the relative mass allocation to gonads (Table III). Contrast analysis indicated that the relative mass allocation to gonads decreased as temperature increased. This trend, however, was non-linear: the reduction was higher between T 2 and T 3 than T 1 and T 2 [Table III and Fig. 3(a)].

Number of eggs

The allometric relationship between N E and M S was significantly more heterogeneous at T 2 and T 3 than at T 1 (LRT for a model allowing respectively variance heterogeneity among T 1 ,T 2 and T 3 , or variance heterogeneity between T 1 and both T 2 and T 3 : χ 2 ,d .f. = 2, P<0•001; χ 2 ,d .f. = 1, P<0•001).

Temperature significantly altered the allometric relationship between N E and M S : N E decreased as temperature increased but this effect decreased as M S increased [Table III and Fig. 3(b)]. This effect partly resulted from higher variance heterogeneity for lighter females according to temperature (i.e. variance heterogeneity resulting from the level of the adjustment covariate which was not taken into account in the model).

Mean egg diameter

Mean D E was significantly more variable at T 3 than at T 1 or T 2 (LRT for a model allowing respectively variance heterogeneity among T 1 ,T 2 and T 3 ,orvariance heterogeneity between T 3 and both T 1 and T 2 : χ 2 ,d.f. = 2, P<0•001; χ 2 ,d.f. = 1, P<0•001).

The analysis indicated that mean D E controlled for the N E was significantly reduced as temperature increased (Table III). Contrast analysis indicated a significant linear trend in the reduction of mean D E with an increasing temperature [Table III and Fig. 3

(c)]. 0⋅6 0⋅4 0⋅2 -0⋅2 -0⋅4 -0⋅8 -0⋅23 -0⋅18 -0⋅13 -0⋅08 -0⋅03 0⋅02 0⋅07 0⋅12 0⋅17 0⋅22 -0⋅6 0⋅0 ln M s In L T 0⋅16 0⋅06 -0⋅04 -0⋅14 -0⋅24 -0⋅34 ln muscle trygliceride content T 1 T 2 T 3
Rearing temperature

Gonad triglyceride content

The gonad triglyceride content varied significantly more at T 2 and T 3 than at T 1 (LRT for a model allowing respectively variance heterogeneity among T 1 ,T 2 and T 3 , or variance heterogeneity between T 1 and both T 2 and T 3 :

χ 2 ,d .f. = 2, P>0•05; χ 2 ,d .f. = 1, P<0•05).
The analysis failed to reveal a relationship between the triglyceride contents of gonads and muscle. Temperature change, however, significantly decreased gonad triglyceride content (Table III). Contrast analysis indicated a linear reduction of the gonad triglyceride content with increasing temperature [Table III and Fig. 3(d)].

EFFECTS OF TEMPERATURE AND POPULATION

The M S , M G , N E and mean D E did not vary according to population. The effects of both temperature and the adjustment covariate were similar to those observed in the analysis with temperature alone, and there were no significant interactions between these variables and the population of origin (Tables IV andV). For clarity, the corresponding illustrations were not represented, but the tendencies are similar to those obtained for temperature alone (Figs 2 and3).

-1 0 0 -1 -1 -2 -2 -2 -3 -3 -3 -4 2 1 0 -1 -2 -3 -4 -5 ln M G ln mean D E ln M S ln M S ln N E ln N E -0⋅8 -0⋅6 -0⋅4 -0⋅2 0 0⋅20 ⋅40 ⋅6 -0⋅8 -0⋅6 -0⋅4 -0⋅2 0⋅0 0⋅20 ⋅40 ⋅6 2 0⋅18 0⋅13 0⋅08 0⋅03 -0⋅02 -0⋅07 -0⋅12 -0⋅17 -0⋅22 ln gonad triglyceride content T 1 T 2 T 3 Rearing temperature (a) (c) (b) (d)
In contrast, muscle triglyceride content and gonad triglyceride content significantly varied with population (Tables IV andV). For muscles, it was higher for PopB than for PopN [Fig. 4(a)]. For gonads, it was higher for PopN than for PopB [Fig. 4(b)].

DISCUSSION

TEMPERATURE EFFECTS ALONE

Normal feeding for C. gobio occurs for a temperature range of 0-26 • C, and the upper lethal temperature has been shown to be >30 • C [START_REF] Elliott | The critical thermal limits for the bullhead, Cottus gobio, from three populations in north-west England[END_REF]. Even if phenotypic plasticity was necessarily present between the populations studied by [START_REF] Elliott | The critical thermal limits for the bullhead, Cottus gobio, from three populations in north-west England[END_REF] and those studied here, there is a high probability that the temperatures applied here were far below the upper limit of normal feeding or death for the species. That being so, the present study showed that an increasing temperature had a strong negative effect on all indicators of somatic and reproductive allocation except for soma fresh mass. The lowest rearing temperature applied was close to the temperature the fish would have experienced in the brooks for the same period (i.e.6•6 • C). Thus, the results agreed with those of [START_REF] Haugen | Population differences in early life-history traits in grayling[END_REF], who showed that populations of grayling Thymallus thymallus (L.) performed better in terms of early life-history indicators (i.e. growth rate and survival during the period of first feeding) at the temperature they normally experience in nature.

The present results also parallel those by [START_REF] Mcmanus | Effects of temperature and salinity on the life history of the sailfin molly (Pisces: Poeciliidae): lipid storage and reproductive allocation[END_REF], who showed that distinct temperatures (i.e.2 3o r2 9 • C) did not affect body mass of male sailfin molly Poecilia latipinna (Lesueur), for four populations raised under controlled laboratory conditions. A causal explanation to this phenomenon remains unclear, but in the present case a hypothesis could be that the duration of the experiment was too short for changes in body mass to become apparent. Numerous studies have examined the effect of temperature on metabolic rates and tissue and organ physiology within lethal limits, showing that performance is usually maximized near the thermal optimum, but drastically decreases with decreasing or increasing temperatures, parallel to an important decrease in enzymatic activity rates [START_REF] Hochachka | Biochemical Adaptation: Mechanism and Process in Physiological Evolution[END_REF]. The relationship between body temperature and metabolism can be described by an asymmetric function, where metabolic rate is maximized at an intermediate temperature: the thermal optimum [START_REF] Angilletta | The evolution of thermal physiology in ectotherms[END_REF]. From an ecological standpoint, this bell-shaped curve can be explained by the principle of allocation [START_REF] Huey | Cost and benefits of lizard thermoregulation[END_REF]: there is a trade-off between the maximum and the range of the performance. An individual that performs well (relatively to other individuals) in its optimal temperature zone should perform relatively poorly at non-optimum temperatures [START_REF] Huey | Is a jack-of-all-temperatures a master of none?[END_REF]. As C. gobio is known to be a coldwater species, it is very probable that the lowest temperature applied here (7 • C) was close to the thermal optimum for fish from these populations, and that the higher temperatures (9 and 12 • C) corresponded with non-optimum temperatures (right side of the bell-shaped curve), reducing metabolism and negatively affecting growth and reproduction.

POPULATION EFFECTS

As far as is known, there are very few studies that have focused on the effects of both temperature and population on ectotherms' life-history strategies, and that include biochemical indicators. One of them is the study by [START_REF] Mcmanus | Effects of temperature and salinity on the life history of the sailfin molly (Pisces: Poeciliidae): lipid storage and reproductive allocation[END_REF], who investigated the effects of temperature and salinity on the lipid storage and reproductive allocation of four natural populations of P. latipinna (male individuals). They found the amount of stored lipids to be altered by temperature variation but the direction and magnitude of the effect substantially varying among populations (a significant interaction between temperature and population was detected). In the present study, no interaction between temperature and population was observed, implying that the two populations reacted similarly to an increase in temperature.

The muscle triglyceride content, however, was higher for the individuals originating from the population of Borne than for those originating from the population of Nonnières, while the opposite was true for gonads (significant effect of the population term but not of the interaction term). In natural conditions, spatial variation in life-history strategy of riverine fish populations has previously been reported for length at maturity [START_REF] Lobón-Cerviá | Environmentally induced spatio-temporal variations in the fecundity of brown trout Salmo trutta L.: trade-offs between egg size and number[END_REF] for S. trutta), fecundity [START_REF] Jonsson | Trade-off between egg mass and egg number in brown trout[END_REF] for S. trutta) and egg diameter [START_REF] Heins | An analysis of factors associated with intraspecific variation in propagule size of the blacktail shiner, Notropis venustus.I nCommunity and Evolutionary Ecology of North American Stream Fishes[END_REF] for blacktail shiner Notropis venustus Girard), but not for lipid contents. As shown in the methods section, the main difference between the habitat characteristics of Borne and Nonnières was the higher variability of the thermal regime in Borne. As fish from Borne live in a more unstable thermal environment than fish from Nonnières, this result agreed with life-history predictions that fishes living in unpredictable environments should be characterized by reproductive traits corresponding to 'pioneering' species, among which low fecundity, small eggs [START_REF] Winemiller | Patterns of life-history diversification in North American fishes: implications for population regulation[END_REF], and, plausibly, low egg quality (i.e. with low energy stores). As far as is known, there are very few studies that have focused on the effects of temperature stability on energy (mainly lipids) allocation. Nevertheless, [START_REF] Hervant | Behavioural, physiological and metabolic responses to long-term starvation and refeeding in a blind cave-dwelling Proteus anguinus and a surface-dwelling Euproctus asper salamander[END_REF] and [START_REF] Hervant | Long-term fasting and realimentation in hypogean and epigean isopods: a proposed adaptive strategy for groundwater organisms[END_REF] proposed a general adaptive strategy for groundwater organisms, including a reduced metabolism, low energetic requirements, low reproductive allocation and large body stores (energy allocation in organs such as muscles, liver in vertebrates). These adaptive responses, favouring investment in somatic energy stores instead of gonadic allocation, might be proposed for numerous organisms as an efficient strategy in harsh and unpredictable environments [START_REF] Hochachka | Biochemical Adaptation: Mechanism and Process in Physiological Evolution[END_REF], and could be considered as a security trade-off allowing to withstand hazardous conditions. The present study clearly showed that muscle and gonad triglyceride contents significantly varied with population and temperature stability. As fish from Borne live in a more unstable thermal environment than fish from Nonnières, the results agreed with this metabolic hypothesis: C. gobio living in unpredictable environments (Borne) should be characterized by both large somatic triglyceride stores (mainly in muscles) and low gonad triglyceride contents.

ABSENCE OF TRADE-OFFS

No trade-off between growth (as measured by soma fresh mass or muscle triglyceride content) and reproduction (as measured by gonad fresh mass, fecundity, mean egg diameter or gonad triglyceride content) was observed. It is likely that under particularly favourable environmental conditions (e.g. ad libitum food and no competition for space or food, as in the present study), trade-offs disappear, because a major limiting factor is no longer present [START_REF] Stearns | The Evolution of Life Histories[END_REF]. Besides, even though negative correlation between growth and fecundity has been repeatedly documented, this is usually in reference to situations in which the effects of reproduction on later growth are measured, a different form of trade-off to the one studied here. This study can be termed an analysis of income trade-off because it balances investment in growth and reproduction during the same season, while other studies can more aptly be termed capital trade-off, balancing investment in one season against investment in a future season [START_REF] Lardner | Growth or reproduction? Resource allocation by female frogs Rana temporaria[END_REF]. Similarly, no consistent trade-off between egg number and egg size was observed, despite that such trade-offs have often been reported in both natural and experimental studies for various taxa (e.g. snakes, fishes and mammals), and usually serve as reference in life-history theory [START_REF] Stearns | The Evolution of Life Histories[END_REF]. In natural environments, variations in egg number and egg size among and within populations of freshwater fishes have been shown to be regulated by local environmental variables, and by latitude and combinations of its covarying factors [START_REF] Lobón-Cerviá | The reproductive tactics of dace in central Siberia: evidence for temperature regulation of the spatio-temporal variability of its life history[END_REF], 1997). Because a trade-off between egg number and mean egg diameter was observed in the field for the two populations studied, the absence of trade-off in this study is probably due to the absence of limitation of food and space.

IMPLICATIONS FOR COTTUS GOBIO POPULATIONS OF CLIMATE WARMING

The lowest rearing temperature was close to the temperature C. gobio would have experienced in the brooks for the same period. Therefore, the results suggested that an increase in water temperature in response to climate warming may negatively affect C. gobio populations of Borne and Nonnières Brooks, as well as other populations from the Bez Basin experiencing similar thermal conditions. For a mean air warming of 3•7 • C in France (Scenario ARPEGE METEO France at 100 years; http://www.cnrm.meteo.fr/gmgec/arpege/arpege.html), a warming in the Drôme catchment (where the two studied brooks are located) of 2•0-2•5 • C during winter and 4•0-4•5 • C during summer has been predicted [START_REF] Pont | Conséquences du changement global sur la distribution et la démographie du chabot (Cottus gobio). Une approche intégrative[END_REF]. According to the present results, this range of increase in air temperature and subsequent corresponding increase in water temperature is enough to negatively affect the C. gobio populations living in these brooks, notably for the Borne Brook where temperature regulation via karstic water inputs does not occur. This phenomenon should be all the more pronounced as C. gobio has restricted movement capacities [START_REF] Knaepkens | Assessment of the movement behaviour of the bullhead Cottus gobio, an endangered European freshwater fish[END_REF] and great difficulties moving upstream of river obstructions such as dams or weirs, even with a low heights (18-20 cm) [START_REF] Utzinger | Effects of environmental parameters on the distribution of bullhead Cottus gobio with particular consideration of the effects of obstruction[END_REF]. The examination of biochemical indicators, however, suggested that populations living in a more variable thermal environment (e.g. the Borne Brook) may be better fitted to cope with increasing temperatures. A monitoring programme is now needed to determine whether or not these slight differences in life-history strategy will allow local populations to maintain themselves.

Thanks are due to Cl. Drevon, A. Abdoli and P. Crouzet for their technical assistance, in the field or in the laboratory, and to C. Rogers for having read a preliminary English version of this manuscript. The authors would also want to gratefully thank G. Ciuperca (Laboratory of Applied Mathematics, Université Claude Bernard-Lyon 1, France) and A. Bertolo (Research Group on Aquatic Ecosystems, Universitéd uQ u ébec à Trois-Rivières, Québec, Canada) for their fruitful and very acute statistical advices with regard to this work, as well as two anonymous reviewers for their comments. This work is part of the Research project 'Potential impact of global warming on distribution and life history strategies of the Bullhead (Cottus gobio L.)', funded by the French Institute for Biodiversity. The experiments carried out for this research work comply with the current French laws.

  Fig. 1. (a) Study area (( ), sampling sites). (b) Illustration of the thermal regime time in the Borne ( )a n d Nonnières ( ) Brooks.

  . 2(a)]. Contrast analyses revealed that as temperature increased significantly, female M S decreased proportionally as the size of the female increased (tests at the first, second and third quartile of the adjustment covariate for a linear component and a quadratic one: F 1,49 , P>0•05; F 1,49 , P<0•01; F 1,49 , P 0•01 and F 1,49 , P>0•05; F 1,49 , P>0•05; F 1,49 , P>0•05).

Fig. 2 .

 2 Fig. 2. Illustrations of the effect of rearing temperature [ ,T 1 (7 • C); ,T 2 (9 • C); ,T 3 (12 • C)] alone on the indicators of relative somatic allocation: (a) fresh soma mass (M S ) and (b) muscle triglyceride content (mean ± s.d.) for the Nonnières Brook population. The dependent variable is either directly represented as a function of rearing temperature (b) or as a function of the adjustment covariate if significant (a) [allometric relationship between M S and total length (L T )] (( ), T 1 ;( ), T 2 ;( )).

Fig. 3 .

 3 Fig. 3. Illustrations of the effect of rearing temperature [T 1 (7 • C), T 2 (9 • C) and T 3 (12 • C)] alone on the indicators of relative reproductive allocation: (a) gonad mass (M G ), (b) number of eggs (N E ), (c) mean egg diameter (D E ) and (d) gonad triglyceride content (mean ± s.d.) for the Nonnières Brook population of Cottus gobio. The dependent variable is either directly represented as a function of rearing temperature (d) or as a function of the adjustment covariate if significant [allometric relationship between M G and fresh soma mass (M S )( a ) ,N E and M S (b) and D E and N E (c)].

Fig. 4 .

 4 Fig. 4. Illustrations of the effect on Cottus gobio of both temperature [T 1 (7 • C) and T 3 (12 • C)] and population [Borne (PopB) and Nonnières (PopN)] on biochemical indicators: triglyceride content of (a) muscle and (b) gonad (mean ± s.d.). The dependent variable is either directly represented as a function of rearing temperature (b) or as a function of the adjustment covariate if significant [allometric relationship between muscle triglyceride content and total length (L T )] (a). (( ), PopN at T 1 ;( ), PopN at T 3 ;( ), PopB at T 1 ;( ), PopB at T 3 ).

Table I .

 I Summary of the main habitat variables for Cottus gobio from the two sampling sites: Borne and Nonnières Brooks Mean river slope (m km -1 )4 •33 •4 Mean water velocity at low flow (cm s -1 )4 3 •34 3 •1 Food resource (benthos) (individuals m -2 ) 8600 9100 * From 1 July 2001 to 30 October 2003 (before fish sampling). c.v., coefficient of variation; D 10 ,D 50 and D 90 , the 10, 50 and 90% quantiles of the grain size distribution.

	Habitat variable	Borne	Nonnières
	Site length (km)	2•31 •5
	Mean elevation (m) Temperature [mean ± s.d.( • C)]*	862 8•78 ± 4•21	679 8•46 ± 1•79
	Temperature (c.v.)0 •48	0•21
	Substratum (D 10 )2 1 •41 9 •9
	Substratum (D 50 )6 6 •15 8 •7
	Substratum (D 90 )	208•3	188•8
	Mean water depth at low flow (cm)	10•21 2 •5
	Mean river width at low flow (m)	3•53 •7

period, i.e.6 •6 • C. The water temperature was checked every hour using MINILOG ® data loggers (VEMCO

Table II .

 II Summary of the models predicting different indicators of the relative somatic investment by Cottus gobio as a function of rearing temperature alone

			M S	Muscle TG
		d.f. n ,d . f . d	FP	d.f. n ,d . f . d	FP
	BIC		-112•1		-40•4
	Covariate	1, 49	722•01 <0•001
	Temperature	2, 49	1•27 >0•05	2, 12•23 6 •63 <0•001
	Covariate × temperature	2, 49	4•23 <0•05
	Linear polynomial	1, 49	1•55 >0•05	1, 16	33•21 <0•001
	Quadratic polynomial	1, 49	0•82 >0•05	1, 12•21 9 •93 <0•001
	M S , fresh soma mass; BIC, Bayesian information criterion; TG, triglycerides; d.f. n , d.f. for the numerator;
	d.f.

d , d.f. for the denominator.

Table III .

 III Summary of the models predicting different indicators of the relative reproductive investment by Cottus gobio , d.f. for the numerator; d.f.

	as a function of rearing	temperature alone	M G N E D E Gonad TG	d.f. FP	BIC 110•1 134 5•6-2 4 •7	Covariate 1, 21•26 9 •7 <0•001 1, 33•51 8 •52 <0•001 1, 33•17 3 •5 <0•001	<0•001 2, 39•11 1 •84 Temperature 2, 41•93 3 •76 <0•001 2, 35•28 •41 0•001 2, 24•96 •85 <0•01	Covariate × temperature 2, 34•24 •13 <0•05	<0•001 1, 45•82 3 •82 Linear polynomial 1, 44•36 6 •58 <0•001 1, 34•71 4 •93 <0•001 7•27 <0•05	>0•05 1, 39•70 •8 Quadratic polynomial 1, 36•65 •07 <0•05 1, 33•50 •76 >0•05 2•27 >0•05	M G , gonad mass; N E , number of eggs; D E , mean egg diameter; BIC, Bayesian information criterion; TG, triglycerides; d.f.

n ,d . f . d FP d.f. n ,d . f . d FP d.f. n ,d . f . d FP d.f. n ,d . f . d n d ,d . f .f o rt h e denominator.

Table IV .

 IV Summary of the models predicting different indicators of the relative somatic investment of Cottus gobio as a function of both rearing temperature and population of origin

			M S		Muscle TG
		d.f. n ,d . f . d	FP	d.f. n ,d . f . d	FP
	BIC		-112•1-2 2 •9
	Covariate	1, 53	733•46 <0•001	1, 50	4•54	<0•05
	Temperature	2, 53	2•75 >0•05	1, 50	104•53	<0•001
	Population			1, 50	32•41	<0•001
	Covariate × temperature	2, 53	6•45 <0•05	1, 50	0•02	>0•05
	Covariate × population			1, 50	5•50	<0•05
	Temperature × population			1, 50	3•03	>0•05
	Covariate × temperature			1, 50	7•41	<0•01
	× population

M S, soma mass; BIC, Bayesian information criterion; TG, triglycerides; d.f. n , d.f. for the numerator; d.f. d, d.f. for the denominator.

Table V .

 V Summary of the models predicting different indicators of the relative reproductive investment of Cottus gobio

	as a function of both	rearing temperature and population of origin	M G N E D E Gonad TG	d.f. n ,d . f . d FP d.f. n ,d . f . d FP	BIC 120•00 93•80 20•80 -33•30	Covariate <0•001 1, 26•83 3 •38 <0•001 1, 27•86 9 •13	82•98 <0•001 1, 26•14 5 •99 <0•001 1, 17•65 •02 <0•05 1, 50 36•89 <0•001 Temperature 1, 28	1, 50 50•08 <0•001 Population	1, 26•11 0 •46 <0•01 1, 33 48•73 <0•001 Covariate × temperature	Covariate × population	Temperature ×	population	Covariate × temperature	× population	d ,d . f .f o rt h e n , d.f. for the numerator; d.f. E , mean egg diameter; BIC, Bayesian information criterion; TG, triglycerides; d.f. E , number of eggs; D M G , gonad mass; N	denominator.

d FP d.f. n ,d . f . d FP d.f. n ,d . f .


J Moreau

Michel Raous
email: raous@lma.cnrs-mrs.fr

Keywords: Unilateral contact, friction, adhesion, non smooth dynamics Friction and Adhesion

The studies carried out on adhesion by the group "Modeling in Contact Mechanics" at the LMA are reviewed in this paper and recent applica-tions are presented. Based on the introduction of the adhesion intensity variable developed by M. Fr´emond, different forms of a model coupling adhesion to unilateral contact and friction have been developed. The formulations are given either under the form of implicit variational in-equalities or the one of complementarity problems. Both quasi-static and dynamic formulations are considered.

The model is non smooth because we do not use any regularization for the unilateral conditions and for the friction, i.e. Signorini conditions and strict Coulomb law are written. In the thermodynamics analysis, the state and the complementarity laws are then written using sub-differentials and differential inclusions because of the non convexity and non differentiability of the potentials. For the dynamics, the formulation is given in term of differential measures in order to deal with the non continuity of the velocities that may occur in the solutions. This work therefore owes much to the theories and the numerical scheme developed by J.

Introduction

In order to describe the smooth transition from a completely adhesive contact to a usual unilateral contact (Signorini conditions) with Coulomb friction, a model based on interface damage has been first developed for quasi-static problems in [START_REF] Raous | Un modèle couplant adhérence et frottement pour le contact unilatéral entre deux solides déformables[END_REF][START_REF] Cangémi | Frottement et adhérence : modèle, traitement numérique et application à l'interface fibre/matrice[END_REF], Raous et al, 1999[START_REF] Raous | Quasistatic Signorini problem with Coulomb friction and coupling to adhesion[END_REF]. Using a dynamic formulation, the model was then extended to account for the brittle behaviour occurring when a crack interacts with fiber-matrix interfaces in composite materials in (Raous-Monerie, 2002[START_REF] Monerie | Fissuration des matériaux composites : rôle de l'interface fibre-matrice[END_REF]. More recently, the model has been used to study metal/concrete interfaces in reinforced concrete in civil engineering (Karray et al, submitted), delamination of coated bodies (Raous et al, 2002), delamination of glued assembling in civil engineering [START_REF] Raous | Modélisation de l'adhésion par collage[END_REF], cohesive masonry [START_REF] Jean | Non-smooth contact dynamics approach of cohesive materials[END_REF][START_REF] Acary | Contribution la modélisation mécanique et numérique des édifices maçonnés[END_REF] and production of wear particles in bio-engineering [START_REF] Baudriller | Approche micro-mécanique et modélisation de l'usure : application aux biomatériaux des prothéses totales de hanche[END_REF].

The quasi-static formulation was extended to deal with hyperelasticity in [START_REF] Bretelle | Unilateral contact with adhesion and friction between two hyperelastic bodies[END_REF]. Mathematical results about the existence of the solutions were given in [START_REF] Cocou | Existence results for unilateral quasistatic contact problems with friction and adhesion[END_REF] without using any regularization on the contact conditions.

1.

The model

The RCC model (Raous-Cangémi-Cocou) has been first given in [START_REF] Raous | Un modèle couplant adhérence et frottement pour le contact unilatéral entre deux solides déformables[END_REF][START_REF] Cangémi | Frottement et adhérence : modèle, traitement numérique et application à l'interface fibre/matrice[END_REF] and then extensively presented in [START_REF] Raous | Consistent model coupling adhesion, friction and unilateral contact[END_REF]. It has been extended to the present form including progressive friction with the term (1 -β) in [START_REF] Monerie | Fissuration des matériaux composites : rôle de l'interface fibre-matrice[END_REF], Raous-Monerie, 2002). Adhesion is characterized in this model by the internal variable β, introduced by Frémond [START_REF] Frémond | Adhérence des solides[END_REF][START_REF] Frémond | Contact with adhesion[END_REF], which denotes the intensity of adhesion. The introduction of a damageable stiffness of the interface ensures a good continuity between the two models during the competition between friction and adhesion. The behaviour of the interface is described by the following relations, where (9.1) gives the unilateral contact with adhesion, (9.2) gives the Coulomb friction with adhesion and (9.3) gives the evolution of the adhesion intensity β. Initially, when the adhesion is complete, the interface is elastic as long as the energy threshold w is not reached. After that, damage of the interface occurs and consequently, on the one hand, the adhesion intensity β and the apparent stiffness β 2 C N and β 2 C T decrease, and on the other hand, friction begins to operate. When the adhesion is completely broken (β = 0), we get the classical Signorini problem with Coulomb friction.

-R r N + β 2 C N u N ≥ 0 ,u N ≥ 0 , -R r N + β 2 C N u N .u N =0, (9.1) R r T = β 2 C T u T ,R r N = R N , R T -R r T ≤µ(1 -β) R N -β 2 C N u N , (9.2) with if R T -R r T <µ(1 -β) R N -β 2 C N u N ⇒ uT =0, if R T -R r T = µ(1 -β) R N -β 2 C N u N ⇒∃ λ ≥ 0, uT = λ(R T -R r T ), β = -w -β (C N u 2 N + C T u T 2 ) -k∆β - b 1/p if β ∈ [0, 1[, β ≤-w -β (C N u 2 N + C T u T 2 ) -k∆β - b 1/p if β
= C T = C, u 0 = ω/C and R 0 = √ ωC) .
It should be noted that the Signorini conditions are strictly imposed when compression occurs. References on other models can be found in [START_REF] Raous | Quasistatic Signorini problem with Coulomb friction and coupling to adhesion[END_REF], Raous et al, 1999) and a comparison between some of them is made in [START_REF] Monerie | Comparaison de lois d'interface fibre/matrice sur la base d'un modèle uni-axial d'expérience de micro-indentation[END_REF] (models developed by Tvergaard-Needleman, Girard-Feyel-Chaboche, Michel-Suquet, Allix-Ladevèze). Using penalization and augmented Lagrangian on a similar model to the RCC one, Talon-Curnier have solved the quasistatic problem using generalized Newton method [START_REF] Talon | A model of adhesion coupled to contact and friction[END_REF].

The thermodynamics

In the framework of continuum thermodynamics, the contact zone is assumed to be a material surface and the local constitutive laws are obtained by choosing two specific forms of the free energy and the dissipation potential associated to the surface. The following thermodynamic variables are introduced: the relative displacements (u N n,u T ) and the adhesion intensity β are chosen as the state variables, and the contact force R and a decohesion force G, as the associated thermodynamic forces. The thermodynamic analysis given for the RCC model in [START_REF] Raous | Consistent model coupling adhesion, friction and unilateral contact[END_REF] has been extended to the present model in [START_REF] Monerie | Fissuration des matériaux composites : rôle de l'interface fibre-matrice[END_REF], Raous-Monerie, 2002) in order to obtain relations (9.1) to (9.3), where friction is progressively introduced in the form of the term (1-β)µ into (9.2).

Expressions (9.4) and (9.5) are adopted for the free energy Ψ(u N ,u T ,β) and the potential of dissipation Φ uT , β . In (9.4), the indicator function I K (where K = {v/v≥ 0} ) imposes the unilateral condition u N ≥ 0 and the indicator function I P (where P = {γ/0 ≤ γ ≤ 1}) imposes the condition β ∈ [0, 1]. In (9.5), the indicator function I C -( β) (where C -= {γ/γ ≤ 0}) imposes that β ≤ 0: the adhesion can only decrease and cannot be regenerated (it is irreversible) in the present model.

Ψ(u N ,u T ,β)= 1 2 β 2 C N u 2 N + 1 2 β 2 C T u T 2 -wβ+ I K (u N )+I P (β) (9.4) Φ uT , β = µ(1-β) R N -β 2 C N u N uT + b p +1 β p+1 + I C -( β) (9.5)
Ψ has a part which is convex but not differentiable and another part which is differentiable but not convex with respect to the pair (u, β). Φ is convex but has a part which is not differentiable. The state laws and the complementarity laws are then written as follows in order to obtain the contact behaviour laws given in section 1 [START_REF] Raous | Consistent model coupling adhesion, friction and unilateral contact[END_REF].

⎧ ⎪ ⎨ ⎪ ⎩ R r T = ∂Ψ d ∂[u T ] R r N -∂Ψ d ∂[u N ] ∈ ∂I K ([u N ]) -G β -∂Ψ d ∂β ∈ ∂I [0,1] (β) (9.6) R N = R r N (R ir T ,G β ) ∈ ∂Φ( [u T ], β) (9.7)

The quasi-static formulation

The formulation and the approximation of quasi-static frictional problems given in [START_REF] Cocou | Formulation and approximation of quasistatic frictional contact[END_REF] have been extended to adhesion problems in [START_REF] Raous | Un modèle couplant adhérence et frottement pour le contact unilatéral entre deux solides déformables[END_REF], Raous et al, 1999[START_REF] Raous | Quasistatic Signorini problem with Coulomb friction and coupling to adhesion[END_REF]. The problem can be here set as the coupling between two variational inequalities (one of which is implicit) and a differential equation.

Problem (P 1 ): Find ( u, β) ∈ W 1,2 (0,T; V ) × W 1,2 (0,T; H) such that u(0) = u 0 ∈ K, β(0) = β 0 ∈ H [0, 1[ and for ∀t ∈ [0,T], u(t) ∈ K, and ∀v ∈ Va ( u, v -˙ u)+j(β, u N ,v T ) -j(β, u N , uT )+ Γ C β 2 C T u T . (v T -uT )ds ≥ ( F,v -˙ u) -R N ,v N -uN , (9.8) -R N ,z-u N + Γ C β 2 C N u N .(z -u N )ds ≥ 0 ∀z ∈ K, (9.9) β = -1/b w -(C N u 2 N + C T u T 2 )β -
a.e. on Γ C , (9.10)

where:

-u =(u 1 ,u 2 ) where u 1 and u 2 define the displacements in Ω 1 and Ω 2 , -

V =(V 1 ,V 2 ), V α = v α ∈ H 1 (Ω α ) 3 ; v α = 0 a.e. on Γ α U ,α=1, 2, -H = L ∞ (Γ c ), -K = v =( v 1 ,v 2 ) ∈ V 1 × V 2 ; v N ≥ 0 a.e. on Γ c
, where Γ c is the contact boundary between the two solids Ω 1 and Ω 2 , -a(.,.) is the bilinear form associated to the elasticity mapping,

-j(β, u N ,v T )= Γ C µ(1 -β) R N -β 2 C N u N v T ds , -F =( F 1 ,F 2 )
are the given force densities applied to solid 1 and to solid 2 respectively. By using an incremental approximation, it has been established in [START_REF] Raous | Un modèle couplant adhérence et frottement pour le contact unilatéral entre deux solides déformables[END_REF], Raous et al, 1999) that the numerical solutions can be obtained by adapting the methods that we have developed for dealing with classical frictional unilateral contacts [START_REF] Raous | Quasistatic Signorini problem with Coulomb friction and coupling to adhesion[END_REF]. The solutions have mainly been obtained as follows:

-either by taking a fixed point on the sliding threshold and solving a sequence of minimization problems with the choice between a projected Gauss-Seidel method (accelerated by relaxation or Aitken processes) and a projected conjugate gradient method [START_REF] Raous | Preconditioned conjugate gradient method for a unilateral problem with friction[END_REF], -or by taking a complementarity formulation, using a mathematical programming method (Lemke's method). With the adhesion model, these solvers are coupled with the numerical integration of the differential equation on β by using θ-methods (θ =1 is often chosen). Implementation of the algorithms has been conducted in our finite element code GYPTIS90 (Latil-Raous, 1991).

The dynamics

Depending on the characteristics of the interface, especially when b tends towards zero (no viscosity for the evolution of the intensity of adhesion), brittle behaviour can be obtained and the inertia effects have to be taken into account. A 3D dynamic formulation has been developed (Raous-Monerie, 2002, Monerie-Raous, 1999, Monerie-Acary, 2001). Because of the non smooth character of these interface laws, the dynamics is written in terms of differential measures as follows (where q denotes the displacement and r the contact force):

M (q, t)d q = F (q, q, t)dt -rdν , (9.11)

where d q is a differential measure associated with q(t):

t 2 t 1 d q = q(t + 2 ) -q(t - 1 ) ∀t 2 >t 1 , (9.12)
and hence:

t 2 t 1 M (q, t)d q = t 2 t 1 F (q, q, t)dt - ]t 1 ,t 2 ]
rdν, (9.13)

q(t 2 )=q(t 1 )+ t 2 t 1 q(τ )dτ . (9.14)
The Non Smooth Contact Dynamics method developed by (Moreau, 1998[START_REF] Moreau | Some numerical methods in multi-body dynamics: application to granular media[END_REF][START_REF] Moreau | Some basics of unilateral dynamics[END_REF][START_REF] Jean | The Non Smooth Contact Dynamics method[END_REF] has been extended to the treatment of the RCC model. A solver for dealing with adhesion has been implemented in the finite element code LMGC [START_REF] Jean | The Non Smooth Contact Dynamics method[END_REF], Monerie-Acary, 2001). Another solver based on complementarity formulation (Lemke's method) dealing with non smooth dynamics problems has been implemented in the finite element code SIMEM3 at the LMA [START_REF] Vola | Consistent time discretization for dynamical frictional contact problems and complementarity techniques[END_REF].

Applications Delamination benchmarks

Various benchmarks for simulating delamination have been developed in the framework of a joint project with the Laboratoire Central des Ponts et Chaussées [START_REF] Raous | Modélisation de l'adhésion par collage[END_REF] focusing on adhesion and gluing in civil engineering. Micro-indentation of a single fiber in a composite material This model was first developed in order to describe the behaviour of the fiber-matrix interface of composite materials [START_REF] Cangémi | Frottement et adhérence : modèle, traitement numérique et application à l'interface fibre/matrice[END_REF], Raous et al, 1999). The parameters of the model were identified on micro-indentation experiments performed at the ONERA on single fibers. Physical and mechanical considerations are taken into account in the identification procedure: for example, the values of the initial interface stiffness C N ,C T are taken to be in the range corresponding to the elastic properties of the oxides located at the interfaces. The viscosity parameter b is particularly difficult to identify: experiments with different loading velocities are required for this purpose. The validity of the model was then confirmed by taking various kinds of loadings (especially cyclic loadings) and comparing experimental results with those obtained in the numer- ical simulation whith the same interface parameters. Details of these studies can be found in [START_REF] Cangémi | Frottement et adhérence : modèle, traitement numérique et application à l'interface fibre/matrice[END_REF], Raous et al, 1999) and later in [START_REF] Monerie | Fissuration des matériaux composites : rôle de l'interface fibre-matrice[END_REF].

Interaction between cracks and fiber/matrix interfaces in composite materials

Again in collaboration with ONERA, during the thesis [START_REF] Monerie | Fissuration des matériaux composites : rôle de l'interface fibre-matrice[END_REF], the RCC model has been used to investigate the different ways in which cracks propagate through a composite material and how they depend on the interface properties. The adhesion model has been used to account both the crack propagation (decohesive crack with no viscosity, i.e. b = 0) and the interface behaviour. Crack bridging, crack trapping and fiber breaking can be observed depending on the interface characteristics (Raous-Monerie, 2002[START_REF] Monerie | Fissuration des matériaux composites : rôle de l'interface fibre-matrice[END_REF]. On Fig. 9.5 adhesion is broken in the black zones.

Metal/concrete interfaces in reinforced concrete

In civil engineering, we are now testing the RCC model for the adhesive contact between steel and concrete in reinforced concrete. Pull-out tests of a steel shaft are being simulated in the framework of a joint project between LMA and ENIT in Tunisia (Karray et al, submitted). In that case, a variable friction coefficient (depending on the sliding displacement) has been used in order to take into account the wear of the surface which occurs during the sliding which seems to be quite significant with concrete. Another version of the way used to introduce friction has been used: (1-β 2 ) instead Delamination of a coated body A simplified approach to the delamination of a coated body has been conducted to simulate the indentation of a thin Chrome layer adhering to a metal body (Raous et al, 2002).

Cohesive masonry

The RCC model has been used in [START_REF] Jean | Non-smooth contact dynamics approach of cohesive materials[END_REF][START_REF] Acary | Contribution la modélisation mécanique et numérique des édifices maçonnés[END_REF] to simulate the behaviour of a cohesive dome. In Fig. 9.8, the deformation of a dome let on pillars and submitted to gravity is shown.

Wear in biomechanics

In order to simulate the wear occurring in biomaterials and especially the production of wear particles, H. Baudriller [START_REF] Baudriller | Approche micro-mécanique et modélisation de l'usure : application aux biomatériaux des prothéses totales de hanche[END_REF] has used the RCC model in a volumic sense. The adhesive model is used on all the interfaces between the finite elements.

Fig. 9 .

 9 Fig. 9.1 and Fig. 9.2 give the normal and tangential behaviour of the interface during loading and unloading (C N= C T = C, u 0 = ω/C and R 0 = √ ωC) .It should be noted that the Signorini conditions are strictly imposed when compression occurs. References on other models can be found in[START_REF] Raous | Quasistatic Signorini problem with Coulomb friction and coupling to adhesion[END_REF], Raous et al, 1999) and a comparison between some of them is made in[START_REF] Monerie | Comparaison de lois d'interface fibre/matrice sur la base d'un modèle uni-axial d'expérience de micro-indentation[END_REF] (models developed by Tvergaard-Needleman, Girard-Feyel-Chaboche, Michel-Suquet, Allix-Ladevèze). Using penalization and augmented Lagrangian on a similar model to the RCC one, Talon-Curnier have solved the quasistatic problem using generalized Newton method[START_REF] Talon | A model of adhesion coupled to contact and friction[END_REF].

Figure 9 . 2 .

 92 Figure 9.2. Tangential behaviour of the interface

  Figure 9.3. Delamination of two layers submitted to vertical traction

  Figure 9.4. Shear delamination of a block submitted to horizontal loading

Figure 9 . 5 .Figure 9 . 6 .

 9596 Figure 9.5. Interaction of a crack and a fiber/matrix interface

FigureFigure 9 . 8 .

 98 Figure 9.7. Indentation of a coated surface

  Fig 9.9 gives an example of production of fragments.

Figure 9 . 9 .

 99 Figure 9.9. Fragmentation of a block

Acknowledgments

M. Cocou, J.-C. Latil, I. Rosu, in the group MMC at the LMA, L. Cangémi, Y. Monerie, during the preparation of their thesis, and M.A. Karray (ENIT) have contributed to the works presented here. Thanks to V. Acary and H. Baudrier for providing figures 9.8 and 9.9.


P Denis

Anne Ferlay

Pierre Noziere

Christine Gerard

Philippe Schmidely

Determination of linoleic acid requirements for dairy cows using a meta-analysis approach


Nazzareno Diodato

Gianni Bellocchi

Spatial probability modelling of forest productivity indicator in Italy

Keywords:

The prediction of forest productivity is essential for sustainable forest management, particularly in countries, like Italy, where forest is an important part of many protected areas. A spatial predictive probability model for forest productivity rates in Italy was developed over the period , based on 135 annually-resolved records of site productivity and auxiliary variables measured at 219 stations. Our analysis shows that the probability of finding high (>7.3 m 3 ha -1 yr -1 ) and low (<5.8 m 3 ha -1 yr -1 ) productivity rates changes across different regions of Italy. The generated spatial patterns contribute to a better understanding of the factors structuring the distribution of forest productivity in Italy because they reflect the dependence of temperature and water availability conditions on the latitudinal and altitudinal location of the study areas. We observed that the temperature control dominates forest productivity at high elevations and latitudes, whereas low-elevation sites in central and southern Italy are more sensitive to water availability. The proposed spatial probability modelling should be further assessed for its possible incorporation into forest management plans.

Introduction

The spatial extent and pattern of wood distribution have a great influence on the microclimate, pedogenesis and hydrology, including soil erosion of both aquatic and terrestrial ecosystems (Diodato and Bellocchi, 2008b;[START_REF] Diodato | Historical evolution of slope instability in the Calore River Basin, Southern Italy[END_REF], linking up to regional and global climatic aspects [START_REF] Gaillard | Past land-use and land-cover change: the challenge of qualification at the subcontinental to global scales[END_REF][START_REF] Morin | Long-term response of forest productivity to climate change is mostly driven by change in tree species composition[END_REF].

In recent decades, Earth scientists have become increasingly interested in the spatial patterns of ecological processes, including the environmental factors that govern forest productivity [START_REF] Ollinger | Estimating regional forest productivity and water yield using an ecosystem model linked to a GIS[END_REF][START_REF] Purves | Predictive models of forest dynamics[END_REF]Rahman et al., 2015;[START_REF] Forzieri | Satellites reveal contrasting responses of regional climate to the widespread greening of Earth[END_REF]. The analysis of complex interactions between spatial distribution patterns of site factors and woody vegetation growing is crucial for understanding the functioning of terrestrial ecosystems, as well as their management and planning challenges [START_REF] Hörsch | Modelling the spatial distribution of montane and subalpine forests in the central Alps using digital elevation models[END_REF]. Accurate estimation of forest productivity requires understanding of the topo-climatic drivers of forest productivity patterns [START_REF] Stegen | Variation in above-ground forest biomass across broad climatic gradients[END_REF], which in turn depend on ecological conditions and geographic location [START_REF] Babst | Site-and species-specific responses of forest growth to climate across the European continent[END_REF][START_REF] Martinez Del Castillo | Contrasting patterns of tree growth of Mediterranean pine species in the Iberian peninsula[END_REF].

The analysis of forest productivity patterns in complex areas is challenging, mostly because the complex topography and soil diversity of certain regions may cause considerable variations in site-specific responses but also because of the sparsity of available bio-climatic information that is required to estimate productivity (Benavides et al., 2009). Forest biomass estimations may be needed at unsampled locations in a region, for which a procedure is required for mapping forest productivity. In countries with diverse geographical and climatic condi tions like Italy (the focus of this study), detailed mapping of all areas is practically and economically difficult. In Italy, with about 90000 km 2 , forests are estimated to cover about one third of the national territory, and the enhancement of high-nature value forestry environments is a strategic objective of the European Rural Development Policy [START_REF] Pignatti | High Nature Value forest areas: a proposal for Italy based on national forest inventory data[END_REF]. The Italian National Forest Inventory (https:// www.inventarioforestale.org) provides, at country and regional levels, high-resolution estimates of the increase in volume of the trees measured at the time of the survey. At the level of smaller administrative units (including associations of municipalities, which are typical of mountain regions), increment data are generally not available, or they are only fragmentarily available [START_REF] Gasparini | A stand-level model derived from National Forest Inventory data to predict periodic annual volume increment of forests in Italy[END_REF]. Moreover, simply mapping forest productivity does not capture the drivers of forest productivity or give managers the knowledge they need about productivity-environment relationships (e.g. [START_REF] Bontemps | Predictive approaches to forest site productivity: recent trends, challenges and future perspectives[END_REF]. Aerial photography is an advanced method for mapping and monitoring, but it is often hampered by observation noise when conditions are cloudy [START_REF] Coops | Characterizing forest growth and productivity using remotely sensed data[END_REF]. In fact, only a few applications in Italy have been reported, limited to some sites (e.g. [START_REF] Maselli | Use of remotely sensed and ancillary data for estimating forest gross primary productivity in Italy[END_REF][START_REF] Laurin | Above-ground biomass prediction by Sentinel-1 multitemporal data in central Italy with integration of ALOS2 and Sentinel-2 data[END_REF].

Simulation modelling in forest productivity can be a practical and effective methodology for studying spatial patterns of forest dynamics and production in relation with topographic and climatic factors (e.g. [START_REF] Goreaud | Spatial structure analysis applied to modelling of forest dynamics: a few examples Oeiras[END_REF][START_REF] Kim | Assessing the impacts of topographic and climatic factors on radial growth of major forest forming tree species of South Korea[END_REF]. Modelling approaches include the use of geographical information systems (GIS) and spatial analysis techniques for estimating forest productivity and searching for sensitive areas that should be devoted to careful timber management practices (e.g. [START_REF] Naesset | Geographical information systems in long-term forest management and planning with special reference to preservation of biological diversity: a review[END_REF][START_REF] Payn | Potential for the use of GIS and spatial analysis techniques as tools for monitoring changes in forest productivity and nutrition, a New Zealand example[END_REF][START_REF] Busch | GIS-based tool for regional assessments and planning processes regarding potential environmental effects of poplar SRC[END_REF][START_REF] Calvert | Toward renewable energy geo-information infrastructures: applications of GIScience and remote sensing that build institutional capacity[END_REF][START_REF] Thiffault | Developing and validating indicators of site suitability for forest harvesting residue removal[END_REF]. In this study, an Italian-wide geospatial modelling was performed to map forest productivity rates (FPR) over a 30-year climate period , based on an extensive network of sample points with information about vegetation and topo-climatic factors. The FPR index was designed to indicate the growth potential of plants based on climate factors only, through which high values mean the occurrence of favourable growing conditions. A network of 135 stations providing meteorological and productivity data for the calculation of FPR (based on [START_REF] Paterson | The Forest Area of the World and its Potential Productivity[END_REF], with revisions by Serrada, 1976) was combined with auxiliary precipitation and elevation data extending over 219 stations (including the first 135 stations), being used as input data in a geostatistical analysis of FPR. With the aim of assessing important factors limiting the spatial distribution of forests, an integrated geostatistics approach was implemented incorporating multiple input datasets (coverage, shape files, raster, grids), which were processed to estimate forest productivity at unsampled locations. Based on the distance between data and data estimation locations, and the description of their spatial configuration and spatial structure by a theoretical semivariogram under the assumption of isotropy (i.e. all directions are equal), our objective was to model spatial patterns of forest productivity, whilst also assessing the uncertainty associated with predicted patterns.

Materials and methods

Study area and data collection

Italy is a Southern European country characterized by a mixture of topographic systems (peninsular shape, islands, mountain chains, valleys and plains). Centrally located in the Mediterranean Sea basin (Fig. 1a), it is an area through which humid air enters generating precipitation with complex shape, location and evolution of the pressure systems in the atmosphere [START_REF] Lionello | Mediterranean Climate Variability[END_REF]. Influences of the Mediterranean climate on vegetation is evident on the Italian coasts and, becoming more attenuated on the pre-Alpine chains and the Apennines [START_REF] Combourieu-Nebout | Climate changes in the central Mediterranean and Italian vegetation dynamics since the Pliocene[END_REF].

Along the Apennines chain, changes from a low to a high altitude produce a rapid transition towards cooler and less arid summers, while a gradient of increasing continentality (colder winters and hotter summers, with more pronounced temperature variations) can be found when moving inward the Alps. The continental thermal regime is more likely in Northern parts of the peninsula (e.g. in the Po Valley), which are partly isolated by the Apennines chain. Such heterogeneity is reflected in geographically varied vegetation patterns (Fig. 1b).

The selected stations at which FPR values (135) were estimated, and the set of stations (219) providing auxiliary precipitation and elevation data (Fig. 2a, orange and black dots, respectively), cover different biogeographical regions and climatic settings (Fig. 2b). The number of stations is lower in the islands (Sicily and Sardinia) than in the mainland where, as evidenced by [START_REF] Voronoi | Nouvelles applications des paramètres continus à la théorie des formes quadratiques[END_REF] tessellation (notshown), stations are uniformly distributed. Basic climate data for the estimation of forest productivity (Fig. 2a, orange sites) were provided by Food Agricultural Organization (FAO)-database LocClim (FAO, 2014) for the World Meteorological Organization's climate normal period 1961-1990 (for which these data are freely available). The station-specific climatic dataset used in the analysis was built upon monthly climatological data, which were directly extracted from the FAO database. Yearly values were derived from monthly ones.

Data were available at 135 stations for all the variables needed for model application. However, the ability to identify the true spatial variability of a dataset depends, to a great extent, on ancillary knowledge of the underlying measured phenomenon (Diodato and Bellocchi, 2008a). Auxiliary precipitation, temperature and elevation data were derived from the Hydrological Annals of the Italian Hydrographic and Oceanographic Service (SIMN), which covered 219 stations equipped with operational rain gauges (accessible on-line -http://www. isprambiente.gov.it/it/progetti/acque-interne-e-marino-costiere-1/ progetto-annali/gli-annali-idrologici -and described by [START_REF] Acquaotta | Assessment of parallel precipitation measurements networks in Piedmont, Italy[END_REF]. Weather data included monthly rainfall totals and mean temperatures, and the mean maximum and minimum temperatures of the warmest and the coolest months, respectively. In addition, sunshine durations and growth periods were retrieved from this database.

Long-term forest productivity rate

According to [START_REF] Paterson | The Forest Area of the World and its Potential Productivity[END_REF], the potential productivity of a forest site is mostly controlled by climate-related factors. They include the flux of solar radiation reaching the ground, the amount of water available for life processes and the time during which thermal resources are favourable to growth (Rahman and Akter, 2015). At a global scale, temperature conditions and water availability appear to be the main controlling factors of forest productivity beyond solar radiation [START_REF] Churkina | Contrasting climatic controls on the estimated productivity of global terrestrial biomes[END_REF][START_REF] Wang | Modeling the spatial and temporal variability in climate and primary productivity across the Luquillo Mountains, Puerto Rico[END_REF]. A measure of long-term climate-induced site productivity, the FPR index (m 3 ha -1 yr -1 ), was derived from concepts by [START_REF] Paterson | The Forest Area of the World and its Potential Productivity[END_REF] and Serrada (1976), as follows:

(1) where P is the annual rainfall total (mm yr -1 ); T m(h) is the mean monthly temperature of the hottest month (°C); G is the length of the growing season, here calculated using the latitude (Lat, degrees), the longitude (Long, degrees), the temperature (T m(h) as in Eq. ( 1)) and the elevation (Ele, m a.s.l.) of the study site:

(2) Compared to the conventional number of months with values of De Martonne's aridity index >20 [=p/(t + 10), where p is monthly precipitation total in mm and t is monthly average temperature in °C; De Martonne, 1926], Eq. ( 2) takes into account winter seasons (during which low temperatures are the main limiting factor for vegetation growth) and site-specific topographic characteristics.

The reported coefficients of Eq. ( 2) were determined from a multivariate regression (R 2 = 0.89), performed at 14 stations (covering an altitudinal range from 100 to 2000 m a.s.l., and located between 8.6 and 16.6 °E and 39 to 46 °N), for which the actual values of the dependent variable G were known (FAO-LocClim). The actual mean value of G is equal to 4 months (±1 standard deviation), with minimum and maximum values of 2 and 7 months, respectively. The above equation reflects situations frequently occurring in the Mediterranean region, with hot months generally being associated with water limitations, negatively affecting forest growth and productivity when situations of intense and prolonged droughts occur [START_REF] Lindner | Climate change impacts, adaptive capacity, and vulnerability of European forest ecosystems[END_REF]. Forest growing seasons logically tend to become longer with increasing latitudes and longitudes, in agreement with more unlikely drought anom alies in the northern regions and towards eastern-central Italy (Fig. 2b). With respect to site elevation, [START_REF] Qiu | Spatiotemporal variability of vegetation phenology with reference to altitude and climate in the subtropical mountain and hill region, China[END_REF] found a positive relationship with the length of growing season up to 1500 m a.s.l., that is, over an altitudinal gradient in which forest productivity is mostly driven by water availability. For elevations above 1500 m a.s.l., where low temperatures become the most prominent limiting factor, G was set equal to two months, that is, the minimum length of forest growing season (Rahman and Akter, 2015).

The factor f in Eq. ( 1) is the evapotranspiration reducer suggested by [START_REF] Paterson | The Forest Area of the World and its Potential Productivity[END_REF], calculated by relating the solar radiation at the poles to the radiation at the site of interest. While this is meaningful for differentiating potential plant productivity on a global scale (as originally intended by Paterson), the modified factor f introduced by Gandullo and Serrada (1977) better reflects the relatively small spatial scope of the present study. In this formulation, f considers the real insolation at each site by using the annual number of sunshine hours (n sun ), as derived from LocClim FAO (2014) database:

(3)

In this study, the values of f fluctuated between 0.58 and 1.10, representing a number of sunshine hours of 1270 and 3285 h per year, respectively; ΔT is the difference between the mean maximum temperature of the warmest month and the mean minimum temperature of the coolest month (°C). The FPR was calculated for the whole of the Italian territory using data from 135 main stations and auxiliary variables from 219 stations, and interpolation procedures for the creation of continued surfaces.

The geostatistical approach for assessing FPR spatial uncertainty

Geostatistical techniques consider that near spatial data values are more related to each other than distant data values (after [START_REF] Tobler | A computer movie simulating urban growth in the Detroit region[END_REF]. The ordinary kriging, in particular, computes an unsampled value (z), knowing its coordinates (x, y) and neighbour values. It makes optimal, unbiased estimates of regionalized variables at unsampled points from values of the same variables at surrounding stations, using Diodato and G. Bellocchi Ecological Indicators xxx (xxxx) xxx-xxx structural analysis and the initial set of measured data [START_REF] Journel | Mining Geostatistics[END_REF]. In the ecological domain, Diodato and Ceccarelli ( 2004) and Diodato et al. (2010a) used kriging approaches to provide geoindicators to classify soil degradation and assess evapotranspiration rates, respectively. Auxiliary variables can also be used, as in co-kriging (e.g. [START_REF] Knotters | A comparison of kriging, co-kriging and kriging combined with regression for spatial interpolation of horizon depth with censored observations[END_REF]. Kriging/co-kriging approaches are probabilistic in nature. They generally return a single value as an estimate at a particular location, but the uncertainty associated with a local estimate can also be calculated, which is useful in decision management [START_REF] Burrough | Principles of Geographical Information Systems[END_REF]. The present analysis targets the latter option by using a type of linear co-kriging, called detrended ordinary co-kriging in the form of probability maps (DOCKpm). In a first stage, a Detrended Ordinary Co-Kriging with Box-Cox transformation (a power transformation of non-normally distributed data to a set of approximately normal data, after [START_REF] Box | An analysis of transformations[END_REF] was applied to remove the local character of sampled values [START_REF] Kitanidis | Introduction to Geostatistics: Application to Hydrogeology[END_REF]. This was possible because forest growth data were available for more than one attribute per sampled location. In this case, one set of data (Z 1 = forest productivity rate), expensive to measure, was sampled infrequently while another set of data (Z 2 = annual precipitation), cheap to measure, had more frequent observations. Under the assumption that Z 1 and Z 2 are spatially correlated, the information about the spatial variation of Z 2 can be used to map Z1. A co-kriged estimate of the variable of interest at the unsampled location s o is given by: (4)

U N C O R R E C T E D P R O O F N.
where z 1 is a vector of the observed primary data (in this case they correspond to annual forest productivity rates), and z 2 is a vector of the observed auxiliary variable (in this case they correspond to annual precipitation data) selected in the s o neighborhood observation s j ; λ 1 and λ 2 are weight vector associated with the distances h o(i) (between s o and s i ) and h o(j) (between s o and s j ), respectively, and estimated by solving a system of equations [START_REF] Johnston | Using ArcGis geostatistical analyst[END_REF].

To help mitigate the uncertainties associated with estimates of wood productivity rates, a probabilistic interpretation (probability maps) was given (after [START_REF] Diodato | Assessment of the spatial uncertainty of nitrates in the aquifers of the Campania Plain (Italy)[END_REF][START_REF] Diodato | Spatial pattern probabilities exceeding critical threshold of annual mean storm-erosivity in Euro-Mediterranean areas[END_REF]. Thus, in a second stage, estimates and their variances were used to calculate and model semivariogram and cross-covariance functions at a range of threshold values (z k ). The semivariogram captures the spatial dependence between samples by plotting semivariance (half the variance of the differences between all possible points spaced a constant distance apart) against separation distance. The cross-covariance gives the covariance of the primary and secondary variables at pairs of points. Both functions quantify the assumption that things nearby tend to be more similar than things that are farther apart. We set z k equal to 5.8 and 7.3 m 3 ha -1 yr -1 , corresponding to the first and the third quartile of forest productivity values, respectively. In this way, intervals below the lowest threshold and above the highest threshold are captured.

Results

Exploratory data analysis

Fig. 3a shows the geographical distribution of the trend, which was used for detrending the FPR data. Local environmental conditions (mostly in longitudinal sections) can be recognized as deviations from the global North-South trend. The descriptive statistics of the FPR data show a moderately skewed distribution with a small number of samples with large values (Fig. 3b). The Box-Cox transformation (with parameter equal to 1.6) reduced such skewness towards an approximately normal distribution (Fig. 3c), and then these modified data were used for the DOCKpm.

The estimated mean FPR was equal to 6.37 ± 1.24 m 3 ha -1 yr -1 , with values ranging from 2.57 m 3 ha -1 yr -1 for the Paganella mountain (46.17 N, 11.07 E, 2129 m a.s.l.) to 8.73 m 3 ha -1 yr -1 at Udine (46.07 N, 13.23 E, 116 m a.s.l.).

Auxiliary information integration

The observational scale of topo-climatological variables used as auxiliary information is of great importance. In this context, indices are usually designed to represent the main drivers or effects across particular lands where only a minimum set of primary data is available. The direct use of annual precipitation as secondary information produced satisfactory results (correlation coefficient equal to 0.69, Fig. 4a). The cross-covariance modelling with a spherical function (curve in Fig. 4b and Eq. ( 6)), with range of ~80 km (radius of the black circle in Fig. 4b1), confirms an important positive spatial relationship between primary (FPR) and auxiliary (precipitation) variables as the distance is reduced among sampled station-points.

Structural analysis and modelling

We fitted a theoretical semivariogram to the observed data under the assumption of isotropy, and obtained the cross-covariance function between FPR and precipitation. Empirical semivariogram and cross-covariance are indicated with γ 11 (h) C 12 (h), respectively, where subscripts 1 and 2 designate FPR and precipitation data, respectively. Parameters such as range a (~80 km), number of lags (assumed equal 7) and lag size h (assumed equal to 40 km ≈ minimum distance among samples) were calibrated iteratively throughout the run. Once verified the assumption of isotropy, the semivariogram was modelled for both the FPR thresholds as a combination of two distinct spatial structures (nugget variance and spherical structure):

(5) while for the cross-covariance functions:

(6) where is the spherical correlation function equal to (7)

The neighbourhood of points was controlled by a search circle describing the isotropic geometry of the neighbourhood, with semi-axis of about 80 and 100 km limiting the number and the configuration of the points used for prediction. The elements closest to each other tend to be more similar because as the data positions move from a position where the value is unknown, they may not be as useful in predicting the value in an unsampled position. At a certain distance, the points will have no correlation with the position of prediction, and could even be in a very different area than the unknown position.

Spatial pattern of probabilities mapping

Fig. 5 shows the final co-kriged map estimates of FPR for the thresholds of (lower than) 5.7 and (higher than) 7.2 m 3 ha -1 yr -1 . In particular, Fig. 5a shows the spatial patterns of the probability of low-productivity areas. These are mostly located in Trentino-Alto Adige (north-east), Umbria and Latium (centre), Apulia and Basilicata (south), and Sicily and southern Sardinia (islands). Fig. 5b illustrates spatial changes of the probability of high-productivity areas, mostly located across Piedmont, Lombardy, Liguria (north-west), Friuli-Venezia Giulia (north-east), Marche and (to a lesser extent), Tuscany (centre), Central and southern Campania, western Basilicata and Calabria (south).

The conclusion is that more suitable areas (that is, with high probability of forest productivity) are located in northern Italy. Instead, Latium, Sardinia, Sicily and Apulia are the regions with high probability to have a poor suitability to forest productivity. In these areas, the scarce precipitation regimes during the growing season, combined with high evapotranspiration rates, cause limited water availability for plants growth.

Cross-validation assessment

The error involved on the expansion of the information from point to landscapes through probability co-kriging estimation at fine grid can be assessed through an estimated indicator's standard error and cross-validation [START_REF] Isaaks | An Introduction to Applied Geostatistics[END_REF]). In the cross-validation diagrams (Fig. 6), the actual values of FPR compare well with the predicted probability indicators for both the highest (Fig. 6a) and lowest (Fig. 6b) threshold values. Mean values are -0.002 (for FPR < 5.7 m 3 ha -1 yr -1 ) and 0.0249 (for FPR > 7.2 m 3 ha -1 yr -1 ), showing lack of systematic errors. Root mean squared error values are 0.420 and 0.400, respectively.

Discussion

Detailed studies are necessary to establish the mechanisms underlying unique functional relationships between climate and forest growth in the complex set of physiological processes governing the response of trees to water and temperature stress. This paper presents a geostatistical modelling framework for proper interpretation of spatially-explicit forest productivity rates (FPR) in Italy, which assumes a set of quantitative data for the location of interest (temperature, rainfall, duration of the growth season) and probabilities associated with ranges of productivity above and below critical values. They include the best available data from a set of stations (135) and secondary data from a wider set (219), as provided by the existing FAOclim database. In this way, the use of topo-climate variables as a source of auxiliary information for estimating the FPR for site-specific and other spatial information is an important issue of this paper. Trend analysis has shown the existence of a non-random (deterministic) component in the spatial distribution of data. In particular, the largest gradient of the FPR data occurs along the northwest to southeast direction (Fig. 3a), which can be attributed to Atlantic influences determining more persistent rainfall in the northwestern part and milder temperatures at lower latitudes [START_REF] Pfister | Early modern Europe[END_REF]. In the application of kriging-type analyses, a gradient could mask the stochastic variability in the range of kriging estimates. Kriging with a trend model (detrended kriging) considers that the local mean varies within the local neighbourhood. The local trend can be, for example, a linear relationship between the investigated variable and auxiliary variables. Detrending the data locally (by subtracting a constant value or a best-fit line from data) was thus a part of the prediction process to satisfy stationary assumptions (e.g. [START_REF] Negreiros | Geographical information systems principles of ordinary kriging interpolator[END_REF]. Here, a constant detrending with a local neighbourhood was applied (Fig. 3a) as a practical approximation, while straight line or polynomial detrending are applicable to more complex situations requiring strong assumptions about parameter estimation (e.g. [START_REF] Myers | Statistical methods for interpolation of spatial data[END_REF].

The data resource was decomposed into a deterministic trend component (which usually changes gradually through space) and an autocorrelated random component (which tends to change more quickly). We assumed that observations on any given year show local variations owing to topography, ground cover and other local factors that are not easily predictable. According to the procedure shown in [START_REF] Goovaerts | Geostatistics for Natural Resources Evaluation[END_REF], once the trend was removed, kriging was carried out on the residuals and the trend was added back to the output surface. The results so obtained indicated that probabilities of high and low productiv ity are distributed without any apparent latitudinal/longitudinal or altitudinal gradient. While suggesting the possibility of using geostatistical spatial modelling to determine FPR values, the present study also showed limitations. These findings are purely geostatistical, with no explicit mechanistic basis related to forest productivity. However, mechanistic modelling of forest growth is data demanding and state-of-the-art model estimates can be accompanied by a large amount of uncertainties (e.g. [START_REF] Charney | Observed forest sensitivity to climate implies large changes in 21 st century North American forest growth[END_REF]. There are also two major criticisms against the Paterson index, that is: (i) it estimates growth potential in terms of timber volume, for which it was originally designed [START_REF] Johnston | Using ArcGis geostatistical analyst[END_REF]Hagglund, 1981), without providing dry matter biomass, and (ii) it does not consider any soil factor. Some other limitations are in the use of pre-established topo-climatic inputs, which do not include surface insolation, second-order constraints provided by slope and orientation of the local terrain, and soil water regime (e.g. Rahman and Akter, 2015). This approach also overlooks the influence of short-term weather events (and year-to-year variability) on estimated productivity.

Relying on credible information regarding forest productivity and its environmental drivers, this study shows how geostatistics methods can be practically implemented to create countrywide spatially explicit probabilistic maps, which could be of direct use for forest management.

Fig. 1 .

 1 Fig. 1. Environmental settings of the study area (a), and map of potential vegetation of Italy (b), as derived from the European vegetation chart by Bohn et al. (2005).

Fig. 2 .

 2 Fig. 2. Geographic distribution of 135 Italian sites for which the FPR (forest productivity rate) was estimated (orange dots), and 219 stations providing auxiliary precipitation co-variables (black points) (a), and spatial patterns of precipitation during the forest growing season (April-September) (b). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 3 .Fig. 4 .

 34 Fig. 3. Geographical trends of forest productivity rates (FPR) (a), statistical distribution of baseline FPR data for the original series over the period 1961-1990 (b), and statistical distribution of baseline FPR data after Box-Cox transformation (c). FPR data are expressed as m 3 ha -1 yr -1 .

  Fig. 5. Spatial patterns of probability with low productivity prone areas (probability of FPR < 5.7 m 3 ha -1 yr -1 ) (a), and high productivity prone areas (probability of FPR > 7.2 m 3 ha -1 yr -1 ) (b).

Fig. 6 .

 6 Fig. 6. Scatterplot of cross-validation between the actual forest productivity rate (FPR) and kriged indicator for two thresholds (blue vertical lines): FPR < 5.7 m 3 ha -1 yr -1 (a) and FPR > 7.2 m 3 ha -1 yr -1 (b). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

  and G. Bellocchi Ecological Indicators xxx (xxxx) xxx-xxx

Acknowledgement

This study is an investigator-driven study, which means generation of the study hypothesis, study design, data collection, analysis and interpretation, decision to publish and the preparation of the manuscript were performed without any financial support.