Julien Claude

Emmanuel Paradis
email: paradis@isem.univ-montp2.fr

Analysis of comparative data using generalized estimating equations

come

Introduction

Comparing species is central in many biological issues. For instance, one of [START_REF] Darwin | On the origin of species by means of natural selection[END_REF] main point in his theory of evolution by natural selection was to show that "wide-ranging, much diffused, and common species vary most." More recently, scaling and allometric relationships have been widely used in biology and physiology [START_REF] Calder | Ecological scaling: mammals and birds[END_REF][START_REF] Peters | The ecological implications of body size[END_REF][START_REF] Schmidt-Nielsen | Scaling: why is animal size important?[END_REF].

Traditionally, these relationships were studied using statistical methods (mainly regressions) which assume that species are independent observations. However, species are not independent observations because they are linked by their phylogenetic relationships. Several methods have been proposed to analyse interspecific comparative data taking this difficulty into account [START_REF] Cheverud | The quantitative assessment of phylogenetic constraints in comparative analyses: sexual dimorphism in body weight among primates[END_REF][START_REF] Felsenstein | Phylogenies and the comparative method[END_REF][START_REF] Grafen | The phylogenetic regression[END_REF][START_REF] Gittleman | Adaptation: statistics and a null model for estimating phylogenetic effects[END_REF][START_REF] Lynch | Methods for the analysis of comparative data in evolutionary biology[END_REF] among others).

The method of phylogenetically independent contrasts (referred to as the contrasts method in this paper) substitutes the original data collected on species by a set of contrasts computed between pairs of sister-species, and between pairs of nodes in the reconstructed phylogeny [START_REF] Felsenstein | Phylogenies and the comparative method[END_REF]. The character states at the node are reconstructed under a model of random evolution [START_REF] Felsenstein | Phylogenies and the comparative method[END_REF]. This method is certainly the most widely used to deal with the problem of non-independence among species. Recently, [START_REF] Martins | Phylogenies and the comparative method: a general approach to incorporating phylogenetic information into the analysis of interspecific data[END_REF] presented a method which can be viewed as an extension of [START_REF] Grafen | The phylogenetic regression[END_REF] method, based on the analysis of a linear model with generalized least squares (GLS) where the variances and covariances of the data are defined with respect to the phylogenetic relationships among species.

In this paper, we present a new approach for the analysis of comparative data taking into account the phylogenetic relationships among species. The traits analysed can be either continuous or discrete, and more than two traits can be analysed simultaneously in a generalized linear modelling framework. The model is fitted using generalized estimating equations where the dependence among species is taken into account with a correlation matrix. We run simulations to assess some statistical properties of our approach. We illustrate its use with an analysis of some data on birds [START_REF] Paradis | Dispersal and spatial scale affect synchrony in spatial population dynamics[END_REF]. [START_REF] Felsenstein | Phylogenies and the comparative method[END_REF] presented the statistical problem with comparative data when closely related species are included in a sample. These species should, in fact, be considered as pseudo-replications of the same observations, and considering them as independent inflates the number of degrees of freedom in the analysis, resulting in an increased type I error rate. For instance, if in a sample of eight species there are four pairs of sibling-species, then the number of degrees of freedom in a linear regression model should be two, and not six (the critical values for a t-test with a significance level of 5% are 2.45 for df = 6, and 4.30 for df = 2).

A generalized estimating equations approach

The strength of the dependence among the observations on two species depends on the phylogenetic distance between them, and how the observed characters evolve through time. Thus, a measure of this dependence can be derived from the distances measured on a phylogenetic tree. An appropriate transformation of these distances, depending on the mode of character evolution, gives a matrix of correlations among species: this matrix is symmetric and its diagonal elements equal to unity (Fig. 1). Generalized estimating equations (GEE) are a procedure to fit regression models taking the correlations among the observations into account [START_REF] Liang | Longitudinal data analysis using generalized linear models[END_REF]. The correlation matrix can take different structures, and its parameter(s) can either be estimated from the data, or fully specified; in the latter case, the structure of the correlation matrix is said to be fixed. In all cases, the parameters from the regression model are estimated from the data. This regression model is a generalized linear model (GLM, McCullagh & Nelder, 1989), meaning that the response variable may be non-normal, for instance, binomial or Poisson.

Consider that we have n species. The phylogenetic relationships among the species are known (topology and branch lengths), and let D be the n × n matrix of pairwise phylogenetic distances. Let y y y be a n × 1 vector of a variable (called the response) whose elements are denoted y i (where the subscript i denotes the ith species, i = 1, . . ., n), and X be the n × p matrix of covariates (or predictors).

Let us assume that the y i follow a marginal distribution belonging to the exponential distribution family, making possible a generalized linear regression of y y y with respect to X :

g(E[y i ]) = x x x T i β β β, (1)
where β β β is a p × 1 vector of regression parameters, x x x T i is the transposed vector of the covariate values for the ith species, and g is a link function. Equation ( 1) is

identical to g(E[y i ]) = β 0 + β 1 x i1 + . . . + β p x ip
, where β 0 is an intercept. This does not differ from a standard GLM. Recall that a GLM is a linear model of the expected mean of a variable belonging to the exponential distribution family (e.g., Gaussian, gamma, Poisson, binomial, for the most commonly used distributions), the expected mean being transformed using a link function [START_REF] Mccullagh | Generalized linear models[END_REF]. In a GLM, the variance of the responses is given by:

Var(y i ) = φV (E[y i ]),
where φ is the dispersion parameter, and V (E[y i ]) is the variance function, both of them being defined with respect to the distribution assumed for y. Note that this is the variance expected under the assumption that all species are independent observations, which we do not want to assume here. If the observations are not independent, we can define the variance-covariance matrix among observations with [START_REF] Liang | Longitudinal data analysis using generalized linear models[END_REF]):

V = φA 1/2 RA 1/2 , ( 2
)
where A is a n × n diagonal matrix defined by diag{φV (E[y i ])}, that is a matrix with all its elements null except the diagonal which contains the variances of the n observations expected under the marginal model, and R is the correlation matrix of the elements of y. If the observations are in fact independent, then R is an n × n identity matrix.

Estimating parameters and testing hypotheses

The main objective of the present approach is to look for relationships between y y y and X . An advantage of the GLM setting is that the y i can follow several types of distributions including Gaussian, gamma, Poisson, or binomial. The two latter distributions are appropriate to model, for instance, numbers and frequencies, respectively. The variables included in X could be continuous or categorical, and the model can include additive, interactive, and nested effects among these predictors.

From the assumptions in the previous section, it is possible to define estimating equations which are consistent estimators of the regression parameters β β β. These generalized estimating equations are:

∂µ µ µ ∂β β β T V -1 (y y y -µ µ µ) = 0, (3)
where µ µ µ is the n × 1 vector of the mean expected responses whose elements µ i (= E[y i ]) are given by g -1 (x x x T i β β β), with g -1 as the reciprocal of the link function in equation ( 1). The Appendix describes a procedure to solve equation (3), and several estimators of the standard-errors of β β β.

The estimates of β β β (denoted β β β) can be found with GEE using information from D and the phylogenetic tree as the correlation matrix R (Fig. 1). We will assume that no element from the matrix R needs to be estimated, so the correlation structure is fixed.

GEE can be fitted easily with existing software [START_REF] Horton | Review of software to fit generalized estimating equation regression models[END_REF].

The estimates of the standard-errors (SE) of β β β can be used to test for the significance of each parameter in the model: under the null hypothesis that β = 0, the ratio β/SE( β) follows a t distribution. In the present implementation of GEE, the number of degrees of freedom of this t-test is not the usual residual number of degrees of freedom of the fitted model (that is the sample size minus the number of estimated parameters): since there is a single cluster of observations, this residual number of degrees of freedom actually overestimates the true number of degrees of freedom. In order to correct for this bias, we used a procedure where the number of degrees of freedom is calculated from the phylogenetic tree. The principle of this procedure is to consider such a tree as a representation of the inter-dependence among observations, and thus as the number of degrees of freedom in the data. First consider, for simplicity, that the tree is ultrametric, each tip being equally distant from the root. The number of degrees of freedom is counted by adding, at each dichotomous node in the tree, one unity multiplied by the ratio of the distance from the node to the tips on the distance from the root to the tips. If the node is trichotomous, then two units (multiplied by the same ratio) are added, if the node is tetrachotomous, three units are added, and so on. This is done starting from the root where two degrees of freedom are counted. We denote the final number df P . This procedure can be generalized to a non-ultrametric tree with:

df P = ∑ tree branch length ∑ n i=1 distance from root to tip i × n.
We will further discuss the need for this correction in the discussion.

To test for the significance of the effect(s) of one or several predictors, we can use an ANOVA-like analysis. Let W denote the variance-covariance matrix of β β β.

Under the null hypothesis of no effect, the quantity β

β β T Ŵ -1 β β β follows an F
distribution with df and df Pdf -1 as numbers of degrees of freedom, where df is the number of degrees of freedom associated with β β β (i.e. the number of estimated parameters for the tested effects).

Simulation study

In order to assess some properties of the approach we present here, we conducted simulations of character evolution under different scenarios. We considered three phylogenetic trees (Fig. 2). All three trees have 32 tips and the same time of evolution from the root to the tips (155 time-steps). Tree A is a balanced phylogeny with equal branch lengths, tree B has strongly clustered tips, and tree C has some imbalance in its terminal branches. We simulated the evolution of four continuous characters (denote them z 1 , z 2 , z 3 , and z 4 ) with a Brownian motion model, and eight discrete binary characters taking the value 0 or 1 (y 1 to y 8 ) with a Markovian model. All characters were set equal to zero at the root of the tree. The continuous characters evolved along the branches by adding, at each time-step, a random normal variate. After reaching a node, the characters evolved independently on each daughter-branch. The characters z 1 and z 2 evolved in a similar way and independently according to: z 1t+1 = z 1t + ε, where t is the time-step along the branches, and ε ∼ N(0, 1). On the other hand, z 3 and z 4 evolved in a correlated way according to: z 3t+1 = z 3t + ε, and

z 4t+1 = z 4t + ζ,
where ζ ∼ N(γz 3t , 1). The parameter γ specifies the strength of the "co-evolution" between z 3 and z 4 . We repeated the simulations with four different values of γ: 0.001, 0.005, 0.01, 0.015. The discrete characters had fixed transition probabilities given by a symmetric matrix: the probability of remaining in the same state (either 0 or 1) was p, whereas the probability of changing was 1p. The characters y 1 and y 2 evolved independently and with the same transition probabilities (p = 0.6), and y 3 and y 4 evolved independently too but with p = 0.8. This higher probability is more likely to result in spurious patterns of association between the states of y 3 and y 4 than for y 1 and y 2 [START_REF] Grafen | A new model for discrete character evolution[END_REF]. The four remaining characters evolved in pair: y 5 evolved in the same way than y 3 and y 4 , but y 6 had transition probabilities which depended on the state of y 5 (they were also equal to 0.6 and 0.4). The same model of coevolution was used for y 7 and y 8 , but the transition probabilities were 0.8 and 0.2.

The simulations were replicated 1000 times for each parameter value.

At the end of the evolution process, the values of the characters at the tips of the tree were analysed with six regressions using GEE: z 1 on z 2 , z 4 on z 3 , y 1 on y 2 , y 3 on y 4 , y 6 on y 5 , and y 8 on y 7 . The regressions with the continuous characters were done assuming a Gaussian distribution for the response and with an identity link; the regressions for the discrete characters assumed a binomial distribution for the response, and used a logit link. We expected the test on the significance of the slope to reject the null hypothesis in 5% of the replications (the expected type I error rate since the null hypothesis was true) with the first regression, whereas the rejection rate of the null hypothesis with the second regression will depend on the power of the test (since the null hypothesis was false). The correlation matrix R was derived from the expected variances and covariances of the characters z 1 , z 2 , z 3 , and z 4 at the tips of the trees under a

Brownian motion model of evolution [START_REF] Butler | The relationships between sexual size dimorphism and habitat use in Greater Antillean Anolis lizards[END_REF][START_REF] Garland | Using the past to predict the present: confidence intervals for regression equations in phylogenetic comparative methods[END_REF].

To compare the GEE approach with a standard contrasts method, we also analysed the simulated data with contrasts. The contrasts were computed following [START_REF] Felsenstein | Phylogenies and the comparative method[END_REF]. For the discrete characters, the values 0 and 1 were used untransformed (see [START_REF] Grafen | Statistical tests for discrete cross-species data[END_REF]. A linear correlation coefficient was then computed between the 31 pairs of contrasts, and its statistical significance was assessed using a t-test with 29 degrees of freedom.

Finally, we analysed some of the simulated data with a GLM (i.e., assuming that the observations are independent) in order to assess its type I error rate. We computed for each analysis and each set of parameters, the rate of rejection of the hypothesis of a null slope at a nominal level of 5%. This rate is (i) the type I error rate if the characters were uncorrelated, or (ii) the power of the test if they were correlated. Among a sample of size 1000, a rate of success is significantly greater than 0.05 if 62, or more, successes are observed (this is a one-sided test since we are not interested in the cases where the rejection rate is significantly smaller than 0.05). Thus, the rejection rate was significantly greater than 5% if we observed 62 or more rejections of the null hypothesis. All these simulations were programmed in R [START_REF] Ihaka | R: a language for data analysis and graphics[END_REF]. The phylogenetically independent contrasts were computed using a C code from the package Phylip [START_REF] Felsenstein | Phylip (Phylogeny Inference Package) version 3.5c[END_REF]) called from R.

The results are summarized in Tables 1 and2. For the continuous characters, the performance of GEE and contrasts were very close; however, the type I error rate of GEE was significantly greater than 5% (around 7%), whereas it was not significantly different from 5% for the contrasts (Table 1). With respect to the power of the tests, the performance of GEE was slightly better than contrasts for the smallest values of γ, but this was the opposite for the greatest values of this parameter (Table 1).

For the discrete characters, GEE performed better: the type I error rates were not significantly different from 5%, except for phylogeny C (Table 2). The type I error rates of contrasts were always significantly greater than 5% for these characters. When the discrete characters evolved in association, the performance were not very good for the lowest value of the parameter (p = 0.6). For the greatest value of this parameter (p = 0.8), the power of the tests were very similar for the different phylogenies, around 37% for GEE, and around 53% for the contrasts (Table 2).

We evaluated the type I error rate of GLM in the case of continuous characters it was 25.4%, 39.5%, and 25.3% for phylogeny A, B, and C, respectively.

Application: relationship between dispersal and population synchrony

There has been a large number of recent publications on the relationships between dispersal and population synchrony (e.g., [START_REF] Lande | Spatial scale of population synchrony: environmental correlation versus dispersal and density regulation[END_REF][START_REF] Kendall | Dispersal, environmental correlation, and spatial synchrony in population dynamics[END_REF]. The strength of synchrony in fluctuations of population density is hypothesized to be driven either by environmental factors (such as climate), or by dispersal connecting local populations [START_REF] Haydon | The effects of large-and small-scale random events on the synchrony of metapopulation dynamics: a theoretical analysis[END_REF]. Interspecific

where the term dispersal*habitat means that the effect of dispersal on synchrony depended on the type of habitat (it was stronger for the species nesting in wet habitats), and the long-term national trend is the slope of the linear change of the global population against time between 1962 and 1995 (the most declining species having more synchronized populations). This left open the issue of the eventual influence of phylogenetic relationships on this relationship since such complex models cannot be fitted directly with contrasts.

We re-analysed these results with GEE separately for breeding and natal dispersal (as done in [START_REF] Paradis | Dispersal and spatial scale affect synchrony in spatial population dynamics[END_REF]. We considered the model described by equation ( 4). The phylogenetic relationships among the 53 species were taken from [START_REF] Sibley | Phylogeny and classification of birds: a study in molecular evolution[END_REF]; we computed the correlation matrix R using the ∆T 50 H distances given by [START_REF] Sibley | Phylogeny and classification of birds: a study in molecular evolution[END_REF]. The number of degrees of freedom of this tree was df P = 16. The results of the F-tests are summarized in Table 3. For both natal and breeding dispersal, the model was strongly significant indicating that the above relationship was not influenced by the phylogenetic relationships among species. The parameter estimates with their standard-errors are given in Table 4. Interestingly, in this analysis the effect of long-term trend appears as just close to significance, whereas it was strongly significant in the GLM without correcting for phylogeny [START_REF] Paradis | Dispersal and spatial scale affect synchrony in spatial population dynamics[END_REF]. When looking at the data on population trend, it appears that some of the most closely related species share a common trend: for instance, doves (three species) or tits (five species) are increasing, whereas thrushes (three species) and finches (five species) are declining.

Discussion

The analysis of comparative data taking phylogeny into account is now widespread among evolutionists and comparative biologists. Our purpose was to present a new approach that tries to be as general as possible since it has some features of the methods currently available. The GEE approach can be viewed as an extension of the GLS one [START_REF] Martins | Phylogenies and the comparative method: a general approach to incorporating phylogenetic information into the analysis of interspecific data[END_REF] with two differences: the dependence among observations is specified with a correlation matrix in GEE (a variance-covariance matrix in GLS), and non-normal distribution can easily be modelled with GEE due to a generalized linear modelling framework. GEE permit the use of categorical predictors as well as continuous ones, similarly to the phylogenetic regression [START_REF] Grafen | The phylogenetic regression[END_REF], or the GLS method [START_REF] Martins | Phylogenies and the comparative method: a general approach to incorporating phylogenetic information into the analysis of interspecific data[END_REF]. Another feature that is common to the GEE approach and some extensions of the contrasts method [START_REF] Grafen | The phylogenetic regression[END_REF][START_REF] Martins | Phylogenies and the comparative method: a general approach to incorporating phylogenetic information into the analysis of interspecific data[END_REF] is the possibility to fit complex models that include continuous, categorical predictors, and possible interactions between them.

Our simulation study showed that the approach with GEE has good statistical properties. Unsurprisingly, it had much lower type I error rates than GLM.

However, the type I error rates of GEE were slightly greater then 5%, and this difference was statistically different. A possible explanation is that this comes from the fact that the number of degrees of freedom we used in these tests was actually an approximation of the true number. Overall, the contrasts performed very well for continuous characters. This confirms several simulation studies [START_REF] Martins | Phylogenetic analyses of the correlated evolution of continuous characters: a simulation study[END_REF][START_REF] Díaz-Uriarte | Testing hypotheses of correlated evolution using phylogenetically independent contrasts: sensitivity to deviations from Brownian-motion[END_REF][START_REF] Harvey | Phylogenetic extinction rates and comparative methodology[END_REF]. For discrete characters, the GEE approach performed globally better than the contrasts one. The type I error rates of GEE were kept smaller than 5%, except for phylogeny C. In all situations the contrasts had inflated type I error rates which certainly comes from the invalid assumption of the distributions of the characters, thus the contrasts method is not robust to violation of this assumption. In terms of power of the tests, they performed very poorly when the transition probabilities were {0.6, 0.4}, but this may come simply from the fact these probabilities being too close to 0.5, and thus the co-evolution between the two discrete characters was too weak to be detected. When the transition probabilities were {0.8, 0.2}, the tests performed better, with greater power for the contrasts than for GEE, but this was obviously at the expense of the increased type I error rate when the null hypothesis was true.

Several recent studies have attempted to compare different approaches to the comparative method [START_REF] Garland | Using the past to predict the present: confidence intervals for regression equations in phylogenetic comparative methods[END_REF][START_REF] Rohlf | Comparative methods for the analysis of continuous variables: geometric interpretations[END_REF]. [START_REF] Garland | Using the past to predict the present: confidence intervals for regression equations in phylogenetic comparative methods[END_REF] showed that the contrasts and GLS methods are identical and expected to give the same results. In theory, the GEE assuming normal errors and the GLS are identical. It may thus be surprising that the contrasts and the GEE yielded different results in our simulations of continuous characters, though the differences were small (Table 1). It should be noted that GEE and GLS are solved with quite different algorithms which may partially explain this discrepancy. Furthermore, the difficulties in estimating the coefficient standard-errors (see below) may be important here as well.

Additionally, it has been demonstrated that different computer programs may give remarkably different results for much simpler problems than the one of comparative analysis such as summary statistics, analysis of variance, or linear regression [START_REF] Mccullough | Assessing the reliability of statistical software: Part I[END_REF][START_REF] Mccullough | Assessing the reliability of statistical software: Part II[END_REF]. Whether we should expect similar discrepancies for the different comparative methods and their computer implementations is not clear, and obviously needs further study.

We did not simulate character evolution with phylogenies having multichotomies because the program code to compute contrasts available for our simulations can only deal with dichotomies (programs that can compute contrasts in presence of multichotomies cannot be called from R). Some contrasts methods are known to be valid in the presence of multichotomies [START_REF] Grafen | The phylogenetic regression[END_REF][START_REF] Purvis | Polytomies in comparative analyses of continuous characters[END_REF]. The procedure to compute the and the associated number of degrees of freedom differ, however, depending on whether the multichotomies should be treated as true instantaneous radiations or as a series of unresolved dichotomies [START_REF] Purvis | Polytomies in comparative analyses of continuous characters[END_REF]. Since, in the GEE approach, we associate some degrees of freedom to each node, the same correction can be applied in the method proposed here, that is only one (fraction of) degree of freedom should be counted for each unresolved multichotomy, rather than the number of daughter-species minus one.

The need to correct the number of degrees of freedom with GEE is not obvious since the dependence among observations is already taken into account with the correlation matrix. However, it is known that GEE estimators of the standard-errors of the regression coefficients (SE(β)) perform poorly when the number of clusters is low [START_REF] Horton | Review of software to fit generalized estimating equation regression models[END_REF][START_REF] Mancl | A covariance estimator for GEE with improved small-sample properties[END_REF].

Our approach was thus to correct the number of degrees of freedom in order to keep the type I error rate at an acceptable level. Furthermore, our simulations showed that the power to detect correlated evolution with GEE is close to that of phylogenetically independent contrasts. However, the estimation of coefficient standard-errors with GEE is still under progress [START_REF] Mancl | A covariance estimator for GEE with improved small-sample properties[END_REF], and some attention will be needed in the context of the application of this method to the analysis of comparative data.

According to our simulations, the GEE approach seems to perform well for the analysis of discrete characters. Other methods exist for similar analyses [START_REF] Sillén-Tullberg | The effect of biased inclusion of taxa on the correlation between discrete characters in phylogenetic trees[END_REF][START_REF] Pagel | Detecting correlated evolution on phylogenies: a general method for the comparative analysis of discrete characters[END_REF][START_REF] Read | Inference from binary comparative data[END_REF]among others). [START_REF] Grafen | Statistical tests for discrete cross-species data[END_REF] assessed the type I error rates of four such methods. To simulate some null data, they used three phylogenies each with 256 tips [START_REF] Grafen | Statistical tests for discrete cross-species data[END_REF]; thus their results cannot be really compared with ours. Clearly, an extensive comparison of the error rates of these methods in different situations (sample sizes, tree balance, pattern of multichotomy) is needed.

The contrasts method can deal only with continuous variables, whereas we have shown that GEE can deal with a variety of distributions because of the GLM setting. This points to the difference in philosophy betwen both approaches. The contrasts method uses a null model of evolution through the phylogenetic tree which is suitable for continuous variables (the Brownian motion model), and then focuses on the contrasts between pairs of taxa and/or nodes which are expected to be normally distributed under the assumed model of evolution. Our GEE approach, focuses on the relationships between the variables observed on recent species, and takes into account their phylogenetic relationships in parameter estimation and statistical testing. We think that both approaches have their own merits. Nevertheless, we believe that the GEE approach is particularly well-suited to characterise relationships between ecological, biological, and physiological variables when the non-independence among taxa (or, more generally, observations) is a potential problem for statistical inference. In this respect, the regression coefficients estimated by GEE can be interpreted functionally in the same way than the coefficients estimated by a standard linear regression. This is due to the fact that the focus of GEE is on the mean structure in the data [START_REF] Palmgren | Exponential family models and statistical genetics[END_REF]. For instance, it is possible to compare a theoretically expected value of the coefficient with a 95% confidence interval β ± 1.96 × SE( β).

Like the contrasts method, the GEE approach makes some assumptions on the model of evolution of the characters through the transformation of phylogenetic branch lengths to obtain the correlation matrix. Though we did not consider such a possibility in this paper, we mentioned above that some parameters of the correlation matrix can be estimated from the data. This is possible by assuming a particular structure of this matrix (exchangeable, unstructured, auto-regressive, among others, see [START_REF] Horton | Review of software to fit generalized estimating equation regression models[END_REF] for more details on correlation matrix structures). An interesting potentiality here would be to include parameters of the model of character evolution in the correlation matrix, and estimate them from the data [START_REF] Liang | Multivariate regression analyses for categorical data (with discussion)[END_REF].

An interesting possibility is to apply different transformations to the elements of the matrix D, or to use different measures for D. For instance, instead of using times of divergence, D could be a matrix of genetic or molecular divergence which is justified if it is assumed that the rate of phenotypic evolution is proportional to the rate of molecular evolution. Thus, it is possible to introduce heterogeneous rates of evolution in the analysis.

GEE can be used to analyse comparative data at the intraspecific level as well.

The concern about phylogenetic dependence has traditionally concentrated on the interspecific level, but the problem exists within a species too, and is probably even more important in this case since the shared history among observations (individuals or local populations) is more recent than in the interspecific case.

To conclude, we present in this paper an approach for the analysis of comparative data taking into account phylogeny by the use of generalized estimating equations. The use of a generalized linear modeling framework allows the analysis of a great variety of models. The focus on the relationships among variables observed on recent species makes possible the interpretation of the relations in terms of constraints.

We are grateful to Rob Freckleton, Christophe Thébaud, and two anonymous referees for helpful comments on a previous version of our paper. This is publication 02-031 of the Institut des Sciences de l' Évolution (Unité Mixte de Recherche 5554 du Centre National de la Recherche Scientifique).

(iv) alternate between steps (ii) and (iii) until convergence.

Liang & Zeger (1986) introduced a robust (or empirical) estimator of W the variance-covariance matrix of β β β. It is robust in the sense that it is consistent (i.e.

it convergences to the true value of W when sample size increases) even if the correlation matrix R is misspecified. However, this robust estimator has very poor properties when the number of independent clusters is small (< 20), so it is not appropriate to the analysis of comparative data where there is a single cluster since all species are linked by their phylogenetic relationships. We prefer the use of a naive (or model-based) estimator of W given by:

W naive = ∂µ µ µ ∂β β β T V -1 ∂µ µ µ ∂β β β -1
.

This estimator is consistent if R is correctly specified; however, [START_REF] Horton | Review of software to fit generalized estimating equation regression models[END_REF] report that it has better statistical properties than the robust estimator when the number of independent clusters is small even if R is wrong.

Furthermore, we were able to verify that the naive estimator has good properties for continuous characters assuming a Gaussian response (see our present simulation study). However, it performs poorly for discrete characters assuming a binomial response since the number of independent clusters is too small [START_REF] Mancl | A covariance estimator for GEE with improved small-sample properties[END_REF]. We chose for binomial responses a quasi-likelihood estimator of W given by:

W quasi = - ∂ 2 ln Q ∂β β β 2 -1 ,

  comparisons are obviously valuable to test such hypotheses since population parameters (like dispersal distance) vary much more at the interspecific level than at the intraspecific one.[START_REF] Paradis | Dispersal and spatial scale affect synchrony in spatial population dynamics[END_REF] showed a positive correlation between dispersal distance and population synchrony for 53 species of birds in Britain and Ireland. An examination of the data revealed that other variables obviously influenced the relationships, and GLMs were used to find which among some candidate variables significantly affect population synchrony. The model finally selected has the following form: population synchrony = dispersal*habitat + long-term national trend,

Fig. 2 .

 2 Fig. 2. The three phylogenies used in the simulation study.

Table 1 .

 1 Results of analyses of simulated continuous characters along the phylogenies on Fig.2. The parameter γ specifies the strength of the association between the two correlated characters. The figures in the table are the rates of rejection of the hypothesis of a null slope between the characters. In the case of uncorrelated evolution, the rejection rate is the type I error rate of the test. † In the case of correlated evolution, the rejection rate is the power of the test (= 1-type II error rate).

	phylogeny method	uncorrelated* correlated characters (with γ =) †
			characters	0.001 0.005 0.01	0.015
	A	GEE	0.066	0.083 0.153 0.372	0.568
		contrasts	0.059	0.067 0.152 0.406	0.653
	B	GEE	0.067	0.072 0.167 0.478	0.706
		contrasts	0.047	0.051 0.153 0.493	0.752
	C	GEE	0.071	0.067 0.154 0.362	0.570
		contrasts	0.042	0.061 0.121 0.396	0.659

*

Table 2 .

 2 Results of analyses of simulated discrete characters along the phylogenies on Fig.2. Schematic representation of the comparative analysis method using generalized estimating equations as proposed in the present paper.

	phylogeny method	p = 0.6		p = 0.8
			uncorrelated correlated	uncorrelated correlated
	A	GEE	0.040	0.057	0.049	0.370
		contrasts	0.087	0.095	0.092	0.513
	B	GEE	0.040	0.028	0.041	0.371
		contrasts	0.078	0.068	0.073	0.536
	C	GEE	0.073	0.084	0.088	0.369
		contrasts	0.128	0.133	0.129	0.538

APPENDIX

The generalized estimating equations in equation (3) can be solved through an iterative process which can be summarized as follows:

(i) compute an initial estimate of β β β, for example, with a GLM;

(ii) compute an estimate of the variance-covariance matrix using equation ( 2);

(iii) update β β β with: [START_REF] Wedderburn | Quasi-likelihood functions, generalized linear models, and the Gauss-Newton method[END_REF]:

.

There is in fact a close connection between GEE and quasi-likelihood (Liang & Zeger, 1986, p. 21).

Programs written in R [START_REF] Ihaka | R: a language for data analysis and graphics[END_REF], a freeware dialect of the S language [START_REF] Becker | The New S Language[END_REF]


Dr Lisa Jacquin
email: lisa.jacquin@univ-tlse3.fr

A Gandar

M Aguirre-Smith

A Perrault

M Le Hénaff

L De Jong

S Paris- Palacios

P Laffaille

S Jean

Paul Sabatier

High temperature aggravates the effects of pesticides in goldfish

Keywords: pesticides, multiple stressors, climate change, freshwater fish, synergic effects, physiology

In human-altered rivers, fish are often conjointly exposed to an increase in water temperature due to global warming and to contamination by organic pollutants such as pesticides, but their combined effects are still elusive. Thermal and chemical stressors could potentially interact because high temperature increases metabolism and toxicant uptake, and can alter the ability of organisms to set up adequate stress responses and to maintain homeostasis. These combined stressors could thus potentially result in higher level of molecular and cellular damage, and stronger effects on behavior and physiology, but experimental evidence across biological levels is still scarce. In this study, goldfish Carassius auratus were experimentally exposed to an environmentally realistic cocktail of pesticides (S-metolachlor, isoproturon, linuron, atrazinedesethyl, aclonifen, pendimethalin and tebuconazol) commonly found in rivers of South-West of France at low or high dose in two different thermal conditions: a common summer temperature (22°C) or a high temperature recorded during heat waves (32°C). Results showed that high temperature alone caused behavioral and physiological changes (increased swimming activity, increased hepatosomatic index, decreased reproductive index) but limited cellular damage. However, high temperature aggravated the effects of pesticides at the molecular and cellular level. Indeed, pesticide exposure resulted in higher genotoxic effects (micronuclei rate) and irreversible cellular damage of the gills and liver (apoptosis, inflammation, necrosis) at 32°C compared to 22°C. This suggests potential synergic effects of climate change and pollution, and highlights the need for multiple stress approaches to better predict the impacts of human activities on aquatic wildlife.

Introduction

In the context of global change, aquatic organisms are exposed to a growing number of environmental stressors that create new selective pressures affecting their behavior, physiology and ultimately their fitness [START_REF] Crozier | Plastic and evolutionary responses to climate change in fish[END_REF][START_REF] Schinegger | Configuration of multiple human stressors and their impacts on fish assemblages in Alpine river basins of Austria[END_REF][START_REF] Schulte | What is environmental stress? Insights from fish living in a variable environment[END_REF]. Although a considerable body of literature addressed the isolated effects of different stressors on aquatic organisms, few studies investigated their combined effects, thereby hindering our ability to predict the effects of human-induced perturbations in realistic environmental conditions [START_REF] Holmstrup | Risk assessment of mixtures of chemicals and combinations of chemicals and natural stressors[END_REF][START_REF] Matthaei | Multiple stressor effects on freshwater fish: a review and metaanalysis[END_REF][START_REF] Nõges | Quantified biotic and abiotic responses to multiple stress in freshwater, marine and ground waters[END_REF][START_REF] Schinegger | Configuration of multiple human stressors and their impacts on fish assemblages in Alpine river basins of Austria[END_REF][START_REF] Schinegger | Untangling the effects of multiple human stressors and their impacts on fish assemblages in European running waters[END_REF]. In human-altered rivers, fish are often conjointly exposed to an increase in water temperature due to global warming, especially during heat waves, and to a concomitant contamination by organic pollutants such as pesticides, which are widely used in urban and agricultural areas [START_REF] Page | Environmental monitoring of selected pesticides and organic chemicals in urban stormwater recycling systems using passive sampling techniques[END_REF]Polard et al., 2011a). In most cases, rivers were contaminated by complex mixtures of pollutants, including several pesticides, which can induce different effects [START_REF] Hayes | Pesticide Mixtures, Endocrine Disruption, and Amphibian Declines: Are We Underestimating the Impact?[END_REF][START_REF] Relyea | A cocktail of contaminants: how mixtures of pesticides at low concentrations affect aquatic communities[END_REF][START_REF] Wang | Single and joint toxic effects of five selected pesticides on the early life stages of zebrafish (Denio rerio)[END_REF]. In addition, thermal and chemical stressors can potentially interact and alter the ability of organisms to set up adequate defensive stress responses and to maintain homeostasis (Gandar et al., 2017a;[START_REF] Nadal | Climate change and environmental concentrations of POPs: A review[END_REF][START_REF] Noyes | The toxicology of climate change: Environmental contaminants in a warming world[END_REF][START_REF] Schiedek | Interactions between climate change and contaminants[END_REF]. Indeed, high temperatures can directly increase the intake and toxicity of organic molecules [START_REF] Laetz | Elevated temperatures increase the toxicity of pesticide mixtures to juvenile coho salmon[END_REF][START_REF] Nadal | Climate change and environmental concentrations of POPs: A review[END_REF]. In addition, high temperature triggers costly stress responses, accelerate metabolism and deplete energy reserves, which could indirectly affect ability of organisms to set up efficient defense responses against chemicals and to limit cellular damage (Gandar et al., 2017a[START_REF] Gandar | Multistress effects on goldfish (Carassius auratus) behavior and metabolism[END_REF][START_REF] Chris | Stress syndromes: Heightened bioenergetic costs associated with contaminant exposure at warm temperatures in teleosts[END_REF].

Increased temperature can be particularly deleterious for ectothermic species as temperature affects most of their behavioral and metabolic processes [START_REF] Manciocco | Global warming and environmental contaminants in aquatic organisms: the need of the etho-toxicology approach[END_REF]. As a consequence, thermal stress might amplify or mask the effects of toxic stressors in fish, making it difficult to predict their combined effects at different levels of biological organization, from the cell to the whole organism.

To tackle this question, we used goldfish Carassius auratus as a model species because of its extensive use in ecotoxicological tests and for its responsiveness to chemical and natural stressors [START_REF] Cavas | In vivo genotoxicity evaluation of atrazine and atrazine-based herbicide on fish Carassius auratus using the micronucleus test and the comet assay[END_REF]Gandar et al., 2017b;[START_REF] Maisano | Haemolytic activity and characterization of nematocyst venom from Pelagia noctiluca (Cnidaria: Scyphozoa)[END_REF][START_REF] Parrino | Comparative study of haematology of two teleost fish (Mugil cephalus and Carassius auratus) from different environments and feeding habits[END_REF]. It is non-native species from Asia, that can now be found in a wide range of lentic and lotic systems across Europe, including in the Garonne watershed, South-West of France [START_REF] Keith | Les Poissons d'eau douce de France[END_REF]. Goldfish were experimentally exposed to an environmentally realistic cocktail of common pesticides (S-metolachlor, isoproturon, linuron, atrazine-desethyl, aclonifen, pendimethalin and tebuconazol) at low and high concentrations observed in rivers of South-West of France [START_REF] Polard | Caractérisation des effets génotoxiques sur poissons de produits phytosanitaires en période de crue[END_REF]Polard et al., 2011a). This experiment was conducted with a temperature of 22°C, currently registered in rivers of South-West of France and corresponding to the optimal temperature for goldfish, and with a high temperature (32°C), which is within the upper tolerance limit of goldfish (0.3°C-12.6°C for the lower limit and 30.8°C-43.6°C for the upper limit, [START_REF] Ford | Temperature tolerance in the goldfish, Carassius auratus[END_REF][START_REF] Fry | The relation of temperature to oxygen consumption in the goldfish[END_REF]. The chosen high temperature of 32°C is a temperature observed in the Garonne river (located in South-West France) during extreme weather events such as heat waves, which will increase in intensity and frequency in the future due to climate change (IPCC, 2013). Fish were exposed to a high temperature or pesticides alone or to both stressors in combination.

To understand better the effects of chemical and thermal stress, alone or together, the behavioral and physiological responses of fish were analyzed at different levels of biological organization. Pesticides generally induce pathological damage at the molecular, cellular and individual levels, depending on the stress intensity and the duration of exposure. At the molecular level, genotoxic effects of pesticides are often observed, with alterations of DNA structure that can result in the production of micronuclei and alter the cell cycle [START_REF] Cavalcante | Genotoxic effects of Roundup on the fish Prochilodus lineatus[END_REF]Polard et al., 2011). Previous proteomic studies on Carassius auratus showed that the cocktail of pesticides used in this study affected proteins involved in cellular inflammatory responses, oxidative stress, and in cytoskeleton and cell structure (e.g. Gandar et al., 2017b;[START_REF] Marchand | Molecular Identification and Expression of Differentially Regulated Genes of the European Flounder, Platichthys flesus, Submitted to Pesticide Exposure[END_REF]. However, the consequences at the cellular level are not known yet.

Oxidative stress and proteome changes caused by pesticides are expected to result in apoptosis or accelerated cell cycle, tissue lesions and potentially carcinogenesis, especially in the liver (Gandar et al., 2017b;[START_REF] Ratn | Phorate induced oxidative stress, DNA damage and differential expression of p53, apaf-1 and cat genes in fish, Channa punctatus (Bloch, 1793)[END_REF][START_REF] Weber | Transcriptome Alterations Following Developmental Atrazine Exposure in Zebrafish Are Associated with Disruption of Neuroendocrine and Reproductive System Function, Cell Cycle, and Carcinogenesis[END_REF][START_REF] Wirbisky | Embryonic atrazine exposure alters zebrafish and human miRNAs associated with angiogenesis, cancer, and neurodevelopment[END_REF].

In this study, we thus focused on histopathological modifications of the gills and liver, because they are the two primary target organs affected by organic pollutants (e.g. [START_REF] Nero | Gill and liver histopathological changes in yellow perch (Perca flavescens) and goldfish (Carassius auratus) exposed to oil sands process-affected water[END_REF].

Gills are an important interface between the external and internal environment and can be a route of contaminant uptake in fish [START_REF] Evans | The fish gill: site of action and model for toxic effects of environmental pollutants[END_REF]. Liver is the primary organ for detoxification and proteomic studies showed differential levels of proteins involved in hepatocyte apoptosis and cirrhosis in goldfish exposed to pesticides (Gandar et al., 2017b).

Given previous proteomic studies (Gandar et al., 2017b), we hypothesized that pesticide exposure would alter liver and gill structure and result in cellular lesions (such as necrosis, inflammatory lesions), especially at high temperature. We also hypothesized that pesticides would reduce gonad mass (often measured by Gonadosomatic Index), liver mass (Hepatosomatic Index) and body condition, because stressful conditions often decrease investment in reproduction [START_REF] Schreck | Stress and fish reproduction: The roles of allostasis and hormesis[END_REF] and lead to reallocation of energy reserves [START_REF] Gandar | Multistress effects on goldfish (Carassius auratus) behavior and metabolism[END_REF].

At the individual level, temperature and organic pollution also often trigger immune changes (here measured as changes in innate immunity) [START_REF] Fatima | Combined effects of herbicides on biomarkers reflecting immune-endocrine interactions in goldfish: Immune and antioxidant effects[END_REF] and behavioural changes (activity, feeding, exploration) [START_REF] Gandar | Multistress effects on goldfish (Carassius auratus) behavior and metabolism[END_REF][START_REF] Jacquin | Evolutionary and immediate effects of crude-oil pollution: depression of exploratory behaviour across populations of Trinidadian guppies[END_REF], although their combined effects are still elusive. Based on previous metabolomic approaches showing synergistic effects of temperature and pesticides on goldfish immune proteins, energy allocation and maintenance (Gandar et al., 2017b(Gandar et al., , 2017a[START_REF] Gandar | Multistress effects on goldfish (Carassius auratus) behavior and metabolism[END_REF], we predicted that pesticides would have stronger effects at 32°C on immunity and behaviour compared to 22°C.

Materials and Methods

Model species

The goldfish (Carassius auratus, Cyprinidae family) was chosen as a biological model for this study because of its extensive use in ecotoxicological tests (e.g. Bretaud et al. 2000, Cavas and Konen 2007, Feng et al. 2013). Introduced in Europe in the eighteenth century for ornamentation, it is a non-native species from Asia that can now be found in a wide range of lentic and lotic systems across Europe, including in the studied Garonne watershed, South West of France [START_REF] Keith | Les Poissons d'eau douce de France[END_REF]. The temperature tolerance of the goldfish is large (0.3°C-12.6°C for the lower limit and 30.8°C -43.6°C for the upper limit) and its optimum temperature is 25°C. [START_REF] Ford | Temperature tolerance in the goldfish, Carassius auratus[END_REF][START_REF] Fry | The relation of temperature to oxygen consumption in the goldfish[END_REF]. A total of 156 individuals of 10-12cm in length were purchased from the fish farm Carpio (Consac, France). They were acclimatized for two weeks in opaque tanks under controlled conditions (18°C with a 12:12h light regime). They were then gradually acclimatized to common (22°C) or high (32°C) experimental temperatures. The temperature was increased by 0.5°C every 12 hours until reaching the experimental temperature following previous studies [START_REF] Gandar | Multistress effects on goldfish (Carassius auratus) behavior and metabolism[END_REF].

Water was aerated and dechlorinated prior to fish introduction. Half of the water was renewed every day and fish were fed daily with commercial pellets (Carpio, France).

Experimental treatments

A total of 156 fish (26 fish in each experimental group) were placed in 30L opaque aquaria (22 and 32°C). Temperature regime was controlled using heaters coupled with thermostats.

Temperature was monitored every 24h to ensure constant conditions throughout the exposure period. The 22°C temperature was chosen as a control temperature since it is close to goldfish optimal temperature range, while 32°C corresponds to a stressful temperature In each temperature group, fish were exposed either to control water (C), or to a low (L) or high dose (H) of a cocktail of pesticides during 96h, resulting in 6 experimental groups of 26 fish (C22, L22, H22, C32, L32 and H32). The cocktail of pesticides contains seven most common pesticides: 6 herbicides (S-metolachlor, isoproturon, linuron, atrazine-desethyl, aclonifen and pendimethalin) and one fungicide (tebuconazol), which are mainly used for intensive agriculture (corn, wheat and sunflower) (details in Table 1). This pesticide composition and concentrations were selected according to preexisting data of a typical pesticide contamination in the rivers of the South-West of France (Polard, 1.

Prior to fish introduction, aquaria and all equipment were pre-saturated with the mixture of pesticides. Air pumps were placed in each aquarium to supply oxygen in water.

Half of water and pesticide solutions was renewed every day to limit concentration variations.

Previous studies showed that this protocol successfully maintained pesticide concentrations in the water and that good oxygen and pH conditions across the experiment [START_REF] Gandar | Multistress effects on goldfish (Carassius auratus) behavior and metabolism[END_REF](Gandar et al. , 2017)). No mortality occurred during the exposure. At the beginning of the experiment at T0, groups did not differ in body mass (χ 2 = 1.51, p= 0.91). At the end of the experiment and after the behavioral tests, fish were weighed (to the nearest 10 mg) and measured (fork length to the nearest mm). Fish were anesthetized and euthanized with benzocaine (80 -120 mg.L -1 ) and concussion.

Molecular level: DNA alterations DNA alterations were estimated on blood smears by counting the Micronuclei (MN) rate in erythrocytes following the protocol described in Polard et al. (2011b) on half of the individuals of each treatment group (N=12 to 14 randomly chosen slides per group). Micronuclei rate is used to quantify the genotoxicity of treatments, which might arise through mutagenic of epigenetic processes [START_REF] Al-Sabti | Fish micronuclei for assessing genotoxicity in water[END_REF][START_REF] Luzhna | Micronuclei in genotoxicity assessment: from genetics to epigenetics and beyond[END_REF][START_REF] Torres De Lemos | River water genotoxicity evaluation using micronucleus assay in fish erythrocytes[END_REF].

Blood was sampled in the caudal vein using 1mL Terumo syringe with a 26G needle and one drop was fixed on slides. Blood slides were fixed in pure ethanol for 15 min and dried for 24 hours. All slides were then stained with fluorescent acridine orange AO (0.003% in Dulbecco PBS) (Polard et al. 2011). Micronuclei yellow-green fluorescence was counted under epi-fluorescence with a U-MWB2 filter at 1000 magnification. Micronuclei rates were calculated as the number of micronuclei per 5000 blood cells.

Cellular level: histopathological alterations of the gills and liver

Histopathological lesions were evaluated in the liver (26 fish per experimental group, 3 observations per individual) and in the gills (3 females per experimental group, 30 observations per individual, 10 measurements per filament).

Gill and liver damage was assessed using standard histological techniques (Fernandes andMazon, 2003, Bernet et al., 1999;[START_REF] Costa | Estuarine ecological risk based on hepatic histopathological indices from laboratory and in situ tested fish[END_REF][START_REF] Cuevas | Histopathological baseline levels and confounding factors in common sole (Solea solea) for marine environmental risk assessment[END_REF][START_REF] Van Dyk | Histological changes in the liver of Oreochromis mossambicus (Cichlidae) after exposure to cadmium and zinc[END_REF].

Organs were dissected and immediately fixed in Bouin's fluid (picric acid 3V, formaldehyde 1V

and acetic acid 2V) for 24h. They were dehydrated through a graded ethanol series (70-100 %),

following by Butanol-1ol bath, cleared in toluene, and embedded in paraffin. Longitudinal serial sections of 10 μm (gills) or 6 µm (liver) were mounted on glass slides, stained with hematoxylineosin HE (gills) or with nuclear fast red NR and picro-indigo-carmin PIC (liver), and observed under Nikon microscope (10X, 40X and 60X). Microphotographies were taken with a light microscope (Laborlux S -Leitz wetzlar Germany) equipped with a digital camera (Imagingsource) and analyzed using Archimed software (Version 5.6.0.).

Gill remodelling was assessed using the percentage of exchange surface area (ESA), corresponding to protruding length of secondary lamellae (10 measurements per filament in the middle of the arch, i.e. 30 measurements per sample, hence, 90 measures per conditions) and was calculated as: ESA (%) = 100*(mean SL -mean ILCM) / mean SL) where: SL= length of secondary lamellae; ILCM = length of interlamellar cell mass. In addition, five filaments located in the middle of the arch (i.e. 5 measures per fish, namely 15 measures per condition) were analysed following [START_REF] Bernet | Histopathology in fish: proposal for a protocol to assess aquatic pollution[END_REF] to calculate the degree of gill damage. Six different types of alterations were recorded (details in Table S1). A gill pathological index (I Gill ) was calculated as: I Gill = ∑ alt (a alt x w alt ) where: a = score value (depending on occurrence) for a specific alteration (i.e. type of lesion), ranging from 0 to 6; w = importance factor for a specific alteration

(1: reversible lesions with low effects; 2. modified cellular architecture: 3: degenerative irreversible alterations leading to a loss of gill function),; alt = specific alteration (see Table S1).

Six main hepatic alterations were recorded in the liver (on 3 random areas of 25 mm²):

lysis-necrosis, mitosis, melano-macrophage aggregates, fibrosis, nuclear alterations, and lipidosis. A global histopathological score was calculated as the mean value of the different scores (details in Table S1).

Organ level: GSI, HSI and body condition

Sex and maturity were determined by gonadal examination. The Hepatosomatic Index HSI representing the energy content of the liver and the Gonadosomatic Index GSI representing the gonad growth status were calculated as liver or gonad mass/total body mass. Peroxidase activity was measured following protocols from [START_REF] Quade | A rapid, direct assay to measure degranulation of bovine neutrophil primary granules[END_REF].

Tetramethylbenzidine (TMB) and hydrogen peroxyde H 2 O 2 were used as enzymatic substrates.

A total of 12.5 µL of plasma was diluted in 0.5 mL of Hank buffer (HBSS). A total of 60 µL of this dilution were placed in a 96-microplate with a 100 µL solution of TMB, DMSO, Na 2 HPO 4 /citric acid 50 mM (pH = 5.0) buffer. After 2 minutes of incubation, 50 µL of 2 M sulfuric acid were added to stop the reaction. The absorbance was measured at 450 nm using a spectrophotometer, and expressed as UA with: UA = A 450 « sample » -A 450 « control » and used as a measure of peroxidase activity.

Because immunity can differ among individuals before the experiment depending on past exposure of fish to disease, we also measured their lysozyme and peroxidase levels in the plasma before the experiment to control for individual immune differences before experimental treatments. Lysozyme or peroxidase variations across the experiment were then calculated as (value after treatment-value before treatment)/value before treatment.

Individual level: behavior

We measured activity, mean speed, trajectory changes, exploration and feeding behaviour at T96

(96 hours post pesticide exposure). These behavioural traits were chosen because they are often affected by stressors and are tightly linked to fish fitness and functional effects (e.g.

Scott and Sloman, 2004, Gandar et al. 2015).

Fish were not fed during the exposure to pesticides to standardize their physiological state and appetite.

First, we recorded the "activity" of fish in the home tank as the total time spent swimming during 5 minutes (i.e., in directed motion and covering a distance > 10% of body length). Activity was recorded by direct observation in the home tank with a companion fish to limit isolation stress and reliably reflects the intrinsic activity of the fish. Time of the day of observation was randomized for each group to avoid any confounding effect of diurnal rhythm.

Second, we recorded more specific features of the trajectory of the fish such as "mean speed", "trajectory changes", and "exploration", because contaminants can have neurotoxic effects that can result in erratic behaviours, affect trajectories and exploration patterns (e.g. [START_REF] Jacquin | Evolutionary and immediate effects of crude-oil pollution: depression of exploratory behaviour across populations of Trinidadian guppies[END_REF]. For this, we transferred each individual in a test tank (50x30cm with 10cm of water) for 10 minutes with the same conditions as the exposure tank. Each focal fish was placed in a transparent cylinder in the middle of the tank. After 2 minutes of acclimation the cylinder was gently removed, using a string on a reel to minimize disturbance. Movements were then recorded precisely for 5 minutes using a camera (Logitech HD Webcam C270) the VirtualDub software and the Tracker software (Tracker 5.0, Open Source Physics 2015) to obtain precise coordinates of the position of the fish at each time point (x, y). This allowed us to obtain the mean speed of the fish (cm/s). The number of trajectory changes was also recorded as the number of rapid (in less than 5 sec) and abrupt (180° and > 100 % of body length) changes of direction.

Lines were drawn on the bottom of the tank to delimit 15 zones of 10 x 10cm.

Exploration was recorded as the number of different zones reached by the fish during 5 minutes [START_REF] Jacquin | Evolutionary and immediate effects of crude-oil pollution: depression of exploratory behaviour across populations of Trinidadian guppies[END_REF]. The position of the different test tanks was changed every 5 fish to avoid any confounding effect of a spatial preference. At the end of the behavioural test, we introduced 10 food pellets in the test tank and recorded the number of pellets consumed during one minute and the latency before feeding to compare the feeding behaviour between groups.

Compliance with ethical standards

All experimental procedures were conducted in compliance with French and European legislation for animal experimentation (European directive 2010/63/UE). Experimental procedures were conducted under the French animal handler's certificate (N° 31-103), the establishment approval for vertebrate experimentation N°A3113002, and were approved by the ethical committee 073 (authorization N°8538).

Statistical analyses

Effects of stressors (temperature or pesticides or both) were evaluated by Generalized Linear Models for parametric data with temperature, pesticide dose, and their interaction as explaining factors. Sex was also added in the initial models but was not retained in any model (all p values > 0.12), except for Gonadosomatic Index where separate analyses were conducted on males and females. Best models were chosen using a backward selection and successively dropping non-significant terms (p > 0.10). Posthoc tests were conducted using Wilcoxon tests.

For non-parametric data, Kruskall-Wallis tests were conducted to compare groups. For feeding behavior, an additional analysis combining feeding latency and feeding status (at least one pellet eaten or no pellet eaten) was analyzed using survival curve differences (survdiff function, Kaplan-Meier method).

Results

DNA alterations and micronuclei

Micronuclei rate significantly differed between experimental groups (Kruskall-Wallis test:

Chisq=14.9, P=0.011). Fish at 22°C exposed to a high dose of pesticides had a marginally higher micronuclei rate than the control group (C22 vs H22: Wilcoxon test: W=47, P=0.057) (Fig 1). Fish at 32°C exposed to a high dose of pesticides had a significantly higher rate of micronuclei compared to the C22 control group (Wilcoxon test, C22 vs H32: W=47, P=0.029) and marginally higher compared to the C32 group (H32 vs C32: W=56, P=0.094) (Fig 1).

Histopathological alterations of the gills

Histopathological alterations of the gills significantly differed between groups (details in Table S2). S2). At 22°C with pesticides in the L22 group, ILCM length decreased, and partial fusion of secondary lamellae (Fu) and abundant mucus occurred (L22, Fig. 2). In the H22 group, the ILCM was very reduced (H22, Fig 2 ), resulting in a significant increase of exchange surface area ESA, that reached 82 % for H22 fish (Table S2).

At high temperature in the C32 group, the elevation of temperature to 32°C also induced a decrease of ILCM in C32 fish (C32, Fig 2) resulting in a significant increase of exchange surface area ESA (77.5 ± 4.4 % in C32 versus 5.0 ± 11.5% in C22, Table S2), and protective mucus was present. Protective mucus was absent in fish exposed both to 32°C and to pesticide mixtures (L32, H32 Fig 2). Lamellar hyperplasia (Hy) was present and pesticides induced lifting of lamellar epithelium (Li), and oedema (Oe) in H22 fish (Fig 2; Table S2). At 32°C with high dose of pesticides (H32), irreversible lesions appearing frequently, especially necrosis (Ne), and oedema (Oe) (Fig 2; Table S2).

As a result, the integrative gill pathological index (I Gill ) thus increased in 32°C group compared to 22°C group and in fish exposed to a high dose of pesticides at 22°C and 32°C (Fig 3A). Fish exposed to both stressors (H32) exhibited a higher level of gill pathological index as compared to all other groups of fish, reflecting a highly degraded cellular structure of the gills (Fig 3A ).

Histopathological alterations of the liver

High temperature alone had no significant effect on fish hepatic structure (C32 vs C22) (Fig 3B).

In contrast, exposure to pesticides induced fibrosis, an increase of the immune response and hepatocyte nuclear alterations typical of an inflammatory syndrome (details in Supplementary Figure S3 and S4), resulting in a higher histopathological score in L22 and H22 fish compared to C22 (Fig 3B). At high temperature, most fish exposed to pesticides (L32 and H32) showed disorganization of some hepatic areas, i.e. an increase in the number of hepatocytes involved in mitosis and inflammation, lipid droplets and a proliferation of blood vessels (details in Supplementary Figure S4). A total of 80% of fish exposed to the high concentration of pesticides and to 32°C (H32) had a significant inflammatory response in the liver (Fig S4). The global histopathologic index shows that perturbation of the liver due to pesticides was more severe at 32°C compared to 22°C (Fig 3B ).

HSI, GSI, body mass loss

The Hepatosomatic Index (HSI) was affected by temperature only (Table 2), with fish exposed to 32°C having a higher HSI than fish exposed to 22°C (t= 2.10, P= 0.042). The Gonadosomatic Index (GSI) of females (but not in males) was significantly affected by temperature and pesticides (Table 2). More precisely, at 22°C, pesticides decreased female GSI at low dose (W=136, P=0.026), and marginally at high dose (W=108, P=0.10), whereas no effect of pesticides was detected at 32°C (all p values >0.12) (Fig S5). Body mass loss was affected by temperature only (Table 2), with fish exposed to 32°C losing more body mass than fish exposed to 22°C (t= -2.63, P= 0.0093).

Immunity

Before exposure to pesticides at T0, there were no differences in lysozyme (χ 2 = 3.69, P= 0.59)

or peroxidase activity (χ 2 = 3.37, P= 0.64) among groups. Lysozyme variations over the experiment were not affected by treatments (Table 2). In contrast, peroxidase activity variations were affected by both temperature and pesticides and their interaction (Table 2).

More precisely, posthoc tests show that 32°C fish exposed to a low dose of pesticides had a higher increase in peroxidase compared to 22°C and 32°C control groups (Fig S6 ).

Behaviour

Activity at T96 was affected by pesticide exposure and by temperature (Table 2), with fish exposed to a high dose of pesticides being less active than control fish whatever the temperature (estimate = -46.63 ± 16.45, t= 2.83, P= 0.0052), and fish from the 32°C groups being more active than fish from the 22°C groups (estimate = 8.34 ± 1.35, t= 6.19, P< 0.001). More specifically, posthoc tests show that at 22°C, fish exposed to a high dose of pesticides were marginally less active than control fish (W= 439, P= 0.064), whereas at 32°C, only fish exposed to a low dose of pesticides were significantly less active than control fish (W= 479, P= 0.0038) (Fig 4A). There was no difference in mean speed, trajectory changes, and exploration between pesticide groups and these behaviours were only affected by temperature (Table 2). More specifically, fish exposed to a high temperature swam faster (estimate = 0.054 ± 0.27, t= 2.00, P= 0.047), changed trajectory more often (estimate = 0.52 ± 0.21, t= 2.42, P= 0.017), and explored more zones compared to fish acclimated to 22°C (estimate = 0.11 ± 0.040, t= 2.59, P= 0.010).

Feeding behavior (number of pellets eaten) was affected by both temperature and pesticides (Table 2). Fish at 32°C ate more pellets than 22°C fish (estimate = 0.029 ± 0.012, P=0.017) (Fig 4B). Fish exposed to a high dose of pesticides ate marginally more than control fish (estimate = 0.25 ± 0.14, P=0.074) (Fig 4B ).

A complementary analysis integrating both the latency to feed and the feeding status (fish eating at least one pellet or no pellet) confirmed that feeding differed among treatment groups (survival analysis: χ²=11.5, P=0.043): a higher proportion of fish exposed to a high temperature only (C32) and to pesticides (H22, L32 and H32) ate at least one pellet and fed more rapidly than control fish C22 (Supplementary Figure S7).

Discussion

The trend of results obtained in this study indicates that environmentally relevant concentrations of pesticides affected goldfish physiology and behaviour, especially at the molecular and cellular levels. Cellular damage induced by pesticides was more severe at 32°C compared to 22°C (summary of results in Fig 5 ) suggesting high deleterious effects of multiple thermal and toxic stressors.

Molecular effects of multiple stressors

At the molecular level, DNA damage reflected by micronuclei rate was higher in fish exposed to both pesticides at high dose and high temperature (H32) compared to control groups. This confirms previous studies showing that pesticides have genotoxic effects and induce DNA alterations in various fish species [START_REF] Cavalcante | Genotoxic effects of Roundup on the fish Prochilodus lineatus[END_REF][START_REF] Cavas | In vivo genotoxicity evaluation of atrazine and atrazine-based herbicide on fish Carassius auratus using the micronucleus test and the comet assay[END_REF]Polard et al., 2011a).

Our results further suggest that genotoxic effects of pesticides could be amplified by high temperatures. In general, high temperature itself can increase DNA damage (e.g. [START_REF] Anitha | Genotoxicity evaluation of heat shock in gold fish (Carassius auratus)[END_REF], but Heat Shock Proteins (HSP) and DNA repair mechanisms can often limit such damage depending on the duration and intensity of stressors [START_REF] Iwama | Heat shock protein expression in fish[END_REF].

Accordingly, micronuclei rate did not differ between control fish at 22°C and 32°C in our study. We also observed that a combined exposure to thermal and toxic stressors led to important DNA alterations and high micronuclei rate during cell division. This suggests that fish exposed to multiple stressors probably entered in the exhausted phase of the stress and became unable to set up DNA protection and/or repair mechanisms, or that these mechanisms might be not sufficient to repair all damage. However, the underlying mechanisms remain to be explored.

Histopathological effects of multiple stressors on gills

At the cellular level, several histopathological alterations were also visible in fish exposed to stressors. Increased temperature only had limited effects on gill structure, and mostly resulted in an increase in exchange surface as expected in this species [START_REF] Sollid | Temperature alters the respiratory surface area of crucian carp Carassius carassius and goldfish Carassius auratus[END_REF]. Indeed, goldfish at low temperature have specific interlamellar cell mass (ILCM) between their gill lamellae. This reduces the exchange surface, limits osmoregulation costs and limits the exposure to pathogens and pollutants [START_REF] Sollid | Temperature alters the respiratory surface area of crucian carp Carassius carassius and goldfish Carassius auratus[END_REF][START_REF] Sollid | Plasticity of respiratory structures -Adaptive remodeling of fish gills induced by ambient oxygen and temperature[END_REF]. In addition, goldfish generally show gill remodelling under high temperatures: the ILCM generally decreases at warmer water temperatures (≥ 25° C), thereby increasing the exchange surface area (ESA) [START_REF] Sollid | Temperature alters the respiratory surface area of crucian carp Carassius carassius and goldfish Carassius auratus[END_REF][START_REF] Sollid | Plasticity of respiratory structures -Adaptive remodeling of fish gills induced by ambient oxygen and temperature[END_REF], as observed in our study at 32°C (C32 vs. C22 group). This increase in gill exchange area could in turn increase their exposure to pollutants. Accordingly, in our study, more severe gill pathological lesions due to pesticides were observed in fish exposed to a high pesticide dose at 32°C (H32) compared to fish exposed to the same pesticide dose at 22°C (H22). Protective mucus that helps preserving the integrity of the gill epithelium (Fernandes and Mazon, 2003) was also reduced in fish exposed to pesticides (H22 and H32), which could participate in the lesions observed.

Several lesions were observed, including lamellar hyperplasia, partial lamellar fusion, lifting of lamellar epithelium and oedema that were very frequent in H32 fish exposed to both stressors, and that can ultimately lead to necrosis. This is consistent with previous studies in goldfish exposed to carbaryl and Polycyclic Aromatic Hydrocarbons (PAHs), which also reported partial lamellar fusion and gill epithelial lifting [START_REF] Nero | Gill and liver histopathological changes in yellow perch (Perca flavescens) and goldfish (Carassius auratus) exposed to oil sands process-affected water[END_REF][START_REF] Pfeiffer | Electron microscopic perspectives of gill pathology induced by 1-naphthyl-N-methylcarbamate in the goldfish (Carassius auratus Linnaeus)[END_REF]. Interestingly, in our study, L22 and L32 groups had a similar gill pathological index (I Gill ) compared to controls, suggesting that fish could maintain their gill functions at the lowest pesticide concentrations (for 96 hours at least), although results should be interpreted with caution given the low sample size [START_REF] Wolf | Evaluating the credibility of histopathology data in environmental endocrine toxicity studies: Credibility of histopathology data in environmental studies[END_REF]. However, at the higher pesticide concentrations (H22 and H32), the gill pathological index (I Gill ) was higher, indicating a higher number and severity of lesions. Fish maintained at 32°C and exposed to a high concentrations of pesticides (H32) exhibited the highest I Gill (especially a high occurrence of necrosis), indicating irreversible damage that will impair ionic regulation, excretory and respiratory functions of the gills. Taken together, these results indicate a limited impact of single stressors (high temperature or pesticides alone) but a high impact of multiple stressors (pesticides at high temperature) on gill cellular structure and integrity, with potential deleterious consequences for gill function and fish survival.

Histopathological effects of multiple stressors on the liver

Several histopathological alterations were also visible in the liver of fish exposed to both stressors. The increase of temperature alone had low effects on the hepatic structure of goldfish, while the exposure to the pesticide induced a pathology characterized by a fibrosis, melano-macrophage centers and hepatocyte nuclear alterations suggesting an inflammatory syndrome, as well as an accumulation of lipids in hepatocytes. Most fish also showed disorganization of some hepatic areas with an increase in the number of hepatocytes involved in mitosis and a proliferation in blood vessels, which may indicate a carcinogenic process.

The pathology was however moderate probably due to the short duration of the exposure (96h). The pathology described was more frequent and acute in fish exposed to the highest pesticide concentration, which indicated a dose-dependent effect. These results are consistent with previous studies indicating that pesticides can trigger hepatic perturbations [START_REF] Paris-Palacios | Metallothionein induction related to hepatic structural perturbations and antioxidative defences in roach (Rutilus rutilus) exposed to the fungicide procymidone[END_REF], 2000;[START_REF] Paris-Palacios | Hepatocyte nuclear structure and subcellular distribution of copper in zebrafish Brachydanio rerio and roach Rutilus rutilus (Teleostei, Cyprinidae) exposed to copper sulphate[END_REF] and cause disorders of fatty acid metabolism in hepatocytes [START_REF] Liu | Organochloride pesticides impaired mitochondrial function in hepatocytes and aggravated disorders of fatty acid metabolism[END_REF]. This is also consistent with a previous study reporting proteome disturbances reflecting liver oxidative stress, hepatic steatosis and necrosis in goldfish exposed to the same cocktail of pesticides (Gandar et al., 2017b). The sequestration of lipophilic pollutants and the perturbations of liver metabolism could lead to the increased number of intracellular lipid droplets observed in hepatocytes of H32 fish. This could indicate hepatosteatosis and is consistent with the hepatotoxicity and non-alcoholic "fatty liver" syndrome observed in vertebrates exposed to pesticides (e.g. (Gandar et al., 2017b;[START_REF] Closs | Effects of the herbicide atrazine in neotropical catfish (Rhamdia quelen)[END_REF][START_REF] Mesnage | Multiomics reveal nonalcoholic fatty liver disease in rats following chronic exposure to an ultra-low dose of Roundup herbicide[END_REF][START_REF] Xing | Oxidative stress response and histopathological changes due to atrazine and chlorpyrifos exposure in common carp[END_REF].

Interestingly, this study shows that the histopathologic perturbation of the liver was more severe when the pesticide exposure was combined to a rise of temperature, suggesting that high temperature aggravates the effects of pesticides on hepatocytes and liver integrity.

Although histopathological data need to be analyzed with caution [START_REF] Wolf | A critical review of histopathological findings associated with endocrine and non-endocrine hepatic toxicity in fish models[END_REF], our results are consistent with a previous study showing different liver proteome profiles in goldfish under multiple stressors compared to a single stressor (Gandar et al., 2017b). This previous study outlined that pesticides at high temperature modified the amount of proteins involved in regulation of oxidative stress, cell cycle and cell death control (Gandar et al., 2017b), which were predicted to cause the observed liver lesions observed in this study.

Multiple stressors effects at the individual level

At higher levels of biological organization (organ and individual level), gonad mass was reduced in females exposed to a high temperature or pesticides, which is consistent with a reduction in reproductive investment in stressed fish [START_REF] Matthaei | Multiple stressor effects on freshwater fish: a review and metaanalysis[END_REF][START_REF] Pankhurst | Effect of stress on reproduction and growth of fis[END_REF][START_REF] Schreck | Stress and fish reproduction: The roles of allostasis and hormesis[END_REF]. Goldfish exposed to both stressors (L32 and H32)

did not have a higher reduction in gonad mass than controls C32, likely because temperature already caused a high reduction of gonad mass. At the individual level, behavioral and immune changes were observed. Activity (as well as speed, trajectory changes and exploration) were increased by temperature as expected by an increase in metabolism and behavioral changes caused by enhanced temperatures [START_REF] Canal | Influence of temperature on surface sediment disturbance by freshwater fish: a microcosm experiment[END_REF][START_REF] López-Olmeda | Thermal biology of zebrafish (Danio rerio)[END_REF]. This is likely due to increased metabolism and higher energy mobilization at high temperature, as expected from previous studies [START_REF] Gandar | Multistress effects on goldfish (Carassius auratus) behavior and metabolism[END_REF].

Accordingly, increases in feeding behavior, body condition and HSI were observed in fish exposed to a high temperature.

Interestingly, activity was reduced by pesticides, which is consistent with the hypoactivity syndrome that is often observed in freshwater fish exposed to herbicides ( [START_REF] Scott | The effects of environmental pollutants on complex fish behaviour: integrating behavioural and physiological indicators of toxicity[END_REF][START_REF] Shinn | Behavioral response of juvenile rainbow trout exposed to an herbicide mixture[END_REF], and could explain the reduction in sediment surface reworking observed in goldfish exposed to the same cocktail of pesticides in a previous study [START_REF] Gandar | Multistress effects on goldfish (Carassius auratus) behavior and metabolism[END_REF]. However, no other behavioral changes (speed, trajectory, exploration) were observed in response to pesticides. It is not clear yet whether the environmental concentrations chosen have neurotoxic effects and further studies are needed to test neurotoxic effects. However, no decrease in body condition or Hepatosomatic Index were observed under pesticide exposure, which could be explained by the short duration of the exposure (96h). Alternatively, lipids might have been reallocated from other life-history traits (growth, reproduction) and/or organs (for instance gonads) to maintain liver mass.

In addition, pesticides did not trigger a stronger change in behavior under high temperature than under normal temperature, maybe because temperature effects masked the effects of pesticides, or because different metabolic strategies were set up under different temperatures.

Further studies specifically testing the allocation of energy in different tissues and function would be useful to identify the metabolic strategy used by goldfish exposed to multiple stressors.

Conclusions

To conclude, a single thermal stress (increased temperature alone) caused behavioral and allocation adjustments (increased activity, feeding and reduced reproductive investment) but did not result in deleterious effects at the cellular level (no deleterious effect on DNA, gills or liver structure). However, high temperature aggravated the effects of pesticides on fish gill and liver. Indeed, pesticides at low temperature (22°C) had limited deleterious effects on behavior, gills and liver structure, but stronger deleterious effects at 32°C. More specifically, combined thermal and toxic stressors resulted in high level of irreversible damage of the gills and liver (apoptosis, inflammation, necrosis), with potential consequences for fish osmoregulation, detoxification and survival. This suggests that high water temperatures could enhance the negative effects of pesticides on freshwater fish. Although the chosen model species Carassius auratus is a non-native species in Europe, it is now widely found in Europeans rivers and lakes, and has a higher range of temperature tolerance than several native European fish species. It is thus likely that several other wild freshwater fish species might be severely impacted by combined thermal and pesticide stressors.

Given that climate change will likely increase the duration and intensity of heat waves in the next decades (IPCC 2014), this raises concerns about potential synergic effects of climate change and pollution, and highlights the need for multiple stress approaches to better predict the impacts of human activities on aquatic wildlife.

  within the upper tolerance limit of goldfish. The 22°C temperature is a common summer temperature in rivers in South-West of France. The 32°C temperature corresponds to a high river temperature recorded during heat waves ( data from Agence de l'Eau Adour-Garonne, Croze et al., 2007), which are expected to increase in frequency and intensity in the next decades due to climate change (IPCC 2013).

  2010; Polard et al., 2011a; Taghavi et al. 2011). Two concentrations were chosen: a low dose of 8.4µg.L -1 (L groups) and a high dose of 42µg.L -1 (H groups). The Low Dose corresponds to a typical pesticide concentration found after spring flood in French rivers (Debenest et al., 2010; Garmouma et al., 1998; Polard et al., 2011a; Taghavi et al., 2011). For instance atrazine and related molecules in the Save river after spring floods. Data from the French national database on water pollution show that these concentrations are in the range of what can be found in other agricultural French rivers (IFEN 2007). The high concentration (H: High dose of 42µg.L -1 ) is five times higher and reflects highest levels of contamination measured in streams and reservoirs of Europe and USA (Battaglin et al., 2000; Graymore et al., 2001; Kreuger, 1998). Details of these concentrations for each pesticide are in

  Body condition was measured by the Fulton condition index as (body mass/ body length 3 ). Individual level: immunity After centrifugation (10min at 3000 rpm), plasma was frozen for subsequent immune analyses. Innate immunity was assessed by measuring lysozyme activity (which reflects the ability of the individual to eliminate bacteria), and peroxidase activity of the plasma (which is used as a marker of leucocyte activation, Eslamloo et al., 2014; Fatima et al., 2007; Ghafoori et al., 2014; Magnadóttir, 2006). To measure lysozyme activity, 8 mg of Micrococcus lysodeikticus were suspended in 40 mL of Na 2 HPO 4 buffer/ 50 mM of citric acid(pH = 5.5). A total of 15 µL of plasma were placed in 96 wells microplates, with 150 µL of revealing solution. Optical density was read at 450 nm from 5 to 60 minutes. The absorbance was measured at 450 nm using a spectrophotometer at an interval of 30 s (total measuring time 60 min). The unit of lysozyme activity is defined as the concentration of sample causing a decrease in turbidity of 0.001/ min (1 UA = Δ 0.001 A 450 /min).

  Figure 2 shows the different types of changes and lesions observed in the gills of exposed fish. At 22°C in the control C22 group, no tissue injuries were observed and the lamellae were embedded in a thick ILCM (C22, Fig 2), resulting in a low Exchange Surface Area (ESA = 5.0 ± 11.5; Table

Figure 4 .

 4 Figure 4. Routine behaviours, i.e. A. activity (time spent swimming) and B. feeding (number of pellets eaten) of goldfish exposed for 96h at 22°C (left) or 32°C (right) to control water (dark grey), a low dose of pesticides (medium grey), or a high dose of pesticide (light grey) Different letters above bars indicate significant differences between groups.

Figure 5 :FigureFigure S5 .

 5S5 Figure5: Summary of the main effects of temperature and pesticides at 22°C or 32°C on goldfish across levels of biological organization (0: no effect, -: negative effect, +:positive effect) . Taken together, results suggest that pesticides induced more severe molecular and cellular damage on DNA structure, gills and liver cells at 32°C compared to 22°C, while other traits were mainly affected by temperature only. This suggests that high temperature aggravates the deleterious non reversible effects of pesticides on cellular and organ structure in a multiple stressor context

Table 1 .

 1

	Pesticides were obtained from Sigma-Aldrich (St. Louis, MO, USA): S-metolachlor
	(CAS-No: 87392-12-9, PESTANAL®, analytical standard), Isoproturon (CAS-No: 34123-59-6,
	PESTANAL®, analytical standard), Linuron (CAS-No: 330-55-2, PESTANAL®, analytical
	standard), Atrazine-desethyl	(CAS-No: 6190-65-4, PESTANAL®, analytical standard),
	Aclonifen (CAS-No: 74070-46-5, PESTANAL®, analytical standard), Pendimethalin (CAS-No:
	40487-42-1, PROWL®, analytical standard), Tebuconazol (CAS-No:107534-96-3,
	PESTANAL®, analytical standard). Pesticides were diluted in acetone (0.3mL.L -1 of water) and
	demineralized water to obtain the two concentrations selected. Acetone (CAS: 67-64-1, Fisher
	Chemical, HPLC solvent) was purchased from Fisher Scientific (Illkirch, France). Control
	groups (C22 and C32) received only acetone.

The composition of the mixture, characteristics and concentrations of pesticides are presented in Table

555 556 Table 1 :

 5561 Composition of the pesticide cocktail

	813
	Chemicals	Family	Use	Action mode	Low Dose	High Dose
					µg/L	µg/L
	S-metolachlor	Chloroacetanilide Herbicide Mitosis perturbation	2.4	12.0
	Linuron	Urea	Herbicide Photosystem II inhibitor	2.0	10.0
	Isoproturon	Sustituted urea	Herbicide Photosystem II inhibitor	1.2	6.0
	Tebuconazole	Triazole	Fungicide Sterol synthesis inhibitor 1.2	6.0
	Aclonifen	Diphenyl ether	Herbicide Carotenoid inhibitor	0.8	4.0
	Atrazin-desethyl	Triazine	Herbicide Photosystem II inhibitor	0.4	2.0
	Pendimethalin	Dinitroaniline	Herbicide Mitosis perturbation	0.4	2.0
	Total				8.4 µg/L	42 µg/L
	814

Table 2 .

 2 Best Generalized Linear Models (GLM) explaining physiological and behavioural816 traits in goldfish exposed to temperature increase and/or pesticides. Sex was included as

	817
	818	explanatory variable but was not retained in any model (see results). Separate models on
	819	males and females were run for GSI. Non parametric analyses were run on other traits and
	820	presented in the text.
			Temperature	Pesticides	Temp*Pesticides
			LR χ 2	P	LR χ 2	P	LR χ 2	P
	Organs and body condition
	Hepatosomatix Index HSI	4.41	0.036	-	-	-	-
	Male Gonadosomatic Index GSI	-	-	-	-	-	-
	Female Gonadosomatic Index GSI	0.30	0.58	9.72	0.0078	9.46	0.0088
	Body mass loss	6.91	0.0086 -	-	-	-
	Immunity
	Lysozyme activity	-	-	-	-	-	-
	Peroxidase activity	0.0002	0.98	5.68	0.058	6.67	0.035
	Behavior
	Activity	38.31	<0.001 10.47	0.0053	-	-
	Speed		4.00	0.045	-	-	-	-
	Trajectory changes	5.86	0.015	-	-	-	-
	Exploration	9.29	0.0022 -	-	-	-
	Feeding	5.84	0.016	9.27	0.0097	-	-

concentrations (Control in white, Low Dose in grey, and High Dose in black) and at two temperatures (22°C on the left and 32°C on the right). A high index indicates a high degree of cellular damage.

Data are mean ± SEM. Different letters above bars indicate significant differences between groups

Acknowledgements

We thank three anonymous reviewers for helpful comments. We are thankful to A Aguettaz, R Munier, A Delgizi, F Clergeaud and S Papillon for technical help. We thank Dr J-P Cornec and Dr X. Moreau for help on histological investigation and analysis. A Gandar was supported by a PhD fellowship from the French Ministry of Research. L. Jacquin was supported by a MCF Idex starting grant at Toulouse University (IDEX-V5RJACQ). This work was supported by an AEAG grant (PHYPAT) and a CNRS EC2CO grant (French national program ECODYN). The EDB laboratory is part of the LABEX TULIP (ANR-11-IDEX-0002-02).

Supplementary Material

Click here to download Supplementary Material: Supplementary Material bis.docx Supplementary Table S1

Supplementary